University of Southern California Dissertations and Theses

Simulating electrical stimulation and recording in a multi-scale model of the hippocampus

(USC Thesis Other)

Simulating electrical stimulation and recording in a multi-scale model of the hippocampus

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

SIMULATING ELECTRICAL STIMULATION AND RECORDING IN A MULTI-
SCALE MODEL OF THE HIPPOCAMPUS

Clayton Scott Bingham

A dissertation submitted to the Faculty of the USC Graduate School in partial fulfillment of
the requirements for the degree of Doctor of Philosophy in Biomedical Engineering at the
University of Southern California

December 2019

ii

TABLE OF CONTENTS
1 INTRODUCTION
2 AM-NEURON MODELING
2.1 ABSTRACT
2.2 INTRODUCTION
2.3 METHODS
2.3.1 Defining the anatomical coordinate space
2.3.2 Admittance Method model construction
2.3.3 Tissue and electrode capacitance
2.3.4 Heterogeneous resistivity
2.3.5 Overview of the large-scale NEURON model
2.3.6 Granule cells
2.3.7 Entorhinal projections to the dentate perforant path
2.3.8 Synapses
2.3.9 Interface between NEURON and Admittance models
2.3.10 Model validation through Local Field Potential prediction
2.3.11 Minimizing power required to achieve population spiking
2.3.12 Multi-objective evaluation of stimulating conditions
2.4 RESULTS
2.4.1 General spatiotemporal properties of threshold
response to changing stimulating location
2.4.2 LFP predictions and qualitative comparison with
experimental recordings for model validation
2.4.3 Sensitivity to changes in resistivity and capacitance
2.4.4 Minimizing power required to achieve population spiking
2.4.5 Recruitment of axons and granule cells
2.4.6 Population spike peak amplitude and sharpness
2.4.7 Power efficiency of activity
2.5 DISCUSSION
2.5.1 Model improvements over previous work
2.5.2 Computational resources and considerations

8
10
10
11
16
16
17
18
19
21
21
22
22
23
23
24
25
28
29

30

33
35
36
36
37
37
38
38

iii

2.5.3 Tissue capacitance and current division at the granule cell layer
2.5.4 Power efficient elicitation of population spiking
2.5.5 Model Limitations and Future Developments
2.6 CONCLUSION
2.7 SUPPORTING INFORMATION
3 GRAPH-BASED MODELS OF CORTICAL AXONS
3.1 ABSTRACT
3.2 INTRODUCTION
3.3 METHODS
3.3.1 Volume and Nodal Constraints
3.3.2 Constraining Patterns of Axon Branching
3.3.3 Synaptic Boutons
3.3.4 Arbor Simplification and Computational Complexity
3.3.5 Alternative Generative Models
3.3.6 Strength-Duration Relationship in Response to
Extracellular Electrical Stimuli
3.3.7 Data/Model Sharing
3.4 RESULTS
3.4.1 Algorithm Assessment and Validation
3.4.2 Arbor Simplification and Computational Complexity
3.4.3 Alternative Generative Models
3.4.4 Strength-Duration Relationship in Response to
Extracellular Electrical Stimuli
3.5 DISCUSSION
3.5.1 Limitations and Alternative Methods
3.6 SUPPORTING INFORMATION
4 ADMITTANCE METHOD FOR ESTIMATING LOCAL FIELD
POTENTIALS GENERATED IN A MULTI-SCALE NEURON MODEL OF
THE HIPPOCAMPUS
4.1 ABSTRACT
4.2 INTRODUCTION
39
40
41
42
43
47
47
48
52
53
54
56
57
58
60

61
62
63
65
68
71
73
73
73
78
81

81
82

iv

4.3 METHODS
4.3.1 Construction of the NEURON Model of the Hippocampus
4.3.2 AM-NEURON Feedforward Simulation
4.3.3 Analytical LFP Estimates as a Proxy Ground Truth
4.3.4 New Meshes for LFP Estimation: tetrahedral vs. hexahedral
4.3.5 Fundamentals of Resistor Network Construction
4.3.6 Incorporating incident endogenous current sources: splitting vs. shifting
4.3.7 Model Sharing and Software Distribution
4.4 RESULTS
4.4.1 New Meshes for LFP Estimation: tetrahedral vs. hexahedral
4.4.2 Incorporating incident endogenous current sources: split vs. shift
4.4.3 Additional Comments on Estimation Accuracy
4.5 DISCUSSION
4.5.1 Spatial Resolution of Sources and Meshes
4.5.2 Mesh Quality
4.5.3 Solvers and Model Complexity
4.6 CONCLUSION
5 CONCLUSION
6 ACKNOWLEDGEMENTS
7 REFERENCES

84
85
86
89
91
91
95
97
98
99
100
101
102
103
104
104
105
106
109
110

v

TABLE OF FIGURES
Fig. 1.1
Conceptual depiction of the Admittance Method-NEURON
model

Fig. 1.2 Resistivity values for each section within the hippocampus

Fig. 1.3. Model volumes in the Kjonigsen rat brain atlas

Fig. 1.4 2-D rendering of the detailed neuronal model

Fig. 1.5 Predicted LFPs

Fig. 1.6 Heatmaps of spatiotemporal patterns of activity

Fig. 1.7 Raster plot of neuronal activity

Fig. 1.8 Predicted LFPs as stimulation amplitude is increased

Fig. 1.9 Composite figure showing impact of varying resistivity

Fig. 1.10 Composite figure showing impact of electrode location

Fig. 1.11 Population spike features at different electrode placements

Fig. 1.12 Multi-Objective Optimization Function output

Fig. S1.1 Strategy for implementing tissue-electrode interface

Fig. S1.2 Raster plot of network activity

Fig. S1.3 Descriptive bar chart of results feeding MOF

Fig. S1.4 Descriptive bar chart of values of U in MOF

Fig 2.1 Axon modelling overview

Fig 2.2. Features of granule cell arbors

Fig 2.3 Flow diagram of ROOTS

Fig 2.4 ROOTS operational efficiency

Fig 2.5 NEURON model simulation complexity

Fig 2.6 Ramer-Douglas-Peucker (RDP) line simplification

vi

Fig 2.7 Example graph models of entorhinal cortical axons

Fig 2.8 ROOTS/TREES comparative analysis: overview

Fig 2.9 ROOTS/TREES comparative analysis: branch order

Fig 2.10 Rendering of a perforant path axon with presynaptic boutons

Fig 2.11 Strength-duration curves

Fig. 2.12 Kernel density estimates of compartmental APs

Fig. S2.1 Pseudocode for ROOTS

Fig. S2.2 Clustering via K-Means

Fig. S2.3 Delaunay triangulation to construct a surface

Fig. S2.4 Djikstra’s ‘most likely path’

Fig. S2.5 Example graph models from TREES

Fig. 3.1 Demarcation of recording/stimulating devices

Fig. 3.2 Current source data

Fig. 3.3 AM-NEURON equivalent circuit

Fig. 3.4 Mesh geometries

Fig. 3.5 Example meshes for LFP estimation

Fig. 3.6 Explanation of cotangent rule for tetrahedra

Fig. 3.7 AM-NEURON bidirectional interface flow diagram

Fig. 3.8 Spatial node alignment

Fig. 3.9 Analysis overview diagram

Fig. 3.10 Impact of mesh geometry

Fig. 3.11 Impact of current source handling algorithms

Fig. 3.12 Estimation accuracy of numerical LFP solutions

Fig. 3.13 Example LFP time-series

vii

Fig. 3.14 Bidirectional AM-NEURON communication via splitting

Fig. 3.15 Bidirectional AM-NEURON communication via shifting

8

1 INTRODUCTION
When communicating with neural systems, there are few in vivo techniques more
widely used than electrode-based electrical stimulation and recording. One of the greatest
feats of neural engineering has been the development of implantable medical devices which
use multielectrode arrays to record and stimulate neurons for the therapeutic treatment of
brain disorders, whether through modulation (i.e. deep brain stimulation) or restoration (i.e.
hippocampal memory prosthesis) [1, 2]. Advancements in medical implant technology have
a trajectory toward “less and more”: smaller electronics and higher density arrays [3, 4, 5,
6]. However, the effectiveness of these engineering advances cannot be maximized without
a fundamental understanding of the neural system that is the basis of the brain disorder,
which includes the electrophysiology and the anatomy of the neuronal network and
principle neurons. This knowledge can inform the design and placement of electrodes to
determine the ideal location, orientation, number of electrodes, and stimulation pattern
necessary to elicit the target response [7, 8, 9].

Advancements in computing have opened doors for more detailed and accurate
mechanistic models of neural tissue, and an impressive corpus of anatomical and
biophysical literature is available to constrain them. This dissertation presents the
construction and co-simulation of a large-scale compartmental model of an in vivo
hippocampal neuronal network and a specialized equivalent-circuit model of the tissue

9

volume-conductor to approximate the spread of electrical currents through extracellular
space [10]. This powerful dual model allows the prediction of neural interactions and
responses to a variety of clinically relevant devices and stimulation protocols [11]. Prior
work with multielectrode arrays has been done to predict the hippocampal output given a
known input [2]. The proposed model will aid the design of multielectrode arrays and
stimulation parameters which would enable reproduction of predicted output in a target
neuronal population, thereby closing the loop and taking another step toward effective
cortical neuro-restorative technology in the hippocampus.

In the midst of constructing this model, important limitations were confronted and
ameliorative solutions were designed to bolster the predictive power of such mechanistic
models. This dissertation sets forth, in three chapters the fundamentals of mechanistic and
parametric modelling of dual neuronal and volume-conduction models of complex neural
tissue systems (chapter 2), generative models of highly-branched axon terminal arbors for
use in grey matter stimulation models (chapter 3), and a candidate algorithm for
bidirectional interactions between neuronal and volume-conduction domains (chapter 4).
Chapters 3 and 4 represent independent efforts to address commonly encountered
limitations—that were identified in the work discussed in chapter 1—in widely used
stimulating and recording models of neural tissue. These models are often used to study
deep brain stimulation and other central nervous system electroceuticals as well as more
general hypothesis-driven investigation of neural circuit dynamics [12, 13].

While the generality of these modelling tools is clear from the literature, we use
the hippocampus and potential hippocampal prostheses as a case study. In addition to
gaining a deeper understanding of the hippocampal system features which helps form the
basis of long-term memory formation, we expect this work to yield an efficient method of
evaluating arrays and patterns of stimulating currents. This outcome represents an
important step toward more effective, well-understood, and safe intracortical neural
implants.

10

2 AM-NEURON MODELING
2.1 Abstract
The ideal form of a neural-interfacing device is highly dependent upon the anatomy of
the region with which it is meant to interface. Multiple-electrode arrays provide a system
which can be adapted to various neural geometries. Computational models of stimulating
systems have proven useful for evaluating electrode placement and stimulation protocols
but have yet to be adequately adapted to the unique features of the hippocampus. As an
approach to understanding potential memory restorative devices, an Admittance Method-
NEURON model was constructed to predict the direct and synaptic response of a region of
the rat dentate gyrus to electrical stimulation of the perforant path. A validation of
estimated local field potentials against experimental recordings was performed and the
results of a bi-linear electrode placement and stimulation amplitude parameter search are
presented. The parametric analysis presented herein suggests that stimulating electrodes
placed between the lateral and medial perforant path, near the crest of the dentate gyrus,
yield a larger relative population response to given stimuli. Beyond deepening our
understanding of the hippocampal tissue system, establishment of this model provides a
method to evaluate candidate stimulating devices and protocols.

11

2.2 Introduction
In communicating with neural systems, there are few in vivo techniques more widely
used than electrode-based stimulation and recording. One of the greatest feats of neural
engineering has been the development of implantable medical devices which use multi-
electrode arrays (MEAs) to record and stimulate neurons for the therapeutic treatment of
brain disorders, whether through modulation (i.e. deep brain stimulation) or restoration (i.e.
Fig. 1.1. A conceptual depiction of the Admittance Method-NEURON model. A
compartmental model of a neuron is represented as a cable model with
compartments connected in series or in parallel according to the morphology. Each
compartment is represented by a circuit model. Voltages calculated by the
Admittance Method are applied as extracellular batteries to the appropriate
compartments to stimulate a NEURON model. Experimentally determined resistivity
values are used for subregions of the hippocampus [48].

12

proposed hippocampal memory prosthesis) [1]. Advancements in medical implant
technology have a trajectory toward “less and more”: smaller electronics and higher density
MEAs [4, 6, 5]. However, the effectiveness of these engineering advances cannot be
maximized without a fundamental understanding of the neural system that is the basis of
the brain disorder, which includes the electrophysiology and anatomy of the neurons and
their network topology. This knowledge can aid the design of MEAs by indicating the ideal
location, orientation, number of electrodes, and stimulation pattern necessary to elicit a
target response.

Damage to the hippocampus can result in memory impairment which profoundly affects
the quality of life; to overcome these defects, a hippocampal memory prosthesis is being
developed, and great strides have been made toward the prediction of hippocampal activity
[1]. In this application, thorough quantitative evaluation of electrode arrays in terms of
ability to efficiently elicit a targeted response is of particular importance due to the
intricate, yet well documented, anatomy of the hippocampus. This intricacy may make it
necessary to implement varied stimulating or recording electrode array geometries and
complex pulse-patterns to elicit or detect a specific pattern of response from the excitable
tissue. The ability to identify or generate a desired response can be dramatically enhanced
by electrode arrays that are appropriately placed.

13

The difficulty of successfully identifying and quantifying neural activation by analyzing
data from low-density LFP recordings, along with increased computational ability in recent
years, has led to the construction of multi-scale computational models that include
stimulating and/or recording electrodes. Virtual stimulation experiments of singular or
simple network neuronal models have been used to study interactions between neural tissue
and externally-applied electromagnetic fields. The modeling focus varies from prosthetic
electrode design and effects of external stimulation [14, 15, 16, 17, 18, 19, 20, 21, 22, 23],
to assessing the accuracy of modern methodologies and parameters [16, 24, 25, 26, 27].

While these models have yielded promising electrode and modeling insights, more
hippocampus specific biological realism is needed to satisfactorily evaluate electrode arrays
for their ability to elicit particular spatiotemporal patterns of neural activation in the
complex neural system. A virtual system is being developed with the goal of evaluating the
efficacy of various electrode positions and stimulus patterns using a large-scale, three-
dimensional Admittance Method-NEURON model of the dentate gyrus [14, 15, 28, 29, 30,
31].
Fig. 1.2. The resistivity values used for each section within the hippocampus [48].

14

At a larger scale, the representation of extracellular space is vital to accurate prediction of
the spatial distribution of current throughout the tissue and devices of a stimulating or
recording system; therefore, accurate dielectric properties for discretizing the space into
circuit components and distinct anatomical representations of the tissue are necessary [32,
25].

Grill and McIntyre have made it clear that explicit neural models are advantageous
because Rattay-style activating functions fail to accurately predict sites of action potential
initiation [22, 33]. While there are a range of possible individual neuronal models, with
varying levels of detail being necessary to accurately recreate certain behaviors, it is
unclear whether neural circuit responses to extracellular electrical stimulation can be
accurately predicted with models much less complex than those described in this paper
[34]. Due to preferential activation of neural subunits (dendrites, somas, and axons), where
axonal chronaxie may be up to 40 times smaller than somatic, it seems important to
implement complex and realistic morphologies which allow such spatial and anatomical
specificity of activity [22, 33, 35]. Additionally, current source-sink dynamics between
branching dendrite and cell body layers of laminated tissue create measurable dipoles in
evoked potentials, features of which can be used for model validation against experimental
recordings.

This work presents a systematic exploration of electrode location and stimulus amplitude
in the dentate gyrus, which yields insights on the combined and individual influence of
anatomy and stimulus amplitude on the population response. The availability of the state of
every neuron, currently infeasible outside an in silico setting, allows the development of
new metrics to describe the spatiotemporal network response and provides an initial
surrogate investigation on describing activity in vivo, should the measurement of data at
this scale become possible.

15

Fig. 1.3. The model dimensions and subdomain volumes change from slice to slice
as can be seen in the Kjonigsen rat brain atlas [37].

16

2.3 Methods
This paper presents the study of a model that combines both three-dimensional electric
field calculation, via the Admittance Method (AM), and a NEURON model comprising
50,000 granule cells from Hendrickson et al., and 10,000 simplified perforant path axons
[36]. These components are spatially arranged to approximate a three-dimensional region of
the dentate gyrus. This combination Admittance Method-NEURON model is a powerful
approach for computing neuronal response to extracellular fields because of its potential for
modeling of not only the extracellular fields created by electrodes but also to provide a
high-resolution view to the behavior of neural processes within that extracellular field.

2.3.1 Defining the anatomical coordinate space
The hippocampal model presented herein was built following the segmented digital
atlas of Kjonigsen et al. [37]. Their atlas contains high resolution images (1.4 pixels/µm) of
30 µm thick slices through the coronal plane of the rat brain, with each image digitally
annotated to show the subdivisions of the hippocampus. While there are other established
rodent brain atlases, such as those of Swanson or Paxinos, the present work is constructed
using the Kjonigsen atlas because of the marginally higher resolution of their series of
coronal slices, marked ease of use, and more complete region demarcation in the
hippocampus [38, 39]. A custom Matlab function was written to threshold each image,
extracting the segmentation and converting the atlas to a voxelized format, where each
voxel was assigned an index corresponding to its subdivision. A particularly clear slice
image was taken at -3.34 mm from the bregma and the extracted boundaries were extruded
400 µm in the septotemporal direction, forming the entire 3-dimensional space within
which the AM-NEURON model was constructed.

17

2.3.2 Admittance Method model construction
An in-house multi-resolution variant of the quasi-static AM was used to simulate the
resulting voltage at nodes throughout this model due to given current input(s) [14, 15, 29,
11, 10, 31].

In this method, a hexahedral, structured mesh is constructed in which the original voxel
size is maintained near boundaries or predefined current source locations. As the distance
from boundaries and sources is increased, voxels are merged. A maximum merged element
size of 64 voxels was used. The original resolution of 12 m-wide cubic voxels was
implemented within the dentate gyrus, maintaining fine resolution in the tissue containing
the superimposed NEURON model of the dentate gyrus. Once a mesh is constructed, it is
overlaid with linear lumped circuit representations of resistors and capacitors based on the
dielectric properties of the tissue and electrodes. Current sources are added into this
network, which may be defined between any two nodes in the model. The resulting circuit
was solved with the bi-conjugate gradient method with a residue of 10
-8
as the stop criteria.
Further details are provided in Cela [29].

The AM was chosen over alternative methods, such as the Finite Element Method
(FEM), due to case-specific advantages exploited in this work. First, the method is
relatively simple to implement relative to FEM. While FEM has the advantage of non-
conformal grids, computational complexity grows with the number of independent current
sources. Using a circuit model allows efficient inclusion of multiple current sources. This is
a useful feature for modeling multi-electrode arrays and estimating feedback potentials.
Second, the circuit description of the model space leads to an intuitive method for
describing bulk neural tissue. This creates a more seamless integration with neural models,
also represented by circuit components, by allowing for voltage at nodes within the circuital
network of the volume conductor to be applied as extracellular voltage sources for neuronal
models. This feature produces a coupled circuital description of bulk tissue and underlying
neural networks. Further, arbitrary circuit components may be added, allowing for the
inclusion of anisotropic dielectric properties, intuitive implementation of various electrode

18

features, circuit representations of neurons in this bulk tissue model, or other phenomena
that can be approximated as an additional circuit, such as the electrode-electrolyte interface
[40] (Fig. S1).

2.3.3 Tissue and electrode capacitance
Due to the low-frequency nature of stimuli (i.e. 95% of Fourier spectrum under 1 kHz,
long pulse-width, and low pulse-frequency) in the simulations performed, introducing
tissue capacitance to the model should yield negligible differences in results [23, 8, 26, 33,
41, 42, 43, 44, 45, 46]. To validate this assumption, various models were simulated with
capacitive tissue properties following Gabriel et al. Further, the basic validity of this
assumption was analyzed by calculating the τ RC of a Cole-Cole representation of tissue bio-
impedance [47]. τ RC is the time which a capacitor charges approximately 63% across a
resistor and can be found by multiplying the resistance and capacitance in an RC circuit, or
resistivity and permittivity in a volumetric model. Based on this analysis and simulation
results presented later, it was deemed appropriate to take advantage of the relative
computational simplicity of purely resistive models, and capacitive properties of bulk tissue
were omitted.

A capacitive layer was added between the tissue and electrode, with a double-layer
capacitance of 10 µF/cm
2
. The method for applying a specific capacitance across the
interface at the electrode surface was to distribute the capacitance into parallel capacitors,
with one applied to each node on the surface of the electrode. Each new capacitor was
placed perpendicularly to the electrode surface and in parallel with tissue-electrode
interface resistors. This was done without an additional spatial component, as the thickness
of the electrode-electrolyte layer is much smaller than can be realized in our model with its
current spatial resolution.

19

2.3.4 Heterogeneous resistivity
The voxelized model was discretized based on the resistive properties of the hippocampal
tissue at low-frequency (<10 kHz). Due to the unavailability of impedance measurements in
the dentate gyrus, model resistivities were assigned based on measurements taken in CA1,
using the four-electrode impedance measurement setup, by Lopez-Aguado, Ibarz, and
Herreras [48]. Many measurements were performed by the experimenters in the apical and
basal dendritic regions, cell body layer, stratum lacunosum-moleculare, and the
hippocampal fissure bordering the outer boundary of the dentate gyrus. The measurements
showed remarkably similar resistivity in the molecular regions, around 2.7 Ω m, with
approximately double the resistance in the cell-body regions. This observation led to the
simplification of the computational model, splitting the hippocampus into four regions,
three molecular and one cell body, with each respective region having homogeneous
resistive properties based on these measurements. The values implemented in the model are
those reported by Lopez-Aguado et al. and are denoted in Figs. 1,2. Lopez-Aguado
concluded that in the CA1 it should not be expected that a somatic to dendritic resistivity
ratio be less than 2.22, though there is no clear indication whether a similar ratio can be
expected in the dentate gyrus.

Because changes in this ratio would cause an increase in the gradient of the
potential and may impact the activation of nearby cells, several simulations were performed
with varying levels of heterogeneity. The original ratio of the cell body layer to molecular
layer was 2.28 (control). This was varied in a series of simulations from cell body
layer/molecular layer ratio of 1:1 to cell body layer/molecular layer ratio of 5:1.

20

To simulate electrical stimulation, two insulated microwire electrodes and a
reference electrode were added to the model in replication of the work of Soussou et al. as
shown in Figs. 1,2 [49]. The electrodes had a diameter of 48 µm and were given resistivity
of platinum (10
-8
Ω m). The insulation encasing the electrodes was 24 µm thick, with
dielectric properties of teflon (10
16
Ω m), and was cut flush with the surface of the
electrodes. A bipolar current source was then applied, with the anode being placed at a
node in the center of one of the micro-wires and the cathode in the other. All simulations
were run with a bi-phasic, square-wave pulse, with a width of 1 ms. The results were tri-
linearly interpolated and applied to each section in the large-scale compartmental
(NEURON) model to predict induced neural activity. Potentials charge the capacitive
membrane, resulting in intracellular currents.

Fig. 1.4. 2-D rendering of the volumetric detailed neuronal model with demarcation
of recording and stimulating devices. Notice the laminar, en passant topographical
features of the perforant path; these are a strong determinant of the behavior of the
tissue in response to extracellular electrical stimulation.

21

2.3.5 Overview of the large-scale NEURON model
The neural response to extracellular electric fields is simulated using a biologically
realistic, computational model of a 400 µm coronal region of a rat dentate gyrus (DG).
Assuming a total hippocampal length of 10 mm and 1,200,000 total GCs, the model
comprises 50,000 granule cells (GC) and 10,000 entorhinal cortical (EC) axons, many
features of which have been previously described [36, 50]. The NEURON simulation
environment (7.3) is used to simulate the compartmental models. The model is parallelized
and run on a 4,040-processor computing cluster supported by the High-Performance
Computing and Communications Center at University of Southern California.

2.3.6 Granule cells
The morphologies of the GC models were generated using the L-NEURON tool [51].
Distributions of structural parameters were extracted from a database containing
morphological reconstructions of GCs, and the distributions of parameters were used to
create unique structures for the GC models [52]. For each generated morphology
compartment lengths were kept uniform. The current morphologies include the dendrites
and somas of the GCs but not the axons.

The biophysical parameters used in the GCs were derived from previously
published computational models of the dentate gyrus [53, 54]. These specifications include
the types, densities, and distributions of the ion channels throughout the neuron model. The
electrophysiological parameters and characteristics of the GC models used in this work
have been previously validated [36].

22

2.3.7 Entorhinal projections to the dentate perforant path
While the eventual goal of this work is to model perforant path fibers along with
their actively signaling parent cells, this current iteration models the perforant path as
singly-branched cables with in vivo Hodgkin-Huxley biophysics. The most recent
morphometric description of axons in the molecular layer of the DG report very narrow
axons in spiny stellate EC cells, the most numerous cell type of the perforant path, with
diameters as small as 0.1µm in between “periodic varicosities”, or thickenings, that likely
correspond to presynaptic boutons [50]. The remaining topological and morphological
features of the axons implemented in this model are simplistic approximations of the axon
arbor statistics previously discussed [36]. A few other well-known and prominent features
of the axons are implemented, namely an en passant connective schema, a primary
bifurcation just above the crest of the DG, and the segregation of the medial and lateral
perforant path, respectively, into the middle and outer thirds of the molecular layer of the
DG. Proximity rules of connectivity between the perforant path fibers and GC dendrites are
implemented to ensure appropriate conduction delays.

2.3.8 Synapses
EC axons are connected to the GCs via 6,000,000 synapses in the outer third of the
molecular layer for lateral perforant path axons, and 5,625,000 synapses in the middle third
of the molecular layer for medial perforant path axons. Upon being triggered by an action
potential, the synapse undergoes a change in conductance which allows current to flow into
the post-synaptic compartment. The time course of the conductance (g) is approximated as
the difference of two exponentials (Eq. 1):

g = e
−t
τ
1
− e
−t
τ
2
(1)

23

The postsynaptic response between entorhinal axons and GCs in this model is
mediated exclusively by AMPA receptors. The maximal conductance for the synapse was
chosen so that the amplitude of a single EPSP at the soma is 0.22 mV, with average τ 1,2 of
1.05 and 5.75, respectively [54].

2.3.9 Interface between NEURON and Admittance models
Because the AM and NEURON models share the same coordinate space, the
voltages from the AM are superimposed on the NEURON model to drive the compartments
of the neurons. Through the addition of voltage sources in series with each membrane,
extracellular potentials may be applied to neuronal compartments within the NEURON
model using the following relationship (Eq. 2):

I
m
(t) =
dV
ext
dt
∗ C
m
(2)

Where I
m
is the transmembrane current which results from an extracellular potential (V
ext
)
charging the membrane capacitance (C
m
). Potentials were applied throughout GC models
at the center of each NEURON section (20-30 locations), and each perforant path axon had
voltage applied to 100 different locations evenly distributed throughout the extent of their
structure. The circuit diagram for the application of the external voltage is shown in Fig. 1.

2.3.10 Model validation through Local Field Potential prediction
To evaluate the biological realism of the model, evoked potentials were estimated
and compared with experimental recordings [49]. Transmembrane current densities were
output at each simulated time-step of stimulation applied near the crest of the model DG
(location 2 in Fig. 4). Stimuli used for comparison included biphasic, charge-balanced, and
square-wave impulses of 1ms width and 200μA amplitude. Field potentials were then

24

estimated for virtual point-electrodes located along the green dots in Fig. 4 using a variant
of the line-source method, where resistivity was calculated as the path-length weighted
average of domain resistivities (Eq. 3 and 4) [55, 56].

ϕ(x , y , z) =
1
4π
∑
I
i
K
i
n
i=0
(3)

K = ∑ σ
j
∗ r
j
n
j=0
(4)

Where Φ is the field potential resulting from a current source, I. K is the discrete
integral of the product of tissue conductance and radial distance (sigma and r, respectively),
where r is found by Eq. 5:

r = √(x− x
0
)
2
+ (y− y
0
)
2
+ (z − z
0
)
2
(5)

2.3.11 Minimizing power required to achieve population spiking
A survey of stimulation parameters was conducted to determine which conditions
best minimized power while still eliciting a concert of activity from the tissue. For each of
the nine electrode locations, designated in Fig. 4, biphasic square-wave impulses (1 ms
pulse width) were delivered over a range of amplitudes (50-650 µA) in increments of 50
µA. The network responses were analyzed for several features intending to thoroughly
describe magnitude, time-course, and velocity of the population response. For this
particular portion of the analysis direct activity, or spiking activity observed prior to the
arrival of synaptic events nearly concurrently with the stimulus, was deliberately excluded.
The remaining indirect, or synaptically driven, activity was measured for total proportion of
cells active, maximum instantaneous proportion activated (2 ms time-bins), maximum
slope of the network activity cumulative density function, as well as the half-height width
of the indirect response.

25

The concept of a population spike (PS) has been a popular one with electro-
physiologists and especially those who work with the hippocampus. The architecture of the
DG, characterized by the common orientation and dense packing of GCs, leads to a
summing effect in the extracellular potentials recorded when EC afferents are stimulated
concurrently [57]. Predictions made in the present study provide clues to the underlying
electrical activity which results in PS observed in local field potentials. While this
phenomenon can be observed via predicted evoked potentials (as in Fig. 5), having access
to the membrane potential of each individual compartment in the NEURON model provides
a level of detail only partially captured by experimentally recorded field potentials. For the
purposes of these results, the PS is a synaptically activated and temporally phasic response
expressed as an aggregation of spiking events observed in the model and not due to direct
activation of granule cells in the dentate gyrus by the stimulation event. These two phases
of the response are easily separable because of the delay associated with the axonal
activation, synaptic conductance, and later suprathreshold activation of granule cells in the
tissue. Fig. 7 demonstrates the features of the PS that are measured in this study, which
combines 117 simulations over a range of stimulation amplitudes in 9 different electrode
locations.

2.3.12 Multi-objective evaluation of stimulating conditions
There are many ways a neural response could be characterized, but knowing which
metric correlates with positive clinical outcomes or a successful experiment is impossible
without defining the specific problem. The following metrics represent several important
features that could generally be useful to optimize: a large PS amplitude would be
necessary to ensure the successful transmission of the response to afferent areas; a short
half-height width would allow greater temporal precision and potentially increase the
temporal resolution at which a particular neural response could be encoded; finally, for a
given pulse-width, a greater power efficiency could be important for safety and to increase

26

the battery life of an implant. A simple, multi-objective optimization function was
constructed using the weighted sums method to determine how each stimulation case would
satisfy these constraints (Eq. 6):

U =
w
1
PS
max
+w
2
PS
efficiency
+w
3
[max (HHW)−HHW]
w
1
+w
2
+w
3
(6)

Where PS
max
is the peak proportion of GCs active, PS
efficiency
is the number of
GC spikes per unit stimulus, HHW is the half-height width of a probability density function
fit to the instantaneous activity (Gaussian), and w1,2,3 are the weights associated with each
metric. Each metric is an independently normalized input to the multi-objective function. In
this no-preference formulation, the weights were naively assumed to have a value of one,
such that the average of the normalized metrics was computed. However, if the importance
of each objective were defined by a decision maker or an automated method based on the
situation, weights could be assigned a priori to yield different results. As formulated, high
values of U indicate elicitation of a distinct and efficient population response (high peak
activity, highly concurrent activity, and high activity per unit stimulus amplitude).
Fig. 1.5. Predicted LFPs (magenta signal) exhibit a characteristic PS as observed in
experimental evoked potentials (black signals) elicited by a 200 μA biphasic, square-
wave pulse [43]. “Positive-going” spikes in the dendritic domain and “negative-
going” spikes in the cell body layer indicate the emergence of a dipole. This is due to
the laminar spatial separation of excitatory post-synaptic potential and action
potential sites of origination.

27

Fig. 1.6. Heatmaps of spatiotemporal patterns of PS activity given the different
locations of stimulation along the transverse axis (time calibrated from time of
stimulation). (A, suprapyramidal PP at 100μA; B, crest PP at 150μA; C,
infrapyramidal PP at 150μA). Activity along the transverse axis is strongly
dependent upon the principle en passant topology of the perforant path.

28

2.4 Results
Once the model was constructed, simulations were performed to determine general
features of the response. The results of more specific studies performed after these test
simulations comprise three main parts: (i) prediction of evoked field potentials, (ii) an
analysis of model sensitivity to changing bulk tissue dielectric properties, and (iii) an
improved stimulating parameter search where the AM-NEURON model was tested for
sensitivity to stimulating location and pulse amplitude with the goal of eliciting a power
efficient PS.
Fig. 1.7. 50k cell/10k axon network response to a 300 µA biphasic pulse (-/+)
delivered by electrode two at time 50 ms. Each red dot represents an action potential
at some location along the transverse axis of the tissue (large Y-axis) in time.
Activity is also measured in terms of instantaneous (red line, 0.25ms bins) and
cumulative response (blue). The total response is temporally separated into a direct
activation phase (orange box) and a population spike (PS) phase (green box).

29

2.4.1 General spatiotemporal properties of threshold response to
changing stimulating location
Initial simulations yielded some nearly universal commonalities: once a critical
mass of PP axons reach threshold (always occurring in compartments most proximal to the
leading cathodal source in this model) PS activity propagates in both directions along the
transverse arc (Fig. 6,7). For Fig. 6, direct activation was excluded to emphasize the
propagating PS. The chief difference remaining between perforant path and cell body
stimulating locations is the activation threshold triggering the pattern that can be seen in
these plots. For PP locations, PS propagates successfully along the transverse axis between
amplitudes 100-150μA, where cell body electrodes required amplitudes >300μA.
Importantly, the topography of PP-dentate connections leads to much greater concurrency
and earlier termination of crest initiated PS.
Fig. 1.8. Predicted LFPs develop a dipole as stimulation amplitude is increased [43].

30

2.4.2 LFP predictions and qualitative comparison with
experimental recordings for model validation
Predicted evoked potentials exhibit characteristic “negative-going” signals in the
cell body layer and “positive-going” signals in the molecular domain (Figs. 5, 8). Upon
comparison with recordings from experimental studies, predicted waveforms appear to
capture the dipole behavior, PS amplitude, and time-course expected from the densely-
packed and highly laminar granule cells [49, 57]. The development of a characteristic
dentate evoked potential, with a dipole along the cell body-dendritic arbor axis, occurs
gradually as stimulus amplitude is increased to 200 μA (Fig. 8).

31

Fig. 1.9 (A). The cell body layers are highlighted in yellow. The white bounding box
shows the borders of the plots in Fig. 8.B. The green-dashed-arrow segment is the
base of Fig. 8.C. Cyan and magenta dots show the electrode locations for dielectric
sensitivity analysis in the perforant path (case #2) and the cell body layer (case #5),
respectively.
(B). Current density map for various levels of heterogeneity in a slice 18 µm beneath
the electrode. Increasing the heterogeneity blocks the current from spreading into the
cell body layer and the hilus toward the CA3 region. (Control = 2.28:1).
(C). Shows voltage difference between each ratio of cell body to molecular layer
resistivity experiment (1:1 in red, 5:1 in blue, control in dotted black ~2.28:1) and
the control from left to right across the green dotted arrow in Fig. 8.A. Stimulating in
the perforant path (top) results in smaller differences than does stimulating in the cell
body layer (bottom).
(D). Shows the NEURON model response to varying resistivity under otherwise
controlled stimulating conditions for two electrodes placed at the crest. (case #2, left
When stimulating in the perforant path of the crest, changes in cell body layer
resistivity have a very limited impact on the features of PS. (case #5, right) When
stimulating in the cell body layer of the crest, changes in resistivity can be seen to
have a very large impact not only on PS features but also on the number of directly
activated cells. (Control = 2.28:1).

Time (ms)
% Activated
Time (ms)
Case #2 Case #5
A
D
B
C
Case #2 and #5

32

Fig. 1.10. Spiking activity at
most locations saturates above
stimulus amplitudes of ~500 μA.
(A) When stimulating at cell
body locations, granule cells
activity versus stimulus
amplitude follows a sigmoidal
trend. (B) When stimulating at
PP locations, the trend is more
logarithmic. (C) The differences
between (A) and (B) gives rise to
greater efficiency at low
amplitude when stimulating in
the PP.

33

2.4.3 Sensitivity to changes in resistivity and capacitance
Fig. 9.A can be referenced as an anatomical guide to results of the sensitivity
analysis presented in 9.B, 9.C, and 9.D. Case #2 and Case #5 correspond to perforant path
and cell body layer electrodes, respectively. Changes in the relative heterogeneity of the
resistive properties of each region of the tissue being modeled results in measurable
changes in the time-course, spatial distribution, and magnitude of the resulting response.
Fig. 9.B/C demonstrates the difference in the spread of current as tissue heterogeneity is
varied between 1:1 and 5:1 (control=2.28:1). In the homogeneous experiments, the current
spreads uniformly through the tissue, penetrating well into the hilar region beyond the cell
body layer. As the heterogeneity of the tissue is increased, the cell body layer inhibits
current flow, acting as a conductive barrier, changing the gradient of the spread and
lowering the density of current within the hilus. The normalized root mean squared error
(NRMSE) of the voltage taken within the dotted bounding box shown in Fig. 8.A between
each experimental case and the control was taken at three different depths within the slice.
For the electrode placed in the perforant path, the NRMSE was less than 0.05 for all depths
and was approximately the same for both the 1:1 and 5:1 experiments. With the electrode
placed in the cell body layer the NRMSE increased to 0.17 for the 5:1 experiment and to
0.37 for 1:1. The NRMSE increases for the 1:1 case due to the increased spread of the
Fig. 1.11. For nearly all stimulating conditions the half-height-width (HHW) of PS
was shorter (ms) when stimulating at the cell body layer and PP at the crest relative
to the supra and infrapyramidal transverse locations. While more variable, nearly all
stimulating amplitudes resulted in larger peak PS amplitudes (%) at the crest relative
to supra/infra locations.

34

current channeled through the cell body layer. The changing dielectric properties primarily
result in changes of potential in the cell body layer but changes also occur in the hilar
region (Fig. 9.A in conjunction with Fig. 9.B/C).

Fig. 9.B shows how the magnitude of the difference between each experiment and
the control test varies across the magenta axis in Fig. 9.A. Placing the stimulating electrode
within the molecular layer of the crest led to a maximum difference of approximately 10
mV within the cell body layer for both the cell body layer/molecular layer 1:1 ratio and cell
body layer/molecular layer 5:1 ratio experiments. However, placing the electrodes within
the cell body layer led to an increase to 60 mV difference. The maximum difference in
voltage for the latter experiment was several times larger, indicating that increased
heterogeneity could lead to a larger difference in model results for this electrode placement.

When these cases were simulated in the NEURON model, the differences were
more apparent (Fig. 9.D). When considering the stimulation of the cell body layer of the
crest, the homogeneous case resulted in a net decrease of activity in all regions of the
Fig. 1.12. Output of the multi-objective optimization function for MPP/LPP
(left) and CBL (right) stimulation cases. High values of U indicate strong PS.
The population spike amplitude and power efficiency were maximized and the
half-height width was minimized with equal weighting in this optimization.
(CBL-cell body layer, LPP/MPP-lateral/medial perforant path). Trend-lines,
added to aid interpretation, are created using a 3
rd
order Savitzky-Golay filter.

35

NEURON model while increases in the resistivity of the conductive barrier feature may
ultimately result in increasingly distal sites of direct activation at a given stimulation
amplitude. The 5:1 resistivity ratio case yielded a net increase of total activity.

Regarding the inclusion of capacitance in the AM field estimation step, for the
model parameters outlined in section II.C, differences in the field potentials of less than
0.1% throughout the modeled tissue were observed in comparison to the model with no
tissue capacitance when each model was stimulated with a biphasic current pulse of 100 µA
magnitude and 1 ms in duration. Further, this assumption is intuitive when considering the
τ RC of a Cole-Cole representation of tissue bioimpedance [47]. The τ RC provides a
simplified estimate of what stimulating currents might be suitably approximated by purely
resistive models. Where τRC is much less than the pulse-width of a current-controlled
stimulus, it can be assumed that tissue capacitance is a negligible determinant of the
potentials around the circuit. For a 1 ms current impulse, with a resistivity of 2.5 Ω m and
permittivity of 2.79 x 10
-6
F/m, capacitors in the volume would charge to 95% in under 21
µs. A complete record of these experiments can be found in supplemental materials (Fig.
S2)

2.4.4 Minimizing power required to achieve population spiking
A key goal of this modeling effort was to discover improved electrode and
stimulation protocols for hippocampal tissue. This study presents results of a thorough
bilinear optimal parameter search, where an improved electrode location and stimulating
amplitude is selected with which to elicit a power efficient PS. The critical results of these
experiments are presented in Figs. 10, 11, and 12. A complete record of these experiments
can be found in supplemental materials (Fig. S3-4)

36

2.4.5 Recruitment of axons and granule cells
The response of the tissue to increasing stimulation amplitude for any location was
nonlinear. While it seemed that axonal recruitment increases in a disjointed manner, there
was a sigmoidal shape to response curves when stimulating in the cell body layer (10.A),
where PP location recruitment pattern was more obviously logarithmic with increasing
stimulus amplitude (10.B). The difference arises from lower thresholds for GC recruitment
in PP locations. Despite unremarkable axonal activation, it was notable that a slightly larger
proportion of GC activity was seen in suprapyramidal blade stimulating locations.

2.4.6 Population spike peak amplitude and sharpness
The population spike peak amplitude and half-height width were calculated by
counting the number of spikes generated by the GCs within 2 ms bins. The peak amplitude
was normalized by the total number of GCs and converted into a percentage. The half-
height width corresponded to the total time during which the population activity was above
half of its peak amplitude. The proportion of granule cells active at peak PS was highly
location dependent, with a markedly greater number active at the crest. As stimulation
amplitude increases there is a general increase in peak proportion, but there was some drop-
off in several cases above 500 µA. This drop-off was exaggerated in the suprapyramidal
stimulating locations where it occurs earlier and more aggressively.

Electrodes at the crest were notable for both larger amplitude and shorter half-
height width PS from the tissue model (Fig. 11). There is no clear trend in PS half-height
width (measured in ms) as stimulation amplitude increases across all locations but there are
some electrodes, particularly at the crest and cell body layer in the suprapyramidal blade
where there was a strongly sloping, and supra-linear shortening of PS half-height width

37

2.4.7 Power efficiency of activity
Power efficiency was calculated by normalizing the peak PS amplitude and the
total proportion of GCs activated by the stimulation amplitude. For locations along the
transverse axis, as the stimulating current source is increasingly centered in between the
lateral and medial perforant path the response per unit stimulus exhibits an exponential
trend (10.C). Stimulating at the cell body layer, however, finds a maximum efficiency at
midrange amplitudes rather than high or low. The greatest efficiencies are seen in perforant
path locations below amplitudes of 250 µA, without regard for transverse extent. These
data were then used as inputs to Eq. 6 to calculate U for all stimulus location/amplitude
pairs. Fig. 14 demonstrates the PS efficiency is best at MPP/LPP locations of the crest at
low amplitudes (~50-100 μA).

2.5 Discussion
As a first step toward modeling the entire hippocampus, a model encompassing a
section of the dentate gyrus that also includes the perforant path axons has been completed.
Where other studies have focused on either the proper implementation of bulk-tissue-level
modeling without network dynamics or the results of much smaller networks using multi-
compartmental neuron models, this study demonstrates how a hybrid, multi-scale AM-
NEURON model reveals properties of hippocampal activity which only emerge with both
scale and detail. Studying the sensitivity of models, like those described in this paper, to
variations in dielectric properties and electrical stimuli will reveal fundamental
characteristics of the hippocampal system. This knowledge guides design of more effective
therapeutic interventions and devices.

38

2.5.1 Model improvements over previous work
By scaling the AM-NEURON method up to the tissue level, the model
demonstrates the capability to predict DG response to a variety of configurations of
stimulating electrodes and stimuli. Previous studies based on the NEURON model used in
this work have reported the existence of naturally emerging spatiotemporal clusters of
activity in the DG due to synaptic input [36]. This new model uses a digitized three-
dimensional hippocampal atlas to preserve the overall anatomy of the macrostructure and
the hippocampal subfields. This model adds the ability to provide an additional source of
input, via electrical stimulation, and includes explicit axons for the EC projections. These
additions grant the ability to explore the temporal and spatial components of localized PS
by allowing a precise view of active propagation along the transverse extent of the
structure, an unexplored topic due to the highly divergent nature of the excitatory EC
inputs. The inclusion of actively conducting axons that follow an anatomically correct
topography is an addition that has allowed more powerful prediction of the spatiotemporal
response to extracellular stimulation; where simplistic activating functions might have only
successfully predicted the immediate and local response of directly activated cells (i.e.
axonal sites of action potential initiation) a hybrid, multi-scale networked AM-NEURON
model can also show how patterns of activation spread due to the initial response to the
stimulation and the subsequent synaptic events that are triggered by propagating action
potentials. GC action potentials can be generated along the transverse axis in locations
several hundred microns away, where stimulating potentials approach zero.

2.5.2 Computational resources and considerations
Computational burden is an important challenge to consider when choosing the
modeling approach outlined in this study. Successful simulation of large networks of
detailed neuronal models requires considerable computational resources. To make this
study possible, the model was parallelized to take advantage of a 4,040 processor HPC
cluster owned by the authors. The simulations reported in Figs. 6 and 7 were run on 600
processors and each took 3.5hrs. to complete, producing 45 GB of data per stimulation

39

case. The PS sensitivity analysis consisted of 117 simulations of a scaled down (1:25 GCs)
version of the model presented in Figs. 10, 11, and 12 were run on 60 processors. These
simulations produced approximately 200 GB of data and required nearly 585 hours of
processor time to simulate just under 12 seconds of network activity. Despite the
interpretability of parametric models, ever increasing model complexity poses significant
computational challenges. Future models with more complete topology, longer simulations,
complex stimulation protocols, and incorporation of feedback potentials, will likely require
more resources than can be brought to bear without further optimization. Methods to reduce
the computational burden include nonparametric multiple-input/multiple-output modeling,
algorithms to simplify dendritic and axonal branching complexity, alternative numerical or
analytical volume conductor methods, and GPU reimplementation [58].

2.5.3 Tissue capacitance and current division at the granule cell
layer
Because of the unique anatomy of the hippocampus, it was clear that careful
consideration of the dielectric properties of bulk neural tissue would be critical to
performing accurate prediction of the network response.

While there is previous evidence of the negligible contribution of capacitance at
low-frequencies ( 1 kHz in this study), this assertion was affirmed for more numerously
populated neuronal models than has been published [33].

The ability of the AM to create a complex heterogeneous model is important in the
hippocampal system because of how the granule cell layer alters the path of injected
current, to avoid regions of high resistivity and changes the gradient of the electric field.
Results from a resistivity ratio sensitivity analysis demonstrate the large impact that
changes of this property can have on the network level activity of a complex tissue system.
As the relative resistivity of the cell body layer was raised, the conductive barrier became

40

more apparent. This study illuminated the role the DG cell body layer could play in
shielding hilar structures and how dynamic changes in resistivity could make this region
more or less vulnerable to stimulating current penetration. While this study presents some
exploration of the impact of varying tissue resistivity and capacitance, the dynamic nature
of neural tissue dielectric properties warrants further study [58, 59].

2.5.4 Power efficient elicitation of population spiking
In the 1970’s, Sam Deadwyler and others popularized the measurement of PS as a
metric of hippocampal tissue response to electrical stimulation [57]. This observation in
local field potentials was adapted as a rough approximation of the quantity and concurrency
of EPSP in granule cells. The present model and the analysis performed was designed to
provide a similar, though much higher resolution, quantitative evaluation of the activity
elicited by electrical stimulation. The proxy measurements of PS used in this study
demonstrate how this model could be used to determine an improved and power efficient
device and impulse with which to elicit PS from the principal cells of the DG.

The clearest conclusion to be made from the series of 117 simulations that explored
the sensitivity of the hybrid model to stimulation location and amplitude is that proximity
to dense regions of axons is the strongest physical system determinant of PS threshold. As
stimulating electrodes were moved closer to the zone between the lateral and medial
perforant path, regardless of transverse extent, a dramatic decrease in activation threshold
was seen (Fig. 10.A/B). In this region, the electrode distance to thoroughfare axons in the
tissue system is minimized. Second only to axonal proximity in determining PS size was
transverse location; electrodes placed at the crest saw both a marked increase in peak
amplitude and a decrease in half-height width of PS (Fig. 11). Also significant are the non-
linear relationships between the system response and stimulus amplitude. The magnitude of
GC response to increasing stimulation amplitude follows a roughly sigmoidal/logarithmic
trend where high sloping regions occur at amplitudes where large recruitments of axons
occur and low sloping regions occur where either few axons are near threshold or most

41

axons are already suprathreshold (Fig. 10.A/B). This indicates that there could be regions of
diminishing returns where each unit of stimulus amplitude results in progressively smaller
increases in the number of GC spikes. The decaying exponential trend seen in the perforant
path response (Fig. 10.C) indicates that fewer and fewer spikes are yielded with each
increase in amplitude across the whole range. This result indicates that power can be
minimized for any arbitrary PS size or half-height width.

When combining these observations from the analysis (as achieved by Fig. 12) it is
clear that if an experimentalist desired to elicit an efficient PS of respectable amplitude and
minimal half-height width then such could be achieved by stimulating between the lamina
of the lateral and medial perforant path at the crest of the DG at an amplitude as low as 100
µA, where this location is already achieving a clear and sharp PS, but not greater than
approximately 250 µA when the efficiency drops off considerably.

Stimulating amplitudes in this paper are sometimes greater than those typically used
in clinical applications, especially for very large and myelinated fibers. Activation
thresholds could have been much greater that those presented given that very small
unmyelinated axons are, theoretically, much less excitable [33, 35]. It should be noted that
this study adds to evidence that near-ideal conditions allow large evoked potentials to be
elicited with 1ms pulses lesser than 300μA. However, not all stimulating locations studied
in this paper had such low thresholds of excitation and these response curves did not begin
to saturate until much higher amplitude pulses were used (i.e. cell body layer of the
infrapyramidal blade).

2.5.5 Model Limitations and Future Developments
First, this study presents predictions of evoked potentials using a different volume
conductor method than that used to estimate the field resulting from stimulation. While
each method is valid and appropriate in application, a more compelling implementation

42

would be to use the same framework for both stimulus and response field estimation. In
future work, the links between bulk tissue and neuronal hierarchical scales will be
strengthened by allowing neuronal feedback potentials to flow between neighboring cells
on a time-stepping basis. In addition to increasing methodological coherence, this feature
would allow a thorough exploration of ephaptic coupling and other activity-dependent
extracellular dynamics in various parts of the hippocampus [59, 60].

As a preliminary work, the extent of the dentate system was limited to a 400µm
section of tissue due to the complexity of the hybrid model. The end goal of the work will
be to expand this section to encompass the entire 10-12mm extent of the dentate gyrus and
incorporate other hippocampal layers. At this stage the effect of anatomic intricacies on the
intracellular propagation of electric fields can be fully appreciated.

Finally, there are neural components present in the DG that are yet to be
implemented in the model. Chief among these are dentate interneurons which would also be
activated by electrical stimuli and interact with the granule cells to change the population
dynamics. While feedback activity isn’t likely to occur fast enough to strongly impact the
initial population response of resting tissue to a single 1 ms impulse in the perforant path,
they may play a larger role in shaping responses to stimulation nearer the granule cell
bodies where these cells synapse. Implementing these network components are a high-
priority going forward so responses to sequences of stimuli can be studied.

2.6 Conclusion
This study shows that the model performs the minimum functions necessary to
improve electrical stimulation systems. However, future iterations of the model and
methodology will be improved in several areas. While the list of limitations and proposed
improvements is not exhaustive, incorporating these features will provide a powerful
framework which can aid the design of devices capable of improved interfacing with this

43

region of the brain. As the detailed neuronal componentry of the model is further expanded
to include more cell types and layers of tissue it will provide a deeper understanding of both
the network properties of the hippocampal circuit and better strategies of inducing and
recording activity via arrays of electrodes. The implications of this work are far reaching as
the community seeks better methods of developing and testing early iterations of prostheses
for the treatment of diverse neurological disorders

2.7 Supporting information
The following materials have been provided as a supplement to the primary
manuscript. These materials include: (i) a detailed description of the approach to modeling
the tissue-electrode interface, (ii) raster plots of the NEURON response to changing
electrode depth and cell body layer resistivity near the crest of the dentate, (iii) complete
tabulated results of stimulus location and pulse amplitude sensitivity analysis, and (iv) the
full evaluation of these results via the multi-objective optimization function described in
2.L. These and core elements of the model have also been made available on ModelDB for
reference. Should additional information prove essential to the understanding of a reader,
please contact the corresponding author.
Fig. S1. The two modes of charge injection, faradaic and capacitive (A), are
accurately modeled by parallel resistors and capacitors (B) [146]. Faradaic
impedance and double layer capacitance at the tissue interface of physical electrodes
with complex geometry can be spatially discretized (C).
A
C B

44

Fig. S2. Shows the NEURON model response to varying resistivity under otherwise
controlled stimulating conditions for two electrodes placed at the crest. (case #2, top
row) When stimulating in the perforant path of the crest, changes in cell body layer
resistivity have a very limited impact on the features of PS. (case #5, bottom row)
When stimulating in the cell body layer of the crest, changes in resistivity can be
seen to have a very large impact not only on PS features but also on the number of
directly activated cells. (Control = 2.28:1).

45

Fig. S3. 117 simulations of nine electrodes at amplitudes 50 µA to 650 µA in
increments of 50 µA. (A) Nonlinear increases in GC activity correspond to axons
recruitment and consequent synaptic activity. (B) Nonlinear recruitment of axons as
stimulation amplitude increases is due to laminar topography. The medial/lateral
performant path (MPP, LPP) of the infra-pyramidal blade and crest show a
sigmoidal trend in GC response to increasing stimulation. (C) The proportion of
GCs active at the peak of PS is location dependent with higher amplitudes at the
crest due to greater concurrency of transverse propagation (D). Due to the greater
concurrency of propagating activity, crest stimulation generates a shorter half-height
width. (E) High concurrency of propagation along the transverse extent when
stimulating at the crest results in greater amplitude PS per unit stimulus. Maxima
indicate optimal stimulation amplitudes. (F) Electrode proximity to dense regions of
axons results in activation threshold depression. Electrodes placed between the
lateral and medial perforant path in the molecular layer of the tissue have greater
activity for all amplitudes relative to the cell body layer (CBL).

46

Fig. S4. Output of the multi-objective optimization function for all stimulation
cases. The population spike amplitude and power efficiency were maximized and
the half-height width was minimized with equal weighting in this optimization.
(CBL-cell body layer, LPP/MPP-lateral/medial perforant path).

47

3 GRAPH-BASED MODELS OF CORTICAL AXONS
3.1 Abstract
Advances in computation and neuronal modeling have enabled the study of entire
neural tissue systems with an impressive degree of biological realism. These efforts have
focused largely on modeling dendrites and somas while largely neglecting axons. The need
for biologically realistic explicit axonal models is particularly clear for applications
involving extracellular electrical stimulation because axons are generally more excitable
than other neuroanatomical subunits. While many modeling efforts can rely on existing
repositories of reconstructed dendritic/somatic morphologies to study real cells or to
estimate parameters for a generative model, such datasets for axons are scarce and
incomplete. Those that do exist may still be insufficient to build accurate models because
the increased geometric variability of axons demands a proportional increase in data. To
address this need, a Ruled-Optimum Ordered Tree System (ROOTS) was developed that
extends the capability of neuronal morphology generative methods to include highly
branched cortical axon terminal arbors. Further, this study presents and explores a clear
use-case for such models in the prediction of cortical tissue response to externally applied
electric fields. The results presented herein comprise (i) a quantitative and qualitative
analysis of the generative algorithm proposed, (ii) a comparison of generated fibers with
those observed in histological studies, (iii) a study of the requisite spatial and

48

morphological complexity of axonal arbors for accurate prediction of neuronal response to
extracellular electrical stimulation, and (iv) an extracellular electrical stimulation strength-
duration analysis to explore probable thresholds of excitation of the dentate perforant path
under controlled conditions. ROOTS demonstrates a superior ability to capture biological
realism in model fibers, allowing improved accuracy in predicting the impact that
microscale structures and branching patterns have on spatiotemporal patterns of activity in
the presence of extracellular electric fields.

3.2 Introduction
In the study of extracellular electrical stimulation of neural systems, spatial and
temporal patterns of activity are strongly influenced by tissue geometry. One established
approach to studying this relationship is through morphologically detailed equivalent
circuit models of neurons, including axons. While these models are invaluable for many
different applications, they are especially useful for prediction of tissue response to
extracellular stimulation, where explicit morphologies aid the prediction of activation
thresholds under varying biological and stimulating conditions [35, 61, 62]. For
biologically realistic network models, a common method involves arranging individual
neuron models in virtual space to reconstruct elements of the tissue system being studied
[10, 22, 25, 24, 12]. This approach enables accurate prediction of membrane potentials in
response to changes in electric field geometry and gradient [56].

49

Despite an understanding that geometry and topology influence activity, much of the
biological realism in these studies is reserved for dendritic rather than axonal arbors. This
lack of realism in axonal fibers becomes especially disconcerting when considering both
that (i) central nervous system axon terminal arbors are often highly branched and tortuous
relative to the mostly straight and long nerves of the periphery, and (ii) under most typical
stimulating conditions, axons have shorter chronaxies than somas and dendrites [63, 64,
62]. It follows that suprathreshold stimulation events result in coupled local (driven by the
injected electric field) and distal (synaptically-driven) activity, with substantial realism
being necessary to predict the spatiotemporal pattern of the resulting response in totality.
Deliberate arrangement of neuronal structures is also useful for accurate model-based
prediction of tissue-tissue interactions due to electric fields arising from endogenous
current sources [65]. Accurate estimation of local field potentials (LFPs) and, therefore,
predictions of the region-specific impact of ephaptic coupling are sensitive to the degree of
biological realism implemented in a model system [10]. Lastly, network models that lack
explicit axonal structures may have unrealistic conduction delays between connected
Fig 2.1. Adding biological realism to axon models for the study of extracellular
electrical stimulation allows more accurate analysis of evoked neural network
activity.

50

populations of neurons, leading to potential prediction errors [66]. While delays may be
trivially added to network connections once they are known, biologically appropriate
behaviors must first be estimated. Therefore, geometrical and anatomical realism may also
be necessary to study emergent network activity such as co-oscillatory activity in
hippocampal networks [67, 68].

Despite the long record of hippocampal observation, axonal morphology is not as
scrupulously described as the somas and dendrites of many cell types. With a few
exceptions, studies yielding explicit reconstructions through morphometric analysis of
neuronal branching have focused on dendritic arbors and overlooked their axons [69, 70,
71]. The dearth of robust datasets is exacerbated by the general observation from staining
experiments performed in the hippocampus that axonal structures, even of the same cell
type, may be less stereotyped than dendritic arbors [72]. Perforant path axon terminal
arbors from layer 2/3 entorhinal cortical spiny stellate cells are not constrained to simple
conical, fanned, or star-shaped volumes like so many dendritic arbors [50]. Geometric and
topological heterogeneity make the prospect of using explicit reconstructions unfeasible for
direct use in computational models which require in situ cell density at tissue scale. This
becomes particularly apparent when considering that the unique geometry of the dentate
gyrus, which changes from septal to temporal poles, requires dramatic inter-region variety
in the volume, orientation, and contour of afferent entorhinal cortical axons.

Despite the general absence of morphometrics for axons, the shape of terminal
fields, distribution of synaptic spines, and a rough measure of anatomic domains of axon
terminal fields from either histology or other imaging methods provides data from which
minimally functional axons can be grown. Explicit dendritic reconstructions, spine
counting, and anterograde and retrograde staining experiments provide information
regarding synaptic targets. Distributions of synaptic targets combined with a knowledge of
the general path and origin of a fiber is sufficient to generate a functional structure that
captures the tertiary conformation and local divergence of the axon terminal arbor. When
representing a neural process as a graph, nodes placed at synapses and edges to connect

51

between them and the soma can effectively reconstruct a functional dendritic arbor [73].
Likewise, presynaptic boutons provide nodes that can be connected to each other and the
soma with a deliberate arrangement of edges to form a spatial network or graph. In other
words, a minimally functional axon connects a parent cell to synaptic targets.

Beyond minimally functional morphologies, sophisticated models of realistic
neuronal branching have been proposed with primary application to either generating
unique and artificial dendritic trees or reconstructing them from a series of images. There
are two chief types of generative models: (i) stochastic [74], or (ii) greedy [75, 76, 77, 78].
Stochastic models operate by sampling distributions of branching statistics extracted from
experimentally measured neurons, while greedy graph-based models are much more
frequently used to reconstruct 3-dimensional dendritic trees from stacks of manually
labeled images [73]. Although each of these approaches have been useful under certain
conditions, neither has been properly adapted for use in the generation of virtual
biologically realistic axons.

Consequently, we present in this paper a new graph-based algorithm for generating
biologically realistic tree representations of axon terminal arbors. The proposed model is
inspected for its utility in studying extracellular electrical stimulation of cortical tissue
through analysis of the impact of arbor topography and morphometry on activation
thresholds and, by extension, spatiotemporal patterns of activity in the hippocampus. The
results of this work comprise (i) a quantitative analysis of the generative algorithm
proposed, (ii) presentation and quantitative/qualitative description of generated fibers, (iii)
comparison to leading alternative methods, (iv) demonstration of a method to reduce spatial
complexity of axonal arbors while maintaining accurate prediction of neuronal response to
extracellular electrical stimulation, and (v) an extracellular electrical stimulation strength-
duration study. The value of these studies is twofold: (i) establishing the novelty and utility
that this modeling system yields, and (ii) determining if stimulation-response recruitment
order (i.e. large before small) for straight, long, large-diameter, and myelinated peripheral
fibers is similarly true of small, highly-branched, and unmyelinated cortical fibers.

52

.
3.3 Methods

Feature Reference
Strictly laminar dentate perforant path [72, 79]
En passant; mostly non-terminal boutons [79, 80]
Distribution of bifurcation angles ≈80.3±35.7° [77]
0.1 µm fiber diameter, ≈0.7 µm boutons [50]
Primary bifurcation at/near crest, envelopes entire
transverse of dentate, continues to CA3
[81, 72, 50, 79]
≈17,700 synapses per EC axon [71, 69, 70]
Myelination—mixed, though clearest images show no
myelination below the crest of dentate
[50]

The study presented here focuses on accurately capturing the emergent spatial
features of spiny stellate entorhinal cortical axons within the dentate gyrus in Sprague-
Dawley rats. The model is designed to be flexible to the inclusion of novel morphometric
criteria as new experimental data becomes available but, at present, it is clearly important
that (i) fibers are constrained to 1-1.5mm within the septotemporal axis [50], (ii) laminar
organization along the transverse axis is inviolate, with medial and lateral entorhinal
cortical axons confined to the middle and outer thirds of the dentate molecular layer,
respectively, and (iii) axons synapse with a pronounced en passant connective schema,
where most pre-synaptic boutons are non-terminal. Many more features and their sources
Table 1. Principle features (left col.) of entorhinal cortical axons found in the dentate
gyrus perforant path and the studies which reported them (right col.). These features,
when combined with the topography of spines in the outer and middle thirds of the
dentate gyrus provide guiding parameters for construction of artificial tree models.

53

are detailed in table 1. A greedy, graph-based, or ruled-optimum ordered tree system
(ROOTS) was developed to control these features: nodes and edges were connected to
minimize wiring costs while satisfying user-defined criteria comprising spatial and
morphometric constraints.

3.3.1 Volume and Nodal Constraints
In the construction of a spatial ordered tree, the number and topography of target
nodes strongly determines the emergent features of the resulting graph. Therefore, the
selection of nodes is an important step in the successful generation of biologically
appropriate axonal trees. This process was executed using the Kjoenigsen rat hippocampal
atlas slice at -3.34 from Bregma to segment model boundaries [37]. Past efforts have
elucidated the approximate number and spatial distribution of synaptic targets within the
dentate perforant path [82, 11, 69]. The approximate number and density of granule cells
found within the dentate region of a 1.5 mm extruded slice was calculated based on density
measurements reported in the literature [83, 84]. Spine counts and the laminar topology of
entorhinal-dentate connections were used to create a pool of synaptic targets and the
number of perforant path arbors from which afferent connections might be formed [71, 69,
70]. An abbreviated table derived from Hendrickson et. al. can be found below to
summarize these parameters as they were used in this study:
Table 2. Entorhinal Cortical-Dentate Gyrus Topological Parameters
Granule Cell # spines: middle 1/3: 1050-1200
Granule Cell # spines: outer 1/3: 1100-1300
To summarize: the number of synaptic spines in the outer and middle third of
hippocampal granule cell dendritic arbors, the number/density of granule cells, and the
number/density of entorhinal cortical cells contributing to the perforant path provide the
Table 2. Parameters describing the Entorhinal Cortical-Dentate Gyrus topology used to
design the topography of synaptic targets for axon fiber growing described in later
sections. Medial Entorhinal cortical fibers synapse in the middle 1/3 and Lateral Entorhinal
cortical fibers synapse in the outer third of the Granule Cell dendritic arbors.

54

necessary arithmetic for deducing the number of synapses made within each axon terminal
field. When combined with observational data regarding the septotemporal range of these
axons, this information provided an approximate volume throughout which nodes could be
distributed and a plausible synthetic terminal arbor could be grown. These same parameters
were used to construct much larger, more complex, and previously validated mechanistic
models of a rat dentate hippocampus, therefore additional details can be found in Bingham
et. al. and Hendrickson et al. [10, 82].

3.3.2 Constraining Patterns of Axon Branching
While many of the possible branching patterns of a functional arbor are
constrained by the topography of synaptic targets, there remain as many as 𝑛 (𝑛 −2)
trees that
span a set of targets (n), few of which are biologically plausible [85]. It follows that
encouraging biological realism requires additional non-trivial steps to constrain branching
Fig 2.2. Features of granule cell arbors and afferent connectivity provides guiding
information for both the sampling of target nodes for generation of axon
morphologies and the simplification of generated morphologies [69].

55

features of generated topologies. Fig 3 is an algorithm flow diagram and pseudocode
(provided in S1 Fig.) for proper preprocessing and successful execution of ROOTS. This
process generated the trees analyzed in later sections of this manuscript. The general
proposed approach is to first identify a collection of potential synaptic targets and then
consider each in turn and whether the fiber could make a connection to that target while
satisfying various data-derived rules for branching and bifurcating. These model inputs
dually exert control over the emergent spatial/geometric features of the entire terminal field
and the branching patterns that develop as the fiber is constructed.

While the core of this algorithm is capable of growing axons with highly particular
geometries, it is limited in the cases where axons execute acute turns without forming
connections that cross the resulting sulcus. An additional rule can be used to mitigate this
case-specific flaw: a relative location sensitive dynamic source point and reference angle
(Fig 3-B). This approach requires additional preprocessing—replacing the initial sorting of
points with respect to a single source point with a more sophisticated sort and a dictionary
of relative source-points and reference directions. This new process is performed by (i)
finding spatial clusters of synaptic targets using K-Means, (ii) fitting a mesh to the spatial
clusters using Delaunay triangulation, and (iii) discovering the most likely path to any
Fig 2.3. For greater efficiency, the generative algorithm was designed to operate in
one of two modes: (A) with or (B) without dynamic source-point updating criteria.
This feature enables generated morphologies to successfully conform to acute turns
of varying radius. Disabling this feature when not needed improves efficiency.

56

cluster center from the origin source-point through the constructed mesh [86, 87]. The most
likely path, found using Dijkstra’s shortest path algorithm, is then used as a lookup table
where the edges leading to the cluster within which any new point being considered may be
found are used to constrain angles of branch extension and bifurcation [88]. These three
new steps results in a path directed sorting of targets. It should be noted that K-means is the
least critical new component and is only used to regularize the triangulation and reduce the
complexity for the steps that follow (between 2500-3000 clusters were used in this
implementation to have the desired effect). This new rule allows successful execution of
acute turns, where the fiber bends backward toward the source point without crossing the
resulting cleft or forming any cycles. Critically, these acute turns are executed without
relaxing constraining parameters for branching and bifurcating.

While there are many clustering algorithms and many sophisticated meshing
algorithms, K-Means and Delaunay algorithms were selected because of their speed and
reliability. Because dynamic sourcing (Fig. 3b vs. 3a) adds considerable computational
burden, it is valuable to select preprocessing methods which do not add to this burden
needlessly. Figures S2-S4 show development of this stage of the algorithm and it’s three
constituents: K-Means clustering, Delaunay triangulation, and Dijkstra’s shortest path
methods.

3.3.3 Synaptic Boutons
In recognition that axons have complex surfaces and are not just a series of smooth
and simple pipes, additional functions were written to allow complexities such as synaptic
boutons to be added to an already grown fiber. These boutons, as they are seen in the
dentate perforant path at non-terminal presynaptic densities, are described by Tamamaki as
“periodic varicosities” [50]. Each bouton is ≈5 µm long, ≈0.7 µm in diameter, with varying
distances between depending on the topography of the synaptic targets. Because neither the
exact distribution of inter-bouton distances nor multisynapse formation behavior is known
in this tissue system, it was assumed that they were uniformly distributed throughout the

57

terminal axon arbors every 25 µm, with each bouton being 5 µm in length and 0.7 µm in
diameter. Because these boutons are non-terminal they are likely to be actively conducting.
Lacking experimental evidence to frustrate this assumption, sodium, calcium, and
potassium channel densities and conductances were implemented with the same parameters
in bouton compartments as in non-bouton axonal regions.

Following development and testing of this algorithm, fibers were exported for
simulation in environments such as NEURON [89]. Later in this paper, extracellular
electrical stimulation studies performed with the generated morphologies are presented to
demonstrate the maturity of this analysis pipeline.

3.3.4 Arbor Simplification and Computational Complexity
Simulation of detailed neuronal models is computationally expensive. While this
study was enabled by non-competitive access to a 4,040 processors computing cluster, there
remain concerns about impractical and unnecessary model complexities. With as many as
17,700 possible synaptic connections made by each entorhinal cortical arbor, it became
clear that an approach to morphology simplification would be necessary to reduce the
computational burden of both generating and simulating arbors in NEURON. While the
exact number is not known, because of the en passant nature of the perforant path fibers it
is likely that a passing fiber synapses more than once with any target granule cell. This
provides the opportunity to generate trees using fewer nodes due to the relative co-locality
of synaptic connections between an arbor and any given granule cell. To examine this
assumption, we generated arbors with 8,850 (two synapses per target cell) or 5,900 (three
synapses per target cell) target cells rather than 17,700 target cells and used a single node
from each to guide arbor growth. Perforant path fibers were generated using each of these
node counts and then each was also line simplified. The Ramer-Douglas-Peucker (RDP)
algorithm for line simplification was used on each of these trees to determine the minimum

58

node count required to approximate full complexity fibers [90]. Using a fraction of the
typical dentate granule cell dendritic arbor height and width (e.g. 20% or ~45μm) to set a
maximum RDP-epsilon ensured that path deviations would be much less likely to make
otherwise probable EC-DG connections anatomically impossible.

3.3.5 Alternative Generative Models
To test the functionality of this algorithm in comparison with others that already
exist, the authors attempted to create satisfactory arbors with two commonly used
alternative tools: the TREES Matlab toolbox [78], and L-NEURON [51, 91].

The TREES toolbox, like ROOTS, depends upon carefully selected nodes or
points to grow a graph; therefore, spanning trees generated by this method can use the same
set of synaptic or cellular targets as those utilized by the proposed system. Unlike ROOTS,
however, the only parameter that can be adjusted to improve the performance of the
resulting morphology is the balancing factor (BF). This BF represents the weighting of
priority for path length versus conductive delay in the spanning tree that is generated. To
study this approach, a topology grown by TREES with 2,097 nodes (this number is justified
later in the manuscript). The BF was calibrated by minimizing a multi-objective function
(MOF), formulated as an unweighted sum of independently normalized mean root square
error for each of the following morphometrics: Euclidean distance/path-length (BF), branch
order, bifurcation angle, and total path length. The BF of the arbor generated by TREES
was set at the value that minimized the difference between the MOF scores of the arbors
generated by each system. To maximize similarity between these two morphologies
through calibration of the TREES BF, direct (Euclidean) vs. path length ratios were
calculated for each carrier node, then a histogram of these data was fit with a kernel density
estimate (KDE via gaussian smoothing). This process was repeated for path length,
branching order, branch length, and branching angle. The KDEs for direct vs. path length

59

and branching angle for the TREES axon were subtracted from those for the ROOTS fiber.
These differences were independently normalized and the RMSE was calculated. Each of
these metrics were summed without any weighting. The TREES BF was then calibrated
through minimizing the summed normalized root mean squared error of these differences
(U) according to Eq. 1.

𝑈 = 𝑑𝐵𝐹 + 𝑑𝐵𝑂 + 𝑑𝑃𝐿 + 𝑑𝐵𝐿 + 𝑑𝐵𝐴 (1)

Where each term represents the normalized summed difference of balancing factor
(dBF), branching order (dBO), total path length (dPL), branch length (dBL), and
bifurcation angle (dBA). Unlike TREES and the proposed algorithm, L-NEURON doesn’t
rely on a preselection of target nodes to construct a topology, but rather depends upon
measures of branching structure of a tree, or morphometrics. Fitted distributions to these
measurements are then stochastically sampled to grow a tree. This approach is intended to
capture branching patterns without much regard for the emergent spatial features of a tree.
Because of the lack of an extensive database of EC axons from which to take branching
measurements, we resorted to using the companion tool L-MEASURE to extract
morphometrics from one of our own generated arbors to gauge the feasibility of using
stochastic methods to generate axonal arbors when such a dataset does become available
[92]. Extracted morphometrics were then fed to L-NEURON to generate a
morphometrically equivalent arbor.

Virtual arbors were ultimately evaluated based on their ability to capture known
spatial features of entorhinal cortical (EC) axon terminal fields, including a complex
geometry which conforms to the contours of the dentate gyrus.

60

3.3.6 Strength-Duration Relationship in Response to
Extracellular Electrical Stimuli
To understand how fiber geometry in the hippocampus gives rise to patterns of
activity, fibers with varying patterns of diameter were simulated in response to a range of
current source-arbor distances and stimulus amplitudes. Images of spiny stellate fibers from
Tamamaki & Nojyo and more evidence from Bindocci et al. show continuous fibers with
‘periodic varicosities’, or non-terminal synaptic boutons on the en passant fibers [80, 50].
Because fiber diameter has been shown to influence both conduction velocity and
excitability, it was important to explore the response characteristics of fibers with this sub-
micron variation in diameter [56]. A fiber was generated by ROOTS and instantiated with
one of three patterns of diameter (Fig. 10 A, B, and C) in the simulation environment,
NEURON 7.6.2, so that its behavior could be simulated [89]. Hodgkin-Huxley membrane
biophysics under in vivo temperature conditions were inserted in all compartments and d-
lambda rules were used to determine appropriate space constants for compartmentalization
of the fibers [89]. All other biophysical features were borrowed from nodal biophysics
described in [64]. All morphological features for these fibers, other than the deliberate
variations in diameter, remained as presented in Table 1.

Monopolar point-source stimuli were used in all stimulation experiments presented
in this article. Electrodes were placed in one of two locations near the primary bifurcation
of the perforant path (demarked in Fig. 10) at distances of 100 µm and 500 µm from the
nearest neuronal compartment.

Following setup, two sets of analysis were performed. The first comprised a
strength-duration study of complex arbors. Square-wave pulses with widths between 0.025-
1.4 ms were delivered at each of the two distances. The extracellular voltage throughout the
model space was estimated using an analog of Coulomb’s law with material resistivity of
3.8 Ω-m (Eq 2 and 3) [55, 56, 10].

61

𝜙 (𝑥 , 𝑦 , 𝑧 ) =
𝐼 0
4𝜋 ∗𝜎 𝑖 ∗𝑟 𝑖 (2)

Where Φ is the field potential resulting from a current source, I. The conductance,
σ, is the inverse of resistivity. Radial distance, r, is found by Eq 3:

𝑟 𝑖 = √(𝑥 𝑖 − 𝑥 0
)
2
+ (𝑦 𝑖 − 𝑦 0
)
2
+ (𝑧 𝑖 − 𝑧 0
)
2
(3)

By adding new voltage sources in series with the circuital elements representing
each section of membrane, extracellular potentials calculated via Eqs. 2 and 3 were applied
to neuronal compartments within the NEURON model using the ‘extracellular’ mechanism.
Stimulation was delivered over a range of current amplitudes designed to cover, at a
minimum, rheobase to twice rheobase for each set of stimulating conditions.

The second analysis examined the same models and stimulating setup but focused
on the temporal distribution of the response at rheobase. The time it took each compartment
in the fiber to reach action potentiation was recorded, used to populate a histogram, and
then gaussian smoothed to present a kernel density estimate (KDE). These plots were used
to identify the impact of boutons on conduction velocity throughout the complex arbor.

3.3.7 Data/Model Sharing
It is important to the authors that the method by which functional cortical axons
were generated be highly accessible to other members of the computational neuroscience
community, especially those studying neuromodulation for the treatment of diverse
neurological disorders. Therefore, the graph-based algorithm used in this study has been

62

compiled in a user-friendly manner and distributed to the Python Package Index under the
name Neural Roots (nRoots) along with concise documentation. The model is also being
prepared for submission to ModelDB where it will be available alongside other elements
presented in this paper.

.
3.4 Results
Results comprise three parts: (i) a qualitative and quantitative assessment of the
proposed algorithm, (ii) a comparative analysis with alternative methods, and (iii) the
presentation of strength-duration curves and an analysis of entorhinal cortical fiber
recruitment order and temporal distribution of the response to extracellular stimulation.

Fig 2.4. Each line represents the time-course of algorithm execution to generate
an axonal graph with dentate perforant path topography using different numbers
of carrier-points, ranging from 300-2000. Each iteration of the algorithm includes
one loop for branch extension and another for bifurcation.

63

3.4.1 Algorithm Assessment and Validation
There are two chief loops within the core ROOTS algorithm: branch extension and
bifurcation. As the algorithm proceeds, the complexity of branch extension decreases while
that of bifurcation increases, this can be seen in Fig 4 as a linear/slightly supra-linear shape
of each individual line. Bifurcation criteria become more difficult to satisfy as the axon
graph becomes more complete, resulting in exponential increases in time-to-completion
when larger and larger axons (more nodes to connect) are generated. Many aspects of these
trends are dependent upon the spatial topography of the nodes themselves. For example, if
nodes fall along a straight enough line, a single branch extension loop with no bifurcations
will complete the axon graph and the time-to-completion will be perfectly linear with
respect to the number of nodes.

Fig 2.5. Entorhinal cortical axons may form as many as 17,700 synapses with
granule cells. However, growing axons with this number of nodes is slow and
simulating axons with this complexity is computationally prohibitive when
attempting to simulate in situ scale/density tissue models.

64

The algorithm described in Fig. 3b yielded fibers shown in Fig. 7. These fibers
captured many of the known features of layer 2/3 entorhinal cortical spiny stellate axons
which make up the dentate perforant path. These features include but not limited to:
distribution of bifurcation angles of approximately 80°±34°; a septal-temporal range of
between 1-1.5mm; laminar organization with the MEC and LEC in the outer and middle
thirds of the dentate molecular layers, respectively; saturation of branching order at a
reasonable level to encourage en passant synapsing; presentation of both DG and CA3
terminal fields with the preservation of a fissure between the two fields; and finally, the tree
structure passed through a plausible topography of DG and CA3 arbor domains which
enables in situ levels of connectivity.

Fig 2.6. Ramer-Douglas-Peucker (RDP) line simplification was studied over a range
of epsilon for reduction in node count and change in total path length. Simplification
was performed for the cases where two or three synapses were made with each
targeted cell (~8,850 or 5,900 starting nodes, left and right plots). Above 95% path
length is preserved by setting epsilon to 38 and 42μm, resulting in reduction of
required number of nodes to 2,944 and 2,097, respectively (~60%).

65

3.4.2 Arbor Simplification and Computational Complexity
Consideration of the complexity of the resulting fiber is essential because
simulation of just a few of these fibers at full complexity could be highly taxing on a
single-processor computer. Fig 5 demonstrates the relationship between simulation
efficiency and fiber complexity using the NEURON engine. Simulation efficiency is
reported as the ratio of clock-time to simulated time. As fiber complexity increases, the
amount of processing time (clock-time) required to complete an otherwise controlled
simulation also increases linearly. Data for Fig 5 were collected on an Acer Aspire TC-885-
UR17 desktop computer.

RDP line simplification was performed on axons with two or three synapses per
target granule cell (8,850 or 5,900 nodes). Fig 6 demonstrates that axons with as few as
2,944 and 2,097 target nodes can be used to approximate the behavior of an arbor with
maximum complexity, 8,850 and 5,900 target granule cells, respectively. These
simplifications can be made while allowing, at most, 38-42 μm deviations from original
contours and less than 5% reduction in total path length.

66

Fig 2.7. Example graph models (red, LEC; teal, MEC) of entorhinal cortical axon
terminal fields, generated by the proposed algorithm. These axons have features
expected in in situ fibers, namely: en passage terminal topography, laminar
architecture (i.e. MEC/LEC are spatially segregated within the dentate molecular
layer), and fibers proceed to the ends of each supra/infrapyramidal blade of the
dentate gyrus and proceed to a terminal field in CA3/2/1.

67

Example graph models of entorhinal cortical axon terminal fields which
incorporate these simplifications, generated by the proposed algorithm are presented in Fig
7. These axons maintain the features (Table 1) expected of in situ fibers, namely: en
passant terminal topography, laminar architecture (i.e. medial (MEC) and lateral (LEC) are
spatially segregated within the dentate molecular layer), and fibers proceed to the ends of
each supra/infrapyramidal blade of the dentate gyrus, continuing into the CA3 stratum-
moleculare.
Fig 2.8. (Top) Dorsal aspect of the enclosed blade (suprapyramidal) portion of an
entorhinal cortical axon generated by one of two methods: Trees Matlab toolbox
(bf=5), and the proposed algorithm. Balancing factor, and connection threshold for
the Trees toolbox were selected to align these distributions aligned favorably.
(Bottom) Fundamental statistical distributions comparing the two methods in terms
of nodal path length, direct distance vs. path distance, branching order, branch
lengths, and bifurcation angles. Despite matching many features of the proposed
algorithm, the Trees generated arbor does not capture the en passant nature of the
perforant path (seen in branching order (A), and length (B)) and has lingering
bifurcations which have extreme angles (C).

68

3.4.3 Alternative Generative Models
The authors took a morphology generated via our own graph-based algorithm and
extracted morphometrics via L-Measure and then, using the stochastic system called L-
Neuron, attempted to regenerate a morphometric equivalent. This exercise failed to return a
morphology which could conform to the topography of the molecular layer of the dentate
gyrus. The chief reason being that purely stochastic methods are unable to prevent
excursion of fibers beyond natural boundaries in the volume of tissue system being
modeled. The general stochastic method is, therefore, unsuitable for generation of
anatomically appropriate axonal morphologies without further modification to constrain
emergent spatial characteristics.

Fig 2.9. Plotting path length vs. branch order makes it clear that ROOTS more
accurately captures the en passant nature of the perforant path than does TREES.
99% of total axon path length of entorhinal cortical spiny stellate cells should be
achieved with fiber of branch order no greater than ≈7.

69

Attempts were also made to generate accurate morphologies using the TREES
MatLab toolbox. These efforts yielded trees superior to those generated by the L-Neuron
method, though still deficient in important ways. An axon generated by TREES was
selected for quantitative and qualitative comparison with that yielded from earlier analysis
using the ROOTS system (Arbor Simplification and Computational Complexity). Each tree
was constructed using the same set of target nodes, therefore all differences in patterns of
branching arise from differences between the two algorithms. The TREES balancing factor
was then calibrated through minimizing a multi-objective function, attempting to find the
best possible match between TREES and ROOTS arbors. This process yielded a functional
axon from the TREES algorithm, though differences remained that could not be resolved
through manipulation of the BF alone (Fig 8). Despite a high degree of similarity in direct
(Euclidean) vs. path length, total path length, and branching angle distributions, TREES
resulted in a higher than expected number of terminal branchlets. This shifted the branching
order and branch length distributions to the right and left, respectively. These shifts
demonstrate difficulty for the TREES approach in appropriately capturing the en passant
connection schema which typifies EC axons of the perforant path, where terminal
branchlets are reportedly rare. Fig 9 (Branch Length vs Branch Order) further highlights
this conformational difference between the two fibers and provides a useful comparison to
experimental measurements.
In this figure, the black line denotes the branch order at which 99% of the total
path length was achieved in layer two spiny stellate axons reported by Budd and colleagues
[75]. While the comparison with Budd’s report is useful, it is made with hesitancy because
Fig 2.10. Rendering of a highly-branched perforant path axon with presynaptic
boutons dispersed evenly throughout the terminal field. This arbor with three
different patterns of diameter were then used in a set of stimulation experiments.

70

of the variability of branching that occurs within the long stretch of fiber between the
entorhinal cortex and the extensively branching terminal arbor. The difficulty TREES has
in controlling branch order arises from the paradigm of control implemented in the system.
More specifically, the balancing of total path length versus conduction time (BF) does not
provide for fine control of the shape and branching properties of an arbor within a volume.

An additional challenge not satisfied by the TREES approach is that of performing
acute turns. A portion of the perforant path innervates the CA3/CA2 region of the
hippocampus and therefore must be able to extend from the end of the enclosed blade of the
dentate into these other domains. The algorithm presented in Fig 3b. (dynamic source)
Fig 2.11. Arbors of three patterns of diameter corresponding to stem and bouton
diameters observed by Tamamaki, and a fiber with 5µm boutons spaced 25µm apart,
were simulated in response to cathodal pulses (25 μs – 1.25 ms) and a range of
stimulus amplitudes from distances of 100 and 500 μm (left, right) [50].

71

describes how updating reference angles allows the fiber to bend around complex
anatomies successfully without violating morphometric constraints (Figs 7 & 10). An
attempt to replicate this feat using the TREES Matlab toolbox is presented in supplemental
Fig 5.

3.4.4 Strength-Duration Relationship in Response to
Extracellular Electrical Stimuli
To determine the importance of including boutons when simulating extracellular
electrical stimulation of axon models, strength-duration curves for three diameter patterns
of a ROOTS arbor, where diameters correspond to bouton and inter-bouton sizes, were
estimated and presented according to the details in Fig 11. When stimulating with anodal
pulses, these results agreed with previous reports in observing that larger diameter fibers
having shorter chronaxies than small, and that this difference is exaggerated by large
electrode-fiber distances and longer stimulation pulses. For cathodal impulses of such small
and highly branched fibers, the model yielded no clear pattern for selective activation of
fibers by diameter. This protocol was repeated for a biologically realistic arbor, complete
with boutons. The recruitment pattern of the boutoned fiber was no different in a cathodal
field but under anodal conditions had activation thresholds between 0.7 µm and 0.1 µm
uniform diameter arbors. The differences between the boutoned fiber and the 0.7 µm fiber
decreased with increasing electrode-fiber distance. These differences dissolved as the
pulse-width approached 700 µs.

72

The impact of boutons on temporal dynamics are presented in Fig 12. The time of
first action potential is plotted along the horizontal axis for each neuronal compartment in
the biologically realistic arbor with the vertical axis representing frequency. A fiber with
each of three patterns of fiber diameters (uniform 0.1, 0.7, or boutoned) were simulated in
response to 1 ms stimuli at rheobase amplitude. Activity was initiated near the electrode
and then actively conducted to other compartments, including those beyond the effective
volume of the extracellular electric field where artifact voltage dropped to tens of
microvolts. In Fig 12 are plotted, as a KDE, the time-course of threshold activity of all
compartments in each case. As the average fiber diameter is increased, the shape of the
KDE shifts to the left, indicating that activity in the arbor is occurring with less delay
following stimulation. This shift is explained by the increased conduction velocity
throughout the arbor which results from changing compartment diameters. The boutoned
fiber nearly mirrored temporal patterns of activity of the uniform 0.7 µm fiber, except at
very close electrode-fiber distances.
Fig. 2.12. Kernel density estimates (1.5 ms gaussian kernel) for compartmental
action potential delays in perforant path arbors with a diameter of 0.7 µm, 0.1 µm,
and a boutoned fiber when stimulated at ~rheobase amplitude for 1 ms. From top to
bottom, left to right, showing two electrode-fiber distances and anodal and cathodal
stimulation. Excepting at very close distances, boutoned fibers had temporal
activation features best approximated by a fiber with uniform 0.7µm diameter.

73

3.5 Discussion
Together with increased computing power, more robust repositories of
electrophysiological and histological information have created opportunities for neural
models of greater complexity. Despite these ever-improving sources of data, biologically
realistic neural modeling has outpaced experimental studies in many areas of inquiry,
providing a testbed for new hypotheses as well as a wealth of preliminary data to aid the
design of superior in-vitro or in-vivo studies. However, the mismatch in biological realism
between experimental and computational modeling, and continuing ambitions to advance
computational neuronal modeling, creates a need for sound approaches to generate
functional neuronal models that can be refined as the quantity and quality of experimental
data improves. The modeling approach presented in this paper demonstrates this utility.

3.5.1 Limitations and Alternative Methods
Despite the demonstrated functionality of the model presented herein, it is not
without limitations. First among potential limitations is biological realism, which may still
be limited due to lack of accurate and well described experimental measurements of fibers
from hippocampal tissue samples. The clearest limitation (though simultaneously the chief
motivation for development) of this modeling approach and the study presented in this
paper is the lack of an extensive dataset describing entorhinal cortical axons in the perforant
path, including rich morphometrics. Existing explicit reconstructions have readily
identifiable errors and, therefore, should not be used as solitary sources of branching
morphometrics or be virtualized and used in stimulation models without sophisticated
automated and/or manual repair [50, 75]. The issues include slicing artifacts which distorts
both distances along or across serial slices and bifurcation angles where branches span
multiple slices; angle and distance mismeasurements due to imaging along a single axis
and, therefore, failing to correct for the impact of rotation of neural processes with respect
to observer perspective; and failures of automated algorithms which may result in orphaned

74

sections or even cycles in the final tree [93]. Despite the unavailability or verifiability of
these data, the modeling methodology presented in this study represents the most detailed
and sophisticated functional model of layer 2/3 entorhinal cortical spiny stellate dentate
perforant path axons to date. In general, the utility of this modeling approach is most
apparent for cases in which explicitly reconstructed morphologies are sparse, poorly
described, and the fibers to be modeled have complex geometry and branching structure.
Further, the parametrization of the model is such that, should this data later become
available, the generative morphologies could be updated to reflect new knowledge of in situ
morphometrics.

While the present study espouses a Ruled-Optimum Ordered Tree System
(ROOTS), there are many algorithms that might achieve topological characteristics
reminiscent of neuronal branching.

Though existing stochastic models may be effective for dendritic arbor generation,
they are not naturally adapted to the construction of axonal morphologies because they fail
to reliably conform to pre-determined geometries or volumes that are irregular or highly
nonsymmetrical [74]. The ability to conform to a predetermined volume or geometry is a
particularly important feature for axons in the hippocampus, which often have lamellar
organization between layers and laminar organization within their terminal fields. For
stochastic models to accomplish this, new pre-processing intensive and highly inefficient
volumetric constraints would be required, adding to algorithmic complexity.

When guided by serial histological sections, spanning-tree algorithms provide an
efficient approach to reconstructing spatial trees. However, without some modification,
traditional minimum-spanning-tree algorithms can result in spatial error and distorted
branching patterns [75]. Using higher resolution images forces significant biological
realism onto the resulting graph, making error more manageable. However, minimum-
spanning-trees generated from points sampled at random from a volume representing an

75

axon terminal field are less biologically realistic due to the extreme path lengths that result.
To moderate this outcome, some investigators have proposed a spanning-tree algorithm
which balances the minimization of membrane against path length [75, 77, 76, 78].
However, optimizing this balance is known to be an np-hard problem due to the direct
conflict introduced by including these two variables in the same loss function [94, 95, 96,
97, 98]. Another challenge to using path-length/membrane-minimization balancing
algorithms is that it only allows indirect control over important tree morphometrics such as
bifurcation angle, branch extension angle, inter-bifurcation length, etc. Prominent
implementations of this approach do not allow explicit morphometric thresholds to be set
and, therefore, extreme branching patterns remain possible. Inflexibility and lack of
sufficient parameterization represent significant limitations of these prior efforts because
branching patterns directly impact spatiotemporal patterns of activity.

It is plausible that established stochastic methods, known to be highly proficient at
recapitulating tree morphometrics, could be modified to allow directedness and improved
conformation to terminal field volumes. Despite these feasible and time-worthy
alternatives, the model presented in this study has demonstrated sufficient performance in
terms of accuracy, flexibility, and computational complexity.

While many alternatives to RDP for tree simplification exist, this method was
preferable to other node cluster-and-merge methods (e.g. Kruskal’s algorithm) because
these alternatives have the potential to introduce large differences in branching patterns,
creating additional challenges in evaluating the equivalence of simplified and unsimplified
fibers. Evaluation of the RDP simplification approach could be done in a very
straightforward fashion because line simplification doesn’t result in changes in branching
patterns--simplified arbors were compared to the reference case on the basis of path length
and subjective evaluation of an acceptable RDP-epsilon, or the maximum distance between
original and simplified contours.

76

Analysis in this paper involving extracellular field estimations used analytical
methods that fail to completely account for anisotropy and heterogeneous resistivity of the
tissue volume or full-wave propagation of electric fields. More sophisticated field
estimation techniques (e.g. Finite Element, Finite Volume, Boundary Element, or
Admittance Methods) that account for complex impedance or by filtering point-source
estimations to account for amplitude dampening and phase-shift of stimuli could have been
used but were deemed unnecessarily complex for the principle questions under examination
in this article [99, 26]. It should be noted, however, that rheobase predictions result from
very long stimulus pulses, where >95% of the frequency domain signal falls below 1 kHz;
it follows that tissue capacitance is only a marginal source of error in these estimates of EC
axon chronaxie. With respect to resistivity, this exercise is further justified because
resistivity measurements in the region performed by Lopez-Aguado et al. [100] show
nearly uniform resistivity throughout the molecular layer of the hippocampus; further,
because the current source in these experiments was placed relatively far from resistive
boundaries, current shunting distortions should be minimized in our estimations of electric
fields in the volume occupied by EC fibers [10]. While future work not limited by these
assumptions will be performed by the authors, doing so here falls outside the scope of the
present study.

The study of strength versus duration of stimulus in the measurement of activation
thresholds as performed in this paper does not provide conclusive data which concretely
establishes the connection between fiber size and excitability for hippocampal networks.
This is due to the realistic elements still missing from the models used, including non-
uniform diameters of fibers, irregular bouton geometry and volumes, varying bouton
topography, and an electrophysiological study of any differences in channel density and
dynamics between boutoned and non-boutoned axon regions. Despite these limitations, it
seems important to consider what impact boutons may have on fiber excitability within the
terminal fields of small, highly-branched, and unmyelinated axon fibers. Importantly, these
results imply that action potentials in extracellular anodally stimulated axon terminal fields
are initiated in boutons with lower input resistance than nearby portions of the fiber and

77

that failing to consider the impact of boutons on conduction velocity will likely result in
meaningful temporal errors. However, thorough experimental work is needed to confirm
this possibility. It is further valuable to recognize that for long pulses delivered at more
distant locations uniformly large diameter fibers approximate the vastly more complex fully
boutoned fibers quite well and may provide a viable approach to reducing the
computational demands of models that include detailed axon terminal arbors.

While much remains to be done, this study represents a step forward for detailed
computational modeling of complex neuronal systems. Where previous models were either
focused on peripheral axons with less complex arbors or used sophisticated methods to
generate dendritic trees but neglected axons altogether, the model presented here
demonstrates an approach to constructing functional axonal morphologies that can be used
for diverse applications, including extracellular electrical stimulation of the cortex. It
should be noted that while ROOTS is itself general, expert knowledge of the tissue system
being virtually recapitulated is required: branching criteria and synaptic or cellular targets
must be provided to the tool. Critical morphometric parameters include bifurcation angles,
branch extension angle, and internode length; these must be independently determined by
users. Despite these requirements, ROOTS presents many opportunities to increase the
sophistication of model-based studies of neural tissue system dynamics and central nervous
system neuromodulating devices such as deep brain stimulation or potential hippocampal
memory prostheses [2, 101]. ROOTS could also be used to support models of optical or
pharmaceutical neuromodulation [102, 103]. The model framework could easily be
extended to incorporate high resolution imaging data or combined with more complex
volume conduction (or diffusion, diffraction, etc.) models to study tissue-electrode
interactions at smaller spatial scales [104]. In addition to supporting prosthesis design,
ROOTS facilitates model-based exploration of the effect of diseases which remodel axons
on proper function of hippocampal tissue (e.g. multiple sclerosis) [105]. More generally,
ROOTS supports efforts to create networks of neuronal models for the study of biologically
plausible spatiotemporal patterns of activity.

78

3.6 Supporting information

S2.1 Fig. Pseudocode of the core algorithm as presented in this study. Other
components such as source point reassignment, and topology simplification
procedures are executed in series with (before or after) the steps outlined here.

79

S2.2 Fig. Clustering is performed via K-Means clustering method. This regularizes
the spatial properties of the mesh which will comprise the volume of the arbor and
reduces computational complexity of subsequent steps.
S2.3 Fig. Delaunay triangulation to construct a surface from spatially clustered
synaptic targets. This yields a network from which ideal paths to potential target
zones can be calculated.

80

S2.4 Fig. Plotting of each ‘most likely path’ between nodes in the Delaunay
triangulated network of k-means cluster centers. Paths are constructed via Dijkstra’s
algorithm. These paths are provided as inputs to the ROOTS algorithm to loosely
guide branching behavior.
S2.5 Fig. Example graph models (red, LEC; teal, MEC) of entorhinal cortical axon
terminal fields, generated by the TREES toolbox method. These axons were
generated with the addition of a terminal field in CA3/2/1. For very deep folds with
small clefts, the TREES MST method has difficulty suppressing trans-cleft
connections that short the cortical circuit without inappropriate manipulation of the
number and density of targets, further exacerbating the proliferation of terminal
branchlets.

81

4 ADMITTANCE METHOD FOR ESTIMATING
LOCAL FIELD POTENTIALS GENERATED IN
A MULTI-SCALE NEURON MODEL OF THE
HIPPOCAMPUS
4.1 Abstract
Significant progress has been made toward model-based prediction of
hippocampal activation in response to extracellular electrical stimulation, but challenges
remain in the accurate and efficient estimation of local field potentials (LFP). Analytical
methods of estimating electric fields are a first-order approximation that may be suitable for
model validation but they are computationally expensive and cannot accurately capture
boundary conditions in heterogeneous tissue. While there exist many appropriate numerical
methods of solving electric fields in neural tissue models, there isn’t an established standard
for mesh geometry nor a well-known rule for handling any mismatch in spatial resolution
or alignment between current sources and mesh nodes in a finite-element or resistor-
network method volume conduction model. Therefore, using a previously published and
validated multi-scale model of the hippocampus, the authors have formulated an algorithm
for bidirectional communication between Admittance Method (AM) volume conduction
models and biologically detailed neural circuit models constructed in NEURON.

82

Development of this algorithm required that we assess meshes of (i) unstructured
tetrahedral and grid-based hexahedral geometries as well as (ii) differing approaches for
managing the spatial misalignment of current sources and mesh nodes. The resulting
algorithm is validated through the comparison of AM predicted evoked potentials with
analytically estimated LFPs. Establishing this method is a critical step toward closed-loop
integration of the AM and NEURON models that could lead to substantial improvement of
the predictive power of multi-scale models of cortical tissue. These models may be used to
deepen our understanding of hippocampal pathologies and the identification of efficacious
electroceutical treatments.

4.2 Introduction
Encouraged by advances in computation, new multi-scale neural network and
electromagnetic models have demonstrated the ability to predict spatiotemporal patterns of
activity in complex neural tissue systems when activated by geometrically and
dielectrically detailed models of stimulating electrodes [11, 106, 107, 108, 109]. Linking
neuronal and tissue volume-conductor scales has previously been performed by estimating
feedforward interactions between extracellular electric fields induced by stimulating
impulses and neural processes which lie within affected volumes. But because behaving
neurons are independent sources of current, it is warranted that numerical methods be
extended to incorporate transmembrane currents into dynamic estimates of extracellular
electric fields generated within complex neural tissue systems.

The Admittance Method (AM) has been established as a valid and intuitive
numerical approach to estimating electric fields induced in animal tissue systems through
extracellular electrical stimulation [29, 31]. In short, this method involves construction and
solution of an equivalent RC circuit used to approximate the attenuation of electric fields in
volume conductors. Numerical approaches like AM (including finite-element method and
finite-volume method) have become favored over analytical methods for two reasons:
numerical methods allow (i) more efficient parallel/simultaneous estimation of multiple

83

locations within an electric field and (ii) superior ability to capture various dielectric
heterogeneities of complex animal tissue systems. For very complex and highly segmented
volumes, both heterogeneous resistivity and anisotropy have been demonstrated to create
important boundary conditions during tissue volume conduction [43, 23, 8, 25, 110, 22,
111, 112, 113]. Additional complexities are presented in neural tissue, which conducts
electricity both passively and actively. In other words, neural tissue systems behave
simultaneously as a volume conductor and as a network of spiking neurons and these two
domains of tissue interact bidirectionally. Active conductance through neural networks,
occurring through changes in ion concentrations across cell membranes, results in
measurable currents in the extracellular space that may entrain or reinforce oscillatory
behavior in neural networks [67, 68, 13]. These measurements, called local field potentials
(LFPs), are critical indicators of living tissue system behavior.

Simulation of neural tissue systems, through large-scale computational modeling,
has become possible through dramatic improvements in computation in the past decade,
namely: parallel computation through either graphic-processor or more traditional CPU
cluster computing. When combined with estimations of electromagnetic fields,
computational models of neurons have demonstrated utility as a method of optimizing
neurostimulating or recording arrays of electrodes [10, 40, 114, 17, 20, 64]. Despite
advancements in these models, there remain significant limitations in numerical or
analytical methodologies for incorporation of endogenous currents into estimations of local
field potentials (LFP) or evoked potentials (EP) for use in analyzing the stimulus-responses
of neurological tissue systems. Barriers to be overcome include: proving the unestablished
theoretical protocols for tetrahedral or hexahedral approximation of a line-source or point-
source within a conductive volume and designing a feasible process for tackling the great
computational burden presented by the task. Such a theoretical solution would ideally allow
currents to give rise to roughly toroidal or spherical electric fields generated by line or
point-sources, respectively, in an infinite homogeneous conductive volume. But it is
unclear at which length scales preserving these hypothetical field geometry matters for the
accuracy of field estimations.

84

This paper seeks to distinguish the Admittance Method as a numerical approach to EM
modeling that is capable of smoothly incorporating neuronal current sources in a meshed
model for prediction of evoked potentials in complex neural tissue systems. The proposed
extension of this method yields the following utilities: (i) a granular algorithm and guidance
for adding dynamic current sources to existing meshed volumes where spatial alignment
between nodes and new current sources is variable, (ii) automated application of very large
numbers of dynamic current sources to existing AM meshed volumes, (iii) intuitive
implementation of tissue volume conduction behaviors via equivalent circuits, and (iv)
simultaneous solution of electric fields at many locations. This last contribution represents
a chief advantage of this method (along with other feasible numerical solutions) over
analytic calculations.

4.3 Methods
The method described here seeks to establish a fundamental approach to adding
dynamic current sources to existing meshed volumes where spatial alignment between
mesh nodes and true current sources is variable. Secondly, this paper seeks to automate
application of very large numbers of dynamic current sources to existing AM meshed
volumes where sources and mesh nodes are presumed to not be co-local. Overcoming these
Fig. 3.1. (A) A 2-D rendering of the biologically realistic NEURON model of a
rodent dentate gyrus with demarcation of recording and stimulating devices. The
slice model is populated with 10k entorhinal cortical axons which synapse with 50k
granule cells (B). The entire model comprises ≈12.25 million synapses and ≈15
million neuronal compartments. Details regarding the construction and application of
this model are described at length in Bingham et al. [106].

85

challenges allows us to take advantage of the superior computational efficiency of
numerical methods of solving volume conduction problems. The algorithms proposed
herein were developed through utilization of data previously published and generated from
a multi-scale model of a rat dentate gyrus. The following sections briefly summarize the
construction, feedforward simulation, and validation of a complex multi-scale AM-
NEURON model of the hippocampus [11, 106]. While this model is extensively described
in previous publications, it benefits the reader for some concise explanations to be provided
here.

4.3.1 Construction of the NEURON Model of the Hippocampus
In previously published work, a model of a rat dentate gyrus slice was
implemented in NEURON [107, 11, 106]. Taking an image from the Kjonigsen and Witter
hippocampal atlas, the cell layer boundaries were extruded 400μm in the septotemporal
direction to yield the entire volume within which the model was constructed [37].

The NEURON model is comprised of 50,000 granule cells and 10,000 entorhinal
cortical axon fibers. Previously validated granule cell models, including morphology and
biophysics, were taken from Hendrickson et al. and arranged along the cell body layer of
the image from the atlas [36, 28]. Then 5,500 lateral perforant path axons were arranged
such that they coursed through the outer third of the granule cell arbors and 4,500 medial
perforant path fibers likewise through the middle third. Each granule cell then formed
synapses with these fibers, prioritizing proximity, totaling ≈5.5 million synapses in the
middle third and ≈6 million in the outer third.

86

4.3.2 AM-NEURON Feedforward Simulation
Volume conduction in neural tissue in this study was modeled using a variant of
the heterogeneous Admittance Method as extensively described in [106, 29, 31]. In brief,
this method involves four steps: (i) image-based sculpting of a model volume and
automated extraction of bounding surfaces between sub-volumes of differing conductivity,
(ii) construction of a resistor (and capacitor) network of varying resolution which contains
Fig. 3.2. This visualization presents the source data used to develop the models
presented in this manuscript. Each ball represents the location of a current source
generated by a neuronal compartment. The color and each ball encodes the charge
polarity and the size denotes the normalized amplitude of the voltage each ball
contributes to the potential as measured at the recording electrode. This
normalization was performed for each time-step independently using the line-source
equation [55]. A feature of note in this data set is that the spatial distribution of
important current sources varies dramatically from beginning to end of the simulated
behavior.

87

or spans the defined volume, (iii) assignment of resistances to individual elements in the
mesh based on experimentally measured resistivity (the estimated real component of
impedance), and (iv) solution by conjugate-gradient descent of a sparse matrix
representation of internode resistances and currents and nodal voltages formulated after
Ohm’s law [115]. This established method provides meaningful flexibility in terms of mesh
construction and dielectric properties of the resulting model when compared to common
alternatives such as FEM [116].

Fig. 3.3. This toy model demonstrates the chief architecture of the equivalent circuits
used to estimate the passive and active propagation of neural network activity in an
AM-NEURON model. Volume conduction is modeled via a specialized resistor and
capacitor network called the Admittance Method (yellow). Intracellular and
transmembrane currents are then modeled using a second set of equivalent circuits
setup and solved within the NEURON simulation environment (blue) [yale].
Differences in spatial resolution between these two domains draws attention to the
impact of potential approaches to handling their interface (i.e. node alignment).

88

RC elements were then added to the AM mesh to represent an array of microwire
stimulating electrodes that could be placed in various locations through the model. Because
both neuronal and volume conductor models are electrical circuits, interfacing the two
intuitively involves passing node voltages or currents back and forth.

Once nodal voltages in the Admittance Method model are solved, voltage at
neuronal compartments are found by interpolating from the 8 nearest nodes and applying
this voltage as an extracellular battery—this is done for each compartment in every neuron
(~200-300 locations for each neuron and about 12.25 million locations in an in situ scale
slice of the dentate gyrus)--this causes currents within each cell that are appropriate given
their location within the field estimated by the Admittance Method model.

With the volume conductor and network models in place, we then added additional
circuit elements to represent virtual stimulating devices. The devices we initially modeled
had real geometry and conductive properties of twisted bipolar platinum microwires
insulated with Teflon. This was implemented with equivalent circuits by discretizing the
tissue-electrode interface with parallel rc components, each representing a mode of current
transmission: faradaic and capacitive charge injection [47, 40].

After this step, we placed the devices in one of three transverse locations and
delivered a biphasic pulse at threshold amplitude. Compartmental currents from this
exercise were then used to calibrate the AM for LFP estimation.

The model was simulated on a 4,040-processor high-performance computing
cluster owned by the authors and supported through the University of Southern California
Center for High-Performance Computing.

89

4.3.3 Analytical LFP Estimates as a Proxy Ground Truth
Previous studies, by these authors and others in the community, used simple
analytical estimations of local field potentials to validate their neuronal models. While the
sources of error associated with such an approach are well understood, the line-source
method of LFP estimate provides reasonably accurate estimates of potentials at a given
observation point in the model. Except in models with significant heterogeneity of
dielectric properties, this estimate can be used as a suitable proxy for ground-truth when
calibrating more sophisticated numerical approaches to estimating electric fields [26, 27,
55].

The analytical method used as a proxy ground-truth in this study is the line-source equation
as follows:
Where ϕ is the field potential, σ is the inverse of tissue resistivity (3.8 Ω-m, obtained from
Lopez-Aguado et. al.), r is the path length from source to recording location [48, 56, 55,
10]. The proxy ground-truth is naively estimated by scaling each current source via Eq.1,
then summing them into a field potential. This process was performed for every time step in
a simulation until the whole time-series is reconstructed. LFP estimates using Eqs. 1 and 2
based on the feedforward AM-NEURON model described in this and previous paragraphs
were validated against LFP recordings from MEA studies of rat hippocampal slices and
published in Bingham et al. and Soussou et al. [10, 11, 49]. The following methods
describe a new analysis of this data.

𝜙 (, 𝑦 , 𝑧 ) =
1
4𝜋 ∑
𝐼 𝑖 𝜎 ∗𝑟 𝑖 𝑛 𝑖 =0
(1)

𝑟 𝑖 = √(𝑥 − 𝑥 𝑖 )
2
+ (𝑦 − 𝑦 𝑖 )
2
+ (𝑧 − 𝑧 𝑖 )
2
(2)

90

Fig. 3.5. Here are two example meshes, one from each principle element geometry.
Each mesh entirely contains the field of current sources for which an LFP is to be
estimated. Hexahedral meshes were created by linear segmentation of the volumes to
ensure voxels were always cubical. Tetrahedral meshes were scripted in python
through the gmsh and tetgen apis using Delaunay triangulation to ensure very high
quality [145] . It should be noted that tetrahedral meshes created using this approach
do not have uniformly equifacial tetrahedra.
Fig. 3.4. Two common mesh geometries, grid-based/structured hexahedral and
unstructured tetrahedral were used to fill 3-dimentional volumes to encapsulate all
current source locations described in Fig. 2. The accuracy and computational
performance of these were analyzed as mesh resolution was gradually increased in
both.

91

4.3.4 New Meshes for LFP Estimation: tetrahedral vs.
hexahedral
Early in this process it was observed that the spatial distribution of currents
contributing to LFP signals differed substantially from the geometry of stimulating currents
(i.e. many disperse small charges vs. a few large and nearly co-local charges) (Fig. 2).
Therefore, we determined to start with a fresh and less constrained set of meshes with
which to study LFP estimation. Without any preconceived notions of which mesh volume
geometry would be best, a bounding box was fit to the point cloud that represented all
compartments in the neuron model which would be contributing currents to any LFP
estimate. Some scalar factor was then applied to this volume to ensure current sources near
the edges of the model would be less affected by boundary-driven current shunting. The
scalar factor used in this study was 30% and was set after determining that the relative field
amplitude of edge-most currents are >95% relaxed by the time they reach the nearest
boundary. This volume was then meshed with either grid-based hexahedral or unstructured
tetrahedral (equifacial ideal) meshes with resolution (maximum edge-length) ranging from
180-40 µm edge-lengths.

4.3.5 Fundamentals of Resistor Network Construction
Tetrahedralizations or hexahedralizations are unbroken networks of nodes and
edges that must then be converted to a valid circuit of virtual resistors (and capacitors for
non-quasistatic models). Rules for deriving resistor parameters from experimental
impedance or resistivity measurements are poorly described or completely absent from the
literature and so we will describe our process here: first for hexahedral and then for
tetrahedral mesh geometry. Because of the orthonormality of structured hexahedral meshes
we assumed edge length, L, and edge cross-sectional area, A, to be the only scalars of
importance in converting local tissue resistivity, ⍴, to resistance, R, according to:
R = ⍴
𝐿 𝐴 (3)

92

A was a constant 0.25 µm
2
for all meshes in this study. Because unstructured meshes
are not orthonormal, it stands to reason that there should be different rules for assigning or
scaling distance derived resistances based on the angle and number of adjoining
connections. Following analysis done by Duffin in the late 1950’s, fractional conductance
within and across tetrahedra can be approximated as proportional to the cotangent of
internal angles [117], after the form:

Where resistance along an edge, 𝑅 𝐴𝐵
, across from the dihedral angle, ϑ, is 6 times
resistivity times the tangent of ϑ scaled by the length of the edge where the dihedral angle is
formed (variables described by Fig. 3.6). The dihedral angle can be reliably calculated with
the following equations:

The denominator in Eq. 5 is the triple product of u, v, and w. This treatment of the
tetrahedral meshing problem later requires that R be summed for edges shared by multiple
conjoined elements to create a unified and continuous network without superfluous co-
parallel resistors. This is done by following traditional logic for summing parallel resistors:

1
𝑅 𝑇𝑜𝑡𝑎𝑙 =
1
𝑅 1
+
1
𝑅 2
+ ⋯ +
1
𝑅 𝑛 (7)
(4)
𝑅 𝐴𝐵
=
6 * ⍴ * tan( ϑ)
|𝐶𝐷 |
(4)
(4)
𝜗 = 𝑐𝑜𝑡 −1
(
√(𝑢 ∙𝑢 )(𝑣 ∙𝑤 )−(𝑢 ∙𝑣 )(𝑢 ∙𝑤 )/√(𝑢 ∙𝑢 )
[𝑢 ,𝑣 ,𝑤 ]
) (5)
(4)
𝑢 = 𝐵 − 𝐴 𝑣 = 𝐶 − 𝐴 𝑤 = 𝐷 − 𝐴 (6)
(4)

93

This overall approach (Eqs. 4-7) has the nice feature of 𝑡𝑎𝑛𝑔𝑒𝑛𝑡 (𝜗 ) going to 0 as ϑ
approaches 0 and infinity as ϑ approaches π/2. It is not clear whether this rule remains
useful for obtuse tetrahedra, though none were generated in these models.

Following construction of a primary mesh and conversion to a circuit, new
resistors were added between the hull nodes and a virtual ground. Each new wire connected
to ground from a hull node had a resistance equal to the average of all other wires
connected to the hull node being considered.

Fig. 3.6. Example tetrahedral element explaining the variables used to apply the
cotangent rule for calculating edge-wise conductance. Resistance from A to B is
calculated by Eq. 3 where 𝑅 𝐴𝐵
is proportional to the ratio of the dihedral angle, ϑ,
and the edge length, CD. A proof of this approach is presented at length in [117].
.

94

Fig. 3.7. AM-NEURON bidirectional interface can be reduced to an intuitive and
methodologically consistent algorithm that ensures field geometry distortion is low
and mostly dependent upon mesh resolution at or near current sources. Eventually,
this algorithm would be executed at each time-step in a simulation: first the
Admittance Method model is solved, then the estimated field is applied to
compartments throughout the NEURON model. Then a single time-step of the
NEURON model is solved and compartmental currents are passed back out to the
Admittance mesh and the whole process begins again to prepare to simulate the next
time-step. Spatial alignment can be accomplished with one of two algorithms:
interpolating/splitting to or from a source or shifting a compartment or source to the
nearest node in the AM mesh. The details of these two approaches are graphically
explained in Fig. 8.

95

4.3.6 Incorporating incident endogenous current sources:
splitting vs. shifting
The fundamental paradigm of the AM was then extended to incorporate incident
and dynamic current sources from NEURON compartments into an already existing
meshed model. For the sake of clarity, a single voxel/single incident current-source
scenario is presented in figure 8. This figure describes two algorithms that will be labeled
descriptively as, current ‘splitting’ and ‘shifting’.

Fig. 3.8. Spatial alignment in feedforward (AM to NEURON) and feedback
(NEURON to AM) interfacing can be accomplished with one of two algorithms:
interpolating/splitting to or from a source, or shifting a compartment or source to the
nearest node in the AM mesh. The accuracy and computational performance of
solving meshes with split vs. shift were analyzed across a range of model
complexities.

96

The steps for current splitting include: (i) relative distance-weighted division of
currents into eight or four parts, for hexahedral and tetrahedral elements, according to each
of the nodes in the voxel that contains the source in the AM mesh. These divided currents
are then applied as a new static current source at each node.

This relatively complex algorithm was also compared with a simpler current
shifting algorithm which applies all current sources at the nearest respective node in the
mesh according to fig. 8. This method became interesting because of the relative
complexity of finding a voxel which contains a source location versus just finding the
nearest mesh node to that current source. LFP estimations via these two algorithms were
compared to each other and a line-source estimated LFP for both computational
performance and estimate accuracy.
Once a resistor network was constructed and properly grounded, either the
splitting or shifting algorithms were used to add a set of currents to the model. While the
geometry of the model remained the same for every timestep, a new netlist was written for
each time step to account for changing current source amplitudes. Current source-node
assignments remained the same for every time step for a given mesh geometry and
resolution. These netlists were then loaded in series and nodal voltages were solved using
the matrix solvers described at length by Cela [29].
Fig. 3.9. This flow diagram provides an overview of the analysis performed in this
manuscript. To establish best practices for interfacing volume conduction models
with large-scale NEURON network models, we explored the impact of mesh
geometry and current source handling algorithms on computational performance and
estimation accuracy. Validation consists of comparison of numerical model estimates
with a proxy ground truth in the form of analytical estimates (line-source equation).

97

4.3.7 Model Sharing and Software Distribution
It is our eventual goal to release a software tool that will enable other investigators
to optimize and later estimate LFPs within tissue models of arbitrary geometry and
conductivity. Therefore, code to generate a meshed volume, integrate current sources, solve
nodal voltages, and visualize results will be distributed via a public Github repository.
There remain, however, significant hurdles for adaptation of this method to different tissue
systems, including but not limited to: tissue geometry, heterogeneous dielectric properties,
and the complexities of constructing a biologically realistic NEURON model capable of
generating meaningful virtual extracellular currents. Further, the value of such a tool may
be reduced should we not focus future efforts on solving significant computational
bottlenecks to be raised in the discussion portion of this manuscript.

Fig. 3.9. To examine the impact of geometry on model and computational
complexity we created meshes of unstructured tetrahedral and structured hexahedral
geometry across a range of spatial resolutions and compared the (A) number of
elements required to mesh a set of nodes and (B) the relative RAM utilization of the
generated models. As models grew in complexity the tetrahedral to hexahedral
element to node ratio converged to ≈6x with tetrahedral models utilizing ≈6.5x more
memory.

98

4.4 Results
The following sections describe results of a comprehensive study of the Admittance
Method for LFP estimation. More specifically the outcome of sensitivity analysis
surrounding mesh geometry (tetrahedral vs. hexahedral) and current source handling
algorithms (splitting vs. shifting). Models were cross-compared for both computational
performance and accuracy with error calculated with respect to a line-source LFP estimate
as a proxy ground truth. This estimate was made at a recording location (denoted in Fig. 1)
that is biologically interesting due to the proximity to the cell body layer and the strong
population spike observed in the evoked potential.

Fig. 3.10. To examine the impact of current source handling algorithms (split vs.
shift) and mesh geometry (tetrahedral vs. hexahedral), meshes of a range of
resolutions were solved and solve times were compared by both geometry and
current source handling algorithm. (A) At low resolutions, shifting sources
outperformed splitting for both mesh geometries but the relative performance
converged to parity at greater model resolutions (≈19k nodes). (B) For all mesh
resolutions studied, tetrahedral meshes were solved more quickly than hexahedral,
though this difference was diminished when using the shifting current handling
algorithm and performance parity was achieved at higher model resolutions (≈15k
mesh nodes.)

99

4.4.1 New Meshes for LFP Estimation: tetrahedral vs.
hexahedral
Neither tetrahedral nor hexahedral meshes were universally superior in the
analyses performed. For very coarse meshes, tetrahedral models were markedly faster to
solve than hexahedral but not without substantial cost in both dynamic memory utilization,
model complexity, and accuracy (Figs. 10, 11, and 12). As models increased in resolution
memory and accuracy differences between the two geometries diminished but tetrahedral
models maintained ≈6-6.5x greater complexity and memory requirements due to higher
element to node ratios.

Fig. 3.11. To examine the estimation accuracy of numerical LFP solutions under
differing current source handling methodologies (shift vs. split) and mesh geometries
(tetrahedral vs. hexahedral), a range of models with varying resolution were
generated, solved, and quantitatively compared with an analytical estimate via the
line-source equation at a particular recording location. (A) While splitting
outperformed shifting at very low resolutions, as meshes became more complex,
shifting and splitting methods reached comparable levels of error (RMSE) with
current shifting consistently achieving lower error for hexahedral meshes above
≈12k mesh nodes. (B) Across the entire range of resolutions studied, hexahedral
meshes achieved lower error. Hexahedral meshes achieved accuracies comparable to
tetrahedral meshes at around half the number of mesh nodes and 1/12 as many
elements.

100

4.4.2 Incorporating incident endogenous current sources: split
vs. shift
At very coarse resolutions current splitting had the advantage of 11-30% better
accuracy at a computational penalty of 220-350% longer solve times. These coarse model
differences were more exaggerated for hexahedral meshes due to the greater number of
nodes per mesh element. Regardless of mesh geometry, differences in accuracy quickly
dissolved by the time node counts reached ≈4k where current shifting in hexahedral meshes
even managed to achieve better accuracy than splitting at higher resolutions.

Because the number and spatial distribution of current sources remained the same
for every meshed model, the preprocessing step of assigning sources to voxels or nodes was
particularly impactful for very coarse meshed models. While this step involves pairwise
relational distance calculation for millions of points and, therefore, increases in complexity
with increasingly complex meshes, it did not increase nearly as fast as did the time required
to solve these circuits.
Fig. 3.12. Example LFP time-series for current shifted tetrahedral and hexahedral
models plotted alongside the analytical (line-source equation). The tetrahedral model
had xx nodes and the hexahedral had xx nodes. The tetrahedral estimate had 14.1%
RMSE with respect to the analytical signal compared to 9.6% RMSE for the
hexahedral signal.

101

4.4.3 Additional Comments on Estimation Accuracy
While hexahedral meshes obtained the greatest accuracy when considering the
entire LFP, poorer accuracy in estimation of the evoked signal was compensated for by
superior accuracy in estimating the stimulus artifact phase of compartmental currents. As
can be seen from Fig. 13, at the highest resolutions studied, the tetrahedral model obtains a
deeper population spike (+54-63ms) than does the hexahedral estimate. However, it also
overshoots the earlier parts of the signal. Indeed, the tetrahedral model manages greater
amplitudes through the whole signal. This difference in amplitude was a consistent trend
for all resolutions and indicates that unstructured tetrahedral meshes have a greater net
resistance than do structured hexahedral meshes of comparable resolution. This is a unique
feature of unstructured vs. structured tetrahedral meshes which have a lower net resistance
than do structured hexahedral meshes constructed from the same set of nodes. This is,
admittedly, the logical outcome of using Eq. 3 to assign edge resistance for both mesh
geometries, tetrahedral and hexahedral.
Fig. 3.13. Synthesis of feedforward AM-NEURON methods demonstrated in
Bingham et al., [106] and feedback AM-NEURON examined in this manuscript via
current splitting yields an algorithm for bidirectional communication in mechanistic
models of neural tissue systems and volume conduction models.
Fig. 3.14. Synthesis of feedforward AM-NEURON methods demonstrated in
Bingham et al., [106] and feedback AM-NEURON examined in this manuscript via
current shifting yields an algorithm for bidirectional communication in mechanistic
models of neural tissue systems and volume conduction models. Except in very
coarse meshes, this more simplistic algorithm produced results comparable to
current splitting in hexahedral or tetrahedral models.

102

4.5 Discussion
The ability to accurately estimate electric fields generated in neural tissue by
stimulating electrodes, evoked behavior of neural tissue that lies within that volume, and
the electric fields that are generated by that behavior completes the critical components of
an algorithm that couples feedforward and feedback interactions between virtual electrical
stimulating, recording devices, and computational models of neural tissue. An algorithm for
bidirectional communication between AM and NEURON domains could prove a powerful
approach to understanding electrode-based brain-machine-interfaces and electroceutical
devices. However, effective use of this algorithm requires that modelers first construct
meaningful mechanistic models. The general value of this algorithm is limited by the
realism of the interfacing models: (i) volume conduction model design is complicated by
geometric and anatomic variety of electrodes and tissue and (ii) prediction of neural tissue
system behavior is complicated by the nuances of tuning biologically detailed models of
neurons and synapses. Through synthesis of results presented in this paper and the papers
that provided the feedforward AM-NEURON model described previously, we propose a
bidirectional algorithm that deepens the cross-scale links between bulk-tissue dielectric
behavior and cellular and network behavior of complex neural tissue systems (Figs. 13 and
14).
A stable and efficient approach to solving fields arising in tissue from endogenous
current sources allows us to begin new inquiry into several new areas of work. Some of
these include (i) providing new clarity regarding the extent to which volume conduction of
transmembrane currents reinforces population synchrony, (ii) use of the model framework
to validate and optimize spike-sorting algorithms for engineering applications, as well as
more fundamental questions about the biology of specific neural circuits by determining
with greater accuracy which components of a neuron or neurons in a network actually cause
the features that emerge in evoked potentials [3]. Each of these questions represents an
avenue of work enabled and encouraged by the modeling framework presented herein.

103

While models constructed via this proposed algorithm may help investigators
answer many questions related to brain-computer-interfaces and their useful application, it
is not without limitations.

We have demonstrated and quantitatively evaluated an approach to estimating
electric fields that arise from endogenous current sources. This is a general method that can
be applied to the study of any neural tissue system with a few important caveats. As our
study system we extended an AM-NEURON model of the rat hippocampus; many features
of this tissue system proved critical for accurate prediction of spatiotemporal patterns of
activity. It is likewise critical to capture the principle geometry and topology of neural
networks and dielectric properties of various regions of tissue systems being modeled.
When combined, these tissue specific properties sculp the geometry of electric fields
generated by neuronal activity. Further, tissue, behavioral, and environmental variables can
change AM meshing requirements and model parameters. This constitutes a substantial
hurdle for proper implementation of AM for LFP estimation; the need for mesh generation
and refinement is a limitation of this modeling approach though not one that disadvantages
AM relative to alternative numerical methods that must face similar meshing challenges
(e.g. FEM) [118].

4.5.1 Spatial Resolution of Sources and Meshes
One feature of model performance that drew our attention and deserves comment
is that the lion-share of error came from the population-spike phase of the evoked potential
(+54-63 ms) and errors arising from the artifact phase (+0-5 ms) were minimized at much
lower mesh resolutions. Because the recording location was so far from the stimulation
location but population-spike activity is very local, LFPs can be estimated accurately with
very coarse meshes without adequately capturing the activity nearest to the virtual
recording electrode [113]. Minimizing error in both population spike phase and artifact
phase likely requires two meshes of differing resolutions or a mesh with two resolutions,
near and distal to the recording site.

104

4.5.2 Mesh Quality
Mesh design challenges represent one specific limitation for most numerical
approaches to solving electromagnetic fields, the Admittance Method included. For every
new problem and new tissue geometry, new meshes must be designed and there is not
likely to be any one optimal mesh volume, fundamental geometry, or resolution.
Technically speaking, each model must find a new optimum for the following three
variables in addition to geometry: resolution, size/volume, placement of the mesh relative
to observation points/sources.

4.5.3 Solvers and Model Complexity
Computational demands and model complexity represents another critical
challenge that must be faced. Approaches to mitigating these demands are being developed
by the authors and others who have a need for mature EM numerical solvers. Potential
approaches include mesh reduction or refinement and parallel solution of the Admittance
Method. Alternative solvers may be considered which have already implemented massively
parallel solution of RC circuits such as the Xyce environment distributed by Sandia
National Labs [119].

Interestingly, for LFP estimation problems it is quite likely that mesh/circuit
solution is not the only computational bottleneck. Preprocessing is not a trivial step in the
process of numerical solution of LFPs. Because this study involved so many simulations
and at varying resolutions, it required repeated discretization of current sources into their
respective mesh model elements. For very small problems this is trivial but computing the
pairwise relational distance of millions of points is very time-consuming and memory
intensive. Our solution was to parallelize this problem but an ideal tool that maximizes
impact in the neural computing community would not require use of a high-performance

105

computing (HPC) cluster. While the authors had access to HPC resources the
implementation that will eventually be distributed will only require a single-machine and
will be multi-thread enabled. This means that investigators must be especially aware of the
potential performance impact of increased quantites and spatial distribution of current
sources relative to their meshed models.

4.6 Conclusion
The present study demonstrates an intuitive numerical approach to estimating local
field potentials generated by a detailed computational model of cortical tissue. AM-
NEURON LFP estimation, the analysis presented herein, and its synthesis with earlier work
makes three recognizable contributions to model-based analysis of complex neural systems
and electroceutical devices: (i) validation of a fundamental approach to incorporating
incident and dynamic current sources into a pre-existing meshed model, (ii) establishment
of an equivalent-circuit method that unifies intracellular and extracellular neuronal models,
and (iii) clarification of the impact of mesh geometry on computational performance and
model accuracy.

106

5 CONCLUSION
Multi-scale modeling of complex neural systems stands between theoretical and
experimental or clinical investigation of neurophysiological mechanisms and is uniquely
capable of providing a deep level of interpretability. Enhancement of the biological realism
of such mechanistic models deepens the community’s understanding of probable outcomes
of varying biological, device, and stimulus parameters. This unmatched level of model
control is sorely needed in the neural engineering of clinically relevant electroceutical
devices where clinicians remain unsure of the value of applying recommended vs.
alternative stimulus parameters or the effectiveness of one anatomical target over another.
When considering potential devices, such as a proposed hippocampal prosthesis, this doubt
is compounded due to a dearth of clinical trials proving such a therapy and, therefore, the
value of multi-scale models is amplified. This thesis demonstrates progress in designing
and improving such a model and outlines its use in pursuit of unanswered questions related
to efficacious interaction between extracellular electrical stimulation and recording devices
and complex neural tissue systems. Where the models described in this thesis end, others
may continue to provide valuable insights by drawing from lessons learned from the efforts
described here.

107

Multi-scale modeling of the central nervous system that spans molecular, cellular,
network, and tissue system levels is a rapidly expanding approach to elucidating emergent
phenomena and potential mechanisms of external control [120]. These models are being
used to understand varying aspects of cognition [121, 122], interactions between sensory
streams and disease states such as schizophrenia [123], linking imaging and other recording
modalities to behavior, enhancement of non-invasive neuromodulation [124, 125],
integrating models of brain vasculature and neural components [126], support the
development of superior brain activity decoding models [127], elucidating the impact of
varying levels of signaling molecules on network behavior [128], and understanding the
impact of the immune response on the integrity of brain machine interfaces [129].
Alongside our own efforts at building such multi-scale models [107, 130] [106] [108] [36],
these can be reasonably extended to understand how brain behavior can be modulated to
encourage more healthy signaling or to discourage aberrant signaling (a la epilepsy).

As the accuracy of multi-scale models in recapitulating complex behaviors and
human imaging modalities become more sophisticated, the opportunity to use computer-
aided prosthesis design to yield immediately actionable conclusions to aid clinicians is
clear. Further, many modeling tasks, such as volume fitting, neural model population,
implant positioning, and dielectric property assignment can be automated and then
simulations can be used to support the complex task of DBS parameter tuning [9].

While already FDA cleared devices provide a proving ground for such automated
technologies it also provides a useful opportunity to examine potential closed-loop
cognitive prostheses as an approach to minimizing power requirements and the
oppressiveness of a control paradigm (less time actively stimulating) [131].

Further, an underappreciated advantage of computer model-based studies over
traditional experimentation is reduced ethical externalities of biomedical research in
animals and humans: many clinical questions can be examined before the opportunity for

108

harm is introduced in animal or human test subjects. Though some new ethical
considerations must inevitably be made before embarking on the construction and
application of multi-scale models of brain structures [132]. Unlike early neurosurgery
experiments done by Penfield and others, new neurological therapies can now be supported
by more than just animal models but also studied using computational models to further
justify decisions made in the operating room. Simulating brain parts and prosthetics in a
computer reduces the number of exploratory surgical experiments that need to be
performed on animals and humans, it has the potential to reduce the likelihood that costly
mistakes will made both in the operating room and in device design, and it can deepen our
understanding of the underlying mechanics of both disease and brain-machine interfaces.

A dozen brilliant research programs have driven the astounding progress and
breakthroughs over the past few decades that have established what is possible with neural
prosthetics technologies. From cochlear implants, retinal implants, to DBS and other
pacemaker-like implants, the restorative potential of these devices is a new and exciting
frontier. However, like all modern technologies, this technology is bigger than any one
engineer, scientist, or physician. Material science, electrophysiology, neurology, surgical
techniques, circuit design, computer modeling, and imaging have all gone into the
development of efficacious electroceuticals and neural prosthetics. It is our collective goal
that the interfacing of these fields will continue to yield a class of devices that will help us
to see past the neurological disorder and provide new inclusiveness to a historically
neglected segment of society.

109

6 ACKNOWLEDGEMENTS
The author thanks Sydney Nguyen and Daniel Nemeth for testing the models used
in chapter four and suggesting various improvements to the code-base in preparation for
software distribution.

The author also thanks Romil Nathan for testing the ROOTS system and
suggesting various improvements to the code-base in preparation for software distribution.

110

7 REFERENCES

[1] T. Berger, R. Hampson, D. Song, A. Goonawardena, V. Marmarelis and S.
Deadwyler, "A cortical neural prosthesis for restoring and enhancing memory,"
Journal of neural engineering, vol. 8, no. 4, p. 046017, 2011.
[2] R. Hampson, D. Song, B. Robinson, D. Fetterhoff, A. Dakos, B. Roeder, X. She, R.
Wicks, M. Witcher, D. Couture and A. Laxton, "Developing a hippocampal neural
prosthetic to facilitate human memory encoding and recall," Journal of neural
engineering, pp. 15(3), 036014, 2018.
[3] G. Einevoll, F. Franke, E. Hagen, C. Pouzat and K. Harris, "Towards reliable spike-
train recordings from thousands of neurons with multielectrodes," Current opinion in
neurobiology, vol. 22, no. 1, pp. 11-17, 2012.
[4] M. Spira and A. Hai, "Multi-electrode array technologies for neuroscience and
cardiology," Nature nanotechnology, vol. 8, no. 2, p. 83, 2013.
[5] P. Brunner, L. Bianchi, C. Guger, F. Cincotti and G. Schalk, "Current trends in
hardware and software for brain–computer interfaces (BCIs)," Journal of neural
engineering, vol. 8, no. 2, p. 025001, 2011.
[6] S. Cogan, "Neural stimulation and recording electrodes," Annu. Rev. Biomed. Eng.,
no. 10, pp. 275-309, 2008.
[7] T. Foutz and C. McIntyre, "Evaluation of novel stimulus waveforms for deep brain
stimulation," Journal of neural engineering, vol. 7, no. 6, p. 066008, 2010.

111

[8] W. Grill, "Modeling the effects of electric fields on nerve fibers: influence of tissue
electrical properties," IEEE Transactions on Biomedical Engineering, pp. 46(8), 918-
928, 1999.
[9] D. Heldman, C. Pulliam, E. Urrea Mendoza, M. Gartner, J. Giuffrida, E.
Montgomery Jr, A. Espay and F. Revilla, "Computer‐Guided Deep Brain Stimulation
Programming for P arkinson's Disease," Neuromodulation: Technology at the Neural
Interface, vol. 19, no. 2, pp. 127-132, 2016.
[10] C. Bingham, K. Loizos, G. Yu, A. Gilbert, J.-M. Bouteiller, D. Song, G. Lazzi and T.
Berger, "Model-Based Analysis of Electrode Placement and Pulse Amplitude for
Hippocampal Stimulation," IEEE Transactions on Biomedical Engineering, pp. vol.
PP, no. 99, pp. 1-1., 2018.
[11] C. Bingham, K. Loizos, G. Yu, A. Gilbert, J. Bouteiller, D. Song, G. Lazzi and T.
Berger, "A large-scale detailed neuronal model of electrical stimulation of the dentate
gyrus and perforant path as a platform for electrode design and optimization," IEEE
Engineering in Medicine and Biology Society (EMBC), , pp. pp. 2794-2797, 2016.
[12] R. Anderson, A. Farokhniaee, K. Gunalan, B. Howell and C. McIntyre, "Action
potential initiation, propagation, and cortical invasion in the hyperdirect pathway
during subthalamic deep brain stimulation," Brain stimulation, vol. 11, no. 5, pp.
1140-1150, 2018.
[13] C. A. Anastassiou and C. Koch, "Ephaptic coupling to endogenous electric field
activity: why bother?," Current opinion in neurobiology, pp. 31, 95-103, 2015.
[14] K. Loizos, C. Cela, R. Marc and G. Lazzi, "Virtual electrode design for increasing
spatial resolution in retinal prosthesis," Healthcare technology letters, vol. 3, no. 2,
pp. 93-97, 2016.
[15] K. Loizos, A. RamRakhyani, J. Anderson, R. Marc and G. Lazzi, "On the
computation of a retina resistivity profile for applications in multi-scale modeling of
electrical stimulation and absorption," Physics in Medicine & Biology, vol. 61, no.
12, p. 4491, 2016.
[16] B. Howell, S. Lad and W. Grill, "Evaluation of intradural stimulation efficiency and
selectivity in a computational model of spinal cord stimulation," PloS one, vol. 9, no.
12, p. e114938, 2014.
[17] B. Howell and W. Grill, "Evaluation of high-perimeter electrode designs for deep
brain stimulation," Journal of Neural Engineering, pp. 11(4), 46026, 2014.

112

[18] H. Kasi, W. Hasenkamp, G. Cosendai, A. Bertsch and P. Renaud, "Simulation of
epiretinal prostheses-Evaluation of geometrical factors affecting stimulation
thresholds," Journal of neuroengineering and rehabilitation, vol. 8, no. 1, p. 44,
2011.
[19] D. Tsai, S. Chen, D. Protti, J. Morley, G. Suaning and N. Lovell, "Responses of
retinal ganglion cells to extracellular electrical stimulation, from single cell to
population: model-based analysis," PloS one, vol. 7, no. 12, p. e53357, 2012.
[20] B. Howell, B. Huynh and W. M. Grill, "Design and in vivo evaluation of more
efficient and selective deep brain stimulation electrodes," Journal of Neural
Engineering, pp. 12(4), 46030, 2015.
[21] H. Martens, E. Toader, M. Decré, D. Anderson, R. Vetter, D. Kipke, K. Baker, M.
Johnson and J. Vitek, "Spatial steering of deep brain stimulation volumes using a
novel lead design," Clinical neurophysiology, vol. 122, no. 3, pp. 558-5, 2011.
[22] C. McIntyre and W. Grill, "Excitation of Central Nervous System Neurons by
Nonuniform Electric Fields," Biophysical Journal, pp. 76(2), 878–888., 1999.
[23] S. Miocinovic, S. F. Lempka, G. S. Russo, C. B. Maks, C. R. Butson, K. E. Sakaie
and C. C. McIntyre, "Experimental and theoretical characterization of the voltage
distribution generated by deep brain stimulation.," Experimental neurology, p. 21,
2009.
[24] W. Grill, "Modeling the effects of electric fields on nerve fibers: influence of tissue
electrical properties," IEEE Transactions on Biomedical Engineering, pp. 46(8), 918-
928, 1999.
[25] B. Howell and C. McIntyre, "Role of Soft-Tissue Heterogeneity in Computational
Models of Deep Brain Stimulation," Brain Stimulation, 2016.
[26] C. A. Bossetti, M. J. Birdno and W. M. Grill, " Analysis of the quasi-static
approximation for calculating potentials generated by neural stimulation," Journal of
neural engineering, pp. 5(1), 44, 2007.
[27] C. Butson and C. McIntyre, "Tissue and electrode capacitance reduce neural
activation volumes during deep brain stimulation," Clinical neurophysiology, vol.
116, no. 10, pp. 2490-2500, 2005.
[28] M. Hines and N. Carnevale, "The NEURON simulation environment," Neural
computation, vol. 9, no. 6, pp. 1179-1209, 1997.
[29] C. Cela, A Multiresolution Admittance Method for Large-Scale Bioelectromagnetic
Interactions, North Carolina State University, 2013.

113

[30] C. S. Bingham, K. Loizos, G. J. Yu, A. Gilbert, J.-M. Bouteiller, D. Song, G. Lazzi
and T. W. Berger, "A large-scale detailed neuronal model of electrical stimulation of
the dentate gyrus and perforant path as a platform for electrode design and
optimization," in IEEE Engineering in Medicine and Biology Society (EMBC),
Orlando, Fl, 2016.
[31] J. Xie, G. J. Wang, L. Yow, C. J. Cela, M. S. Humayun, J. D. Weiland and H. Jadvar,
"Modeling and percept of transcorneal electrical stimulation in humans," IEEE
Transactions on Biomedical Engineering, pp. 1932-1939, 2011.
[32] W. J. Grill, "Modeling the effects of electric fields on nerve fibers: influence of tissue
electrical properties," Biomedical Engineering, IEEE Transactions on, pp. 46(8),
918-928., 1999.
[33] F. Rattay, D. Jaeger and R. Jung, in Encyclopedia of computational neuroscience,
Springer Publishing Company, Incorporated., 2015.
[34] A. Herz, T. Gollisch, C. Machens and D. Jaeger, "Modeling single-neuron dynamics
and computations: a balance of detail and abstraction," Science, vol. 314, no. 5796,
pp. 80-85, 2006.
[35] L. Nowak and J. Bullier, " Axons, but not cell bodies, are activated by electrical
stimulation in cortical gray matter," Experimental Brain Research, pp. 118(4), 477–
488. , 1998.
[36] P. J. Hendrickson, J. Y. Gene, D. Song and T. W. Berger, "A million-plus neuron
model of the hippocampal dentate gyrus: critical role for topography in determining
spatiotemporal network dynamics," IEEE Transactions on Biomedical Engineering,
p. 63(1, 2016.
[37] L. Kjonigsen, T. Leergaard, J. Kjode, M. Witter and J. Bjaalie, "Interactive
histological atlas system for anatomical parcellation of the rat hippocampus proper,
fasciola and associated parahippocampal cortex," In Front. Neuroinform, 2008.
[38] G. Petrovich, N. Canteras and L. Swanson, "Combinatorial amygdalar inputs to
hippocampal domains and hypothalamic behavior systems," Brain research reviews,
vol. 38, no. 1-2, pp. 247-289, 2001.
[39] B. Vogt and G. Paxinos, "Cytoarchitecture of mouse and rat cingulate cortex with
human homologies," Brain Structure and Function, vol. 219, no. 1, pp. 185-192,
2014.
[40] L. Geddes, "Historical Evolution of Circuit Models for the Electrode-Electrolyte
Interface," Annals of Biomedical Engineering, pp. 25, 1-14, 1997.

114

[41] N. Logothetis, C. Kayser and A. Oeltermann, "In vivo measurement of cortical
impedance spectrum in monkeys: implications for signal propagation," Neuron, vol.
55, no. 5, pp. 809-823, 2007.
[42] S. Joucla, P. Branchereau, D. Cattaert and B. Yvert, "Extracellular neural
microstimulation may activate much larger regions than expected by simulations: a
combined experimental and modeling study," PLoS One, vol. 7, no. 8, p. e41324,
2012.
[43] S. Joucla and B. Yvert, "Modeling extracellular electrical neural stimulation: From
basic understanding to MEA-based applications," Journal of Physiology-Paris, pp.
106(3), 146–158, 2012.
[44] C. Gabriel, S. Gabriel and Y. Corthout, "The dielectric properties of biological
tissues: I. Literature survey," Physics in medicine & biology, vol. 41, no. 11, p. 2231,
1996.
[45] S. Gabriel, R. Lau and C. Gabriel, "The dielectric properties of biological tissues: III.
Parametric models for the dielectric spectrum of tissues," Physics in Medicine &
Biology, vol. 41, no. 11, p. 2271, 1996.
[46] M. Wilson, M. Elbohouty, L. Voss and D. Steyn-Ross, "Electrical impedance of
mouse brain cortex in vitro from 4.7 kHz to 2.0 MHz," Physiological Measurement,
pp. 35, 267–81, 2014.
[47] K. Cole and R. Cole, "Dispersion and Absorption in Dielectrics I. Alternating Current
Characteristics," The Journal of Chemical Physics, pp. 9(4), 341, 1941.
[48] L. López-Aguado, J. Ibarz and O. Herreras, "Activity-dependent changes of tissue
resistivity in the CA1 region in vivo are layer-specific: Modulation of evoked
potentials," Neuroscience, pp. 108(2), 249–262, 2001.
[49] W. Soussou, G. Gholmieh, M. Han, A. Ahuja, D. Song, M. Hsiao, Z. Wang, A.
Tanguay and T. Berger, "Mapping Spatio-Temporal Electrophysiological Activity in
Hippocampal Slices with Conformal Planar Multi-Electrode Arrays," in Advances in
Network Electrophysiology - Using Multi-Electrode Arrays, Boston, MA., Springer,
2006, pp. (4), 127–152.
[50] N. Tamamaki and Y. Nojyo, "Projection of the entorhinal layer II neurons in the rat
as revealed by intracellular pressure-injection of neurobiotin," Hippocampus, pp.
3(4), 471–80, 1993.
[51] G. Ascoli and J. Krichmar, "L-neuron: A modeling tool for the efficient generation
and parsimonious description of dendritic morphology," Neurocomputing, pp. 32,
1003–1011, 2000.

115

[52] R. Williams and S. Matthysse, "Morphometric analysis of granule cell dendrites in
the mouse dentate gyrus," Journal of Comparative Neurology, vol. 215, no. 2, pp.
154-164, 1983.
[53] V. Santhakumar, I. Aradi and I. Soltesz, "Role of mossy fiber sprouting and mossy
cell loss in hyperexcitability: a network model of the dentate gyrus incorporating cell
types and axonal topography," Journal of neurophysiology, vol. 93, no. 1, pp. 437-
453, 2005.
[54] G. Yuen and D. Durand, "Reconstruction of hippocampal granule cell
electrophysiology by computer simulation," Neuroscience, vol. 41, no. 2-3, pp. 411-
423, 1991.
[55] G. Holt and C. Koch, "Electrical interactions via the extracellular potential near cell
bodies," Journal of Computational Neuroscience, pp. 6(2), 169–184, 1999.
[56] J. Clark and R. Plonsey, "A Mathematical Study of Nerve Fiber Interaction,"
Biophysical Journal, pp. 10(10), 937–957, 1970.
[57] S. Deadwyler, J. West, C. Cotman and G. Lynch, "Physiological studies of the
reciprocal connections between the hippocampus and entorhinal cortex,"
Experimental neurology, vol. 49, no. 1, pp. 35-57, 1975.
[58] E. Hu, J. Bouteiller, D. Song, M. Baudry and T. Berger, "Volterra representation
enables modeling of complex synaptic nonlinear dynamics in large-scale
simulations," Frontiers in computational neuroscience, vol. 9, p. 112, 2015.
[59] M. Elbohouty, M. Wilson, L. Voss, D. Steyn-Ross and L. Hunt, "In vitro electrical
conductivity of seizing and non-seizing mouse brain slices at 10 kHz," Physics in
Medicine & Biology, vol. 58, no. 11, p. 3599, 2013.
[60] L. Fritz and A. Gardner-Medwin, "The effect of synaptic activation on the
extracellular potassium concentration in the hippocampal dentate area, in vitro,"
Brain research, vol. 112, no. 1, pp. 183-187, 1976.
[61] L. Nowak and J. Bullier, "Spread of stimulating current in the cortical grey matter of
rat visual cortex studied on a new in vitro slice preparation," Journal of Neuroscience
Methods, pp. 67(2), 237–248, 1996.
[62] J. Ranck, "Which elements are excited in electrical stimulation of mammalian central
nervous system: A review," Brain Research, pp. 98(3), 417–440, 1975.
[63] F. Rattay, L. Paredes and R. Leao, "Strength–duration relationship for intra-versus
extracellular stimulation with microelectrodes," Neuroscience, pp. 214, 1-13, 2012.

116

[64] M. Johnson and C. McIntyre, "Quantifying the neural elements activated and
inhibited by globus pallidus deep brain stimulation," J Neurophysiol , pp.
100(5):2549-2563, 2008.
[65] C. A. Anastassiou and C. Koch, "Ephaptic coupling to endogenous electric field
activity: why bother?," Current opinion in neurobiology, pp. 31, 95-103., 2015.
[66] J. Kim, H. Lee, W. Choi and K. Lee, "Encoding information into autonomously
bursting neural network with pairs of time-delayed pulses," Scientific reports, vol. 9,
no. 1, p. 1394, 2019.
[67] M. Whittington, I. Stanford, S. Colling, J. Jefferys and R. Traub, "Spatiotemporal
patterns of ? frequency oscillations tetanically induced in the rat hippocampal slice,"
The Journal of Physiology, pp. 502(3), 591–607, 1997.
[68] P. Fries, "A mechanism for cognitive dynamics: neuronal communication through
neuronal coherence," Trends in Cognitive Sciences, pp. 9(10), 474–480, 2005.
[69] B. Claiborne, D. Amaral and W. Cowan, "Quantitative, three-dimensional analysis of
granule cell dendrites in the rat dentate gyrus," J. Comp. Neurol, p. 302: 206–219. ,
1990.
[70] N. L. Desmond and W. B. Levy, "Granule cell dendritic spine density in the rat
hippocampus varies with spine shape and location," Neurosci. Lett., pp. vol. 54, nos.
2/3, pp. 219–224, 1985.
[71] K. Hama, T. Arii and T. Kosaka, "Three-dimensional morphometrical study of
dendritic spines of the granule cell in the rat dentate gyrus with HVEM stereo
images," J. Electron Microscopy Tech., pp. vol. 12, no. 2, pp. 80–87, 1989.
[72] A. Hjorth-Simonsen and B. Jeune, "Origin and termination of the hippocampal
perforant path in the rat studied by silver impregnation," Journal of Comparative
Neurology, pp. 144(2), 215-231, 1972.
[73] E. Türetken, G. González, C. Blum and P. Fua, "Automated reconstruction of
dendritic and axonal trees by global optimization with geometric priors,"
Neuroinformatics, pp. 9(2-3), 279-302, 2011.
[74] G. Rozenberg and A. Salomaa, The mathematical theory of L systems, Academic
press, 1980.
[75] J. Budd, K. Kovács, A. Ferecskó, P. Buzás, U. Eysel and Z. Kisvárday, "Neocortical
axon arbors trade-off material and conduction delay conservation," PLoS
computational biology, pp. 6(3), e1000711, 2010.

117

[76] H. Cuntz, A. Borst and I. Segev, "Optimization principles of dendritic structure,"
Theoretical Biology and Medical Modelling, pp. 4(1), 21, 2007.
[77] J. Budd and Z. Kisvarday, "Communication and wiring in the cortical connectome,"
Frontiers in Neuroanatomy, pp. 5(68), 37–49., 2012.
[78] H. Cuntz, F. Forstner, A. Borst, M. Häusser, S. Kuhlman and P. Saggau, "One Rule
to Grow Them All: A General Theory of Neuronal Branching and Its Practical
Application," PLoS Computational Biology, pp. 6(8), e1000877. , 2010.
[79] M. Witter, "The perforant path: projections from the entorhinal cortex to the dentate
gyrus," Prog. Brain Res. , vol. 163, pp. 43-61, 2007.
[80] E. Bindocci, I. Savtchouk, N. Liaudet, D. Becker, G. Carriero and A. Volterra,
"Three-dimensional Ca2+ imaging advances understanding of astrocyte biology,"
Science, vol. , p. 356:6339, 2017.
[81] S. P. Schwartz and P. D. Coleman, "Neurons of origin of the perforant path,"
Experimental neurology, vol. 74, no. 1, pp. 305-312, 1981.
[82] P. J. Hendrickson, G. Yu, D. Song and T. W. Berger, "A million-plus neuron model
of the hippocampal dentate gyrus: critical role for topography in determining
spatiotemporal network dynamics," IEEE Transactions on Biomedical Engineering,
2016.
[83] P. E. Patton and B. McNaughton, "Connection matrix of the hippocampal formation:
I. The dentate gyrus," Hippocampus, pp. vol. 5, no. 4, pp. 245–286, 1995.
[84] F. B. Gaarskjaer, "Organization of the mossy fiber system of the rat studied in
extended hippocampi. I. Terminal area related to number of granule and pyramidal
cells," J. Comparative Neurol., pp. vol. 178, no. 1, pp. 49–72, 1978.
[85] A. Cayley, "A theorem on trees," Quart. J. Pure Appl. Math., p. 23: 376–378, 1889.
[86] J. Hartigan and M. Wong, "Algorithm AS 136: A k-means clustering algorithm,"
Journal of the Royal Statistical Society, pp. 28(1), 100-108., 1979.
[87] L. Chew, "Constrained delaunay triangulations," Algorithmica, pp. 4(1-4), 97-108,
1989.
[88] J. Chen, "Dijkstra’s shortest path algorithm," Journal of Formalized Mathematics, pp.
15(9), 237-247, 2003.

118

[89] M. Hines and N. Carnevale, The NEURON simulation environment. In: The
Handbook of Brain Theory and Neural Networks, 2nd ed, Cambridge, MA: MIT
Press, 2003, pp. 769-773.
[90] A. Saalfeld, "Topologically consistent line simplification with the Douglas-Peucker
algorithm," Cartography and Geographic Information Science, pp. 26(1), 7-18.,
1999.
[91] R. Scorcioni and G. Ascoli, "Algorithmic reconstruction of complete axonal
arborizations in rat hippocampal neurons," Neurocomputing, pp. 65, 15–22, 2005.
[92] R. Scorcioni, S. Polavaram and G. Ascoli, "L-Measure: a web-accessible tool for the
analysis, comparison and search of digital reconstructions of neuronal morphologies,"
Nature protocols, pp. 3(5), 866, 2008.
[93] P. Quilichini, A. Sirota and G. Buzsáki, "Intrinsic circuit organization and theta–
gamma oscillation dynamics in the entorhinal cortex of the rat," Journal of
Neuroscience, vol. 30, no. 33, pp. 11128-11142, 2010.
[94] T. Hu, "Optimum Communication Spanning Trees," SIAM J of Comput., pp. 3, 188–
195, 1974.
[95] C. Alpert, T. Hu, J. Huang, A. Kahng and D. Karger, "Prim-Dijkstra tradeoffs for
improved performance-driven routing tree design.," IEEE Transactions on Computer-
Aided Design of Integrated Circuits and Systems, pp. vol. 14, no. 7, pp. 890-896.,
1995.
[96] S. Khuller, B. Raghavachari and N. Young, "Balancing minimum spanning trees and
shortest-path trees," Algorithmica, p. 14: 305, 1995.
[97] B. Wu, K.-M. Chao and C. Tang, "Approximation algorithms for some optimum
communication spanning tree problems," Discrete Applied Mathematics, pp. Vol.
102:3, pgs. 245-266, 2000.
[98] M. Gastner and M. Newman, "Shape and efficiency in spatial distribution networks,"
J. Stat. Mech., p. P01015. , 2006.
[99] A. Al-Humaidi, Resistor networks and finite element models, The University of
Manchester, 2011.
[100] L. López-Aguado, J. Ibarz and O. Herreras, "Activity-dependent changes of tissue
resistivity in the CA1 region in vivo are layer-specific: Modulation of evoked
potentials," Neuroscience, pp. 108(2), 249–262, 2001.

119

[101] C. McIntyre, "Computational modeling of deep brain stimulation," Neuromodulation,
pp. 171-178, 2009.
[102] T. J. Foutz, R. L. Arlow and C. C. McIntyre, "Theoretical principles underlying
optical stimulation of a channelrhodopsin-2 positive pyramidal neuron," American
Journal of Physiology-Heart and Circulatory Physiology, 2012.
[103] J. M. C. Bouteiller and T. W. Berger, " The Role of Simulations in
Neuropharmacology," in In Computational Neurology and Psychiatry, Springer,
Cham., 2017, pp. (pp. 429-448)..
[104] J. L. Lujan, A. Chaturvedi, K. S. Choi, P. E. Holtzheimer, R. E. Gross, H. S. Mayberg
and C. C. McIntyre, "Tractography-activation models applied to subcallosal cingulate
deep brain stimulation," Brain stimulation, pp. 6(5), 737-739, 2013.
[105] I. Michailidou, J. Willems, E. Kooi, C. van Eden, S. Gold, J. Geurts, F. Baas, I.
Huitinga and V. Ramaglia, "Complement C 1q-C 3–associated synaptic changes in
multiple sclerosis hippocampus," Annals of neurology, pp. 77(6), 1007-1026, 2015.
[106] C. S. Bingham, K. Loizos, G. J. Yu, A. Gilbert, J.-M. Bouteiller, D. Song, G. Lazzi
and T. W. Berger, "Model-Based Analysis of Electrode Placement and Pulse
Amplitude for Hippocampal Stimulation," IEEE Transactions on Biomedical
Engineering, pp. vol. PP, no. 99, pp. 1-1, 2018.
[107] C. S. Bingham, J. M. C. Bouteiller, D. Song and T. W. Berger, "Graph-Based Models
of Cortical Axons for the Prediction of Neuronal Response to Extracellular Electrical
Stimulation," in IEEE Engineering in Medicine and Biology Society (EMBC),
Honolulu, 2018.
[108] A. Gilbert, C. Bingham, K. Loizos, A. RamRakhyani, P. Hendrickson, G. Lazzi and
T. Berger, "A 3-D admittance-level computational model of a rat hippocampus for
improving prosthetic design," IEEE Engineering in Medicine and Biology Society
(EMBC), pp. 2295-2298, 2015.
[109] J. Cline, C. Bingham, K. Loizos, G. Yu, P. Hendrickson, J. Bouteiller, T. Berger and
G. Lazzi, "Estimation of initiated local field potential by neurons in heterogeneous
tissue environment using admittance method," IEEE USNC-URSI, pp. 310-310, 2015.
[110] C. A. Bossetti, M. J. Birdno and W. M. Grill, " Analysis of the quasi-static
approximation for calculating potentials generated by neural stimulation," Journal of
Neural Engineering, pp. 5(1), 44–53, 2008).
[111] C. C. McIntyre, "Computational modeling of deep brain stimulation,"
Neuromodulation, pp. 171-178, 2009.

120

[112] M. Bazhenov, P. Lonjers, S. Skorheim, C. Bedard and A. Destexhe, "Non-
homogeneous extracellular resistivity affects the current-source density profiles of
up-down state oscillations," Philosophical Transactions, pp. 369(1952), 3802–1,
2011.
[113] L. G. Nowak and J. Bullier, "Spread of stimulating current in the cortical grey matter
of rat visual cortex studied on a new in vitro slice preparation," Journal of
Neuroscience Methods, pp. 67(2), 237–248, 1996.
[114] F. Rattay, L. Paredes and R. Leao, "Strength–duration relationship for intra-versus
extracellular stimulation with microelectrodes," Neuroscience, pp. 214, 1-13, 2012.
[115] L. Zhang, W. Zhou and D. Li, "A descent modified Polak–Ribière–Polyak conjugate
gradient method and its global convergence," IMA Journal of Numerical Analysis,
vol. 26, no. 4, pp. 629-640, 2006.
[116] A. Al-Humaidi, Resistor networks and finite element models, United Kingdom: The
University of Manchester, 2001.
[117] R. Duffin, "Distributed and lumped networks," Journal of Mathematics and
Mechanics, pp. 793-826, 1959.
[118] D. Pridmore, G. Hohmann, S. Ward and W. Sill, "An investigation of finite-element
modeling for electrical and electromagnetic data in three dimensions," Geophysics,
vol. 46, no. 7, pp. 1009-1024, 1981.
[119] S. Hutchinson, E. Keiter, R. Hoekstra, H. Watts, A. Waters, T. Russo and C.
BOGDAN, "THE Xyce™ PARALLEL ELECTRONIC SIMULATOR–AN
OVERVIEW," Parallel Computing: Advances and Current Issues, pp. 165-172,
2002.
[120] R. Betzel and D. Bassett, "Multi-scale brain networks," Neuroimage, vol. 160, pp.
73-83, 2017.
[121] F. Turkheimer, R. Leech, P. Expert, L. Lord and A. Vernon, "The brain's code and its
canonical computational motifs. From sensory cortex to the default mode network: a
multi-scale model of brain function in health and disease," Neuroscience &
Biobehavioral Reviews, vol. 55, pp. 211-222, 2015.
[122] X. Liu, Y. Zeng, T. Zhang and B. Xu, "Parallel brain simulator: a multi-scale and
parallel brain-inspired neural network modeling and simulation platform," Cognitive
Computation, vol. 8, no. 5, pp. 967-981, 2016.
[123] M. Schirner, A. McIntosh, V. Jirsa, G. Deco and P. Ritter, "Inferring multi-scale
neural mechanisms with brain network modelling," Elife, vol. 7, p. e28927, 2018.

121

[124] H. Seo and S. Jun, "Multi-scale computational models for electrical brain
stimulation," Frontiers in human neuroscience, vol. 11, p. 515, 2017.
[125] M. Soldati, M. Mikkonen, I. Laakso, T. Murakami, Y. Ugawa and A. Hirata, "A
multi-scale computational approach based on TMS experiments for the assessment of
electro-stimulation thresholds of the brain at intermediate frequencies," Physics in
Medicine & Biology, vol. 63, no. 22, p. 225006, 2018.
[126] A. Marreiros, O. Eschenko and N. Logothetis, "Whole-brain mapping of state-
dependent cortical responses to electrical stimulation," In 22nd Annual Meeting of the
Organization for Human Brain Mapping, 2016.
[127] K. Amunts, C. Ebell, J. Muller, M. Telefont, A. Knoll and T. Lippert, "The human
brain project: creating a European research infrastructure to decode the human brain,"
Neuron, vol. 92, no. 3, pp. 574-581, 2016.
[128] C. Kaya, E. Block, A. Sorkin, J. Faeder and I. Bahar, "Multi-Scale Spatial
Simulations Reveal the Effect of Dopamine Transporter Localization on Dopamine
Neurotransmission," Biophysical Journal, vol. 110, no. 3, p. 632a, 2016.
[129] N. V. A. Michelson, J. Eles, J. Salatino, E. Purcell, J. Williams, X. Cui and T. Kozai,
"Multi-scale, multi-modal analysis uncovers complex relationship at the brain tissue-
implant neural interface: New Emphasis on the Biological Interface," Journal of
neural engineering, vol. 15, no. 3, p. 033001, 2018.
[130] E. Hu, A. Mergenthal, C. Bingham, D. Song, J. Bouteiller and T. Berger, "A
glutamatergic spine model to enable multi-scale modeling of nonlinear calcium
dynamics," Frontiers in Computational Neuroscience, vol. 12, 2018.
[131] M. Lo and A. Widge, "Closed-loop neuromodulation systems: next-generation
treatments for psychiatric illness," International Review of Psychiatry, vol. 29, no. 2,
pp. 191-204, 2017.
[132] M. Christen, N. Biller-Andorno, B. Bringedal, K. Grimes, J. Savulescu and H.
Walter, "Ethical challenges of simulation-driven big neuroscience," AJOB
Neuroscience, vol. 7, no. 1, pp. 5-17, 2016.
[133] B. Howell and C. McIntyre, "Role of Soft-Tissue Heterogeneity in Computational
Models of Deep Brain Stimulation," Brain Stimulation, 2016.
[134] C. McIntyre and W. Grill, "Excitation of Central Nervous System Neurons by
Nonuniform Electric Fields," Biophysical Journal, pp. 76(2), 878–888, 1999.

122

[135] L. Nowak and J. Bullier, "Spread of stimulating current in the cortical grey matter of
rat visual cortex studied on a new in vitro slice preparation," Journal of Neuroscience
Methods, pp. 67(2), 237–248, 1996.
[136] M. Whittington, I. Stanford, S. Colling, J. Jefferys and R. Traub, "Spatiotemporal
patterns of γ frequency oscillations tetanically induced in the rat hippocampal slice,"
The Journal of Physiology, pp. 502(3), 591–607, 1997.
[137] P. Fries, "A mechanism for cognitive dynamics: neuronal communication through
neuronal coherence," Trends in Cognitive Sciences, pp. 9(10), 474–480, 2005.
[138] M. Johnson and C. McIntyre, "Quantifying the neural elements activated and
inhibited by globus pallidus deep brain stimulation," J Neurophysiol, pp.
100(5):2549-2563, 2008.
[139] L. Kjonigsen, T. Leergaard, J. Kjode, M. Witter and J. Bjaalie, "Interactive
histological atlas system for anatomical parcellation of the rat hippocampus proper,
fasciola and associated parahippocampal cortex.," In Front. Neuroinform., 2008.
[140] M. Hines and N. Carnevale, The NEURON simulation environment. In: The
Handbook of Brain Theory and Neural Networks, Cambridge, MA: MIT Press, 2003,
pp. pp. 769-773..
[141] G. Einevoll, H. Lindén, T. Tetzlaff, S. Łeski and K. Pettersen, "Local field
potentials," in Principles of Neural Coding, 2013, p. 37.
[142] G. Holt and C. Koch, "Electrical interactions via the extracellular potential near cell
bodies," Journal of Computational Neuroscience, pp. 6(2), 169–184, 1999.
[143] J. Clark and R. Plonsey, "A Mathematical Study of Nerve Fiber Interaction,"
Biophysical Journal, pp. 10(10), 937–957, 1970.
[144] I. D. Hentall, G. Zorman, S. Kansky and H. L. Fields, "Relations among threshold,
spike height, electrode distance, and conduction velocity in electrical stimulation of
certain medullospinal neurons," Journal of Neurophysiology, p. 51(5), 1984.
[145] L. Chew, "Constrained delaunay triangulations," Algorithmica, pp. 4(1-4), 97-108,
1989.
[146] D. Merrill, M. Bikson and J. Jefferys, "Electrical stimulation of excitable tissue:
design of efficacious and safe protocols," Journal of neuroscience methods, vol. 141,
no. 2, pp. 171-198, 2005.

123

size-fits all problems mesh-generation
• Edge rather than node-centric algorithm because of volume
subdivision – can’t specify arbitrary nodes at which you want the
mesh to align
• Difficult to implement smooth source/sink gradients (focuses on
boundary proximity gradients)

Abstract (if available)

Abstract When communicating with neural systems, there are few in vivo techniques more widely used than electrode-based electrical stimulation and recording. One of the greatest feats of neural engineering has been the development of implantable medical devices which use multielectrode arrays to record and stimulate neurons for the therapeutic treatment of brain disorders, whether through modulation (i.e. deep brain stimulation) or restoration (i.e. hippocampal memory prosthesis). Advancements in medical implant technology have a trajectory toward “less and more”: smaller electronics and higher density arrays. However, the effectiveness of these engineering advances cannot be maximized without a fundamental understanding of the neural system that is the basis of the brain disorder, which includes the electrophysiology and the anatomy of the neuronal network and principle neurons. This knowledge can inform the design and placement of electrodes to determine the ideal location, orientation, number of electrodes, and stimulation pattern necessary to elicit the target response. ❧ Advancements in computing have opened doors for more detailed and accurate mechanistic models of neural tissue, and an impressive corpus of anatomical and biophysical literature is available to constrain them. This dissertation presents the construction and co-simulation of a large-scale compartmental model of an in vivo hippocampal neuronal network and a specialized equivalent-circuit model of the tissue volume-conductor to approximate the spread of electrical currents through extracellular space. This powerful dual model allows the prediction of neural interactions and responses to a variety of clinically relevant devices and stimulation protocols. Prior work with multielectrode arrays has been done to predict the hippocampal output given a known input. The proposed model will aid the design of multielectrode arrays and stimulation parameters which would enable reproduction of predicted output in a target neuronal population, thereby closing the loop and taking another step toward effective cortical neuro-restorative technology in the hippocampus. ❧ In the midst of constructing this model, important limitations were confronted and ameliorative solutions were designed to bolster the predictive power of such mechanistic models. This dissertation sets forth, in three chapters the fundamentals of mechanistic and parametric modelling of dual neuronal and volume-conduction models of complex neural tissue systems, generative models of highly-branched axon terminal arbors for use in grey matter stimulation models, and a candidate algorithm for bidirectional interactions between neuronal and volume-conduction domains. Chapters 3 and 4 represent independent efforts to address commonly encountered limitations—that were identified in the work discussed in chapter 1—in widely used stimulating and recording models of neural tissue. These models are often used to study deep brain stimulation and other central nervous system electroceuticals as well as more general hypothesis-driven investigation of neural circuit dynamics. ❧ While the generality of these modelling tools is clear from the literature, we use the hippocampus and potential hippocampal prostheses as a case study. In addition to gaining a deeper understanding of the hippocampal system features which helps form the basis of long-term memory formation, we expect this work to yield an efficient method of evaluating arrays and patterns of stimulating currents. This outcome represents an important step toward more effective, well-understood, and safe intracortical neural implants.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Functional consequences of network architecture in rat hippocampus: a computational study

PDF

Multi-region recordings from the hippocampus with conformal multi-electrode arrays

PDF

Therapeutic electrical stimulation strategies for neuroregeneration and neuroprotection of retinal neurons

PDF

Nonlinear modeling of causal interrelationships in neuronal ensembles: an application to the rat hippocampus

PDF

Data-driven modeling of the hippocampus & the design of neurostimulation patterns to abate seizures

PDF

Nonlinear models for hippocampal synapses for multi-scale modeling

PDF

Next generation neural interfaces for neural modulation and prosthesis

PDF

Experimental and computational models for seizure prediction

PDF

Decoding memory from spatio-temporal patterns of neuronal spikes

PDF

An in vitro model of a retinal prosthesis

PDF

Optimal electrical stimulation of smooth muscle

PDF

Nonlinear dynamical modeling of single neurons and its application to analysis of long-term potentiation (LTP)

PDF

Computational investigation of glutamatergic synaptic dynamics: role of ionotropic receptor distribution and astrocytic modulation of neuronal spike timing

PDF

Penetrating parylene neural probe array for dense, in vivo, chronic recordings

PDF

Manipulation of RGCs response using different stimulation strategies for retinal prosthesis

PDF

Computational approaches to establish safety and efficacy assessments of electrical stimulation to peripheral nerve

PDF

On the electrophysiology of multielectrode recordings of the basal ganglia and thalamus to improve DBS therapy for children with secondary dystonia

PDF

Developing a neural prosthesis for hippocampus: proof-of-concept using the in vitro slice

PDF

Adaptive octree meshing of the extracellular space in neural simulations

PDF

A biomimetic approach to non-linear signal processing in ultra low power analog circuits

Asset Metadata

Creator Bingham, Clayton Scott (author)

Core Title Simulating electrical stimulation and recording in a multi-scale model of the hippocampus

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Biomedical Engineering

Publication Date 12/09/2019

Defense Date 05/29/2019

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag computational neuroscience,computer modeling,electrical stimulation,electromagnetics,electrophysiology,hippocampus,medical devices,neuroprosthetics,OAI-PMH Harvest

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Berger, Theodore W. (committee chair), Bouteiller, Jean-Marie (committee member), Lazzi, Gianluca (committee member), Marmarelis, Vasilis (committee member)

Creator Email clayton.bingham@gmail.com,csb119@case.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-248178

Unique identifier UC11673174

Identifier etd-BinghamCla-8023.pdf (filename),usctheses-c89-248178 (legacy record id)

Legacy Identifier etd-BinghamCla-8023.pdf

Dmrecord 248178

Document Type Dissertation

Rights Bingham, Clayton Scott

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA