Source code for compartment_net_no_array

# compartment_net_no_array.py ---
#
# Filename: compartment_net_no_array.py
# Description:
# Author:Subhasis Ray
# Maintainer:
# Created: Sat Aug 11 14:30:21 2012 (+0530)
# Version:
# Last-Updated: Tue May  7 18:26:26 2013 (+0530)
#           By: subha
#     Update #: 974
# URL:
# Keywords:
# Compatibility:
#
#

# Commentary:
#
#
#

# Change log:
#
#
#
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING.  If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
# Floor, Boston, MA 02110-1301, USA.
#
#

# Code:
"""
Following, is a demo to create a network of single compartmental neurons connected
via alpha synapses. This is same as compartment_net.py except that
we avoid ematrix and use single melements.
"""

import sys
sys.path.append('../../python')
import os
os.environ['NUMPTHREADS'] = '1'
import numpy as np
import matplotlib.pyplot as plt

import moose
from moose import utils

EREST_ACT = -70e-3

# Gate equations have the form:
#
# y(x) = (A + B * x) / (C + exp((x + D) / F))
#
# where x is membrane voltage and y is the rate constant for gate
# closing or opening

Na_m_params = [1e5 * (25e-3 + EREST_ACT),   # 'A_A':
                -1e5,                       # 'A_B':
                -1.0,                       # 'A_C':
                -25e-3 - EREST_ACT,         # 'A_D':
               -10e-3,                      # 'A_F':
                4e3,                     # 'B_A':
                0.0,                        # 'B_B':
                0.0,                        # 'B_C':
                0.0 - EREST_ACT,            # 'B_D':
                18e-3                       # 'B_F':
               ]
Na_h_params = [ 70.0,                        # 'A_A':
                0.0,                       # 'A_B':
                0.0,                       # 'A_C':
                0.0 - EREST_ACT,           # 'A_D':
                0.02,                     # 'A_F':
                1000.0,                       # 'B_A':
                0.0,                       # 'B_B':
                1.0,                       # 'B_C':
                -30e-3 - EREST_ACT,        # 'B_D':
                -0.01                    # 'B_F':
                ]
K_n_params = [ 1e4 * (10e-3 + EREST_ACT),   #  'A_A':
               -1e4,                      #  'A_B':
               -1.0,                       #  'A_C':
               -10e-3 - EREST_ACT,         #  'A_D':
               -10e-3,                     #  'A_F':
               0.125e3,                   #  'B_A':
               0.0,                        #  'B_B':
               0.0,                        #  'B_C':
               0.0 - EREST_ACT,            #  'B_D':
               80e-3                       #  'B_F':
               ]
VMIN = -40e-3 + EREST_ACT
VMAX = 120e-3 + EREST_ACT
VDIVS = 30000

soma_dia = 30e-6

def create_na_chan(path):
    na = moose.HHChannel('%s/na' % (path))
    na.Xpower = 3
    xGate = moose.HHGate(na.path + '/gateX')
    xGate.setupAlpha(Na_m_params +
                      [VDIVS, VMIN, VMAX])
    na.Ypower = 1
    yGate = moose.HHGate(na.path + '/gateY')
    yGate.setupAlpha(Na_h_params +
                      [VDIVS, VMIN, VMAX])
    na.Ek = 115e-3 + EREST_ACT
    return na

def create_k_chan(path):
    k = moose.HHChannel('%s/k' % (path))
    k.Xpower = 4.0
    xGate = moose.HHGate(k.path + '/gateX')
    xGate.setupAlpha(K_n_params +
                      [VDIVS, VMIN, VMAX])
    k.Ek = -12e-3 + EREST_ACT
    return k

def create_compartment(path):
    comp = moose.Compartment(path)
    comp.diameter = soma_dia
    comp.Em = EREST_ACT + 10.613e-3
    comp.initVm = EREST_ACT
    sarea = np.pi * soma_dia * soma_dia
    comp.Rm = 1/(0.3e-3 * 1e4 * sarea)
    comp.Cm = 1e-6 * 1e4 * sarea
    if moose.exists('/library/na'):
        nachan = moose.element(moose.copy('/library/na', comp, 'na'))
    else:
        nachan = create_na_chan(comp.path)
    nachan.Gbar = 120e-3 * sarea * 1e4
    moose.showfield(nachan)
    moose.connect(nachan, 'channel', comp, 'channel')
    if moose.exists('/library/k'):
        kchan = moose.element(moose.copy('/library/k', comp, 'k'))
    else:
        kchan = create_k_chan(comp.path)
    kchan.Gbar = 36e-3 * sarea * 1e4
    moose.connect(kchan, 'channel', comp, 'channel')
    synchan = moose.SynChan(comp.path + '/synchan')
    synchan.Gbar = 1e-8
    synchan.tau1 = 2e-3
    synchan.tau2 = 2e-3
    synchan.Ek = 0.0
    m = moose.connect(comp, 'channel', synchan, 'channel')
    spikegen = moose.SpikeGen(comp.path + '/spikegen')
    spikegen.threshold = 0.0
    m = moose.connect(comp, 'VmOut', spikegen, 'Vm')
    return comp

def test_compartment():
    n = moose.Neutral('/model')
    lib = moose.Neutral('/library')
    create_na_chan(lib.path)
    create_k_chan(lib.path)
    comp = create_compartment('/model/soma')
    pg = moose.PulseGen('/model/pulse')
    pg.firstDelay = 50e-3
    pg.firstWidth = 40e-3
    pg.firstLevel = 1e-9
    moose.connect(pg, 'output', comp, 'injectMsg')
    d = moose.Neutral('/data')
    vm = moose.Table('/data/Vm')
    moose.connect(vm, 'requestOut', comp, 'getVm')
    gK = moose.Table('/data/gK')
    moose.connect(gK, 'requestOut', moose.element('%s/k' % (comp.path)), 'getGk')
    gNa = moose.Table('/data/gNa')
    moose.connect(gNa, 'requestOut', moose.element('%s/na' % (comp.path)), 'getGk')
    # utils.resetSim(['/model', '/data'], 1e-6, 1e-4, simmethod='ee')
    assign_clocks(['/model'], 1e-6, 1e-4)
    simtime = 100e-3
    moose.start(simtime)
    t = np.linspace(0, simtime, len(vm.vector))
    plt.subplot(221)
    plt.title('Vm')
    plt.plot(t, vm.vector)
    plt.subplot(222)
    plt.title('Conductance')
    plt.plot(t, gK.vector, label='GK')
    plt.plot(t, gNa.vector, label='GNa')
    plt.legend()
    plt.subplot(223)
    ma = moose.element('%s/na/gateX' % (comp.path)).tableA
    mb = moose.element('%s/na/gateX' % (comp.path)).tableB
    ha = moose.element('%s/na/gateY' % (comp.path)).tableA
    hb = moose.element('%s/na/gateY' % (comp.path)).tableB
    na = moose.element('%s/k/gateX' % (comp.path)).tableA
    nb = moose.element('%s/k/gateX' % (comp.path)).tableB
    plt.plot(1/mb, label='tau_m')
    plt.plot(1/hb, label='tau_h')
    plt.plot(1/nb, label='tau_n')
    plt.legend()
    plt.subplot(224)
    plt.plot(ma/mb, label='m_inf')
    plt.plot(ha/hb, label='h_inf')
    plt.plot(na/nb, label='n_inf')
    plt.legend()
    plt.show()
    plt.close()

[docs]def create_population(container, size): """Create a population of `size` single compartmental neurons with Na and K channels. Also create SpikeGen objects and SynChan objects connected to these which can act as plug points for setting up synapses later.""" path = container.path # Contrast this with # comps = moose.vec(path+'/soma', size, 'Compartment') comps = [create_compartment(path+'/soma_%d' % (ii)) for ii in range(size)] spikegens = [] synchans = [] Em = EREST_ACT + 10.613e-3 initVm_array = [EREST_ACT] * size Em_array = [Em] * size # initVm_array = np.random.normal(EREST_ACT, np.abs(EREST_ACT) * 0.1, size) # Em_array = np.random.normal(Em, np.abs(Em) * 0.1, size) for comp, initVm, Em in zip(comps, initVm_array, Em_array): comp.Em = Em comp.initVm = initVm synchan = moose.element(comp.path + '/synchan') synchans.append(synchan) spikegen = moose.element(comp.path + '/spikegen') spikegens.append(spikegen) return {'compartment': comps, 'spikegen': spikegens, 'synchan': synchans}
[docs]def make_synapses(spikegen, synchan, delay=5e-3): """ Create synapses from spikegens to synchans in a manner similar to OneToAll connection. spikegen: list of spikegen objects These are sources of synaptic event messages. synchan: list of synchan objects These are the targets of the synaptic event messages. delay: mean delay of synaptic transmission. Individual delays are normally distributed with sd=0.1*mean. """ scount = len(spikegen) for ii, sid in enumerate(synchan): s = moose.SynChan(sid) sh = moose.SimpleSynHandler( sid.path + "/synh" ) moose.connect( sh, "activationOut", s, "activation" ) sh.synapse.num = scount delay_list = np.random.normal(delay, delay*0.1, scount) # print delay_list for jj in range(scount): sh.synapse[jj].delay = delay_list[jj] # Connect all spikegens to this synchan except that from # same compartment - we assume if parents are same the two belong to the same compartment if s.parent.path != spikegen[jj].parent.path: m = moose.connect(spikegen[jj], 'spikeOut', moose.element(sh.path + '/synapse'), 'addSpike')
[docs]def two_populations(size=2): """An example with two population connected via synapses.""" net = moose.Neutral('network2') pop_a = create_population(moose.Neutral('/network2/pop_A'), size) pop_b = create_population(moose.Neutral('/network2/pop_B'), size) make_synapses(pop_a['spikegen'], pop_b['synchan']) make_synapses(pop_b['spikegen'], pop_a['synchan']) pulse = moose.PulseGen('/network2/net2_pulse') pulse.firstLevel = 1e-9 pulse.firstDelay = 0.05 # disable the pulsegen pulse.firstWidth = 0.02 moose.connect(pulse, 'output', pop_a['compartment'][0], 'injectMsg') data = moose.Neutral('/data') vm_a = [moose.Table('/data/net2_Vm_A_%d' % (ii)) for ii in range(size)] for tab, comp in zip(vm_a, pop_a['compartment']): moose.connect(tab, 'requestOut', comp, 'getVm') vm_b = [moose.Table('/data/net2_Vm_B_%d' % (ii)) for ii in range(size)] for tab, comp in zip(vm_b, pop_b['compartment']): moose.connect(tab, 'requestOut', comp, 'getVm') gksyn_a = [moose.Table('/data/net2_Gk_syn_a_%d' % (ii)) for ii in range(size)] for tab, synchan in zip(gksyn_a, pop_a['synchan']): moose.connect(tab, 'requestOut', synchan, 'getGk') gksyn_b = [moose.Table('/data/net2_Gk_syn_b_%d' % (ii)) for ii in range(size)] for tab, synchan in zip(gksyn_b, pop_b['synchan']): moose.connect(tab, 'requestOut', synchan, 'getGk') pulsetable = moose.Table('/data/net2_pulse') pulsetable.connect('requestOut', pulse, 'getOutputValue') return {'vm_a': vm_a, 'vm_b': vm_b, 'gksyn_a': gksyn_a, 'gksyn_b': gksyn_b, 'pulse': pulsetable,}
[docs]def single_population(size=2): """Example of a single population where each cell is connected to every other cell. Creates a network of single compartmental cells under /network1 and a pulse generaor """ net = moose.Neutral('network1') pop = create_population(moose.Neutral('/network1'), size) make_synapses(pop['spikegen'], pop['synchan']) pulse = moose.PulseGen('/network1/net1_pulse') pulse.firstLevel = 1e-9 pulse.firstDelay = 0.05 pulse.firstWidth = 0.02 moose.connect(pulse, 'output', pop['compartment'][0], 'injectMsg') data = moose.Neutral('/data') vm = [moose.Table('/data/net1_Vm_%d' % (ii)) for ii in range(size)] for tab, comp in zip(vm, pop['compartment']): moose.connect(tab, 'requestOut', comp, 'getVm') gksyn = [moose.Table('/data/net1_Gk_syn_%d' % (ii)) for ii in range(size)] for tab, synchan in zip(gksyn, pop['synchan']): moose.connect(tab, 'requestOut', synchan, 'getGk') pulsetable = moose.Table('/data/net1_pulse') pulsetable.connect('requestOut', pulse, 'getOutputValue') return {'vm': vm, 'gksyn': gksyn, 'pulse': pulsetable,}
inited = False
[docs]def assign_clocks(model_container_list, simdt, plotdt): """ Assign clocks to elements under the listed paths. This should be called only after all model components have been created. Anything created after this will not be scheduled. """ global inited # `inited` is for avoiding double scheduling of the same object if not inited: print(('SimDt=%g, PlotDt=%g' % (simdt, plotdt))) moose.setClock(0, simdt) moose.setClock(1, simdt) moose.setClock(2, simdt) moose.setClock(3, simdt) moose.setClock(4, plotdt) for path in model_container_list: print(('Scheduling elements under:', path)) moose.useClock(0, '%s/##[TYPE=Compartment]' % (path), 'init') moose.useClock(1, '%s/##[TYPE=Compartment]' % (path), 'process') moose.useClock(2, '%s/##[TYPE=SynChan],%s/##[TYPE=HHChannel]' % (path, path), 'process') moose.useClock(3, '%s/##[TYPE=SpikeGen],%s/##[TYPE=PulseGen]' % (path, path), 'process') moose.useClock(4, '/data/##[TYPE=Table]', 'process') inited = True moose.reinit()
[docs]def main(): """ A demo to create a network of single compartmental neurons connected via alpha synapses. This is same as compartment_net.py except that we avoid ematrix and use single melements. """ # test_compartment() # this calls assign_clocks - after which nothing else will be scheduled. simtime = 0.1 simdt = 0.25e-5 plotdt = 0.25e-3 size = 2 data1 = single_population(size=size) data2 = two_populations(size=size) assign_clocks(['/network1', '/network2'], simdt, plotdt) # assign_clocks(['/network1'], simdt, plotdt) moose.start(simtime) plt.figure(1) plt.suptitle('Single population') plt.subplot(211) for vm in data1['vm']: t = np.linspace(0, simtime, len(vm.vector)) plt.plot(t, vm.vector, label=vm.path) plt.plot(np.linspace(0, simtime, len(data1['pulse'].vector)), data1['pulse'].vector * 1e6, label='Inject(uA)') plt.legend() plt.subplot(212) for gk in data1['gksyn']: t = np.linspace(0, simtime, len(gk.vector)) plt.plot(t, gk.vector, label=gk.path) plt.legend() plt.figure(2) plt.suptitle('Two populations') plt.subplot(221) for vm in data2['vm_a']: t = np.linspace(0, simtime, len(vm.vector)) plt.plot(t, vm.vector, label=vm.path) plt.plot(np.linspace(0, simtime, len(data2['pulse'].vector)), data2['pulse'].vector*1e6, label='Inject(uA)') plt.legend() plt.subplot(223) for vm in data2['vm_b']: t = np.linspace(0, simtime, len(vm.vector)) plt.plot(t, vm.vector, label=vm.path) plt.legend() plt.subplot(222) for gk in data2['gksyn_a']: t = np.linspace(0, simtime, len(gk.vector)) plt.plot(t, gk.vector, label=gk.path) plt.legend() plt.subplot(224) for gk in data2['gksyn_b']: t = np.linspace(0, simtime, len(gk.vector)) plt.plot(t, gk.vector, label=gk.path) plt.legend() plt.show()
# # compartment_net_no_array.py ends here