A single cell model from the Allen Brain Atlas#
In this tutorial we’ll see how we can take a model from the Allen Brain Atlas and run it with Arbor.
Note
Concepts covered in this example:
Take a model from the Allen Brain Atlas.
Load a morphology from an
swc
file.Load a parameter fit file and apply it to a
arbor.decor
.Building a
arbor.cable_cell
representative of the cell in the model.Building a
arbor.recipe
reflective of the cell stimulation in the model.Running a simulation and visualising the results.
Obtaining the model#
We need a model for which morphological data is available. We’ll take
a cell from the mouse’s visual cortex.
As the README
in the “Biophysical - all active” model (zip) file explains, the data set is created for use with the Allen SDK and the
Neuron simulator. Instructions about running the model with the
Allen SDK and Neuron can be found here.
For your convenience, this tutorial comes with a reference trace pre-generated.
In the “Biophysical - all active” model (zip) file you’ll find:
a directory with mechanisms (modfiles). For your convenience, Arbor already comes with the
arbor.allen_catalogue()
, containing mechanisms used by the Allen Brain Atlas.an
swc
file describing the morphology of the cell.fit_parameters.json
describing a set of optimized mechanism parameters. This information will need to be loaded into the Arbor simulation.
We will replicate the “Sweep 35” experiment, which applies a current of 150 nA for a duration of 1 s.
Morphology and Labels#
morphology = A.load_swc_neuron(base / swc).morphology
Step (1) loads the swc
file using arbor.load_swc_neuron()
. Since
the swc
specification is informal, a few different interpretations exist,
and we use the appropriate one. The interpretations are described here.
# In addition, label the midpoint of the soma.
labels = A.label_dict().add_swc_tags()
labels["midpoint"] = "(location 0 0.5)"
Step (2) sets the labels to the defaults of the swc
specification,
plus a label for the midpoint of the soma. (You can verify in the swc
file,
the first branch is the soma.)
The parameter fit#
The most complicated part is transferring the values for the appropriate
parameters in parameter fit file to an arbor.decor
. The file file is a
json
file, which is fortunate; Python comes with a json
package in its
standard library. The passive and conditions block contains cell-wide
defaults, while the genome section contains the parameters for all the
mechanism properties. In certain cases, parameters names include the mechanism
name, so some processing needs to take place.
Step (3) shows the precise steps needed to load the fit parameter file into
a list of global properties, region specific properties, reversal potentials,
and mechanism parameters. This is not a generic function that will successfully
load any Allen model, but it can be used as a starting point. The function
composes 4 components out of the json
file:
global electro-physiological parameters
a set of electro-physiological parameters per region,
a set of reversal potentials per ion species and region,
a set of mechanisms with parameters per region.
def load_allen_fit(fit):
with open(fit) as fd:
fit = json.load(fd)
# cable parameters convenience class
@dataclass
class parameters:
cm: Optional[U.quantity] = None
temp: Optional[U.quantity] = None
Vm: Optional[U.quantity] = None
rL: Optional[U.quantity] = None
param = defaultdict(parameters)
mechs = defaultdict(dict)
for block in fit["genome"]:
mech = block["mechanism"] or "pas"
region = block["section"]
name = block["name"]
value = float(block["value"])
if name.endswith("_" + mech):
name = name[: -(len(mech) + 1)]
elif mech == "pas":
# transform names and values
if name == "cm":
param[region].cm = value * U.uF / U.cm2
elif name == "Ra":
param[region].rL = value * U.Ohm * U.cm
elif name == "Vm":
param[region].Vm = value * U.mV
elif name == "celsius":
param[region].temp = value * U.Celsius
else:
raise Exception(f"Unknown key: {name}")
continue
else:
raise Exception(f"Illegal combination {mech} {name}")
mechs[(region, mech)][name] = value
regs = [(r, vs) for r, vs in param.items()]
mechs = [(r, m, vs) for (r, m), vs in mechs.items()]
default = parameters(
temp=float(fit["conditions"][0]["celsius"]) * U.Celsius,
Vm=float(fit["conditions"][0]["v_init"]) * U.mV,
rL=float(fit["passive"][0]["ra"]) * U.Ohm * U.cm,
)
ions = []
for kv in fit["conditions"][0]["erev"]:
region = kv["section"]
for k, v in kv.items():
if k == "section":
continue
ion = k[1:]
ions.append((region, ion, float(v) * U.mV))
return default, regs, ions, mechs, fit["fitting"][0]["junction_potential"]
The decor#
With the ingredients for the arbor.decor
extracted, we continue with
the function that will return the cable cell from the model as an
arbor.cable_cell
.
# (1) Load the swc file passed into this function
morphology = A.load_swc_neuron(base / swc).morphology
# (2) Label the region tags found in the swc with the names used in the parameter fit file.
# In addition, label the midpoint of the soma.
labels = A.label_dict().add_swc_tags()
labels["midpoint"] = "(location 0 0.5)"
# (3) A function that parses the Allen parameter fit file into components
dflt, regions, ions, mechanisms, offset = load_allen_fit(base / fit)
# (4) Instantiate an empty decor.
decor = A.decor()
# (5) assign global electro-physiology parameters
decor.set_property(
tempK=dflt.temp,
Vm=dflt.Vm,
cm=dflt.cm,
rL=dflt.rL,
)
# (6) override regional electro-physiology parameters
for region, vs in regions:
decor.paint(
f'"{region}"',
tempK=vs.temp,
Vm=vs.Vm,
cm=vs.cm,
rL=vs.rL,
)
# (7) set reversal potentials
for region, ion, e in ions:
decor.paint(f'"{region}"', ion=ion, rev_pot=e)
decor.set_ion("ca", int_con=5e-5 * U.mM, ext_con=2.0 * U.mM, method="nernst/x=ca")
# (8) assign ion dynamics
for region, mech, values in mechanisms:
nm = mech
vs = {}
sp = "/"
for k, v in values.items():
if mech == "pas" and k == "e":
nm = f"{nm}{sp}{k}={v}"
sp = ","
else:
vs[k] = v
decor.paint(f'"{region}"', A.density(A.mechanism(nm, vs)))
# (9) attach stimulus and detector
decor.place('"midpoint"', A.iclamp(0.2 * U.s, 1 * U.s, 150 * U.pA), "ic")
decor.place('"midpoint"', A.threshold_detector(-40 * U.mV), "sd")
# (10) discretisation strategy: max compartment length
decor.discretization(A.cv_policy_max_extent(20))
# (11) Create cell
return A.cable_cell(morphology, decor, labels), offset
Step (4) creates an empty arbor.decor
.
Step (5) assigns global (cell-wide) properties using
arbor.decor.set_property()
. In addition, initial internal and external
calcium concentrations are set, and configured to be determined by the Nernst
equation.
Note
- Setting the calcium reversal potential to be determined by the Nernst
equation has to be done manually, in order to mirror an implicit Neuron behavior, for which the fit parameters were obtained. This behavior can be stated as the following rule:
If the internal or external concentration of an ion is written, and its reversal potential is read but not written, then the Nernst equation is used continuously during the simulation to update the reversal potential of the ion according to the Nernst equation
Step (6) overrides the global properties for all regions for which the fit
parameters file specifies adapted values. Regional properties are painted
, and are painted over (e.g. replacing) the defaults.
Step (7) sets the regional reversal potentials.
Step (8) assigns the regional mechanisms.
Now that the electro-physiology is all set up, let’s move on to the experimental setup.
Step (9) configures the stimulus
of 150 nA for a
duration of 1 s, starting after 200 ms of the start of the simulation. We’ll
also install a arbor.threshold_detector
that triggers at -40 mV. (The
location is usually the soma, as is confirmed by coordinates found in the
experimental dataset at 488683423.nwb/general/intracellular_ephys/Electrode
1/location
)
Step (10) specifies a maximum control volume length of 20 μm.
Step (11) returns the arbor.cable_cell
.
The model#
cell, offset = make_cell(here, "single_cell_allen.swc", "single_cell_allen_fit.json")
model = A.single_cell_model(cell)
# (13) Set the probe
model.probe("voltage", '"midpoint"', "Um", frequency=1 / (5 * U.us))
# (14) Install the Allen mechanism catalogue.
model.properties.catalogue.extend(A.allen_catalogue(), "")
# (15) Run simulation
model.run(tfinal=1.4 * U.s, dt=5 * U.us)
Step (12) instantiates the arbor.cable_cell
and an
arbor.single_cell_model
.
Step (13) shows how to install a probe to the "midpoint"
, with a
sampling frequency of 200 kHz.
Step (14) installs the arbor.allen_catalogue
, thereby making its
mechanisms available to the definitions added to the decor.
Step (15) starts the simulation for a duration of 1.4 s and a timestep of 5 ms.
The result#
Let’s look at the result! In step (16) we first load the reference generated
with Neuron and the AllenSDK. Then, we extract Arbor’s output, accessible after
the simulation ran at arbor.single_cell_model.traces
. Then, we plot
them, together with the arbor.single_cell_model.spikes
in step
(17).
reference = (
1e3 * pd.read_csv(here / "single_cell_allen_neuron_ref.csv")["U/mV"][:-1] + offset
)
# (17) Plot
df = pd.concat(
[
pd.DataFrame(
{
"Simulator": "Arbor",
"t/ms": model.traces[0].time,
"U/mV": model.traces[0].value,
}
),
pd.DataFrame(
{
"Simulator": "Neuron",
"t/ms": model.traces[0].time,
"U/mV": reference.values,
}
),
],
ignore_index=True,
)
sns.relplot(data=df, kind="line", x="t/ms", y="U/mV", hue="Simulator", errorbar=None)
plt.scatter(
model.spikes, [-40] * len(model.spikes), color=sns.color_palette()[2], zorder=20
)
plt.bar(
200,
max(reference) - min(reference),
1000,
min(reference),
align="edge",
label="Stimulus",
color="0.9",
)
plt.savefig("single_cell_allen_result.pdf")
Note
The careful observer notices that this trace does not match the experimental
data shown on the Allen website (or in the 488683423.nwb
file). Sweep 35
clearly has 5 spikes, not 4. That is because in the Allen SDK, the axon in the
swc
file is replaced with a stub, see this paper and this
AllenSDK Issue.
However, that adapted morphology is not exportable back to a modified swc
file. When we tried to mimic the procedure, we did not obtain the experimental
trace.
Therefore, we used the unmodified morphology in Arbor and the Neuron reference (by commenting out the changes the Allen SDK makes to the morphology) in order to make a 1:1 comparison possible.
The full code#
You can find the source code for this example in full at python/examples/single_cell_allen.py
.