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:

  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.

The morphology

def make_cell(swc, fit):
    # (1) Load the swc file passed into this function
    morphology = arbor.load_swc_neuron(swc)
    # (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 somarbor.
    labels = arbor.label_dict().add_swc_tags()
    labels["midpoint"] = "(location 0 0.5)"

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.

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.

# (3) A function that parses the Allen parameter fit file into components for an arbor.decor
# NB. Needs to be adjusted when using a different model
def load_allen_fit(fit):
    with open(fit) as fd:
        fit = json.load(fd)

    # cable parameters convenience class
    @dataclass
    class parameters:
        cm: float = None
        tempK: float = None
        Vm: float = None
        rL: float = 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":
                # scaling factor NEURON -> Arbor
                param[region].cm = value / 100.0
            elif name == "Ra":
                param[region].rL = value
            elif name == "Vm":
                param[region].Vm = value
            elif name == "celsius":
                param[region].tempK = value + 273.15
            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(
        tempK=float(fit["conditions"][0]["celsius"]) + 273.15,
        Vm=float(fit["conditions"][0]["v_init"]),
        rL=float(fit["passive"][0]["ra"]),
    )

    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)))

    return default, regs, ions, mechs, fit["fitting"][0]["junction_potential"]

    # (3) A function that parses the Allen parameter fit file into components for an arbor.decor
    dflt, regions, ions, mechanisms, offset = load_allen_fit(fit)

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.

# (4) Instantiate an empty decor.
decor = arbor.decor()
# (5) assign global electro-physiology parameters
decor.set_property(tempK=dflt.tempK, 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.tempK, Vm=vs.Vm, cm=vs.cm, rL=vs.rL)
# (7) set reversal potentials
for region, ion, e in ions:
    decor.paint(f'"{region}"', ion_name=ion, rev_pot=e)
decor.set_ion(
    "ca", int_con=5e-5, ext_con=2.0, method=arbor.mechanism("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}"', arbor.density(arbor.mechanism(nm, vs)))
# (9) attach stimulus and detector
decor.place('"midpoint"', arbor.iclamp(200, 1000, 0.15), "ic")
decor.place('"midpoint"', arbor.threshold_detector(-40), "sd")
# (10) discretisation strategy: max compartment length
decor.discretization(arbor.cv_policy_max_extent(20))

# (11) Create cell
return arbor.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

# (12) Create cell, model
cell, offset = make_cell("single_cell_allen.swc", "single_cell_allen_fit.json")
model = arbor.single_cell_model(cell)

# (13) Set the probe
model.probe("voltage", '"midpoint"', frequency=200)

# (14) Install the Allen mechanism catalogue.
model.properties.catalogue.extend(arbor.allen_catalogue(), "")

# (15) Run simulation
model.run(tfinal=1400, dt=0.005)

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).

# (16) Load and scale reference
reference = (
    1000.0 * pandas.read_csv("single_cell_allen_neuron_ref.csv")["U/mV"].values[:-1]
    + offset
)

# (17) Plot
df_list = []
df_list.append(
    pandas.DataFrame(
        {
            "t/ms": model.traces[0].time,
            "U/mV": model.traces[0].value,
            "Simulator": "Arbor",
        }
    )
)
df_list.append(
    pandas.DataFrame(
        {"t/ms": model.traces[0].time, "U/mV": reference, "Simulator": "Neuron"}
    )
)
df = pandas.concat(df_list, ignore_index=True)
seaborn.relplot(
    data=df, kind="line", x="t/ms", y="U/mV", hue="Simulator", errorbar=None
)
plt.scatter(
    model.spikes, [-40] * len(model.spikes), color=seaborn.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.svg")
../_images/single_cell_allen_result.svg

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.

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 here and here. 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.