Minimal parallel ring network

Moderator: wwlytton

Post Reply
pascal
Posts: 86
Joined: Thu Apr 26, 2012 11:51 am

Minimal parallel ring network

Post by pascal » Fri Mar 27, 2020 6:14 pm

I'm trying to implement the code for a minimal parallel ring network, from the 2019 NEURON Course (I could not find the code anywhere online, so I adapted it from the course book).

When I run the program (in serial, for now), for some reason the code to stimulate the 0th cell has no effect. When I plot the voltage traces of all 10 cells, they are identical (just leveling off to their resting membrane potential). I've tried increasing the weight (setting ncstim.weight[0] to large values, such as 1000), but this has no effect. I also tried placing the stimulation code within the 'Network' object, but that didn't work, either. Thanks for the help, and here's the code:

Code: Select all

from neuron import h
from matplotlib import pyplot as plt

h.load_file('stdrun.hoc')

class Pyramidal:
    def __init__(self, gid):
        self._gid = gid
        self._createSections()
        self._defineBiophysics()
        self._register_netcon()
        
    def __repr__(self):
        return "Pyramidal[%d]" % self._gid
    
    def _createSections(self):
        self.soma=h.Section(name='soma', cell=self) 
        self.dend=h.Section(name='dend', cell=self)
        
    def _defineBiophysics(self):
        self.soma.insert('hh')
        self.dend.insert('pas')
        
    def _register_netcon(self):
        self.nc = h.NetCon(self.soma(0.5)._ref_v, None, sec=self.soma)
        pc.set_gid2node(self._gid, idhost)
        pc.cell(self._gid, self.nc)
        self.spike_times = h.Vector()
        self.nc.record(self.spike_times)

class Network:
    def __init__(self,num):
        mygids= list(range(idhost,num,nhost))
        self.cells = [Pyramidal(i) for i in mygids]
        #set up an excitable ExpSyn on each cell's dendrites
        self.syns = [h.ExpSyn(cell.dend(0.5)) for cell in self.cells]
        for syn in self.syns:
            syn.e = 0.0
        # connect cell (i-1)%num to cell i
        self.ncs = []
        for i, syn in zip(mygids,self.syns):
            nc = pc.gid_connect((i-1) % num, syn)
            nc.delay = 1.0
            nc.weight[0] = 1.0
            self.ncs.append(nc)
           
pc = h.ParallelContext()
pc.set_maxstep(10) 
idhost = int(pc.id())
nhost = int(pc.nhost())

n = Network(10)

#drive the 0th cell
if pc.gid_exists(0):
    stim = h.NetStim()
    stim.number = 1
    stim.start = 3.0
    ncstim = h.NetCon(stim,n.syns[0])
    ncstim.delay = 1.0
    ncstim.weight[0] = 1000.0

t = h.Vector().record(h._ref_t)
v = [h.Vector().record(cell.soma(0.5)._ref_v) for cell in n.cells]

h.stdinit()
pc.psolve(500)

for myv in v:
    plt.plot(t,myv)
plt.xlabel('Time (ms)')
plt.ylabel('Voltage (mV)')

pascal
Posts: 86
Joined: Thu Apr 26, 2012 11:51 am

Re: Minimal parallel ring network

Post by pascal » Fri Mar 27, 2020 10:41 pm

Never mind, I found my stupid mistake: I forgot to connect the dendrite to the soma. I was tempted to take this post down to save the embarrassment, but I figured I'd leave it up for posterity (since I couldn't find the code for this Python parallelized ring network anywhere online).

Here's the working code:

Code: Select all

from neuron import h
from matplotlib import pyplot as plt

h.load_file('stdrun.hoc')

class Pyramidal:
    def __init__(self, gid):
        self._gid = gid
        self._createSections()
        self._defineBiophysics()
        self._register_netcon()
        
    def __repr__(self):
        return "Pyramidal[%d]" % self._gid
    
    def _createSections(self):
        self.soma=h.Section(name='soma', cell=self) 
        self.dend=h.Section(name='dend', cell=self)
        
    def _defineBiophysics(self):
        self.soma.insert('hh')
        self.dend.insert('pas')
        self.dend.connect(self.soma(1))
        
    def _register_netcon(self):
        self.nc = h.NetCon(self.soma(0.5)._ref_v, None, sec=self.soma)
        pc.set_gid2node(self._gid, idhost)
        pc.cell(self._gid, self.nc)
        self.spike_times = h.Vector()
        self.nc.record(self.spike_times)

class Network:
    def __init__(self,num):
        mygids= list(range(idhost,num,nhost))
        self.cells = [Pyramidal(i) for i in mygids]
        #set up an excitable ExpSyn on each cell's dendrites
        self.syns = [h.ExpSyn(cell.dend(0.5)) for cell in self.cells]
        for syn in self.syns:
            syn.e = 0.0
        # connect cell (i-1)%num to cell i
        self.ncs = []
        for i, syn in zip(mygids,self.syns):
            nc = pc.gid_connect((i-1) % num, syn)
            nc.delay = 1.0
            nc.weight[0] = 1.0
            self.ncs.append(nc)
           
pc = h.ParallelContext()
pc.set_maxstep(10) 
idhost = int(pc.id())
nhost = int(pc.nhost())

n = Network(10)

#drive the 0th cell
if pc.gid_exists(0):
    stim = h.NetStim()
    stim.number = 1
    stim.start = 3.0
    ncstim = h.NetCon(stim,n.syns[0])
    ncstim.delay = 1.0
    ncstim.weight[0] = 1000.0

t = h.Vector().record(h._ref_t)
v = [h.Vector().record(cell.soma(0.5)._ref_v) for cell in n.cells]

h.stdinit()
pc.psolve(500)

for myv in v:
    plt.plot(t,myv)
plt.xlabel('Time (ms)')
plt.ylabel('Voltage (mV)')

ted
Site Admin
Posts: 5724
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Minimal parallel ring network

Post by ted » Sat Mar 28, 2020 11:05 am

Thank you for posting your question, and for finding the problem and posting the solution so others who read this thread can learn.

Bugs are common in non-trivial software projects. They can be minimized and often prevented by testing at every stage of development. In the context of implementing a computational model, the algorithm is

Code: Select all

given any X that will be used in the model
  make sure that X works properly
This applies to every
voltage- or ligand-gated ion channel
ion transport or accumulation mechanism
synaptic mechanism or gap junction
cell class
and anything else that will be involved.

And even then bugs can still occur--but debugging will be easier because the individual components will have been tested.

Post Reply