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.


Concepts covered in this example:

  1. Take a model from the Allen Brain Atlas.

  2. Load a morphology from an swc file.

  3. Load a parameter fit file and apply it to a arbor.decor.

  4. Building a arbor.cable_cell representative of the cell in the model.

  5. Building a arbor.recipe reflective of the cell stimulation in the model.

  6. 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:

  1. global electro-physiological parameters

  2. a set of electro-physiological parameters per region,

  3. a set of reversal potentials per ion species and region,

  4. 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
    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
                raise Exception(f"Unknown key: {name}")
            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":
            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

# (6) override regional electro-physiology parameters
for region, vs in regions:

# (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 = ","
            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

# (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.


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(
                "Simulator": "Arbor",
                "t/ms": model.traces[0].time,
                "U/mV": model.traces[0].value,
                "Simulator": "Neuron",
                "t/ms": model.traces[0].time,
                "U/mV": reference.values,

sns.relplot(data=df, kind="line", x="t/ms", y="U/mV", hue="Simulator", errorbar=None)
    model.spikes, [-40] * len(model.spikes), color=sns.color_palette()[2], zorder=20
    max(reference) - min(reference),

Plot of experiment 35 of the Allen model, compared to the reference generated by the AllenSDK. In green: the threshold detector output; in shaded grey: the stimulus.#


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.