A simple dendrite

In this example, we will set up a single dendrite represented by a passive cable. On one end the dendrite is stimulated by a short current pulse. We will investigate how this pulse travels along the dendrite and calculate the conduction velocity.

Note

Concepts covered in this example:

  1. Creating a simulation recipe of a single dendrite.

  2. Placing probes on the morphology.

  3. Running the simulation and extracting the results.

  4. Investigating the influence of control volume policies.

The full code

You can find the full code of the example at python/examples/single_cell_cable.py

Executing the script will run the simulation with default parameters.

Walk-through

We set up a recipe for the simulation of a single dendrite with some parameters.

class Cable(arbor.recipe):
    def __init__(self, probes,
                 Vm, length, radius, cm, rL, g,
                 stimulus_start, stimulus_duration, stimulus_amplitude,
                 cv_policy_max_extent):
        """
        probes -- list of probes

        Vm -- membrane leak potential
        length -- length of cable in μm
        radius -- radius of cable in μm
        cm -- membrane capacitance in F/m^2
        rL -- axial resistivity in Ω·cm
        g -- membrane conductivity in S/cm^2

        stimulus_start -- start of stimulus in ms
        stimulus_duration -- duration of stimulus in ms
        stimulus_amplitude -- amplitude of stimulus in nA

        cv_policy_max_extent -- maximum extent of control volume in μm
        """

Implementing the cell_description member function constructs the morphology and sets the properties of the dendrite as well as the current stimulus and the discretization policy.

    def cell_kind(self, gid):
        assert gid == 0
        return arbor.cell_kind.cable

    def cell_description(self, gid):
        """A high level description of the cell with global identifier gid.

        For example the morphology, synapses and ion channels required
        to build a multi-compartment neuron.
        """
        assert gid == 0

        tree = arbor.segment_tree()

        tree.append(arbor.mnpos,
                    arbor.mpoint(0, 0, 0, self.radius),
                    arbor.mpoint(self.length, 0, 0, self.radius),
                    tag=1)

        labels = arbor.label_dict({'cable': '(tag 1)',
                                   'start': '(location 0 0)'})

        decor = arbor.decor()
        decor.set_property(Vm=self.Vm)
        decor.set_property(cm=self.cm)
        decor.set_property(rL=self.rL)

        decor.paint('"cable"',
                    arbor.mechanism(f'pas/e={self.Vm}', {'g': self.g}))

        decor.place('"start"', arbor.iclamp(self.stimulus_start, self.stimulus_duration, self.stimulus_amplitude), "iclamp")

        policy = arbor.cv_policy_max_extent(self.cv_policy_max_extent)
        decor.discretization(policy)

        return arbor.cable_cell(tree, labels, decor)

We parse the command line arguments, instantiate the recipe, run the simulation, extract results and plot:

if __name__ == "__main__":

    parser = argparse.ArgumentParser(description='Cable')

    parser.add_argument(
        '--Vm', help="membrane leak potential in mV", type=float, default=-65)
    parser.add_argument(
        '--length', help="cable length in μm", type=float, default=1000)
    parser.add_argument(
        '--radius', help="cable radius in μm", type=float, default=1)
    parser.add_argument(
        '--cm', help="membrane capacitance in F/m^2", type=float, default=0.01)
    parser.add_argument(
        '--rL', help="axial resistivity in Ω·cm", type=float, default=90)
    parser.add_argument(
        '--g', help="membrane conductivity in S/cm^2", type=float, default=0.001)

    parser.add_argument('--stimulus_start',
                        help="start of stimulus in ms", type=float, default=10)
    parser.add_argument('--stimulus_duration',
                        help="duration of stimulus in ms", type=float, default=0.1)
    parser.add_argument('--stimulus_amplitude',
                        help="amplitude of stimulus in nA", type=float, default=1)

    parser.add_argument('--cv_policy_max_extent',
                        help="maximum extent of control volume in μm", type=float,
                        default=10)

    # parse the command line arguments
    args = parser.parse_args()

    # set up membrane voltage probes equidistantly along the dendrites
    probes = [arbor.cable_probe_membrane_voltage(
        f'(location 0 {r})') for r in np.linspace(0, 1, 11)]
    recipe = Cable(probes, **vars(args))

    # create a default execution context and a default domain decomposition
    context = arbor.context()
    domains = arbor.partition_load_balance(recipe, context)

    # configure the simulation and handles for the probes
    sim = arbor.simulation(recipe, domains, context)
    dt = 0.001
    handles = [sim.sample((0, i), arbor.regular_schedule(dt))
               for i in range(len(probes))]

    # run the simulation for 30 ms
    sim.run(tfinal=30, dt=dt)

    # retrieve the sampled membrane voltages and convert to a pandas DataFrame
    print("Plotting results ...")
    df_list = []
    for probe in range(len(handles)):
        samples, meta = sim.samples(handles[probe])[0]
        df_list.append(pandas.DataFrame({'t/ms': samples[:, 0], 'U/mV': samples[:, 1], 'Probe': f"{probe}"}))

    df = pandas.concat(df_list,ignore_index=True)
    seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Probe",ci=None).set(xlim=(9,14)).savefig('single_cell_cable_result.svg')

The output plot below shows how the current pulse is attenuated and gets broader the further it travels along the dendrite from the current clamp.

../_images/single_cell_cable_result.svg

The conduction velocity in simulation is calculated from the time of the maximum of the membrane potential.

    # calculcate the idealized conduction velocity, cf. cable equation
    data = [sim.samples(handle)[0][0] for handle in handles]
    rm = get_rm(args.g*1/(0.01*0.01))
    taum = get_taum(args.cm, rm)
    lambdam = get_lambdam(args.radius*1e-6, rm, args.rL*0.01)
    vcond_ideal = get_vcond(lambdam, taum)

    # take the last and second probe
    delta_t = (get_tmax(data[-1]) - get_tmax(data[1]))

    # 90% because we took the second probe
    delta_x = args.length*0.9
    vcond = delta_x*1e-6/(delta_t*1e-3)

    print(f"calculated conduction velocity: {vcond_ideal:.2f} m/s")
    print(f"simulated conduction velocity:  {vcond:.2f} m/s")

Please have in mind that the calculated (idealized) conduction velocity is only correct for an infinite cable.

calculated conduction velocity: 0.47 m/s
simulated conduction velocity:  0.50 m/s

When we set the control volume policy to cover the full dendrite, all probes will see the same voltage.

./single_cell_cable.py --length 1000 --cv_policy_max_extent 1000

Output:

../_images/single_cell_cable_result_cv_policy.svg