Ball-and-stick: 4 - Ring Network Classes

This is the fourth part of a tutorial series where we build a multicompartment cell and evolve it into a network of cells running on a parallel machine. In this part, we take the functionality of the ring network we constructed in the previous page and encapsulate it into various classes so that the network is more extensible. We also begin parameterizing the model so that particular values are not hard-coded, but remain variable so that the model is flexible.

The first thing we will do is pull in our necessary imports. was defined in the previous part of this tutorial.

import numpy
import simrun
from neuron import h, gui
from math import sin, cos, pi
from matplotlib import pyplot

Create a generic cell class

We have evolved our BallAndStick class into a fairly nice class, but it can still be enhanced. For starters, it would be convenient for the cell to know about and maintain a list of its synapses. Additionally, should we ever need to model more than one cell type, it would be nice to specify a template and generic functionality once for all cell classes. Object-oriented capabilities of Python permit this. The following script block defines a generic Cell class that other cells can inherit from. As a template, it mandates that child classes implement methods called in __init__(). It also permits 3D placement, gives cells an ‘all’ SectionList, and also uses a list of synapses and a connection routine.

class Cell:
    """Generic cell template."""
    def __init__(self):
        self.x = 0; self.y = 0; self.z = 0
        self.soma = None
        self.synlist = [] #### NEW CONSTRUCT IN THIS WORKSHEET
    def create_sections(self):
        """Create the sections of the cell. Remember to do this
        in the form::

            h.Section(name='soma', cell=self)

        raise NotImplementedError("create_sections() is not implemented.")
    def build_topology(self):
        """Connect the sections of the cell to build a tree."""
        raise NotImplementedError("build_topology() is not implemented.")
    def define_geometry(self):
        """Set the 3D geometry of the cell."""
        raise NotImplementedError("define_geometry() is not implemented.")
    def define_biophysics(self):
        """Assign the membrane properties across the cell."""
        raise NotImplementedError("define_biophysics() is not implemented.")
    def create_synapses(self):
        """Subclasses should create synapses (such as ExpSyn) at various
        segments and add them to self.synlist."""
        pass # Ignore if child does not implement.
    def build_subsets(self):
        """Build subset lists. This defines 'all', but subclasses may
        want to define others. If overriden, call super() to include 'all'."""
        self.all = h.SectionList()
    def connect2target(self, target, thresh=10):
        """Make a new NetCon with this cell's membrane
        potential at the soma as the source (i.e. the spike detector)
        onto the target passed in (i.e. a synapse on a cell).
        Subclasses may override with other spike detectors."""
        nc = h.NetCon(self.soma(1)._ref_v, target, sec = self.soma)
        nc.threshold = thresh
        return nc
    def is_art(self):
        """Flag to check if we are an integrate-and-fire artificial cell."""
        return False
    def set_position(self, x, y, z):
        Set the base location in 3D and move all other
        parts of the cell relative to that location.
        for sec in self.all:
            for i in range(sec.n3d()):
                        sec.diam3d(i), sec=sec)
        self.x = x; self.y = y; self.z = z
    def rotateZ(self, theta):
        """Rotate the cell about the Z axis."""
        rot_m = numpy.array([[numpy.sin(theta), numpy.cos(theta)],
                [numpy.cos(theta), -numpy.sin(theta)]])
        for sec in self.all:
            for i in range(sec.n3d()):
                xy =[sec.x3d(i), sec.y3d(i)], rot_m)
                h.pt3dchange(i, float(xy[0]), float(xy[1]), sec.z3d(i),

Test this:

class ChildCell(Cell):

cell = ChildCell()

By design, an exception is raised letting us know that we need to at least override create_sections() before proceeding to have a valid subclass object of Cell.

Define BallAndStick as an extension to the base Cell class

This Cell object serves as a template and reduces the coding in inherited objects. Look how much leaner we can make BallAndStick.

class BallAndStick(Cell):  #### Inherits from Cell
    """Two-section cell: A soma with active channels and
    a dendrite with passive properties."""
    #### __init__ is gone and handled in Cell.
    #### We can override __init__ completely, or do some of
    #### our own initialization first, and then let Cell do its
    #### thing, and then do a bit more ourselves with "super".
    #### def __init__(self):
    ####     # Do some stuff
    ####     super(Cell, self).__init__()
    ####     # Do some more stuff
def create_sections(self):
    """Create the sections of the cell."""
    self.soma = h.Section(name='soma', cell=self)
    self.dend = h.Section(name='dend', cell=self)
def build_topology(self):
    """Connect the sections of the cell to build a tree."""
def define_geometry(self):
    """Set the 3D geometry of the cell."""
    self.soma.L = self.soma.diam = 12.6157 # microns
    self.dend.L = 200                      # microns
    self.dend.diam = 1                     # microns
    self.dend.nseg = 5
def define_biophysics(self):
    """Assign the membrane properties across the cell."""
    for sec in self.all: # 'all' exists in parent object.
        sec.Ra = 100     # Axial resistance in Ohm * cm = 1       # Membrane capacitance in micro Farads / cm^2
        # Insert active Hodgkin-Huxley current in the soma
        for seg in self.soma:
            seg.hh.gnabar = 0.12  # Sodium conductance in S/cm2
            seg.hh.gkbar = 0.036  # Potassium conductance in S/cm2
   = 0.0003    # Leak conductance in S/cm2
            seg.hh.el = -54.3     # Reversal potential in mV
        # Insert passive current in the dendrite
        for seg in self.dend:
            seg.pas.g = 0.001  # Passive conductance in S/cm2
            seg.pas.e = -65    # Leak reversal potential mV
def shape_3D(self):
    Set the default shape of the cell in 3D coordinates.
    Set soma(0) to the origin (0,0,0) and dend extending along
    the X-axis.
    len1 = self.soma.L
    h.pt3dadd(0, 0, 0, self.soma.diam, sec=self.soma)
    h.pt3dadd(len1, 0, 0, self.soma.diam, sec=self.soma)
    len2 = self.dend.L
    h.pt3dadd(len1, 0, 0, self.dend.diam, sec=self.dend)
    h.pt3dadd(len1 + len2, 0, 0, self.dend.diam, sec=self.dend)
#### build_subsets, rotateZ, and set_location are now in cell object. ####
#### NEW STUFF ####
def create_synapses(self):
    """Add an exponentially decaying synapse in the middle
    of the dendrite. Set its tau to 2ms, and append this
    synapse to the synlist of the cell."""
    syn = h.ExpSyn(self.dend(0.5))
    syn.tau = 2
    self.synlist.append(syn) # synlist is defined in Cell

Make a Ring class

Encapsulating code into discrete objects is not only conceptually useful for code management, but as we know with cell objects, it lets us make several instances of the object for use in a network. Thinking ahead, we may very well need several networks – each network configured differently. This allows scripting of several simulations en masse, either in a for loop that sequentially processes the networks, or it can be used with NEURON’s subworlds architecture in a parallel context.

class Ring:
    """A network of *N* ball-and-stick cells where cell n makes an
    excitatory synapse onto cell n + 1 and the last, Nth cell in the
    network projects to the first cell.
    def __init__(self, N=5, stim_w=0.004, stim_number=1,
            syn_w=0.01, syn_delay=5):
        :param N: Number of cells.
        :param stim_w: Weight of the stimulus
        :param stim_number: Number of spikes in the stimulus
        :param syn_w: Synaptic weight
        :param syn_delay: Delay of the synapse
        self._N = N              # Total number of cells in the net
        self.cells = []          # Cells in the net
        self.nclist = []         # NetCon list
        self.stim = None         # Stimulator
        self.stim_w = stim_w     # Weight of stim
        self.stim_number = stim_number  # Number of stim spikes
        self.syn_w = syn_w       # Synaptic weight
        self.syn_delay = syn_delay  # Synaptic delay
        self.t_vec = h.Vector()   # Spike time of all cells
        self.id_vec = h.Vector()  # Ids of spike times
        self.set_numcells(N)  # Actually build the net.
    def set_numcells(self, N, radius=50):
        """Create, layout, and connect N cells."""
        self._N = N
    def create_cells(self, N):
        """Create and layout N cells in the network."""
        self.cells = []
        r = 50 # Radius of cell locations from origin (0,0,0) in microns
        N = self._N
        for i in range(N):
            cell = BallAndStick()
            # When cells are created, the soma location is at (0,0,0) and
            # the dendrite extends along the X-axis.
            # First, at the origin, rotate about Z.
            cell.rotateZ(i * 2 * pi / N)
            # Then reposition
            x_loc = cos(i * 2 * pi / N) * r
            y_loc = sin(i * 2 * pi / N) * r
            cell.set_position(x_loc, y_loc, 0)
    def connect_cells(self):
        """Connect cell n to cell n + 1."""
        self.nclist = []
        self.spike_times = []
        N = self._N
        for i in range(N):
            src = self.cells[i]
            tgt_syn = self.cells[(i+1)%N].synlist[0]
            nc = src.connect2target(tgt_syn)
            nc.weight[0] = self.syn_w
            nc.delay = self.syn_delay
            spike_times = h.Vector()
    def connect_stim(self):
        """Connect a spiking generator to the first cell to get
        the network going."""
        self.stim = h.NetStim()
        self.stim.number = self.stim_number
        self.stim.start = 9
        self.ncstim = h.NetCon(self.stim, self.cells[0].synlist[0])
        self.ncstim.delay = 1
        self.ncstim.weight[0] = self.stim_w # NetCon weight is a vector.

Test the network

Let’s make a ring object, render it, and run a simulation.

ring = Ring()
shape_window = h.PlotShape()
shape_window.exec_menu('Show Diam')
soma_v_vec, dend_v_vec, t_vec = simrun.set_recording_vectors(ring.cells[0])
simrun.show_output(soma_v_vec, dend_v_vec, t_vec)

Let’s see a spike plot.

spikes = ring.spike_times
for i, spike_times in enumerate(spikes):
    pyplot.vlines(spike_times, i + 0.5, i + 1.5)

Run a few networks

Let’s run other instances of the net. The code below keeps the variable spikes from our default run, but replaces the net with a new instance with arguments that we pass in, drawing a second set of spikes in red.

ring = Ring(syn_w=.005) # Try different weights, for example.
spikes2 = ring.spike_times
for i, spike_times in enumerate(spikes):
    pyplot.vlines(spike_times, i + 0.5, i + 1.5, color='black')
for i, spike_times in enumerate(spikes2):
    pyplot.vlines(spike_times, i + 0.5, i + 1.5, color='red')

In both simulations, the first spike occurs at 12.625 ms. After that, the red spikes lag the black ones by steadily increasing amounts.

This concludes this part of the tutorial. The next part translates this serial implementation into a parallel model.

Here we specify different labels for the two sets of spikes because if a label was reused (or both were omitted), then the first raster would be replaced with the second.