of 33
1
Multi
-
modal characterization and simulation of human epileptic circuitry
Anatoly Buchin
1
*
,
Rebecca de Frates
1,x
,
Anirban Nandi
1,x
,
Rusty Mann
1
,
Peter Chong
1
,
Lindsay
Ng
1
, Jeremy Miller
1
, Rebecca Hodge
1
,
Brian Kalmbach
1,2
,
Soumita
Bose
1,
3
,
Ueli Rutishauser
4
,
5
,
Stephen McConoughey
1
,
Ed Lein
1
,7
,
Jim Berg
1
,
Staci Sorensen
1
,
Ryder Gwinn
6
,
Christof
Koch
1
,
Jonathan Ting
1
,7
, Costas A. Anastassiou
1,
8
*
Affiliations
:
1
Allen Institute for Brain Science, Seattle, WA, USA*
2
University of Washington, Seattle, USA
3
CiperHealth
,
San Francisco
, USA
4
Cedars
-
Sinai Medical Center, Los Angeles CA, USA
5
California Institute of Technology, Pasadena, CA, USA
6
Swedish Medical C
enter, Seattle, WA, USA
7
University of Washington, Seattle, WA, USA
8
University of British Columbia, Vancouver, BC, CA
x
Equal author contribution
Current address: Institute for Advanced Clinical Trials for Children, WA, USA
*
Corresponding authors:
CAA (
costasa@alleninstitute.org
,
costas.anastassiou@gmail.com
)
and
AB
(
anatolyb@
alleninstitute.org
)
2
S
u
pporting
m
aterial
and methods
Tissue preparation
Human brain tissue from neurosurgical origin w
as
made available through the generosity of
tissue donors
(Table S1)
.
Excised human brain
tissue was collected from the brain within 1
-
3
minutes of resection (up to 5 minutes in rare cases), and transported in chilled, oxygenated
ACSF.VII from the hospital to the Allen Institute laboratories
within
15
-
30 minutes. Tissue was
then mounted for sl
ice preparation on the chuck of a Compresstome VF
-
200 or VF
-
300 vibrating
microtome (Precisionary Instruments) to be sliced perpendicular to pial surface.
Each human
tissue slice was mounted in the recording chamber and inspected to ensure that the entire
hippocampal depth was intact. Regions of the hippocampus were
identified
visually.
Dentate
gyrus granule cells along the g
ranule cell
layer
were
targeted for the full dataset.
More details on
tissue preparation can be found in
(
43
)
under section “Electrophysiology overview”.
Single nucleus RNA sequencing analysis
I
ntact nuclei from
the hippocampus of four
human tissue donors
H16.0
6.008 (WG1),
H17.06.015 (WG1), H16.06.009 (WG4) and H16.06.010 (WG4)
were collected and
processed for sequencing using 10X Genomics as follows. To generate single nuclei, flash frozen
hippocampal samples were Dounce homogenized, stained with mouse anti
-
N
euN antibody
(Millipore,
FCMAB317PE
), and NeuN+ (~70%) and NeuN
-
(~30%) nuclei were collected using
fluorescence activated cell sorting (FACS). Sorted nuclei suspensions were concentrated to
~1000 nuclei/μl and loaded onto a 10X Genomics version 3 single
-
c
ell RNA
-
sequencing chip
with a targeted capture of ~8000 nuclei per specimen. Nuclei (n = 29,212) were sequenced at a
median read depth of 83,474 ± 24,111 reads/nucleus (median ± standard deviation). Median gene
3
detection was 5295 ± 240 genes/nucleus and t
he median number of UMIs/nucleus was 16,088 ±
1790.
For more information on RNA
-
seq technology applied in human brain tissue, see
(
44
,
45
)
.
RNA
-
seq data was curated to remove nuclei with less
than 1’000 genes detected, leaving 27,242
nuclei remaining. Nuclei were then assigned to broad classes (GCs, other excitatory, inhibitory,
and non
-
neuronal) via supervised clustering, using the 100 most selective genes for each class.
GC genes were defi
ned as genes with the most significant differential expression between “DG”
and “HiF” in the “Differential Search” in the Allen Human Brain Atlas (
http://human.brain
-
map.org/
(
24
)
). Similarly, markers for the other classes were selected using the RNA
-
seq Data
Navigator: Human
a part of the Allen Cell Types
Database (
http://celltypes.brain
-
map.org
(
45
)
)
by running “Find Marker Genes” for “Ce
ll Category” against all other cells. This identified
genes specific for each class based on single nucleus RNA
-
seq in human middle temporal gyrus.
The RNA
-
seq data was normalized by taking the unique molecular identifier (UMI) count for
each cell, calcu
lating counts per million (CPM), and then converting to logarithmic space:
EXPR_norm = log2(CPM(UMI)+1). We then performed supervised
k
-
means clustering with
k
=4
on all normalized data for nuclei using only these differential genes to identify cells most l
ikely
corresponding to each cell class. 12’011,
6
386
,
5
’302, and
3
543
nuclei were found
corresponding to GCs, other excitatory, inhibitory, and non
-
neuronal cells, respectively.
Expression of select ion channels across broad classes is shown in Figure
4. Following analysis
of the computational models (see below), gene expression in GCs was analyzed for genes coding
for the following ion channels:
KCNMA1
(BK),
CACNA1b
(Cav2.2) and
KCNJ2
(Kir 2.1). Box
plots showing gene expression for each donor are pr
esented in Figure 4 without a statistical
assessment due to the small sample size for this type of analysis (N=4 donors).
4
Electrophysiological recordings
For electrophysiological recordings,
human brain slices
w
ere
mounted in a custom designed
chamber and held in place with a slice anchor (Warner Instruments). The slice was bathed in
ACSF at 34
±
1
o
C (warmed by an npi hpt
-
2 flow
-
through heater and thermafoil, controlled by an
npi TC
-
20) at a rate of 2 ml per minute (G
ilson Minipuls 3 pump). As a target temperature, 3
4
o
C
was chosen to approximate physiological conditions, but with a safety buffer so as to not exceed
3
7
o
C.
B
ath temperature was continuously monitored. ACSF oxygenation was maintained by
bubbling 95% O2, 5%
CO2 gas in the specimen reservoir as well as delivering across the surface
of the incubation bath via the custom designed specimen chamber. Solution oxygen levels,
temperature and flow rate were monitored; values were measured, documented, and calibrated
weekly to verify system quality.
Thick walled borosilicate glass (Sutter BF150
-
86
-
10) electrodes were manufactured (Sutter
P1000 electrode puller) with a resistance of 3 to 7 MΩ. Prior to recording, electrodes were filled
with 20 μl of
i
nternal
s
olution
S
OP consisting of potassium gluconate, HEPES and other
components
(
43
)
with
b
iocytin, which was thawed fresh every day and kept o
n ice. The pipette
was mounted on a Multiclamp 700B amplifier headstage (Molecular Devices) fixed to a
micromanipulator (PatchStar, Scientifica).
Electrophysiology signals were recorded using an ITC
-
18 Data Acquisition Interface (HEKA).
Commands were gene
rated, signals processed, and amplifier metadata was acquired using a
custom acquisition software program, written in Igor Pro (Wavemetrics). Data were filtered
5
(Bessel) at 10 kHz and digitized at 50
k
Hz. Data were reported uncorrected for the measured
-
14
mV liquid junction potential between the electrode and bath solutions.
Upon break
-
in and formation of a stable seal (typically within the first minute and not more than
3 minutes after break
-
in), the resting membrane potential of the neuron was recorded.
All
recordings were bridge balanced and systematically checked for access resistance matching the
quality control criteria described in
(
46
)
under section “Electrophysiology overview”.
Morphology rec
onstruction and feature extraction
Morphological reconstructions were generated using previously described methods
(
43
)
under
se
ction “Morphology and Histology Overview”.
Briefly, biocytin
-
filled
neurons
were stained via
diaminobenzidine reaction and imaged at 63x magnification on a Zeiss Axio
Imager 2.
Individual cells were digitally reconstructed
with custom
-
written
software (Vaa3D and Mozak)
to create accurate
,
whole
-
neuron representations saved in the SWC format.
A feature extraction
suite was adopted and customized enabling the analysis of
dendritic and somatic reconstructions
and extraction of a number of features related to branching pattern, size, density, soma position,
etc. resulting in 59 morphology features (Table S4).
Dendritic spine densities were assessed using Neurolucida 360
(v2017.01.1)
and Imaris (v9.3).
Initial and terminal 100
μm
sections
of
the
dendritic region outside the granule cell layer were
visualized for manual spine density estimation
along
10
μ
m
-
long segments
per dendrit
ic branch
.
6
Data analysis
We used the
standard Python libraries for analyzing the electrophysiological and morphological
data: Mann
-
Whitney U
-
testing (
mannwhitneyu
from statistics library), random forest
classification (
RandomForestClassifier
with 200 decision trees from scikit
-
learn lirary) a
nd
regression analysis (
linalg.lstsq
from numpy). For dimensionality reduction to perform tSNE we
used the sklearn.manifold package.
Electrophysiology feature extraction and single
-
cell model setup
From whole
-
cell patch
-
clamping experiments, e
lectrophysio
logical responses
to a battery of
standardized current stimuli
(1 s
-
long dc current injections of increasing amplitude) were
analyzed resulting in a set
of
subthreshold and spiking
features for each experiment (Table S2)
.
Electrophysiology features such as spike timing, amplitude, width, etc. were obtained for each
experiment. For a subset of features, regression was used to assess how a particular
electrophysiology feature changes with increasing intracellular stimulation
amplitude. In total, 30
electrophysiology features were extracted from our
in vitro
experiments (Table S2) and analyzed
with a number of statistical techniques for WG
-
dependent differences. Details on the
implementation of the feature extraction analysis a
nd relevant code is publicly available through
(
47
)
.
The
model
optimization procedure
for generating
biophysically realistic
single
-
neuron models is
feature
-
based and
attempts to
set
the somatic, axon initial segment, and dendritic properties in
such manner to capture features of intracellular somatic responses
for
a number of standardized
current stimulation waveforms.
For model generation, electrophysiology
features w
ere extracted
7
using
a f
eature
e
xtraction
l
ibrary (eFEL) developed in
(
48
)
.
Specifically,
11 electrophysiolog
y
features
were
extracted for each experiment
(Fig.
S
4)
and their mean and standard deviation (std)
was
computed for a particular stimulation waveform.
I
f more than
one
sweep
of the same
stimulation waveform exists (majority of the experiments), then the std of that particular
waveform
was
used. If only a single
sweep
of the experiment exist
s
,
a
default value of 10%
was
used for the std.
The compartments of human granule
and basket cells were separated into three zones: the axon
initial segment (AIS), the soma and dendrites.
In our models, t
he AIS was represented by one
fixed
-
length section with a total length of 30 μm and 1 μm diameter.
In all models, passive and active
properties were optimized using the same fitting procedure. For
passive properties, one value for the specific capacitance (
cm
), passive conductance (
g_pas
),
passive reversal potential (
e_pas
), and cytoplasmic resistivity (
Ra
) was uniformly set across all
compartments. Notably, the values of these parameters were part of the genetic optimization
procedure (see below). Active channel mechanisms were spatially uniformly distributed in the
AIS, soma, and dendrites with every zone receiving a separate set of c
hannels (Table S5). The
same approach was used for generating biophysically detailed basket cell models from human
medial temporal gyrus with slice electrophysiology and morphology data originating from
(
46
)
(specimen IDs: 529807751, 541536216; Table S6). The final excitatory and inhibitory models
will be made available on
https://github.com/AllenInstitute/epilepsy_human_dg
.
8
Parameter
optimization and single
-
cell model generatio
n workflow
For the single
-
cell model optimization and generation, the BluePyOpt framework was adopted
(
https://github.com/AllenInstitute/All
-
active
-
Workflow
,
https://github.com/AllenInstitute/All
-
active
-
Manuscript
)
(
48
)
which utilizes multi
-
obj
ective optimization relying on an evolutionary
optimization algorithm. Briefly, for every electrophysiological feature out of the 11 used to build
computational models (Fig.
S
4), an absolute standard score was calculated
Z
i
=|f
i
-
μ
i
|/
s
i
with the
feature valu
e (
f
i
) measured from the output traces of the models and
!
,
!
being the
experimentally measured mean and standard deviation,
respectively
(
48
)
.
The development
of every biophysically detailed single
-
cell model involved a 3
-
stage
optimization workflow that iteratively focused on a set of conductances. In the first stage,
passive model properties like membrane resistance and capacitance were set by fitting
subthreh
old responses and subthreshold features such as
voltage_sag
,
voltage_base
and
voltage_after_stim
. In the second stage, 15 active ionic conductances were included while
keeping the membrane capacitance and passive reversal potential fixed. Notably, the init
ial
parameter ranges for all GC models (both for passive and active ion channels) were chosen
irrespective of WG (Table S5) and optimized for 100 generations. In this stage, ten spiking traces
were used for model training along with subthreshold responses.
In the final stage, the best
model from the previous stage was chosen, its parameter range extended by two
-
fold for all
conductances and optimization was re
-
initialized for at least 50 generations. As part of the
BluePyOpt framework, multiple models were
generated for every cell with the best individual
considered as the one with the smallest sum of objective values. The best model for each
9
experiment was used in the pool of models selected for the network simulations (see below).
Concurrently, for every c
ell, the 10 best models (“hall of fame” models according to
(
25
,
48
)
)
were selected and used for pairwise conductance comparisons between WG1 and WG4 (Fig.
S
6).
For the evolutionary algorithm,
a population size of typically 320 individuals on 320 cores of a
BlueGene/P or Cori Haswell (Cray XC40) supercomputer was used. Generation of a single,
biophysically detailed, single
-
neuron model required approximately 200’000 core
-
hours of
computational t
ime on the aforementioned architectures.
Hippocampal network model setup
A bio
-
plausible DG network model was constructed consisting of 506 single
-
cell models with
t
he
details of the network
relying on
(
18
)
. Briefly, GC and BC cells in our network were
positioned via a topographic approach with somata of the models positioned
in
3
D
space
along
a
ring
with a radius of
15
00
μ
m for
GCs
and 750
μ
m for
BC
s. The ring was open
ed
to form the
arch roughly correspond
i
ng
to
the physical extent of human DG (Fig
.
5A).
S
ynaptic connectivity in the network
depends
on the location of
a
neuron
along
the ring
(
characterized by the
arch
-
angle
)
.
The overall neural network structure was adapted from
(
18
)
.
We have modified the number of connections per synapse to make
it more realistic
and
in
agreement
to
(
49
)
.
The following connectivity structure was chosen to represent the
DG
circuit
with periodic
boundary
conditions
(
compensat
ing
for border
artifacts)
. Every
GC
formed
50
connections
with its
100 closest
GCs
along
the ring
with the
e
xact connecti
vity between the
postsynaptic target pool
chosen randomly. Every
GC
also project
ed
to 3
of
its
closest
BCs
along
10
the ring.
Recurrent GC
-
GC connections formed
2 to 5 syn
apses w
ith
the exact number chosen
randomly. The location of
GC
synapses
was
chosen
along dendrites within 150
μ
m
from the
soma
. All excitatory connections of
GCs
correspond
ed
to AMPA synapse kinetics
(
18
)
.
Every
BC
created
100 out of 140 possible connections with
its topographically
closest
GCs with
synapses located
on
the soma of
GCs
.
Recurrent
B
C connectivity was also implemented
s
o
that
every
BC
is connected to
its two
closest neighbors
with t
he exact
location of inhibitory synapses
chosen randomly
within 150
μ
m
from the BC soma.
All inhibitory connections correspond to
GABA
-
A synapse kinetics
(
18
)
. Synaptic kinetics were approximated using biexponential
synapse models (`exp2syn’ in
(
50
)
) involving two time constants (
!
a
nd
"
) as well as the
reversal potential (
V
rev
). The values used for these parameters are the following (adopted from
(
18
)
):
synapse type
!
/ ms
"
/ ms
V
rev
/ mV
GC
à
GC (AMPA)
1.5
1.5
0.0
GC
à
BC (AMPA)
0.3
0.6
0.0
Perforant path
à
GC
1.5
5.5
0.0
BC
à
BC (GABA
-
A)
0.16
1.8
-
70.0
BC
à
GC (GABA
-
A)
0.26
5.5
-
70.0
Perforant path
à
B
C
2.0
6.3
0.0
To activate the network, GCs
and
BCs
receive
d
background synaptic input
along
their dendrites.
G
Cs
receive
d
5
external
connections per cell with 5 to 15 synapses
each
the exact number of
11
synapses
was
chosen randomly
from a homogenous distribution
. B
Cs
receiv
ed
the same amount
of background synaptic input. Every background connection
was
activated by
a
Poisson
process
(rate: 3 Hz)
emulating
perforant path
input
.
B
ackground synaptic input
corresponded
to AMPA
synapse
kinetics
both for
GC
s and
BC
s
(
18
)
.
All network simulations were
setup and performed
using
inhouse developed
soft
ware based on
(
51
)
.
The versions of the DG network presented in Fig. 5 will be made publicly available through
https://github.com/AllenInstitute/epilepsy_human_dg
12
Supplementary
figure 1
Supplementary
f
igure
1
: Patient
-
specific electrophysiological properties.
(A)
tSNE
5
visualization of 30 electrophysiology features per cell for 112 DG GCs from all subjects (blue:
WG1; red: WG4; 4 WG1, 3 WG4 subjects; Table S1).
(B)
Same as panel (A) with colors
corresponding to different patients.
13
Supplementary
figure 2
Supplementary figure
2
: Location of sampled granule cells within human dentate gyrus.
5
(A)
Distance between the soma of whole
-
cell patch
-
clamped and morphologically reconstructed
DGs from the inner side of DG. Statistical significance testing (Mann
-
Whitn
ey U
-
test,
significance level
=
0.05) shows no difference in soma depth location of sampled neurons across
DG
.
(B)
Number of sampled GCs located in the upper vs. lower DG blade indicate the absence
of a location
-
bias. Data are from all reconstructed cells
(102) for the entire patient cohort (7
10
patients, Table S1).
(C), (D)
Raw data for the depth and soma location (upper/lower blade) of
14
sampled GCs within the GC
-
layer for WG1 and WG4. The separating line indicates the fictitious
landmark separating upper and
lower blade in our analysis.
15
Supplementary
figure 3
Supplementary f
igure
3
. Classification of
both
electrophysiological and morphological
features with degree of sclerosis.
(
A
) tSNE
visualization of electrophysiological and
morphological features (blue: WG1; red: WG4) of GCs with both data modalities present (
77
5
GC
s
out of 112 electrophysiologically recorded cells and 102 morphologically reconstructed had
both data modalities, i.
e
. 6
8%
out of all possible electrophysiology experiments
and 75% out
of
all morphology reconstructions).
(B)
Pairwise comparison of WG1 and WG4 between 79
morpho
-
electric features
(o
range line
:
significance level
p
-
value
= 0.05 / 79 = 0.0006
; stars
indicate fe
atures of statistical signifance)
.
(C)
9 most significant
morpho
-
electric feature
s
10
16
between WG1 and WG4 cells
(order determined per the
p
-
value from the pairwise comparisons
in panel B)
.
(D)
Feature weights of random forest classifier trained on
the same set of
morpho
-
electric features
as in panel B (green boxes:
feature
s
shared
between
pairwise comparison
, panel
B, and
the
random forest
analysis, panel D)
.
5