Source code for cubeMeshSigNeur

# HsolveInstability.py ---

# Commentary:
#
# A toy compartmental neuronal + chemical model in just a cubic volume
#
# 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:

import sys
sys.path.append('../../python')
import os
os.environ['NUMPTHREADS'] = '1'
import math
import numpy
import pylab

import moose

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 = -30e-3 + EREST_ACT
VMAX = 120e-3 + EREST_ACT
VDIVS = 3000

[docs]def createSquid(): """Create a single compartment squid model.""" parent = moose.Neutral ('/n' ) elec = moose.Neutral ('/n/elec' ) compt = moose.SymCompartment( '/n/elec/compt' ) Em = EREST_ACT + 10.613e-3 compt.Em = Em compt.initVm = EREST_ACT compt.Cm = 7.85e-9 * 0.5 compt.Rm = 4.2e5 * 5.0 compt.Ra = 7639.44e3 compt.length = 100e-6 compt.diameter = 4e-6 nachan = moose.HHChannel( '/n/elec/compt/Na' ) nachan.Xpower = 3 xGate = moose.HHGate(nachan.path + '/gateX') xGate.setupAlpha(Na_m_params + [VDIVS, VMIN, VMAX]) xGate.useInterpolation = 1 nachan.Ypower = 1 yGate = moose.HHGate(nachan.path + '/gateY') yGate.setupAlpha(Na_h_params + [VDIVS, VMIN, VMAX]) yGate.useInterpolation = 1 nachan.Gbar = 0.942e-3 nachan.Ek = 115e-3+EREST_ACT moose.connect(nachan, 'channel', compt, 'channel', 'OneToOne') kchan = moose.HHChannel( '/n/elec/compt/K' ) kchan.Xpower = 4.0 xGate = moose.HHGate(kchan.path + '/gateX') xGate.setupAlpha(K_n_params + [VDIVS, VMIN, VMAX]) xGate.useInterpolation = 1 kchan.Gbar = 0.2836e-3 kchan.Ek = -12e-3+EREST_ACT moose.connect(kchan, 'channel', compt, 'channel', 'OneToOne') return compt
def createSynapseOnCompartment( compt ): FaradayConst = 96485.3415 # s A / mol length = compt.length dia = compt.diameter gluR = moose.SynChan( compt.path + '/gluR' ) gluR.tau1 = 4e-3 gluR.tau2 = 4e-3 gluR.Gbar = 1e-6 gluR.Ek= 10.0e-3 moose.connect( compt, 'channel', gluR, 'channel', 'Single' ) gluSyn = moose.SimpleSynHandler( compt.path + '/gluR/sh' ) moose.connect( gluSyn, 'activationOut', gluR, 'activation' ) gluSyn.synapse.num = 1 # Ca comes in through this channel, at least for this example. caPool = moose.CaConc( compt.path + '/ca' ) caPool.CaBasal = 1e-4 # 0.1 micromolar caPool.tau = 0.01 B = 1.0 / ( FaradayConst * length * dia * dia * math.pi / 4) B = B / 20.0 # scaling factor for Ca buffering caPool.B = B moose.connect( gluR, 'IkOut', caPool, 'current', 'Single' ) # Provide a regular synaptic input. synInput = moose.SpikeGen( '/n/elec/compt/synInput' ) synInput.refractT = 47e-3 synInput.threshold = -1.0 synInput.edgeTriggered = 0 synInput.Vm( 0 ) syn = moose.element( gluSyn.path + '/synapse' ) moose.connect( synInput, 'spikeOut', syn, 'addSpike', 'Single' ) syn.weight = 0.2 syn.delay = 1.0e-3 return gluR def createPool( compt, name, concInit ): pool = moose.Pool( compt.path + '/' + name ) pool.concInit = concInit pool.diffConst = 5e-11 return pool # This is a Ca-activated enzyme that phosphorylates and inactivates kChan # as per the following scheme: # Ca + inact_kinase <===> Ca.kinase # kChan ----- Ca.kinase -----> kChan_p # kChan_p -------> kChan def createChemModel( neuroCompt ): dendCa = createPool( neuroCompt, 'Ca', 1e-4 ) dendKinaseInact = createPool( neuroCompt, 'inact_kinase', 1e-4 ) dendKinase = createPool( neuroCompt, 'Ca.kinase', 0.0 ) dendTurnOnKinase = moose.Reac( neuroCompt.path + '/turnOnKinase' ) moose.connect( dendTurnOnKinase, 'sub', dendCa, 'reac' ) moose.connect( dendTurnOnKinase, 'sub', dendKinaseInact, 'reac' ) moose.connect( dendTurnOnKinase, 'prd', dendKinase, 'reac' ) dendTurnOnKinase.Kf = 50000 dendTurnOnKinase.Kb = 1 dendKinaseEnz = moose.Enz( dendKinase.path + '/enz' ) dendKinaseEnzCplx = moose.Pool( dendKinase.path + '/enz/cplx' ) kChan = createPool( neuroCompt, 'kChan', 1e-3 ) kChan_p = createPool( neuroCompt, 'kChan_p', 0.0 ) moose.connect( dendKinaseEnz, 'enz', dendKinase, 'reac', 'OneToOne' ) moose.connect( dendKinaseEnz, 'sub', kChan, 'reac', 'OneToOne' ) moose.connect( dendKinaseEnz, 'prd', kChan_p, 'reac', 'OneToOne' ) moose.connect( dendKinaseEnz, 'cplx', dendKinaseEnzCplx, 'reac', 'OneToOne' ) dendKinaseEnz.Km = 1e-4 dendKinaseEnz.kcat = 20 dendPhosphatase = moose.Reac( neuroCompt.path + '/phosphatase' ) moose.connect( dendPhosphatase, 'sub', kChan_p, 'reac' ) moose.connect( dendPhosphatase, 'prd', kChan, 'reac' ) dendPhosphatase.Kf = 1 dendPhosphatase.Kb = 0.0 def makeModelInCubeMesh(): compt = createSquid() createSynapseOnCompartment( compt ) chem = moose.Neutral( '/n/chem' ) neuroMesh = moose.CubeMesh( '/n/chem/neuroMesh' ) coords = [0] * 9 coords[3] = compt.length coords[4] = compt.diameter coords[5] = compt.diameter coords[6] = compt.length coords[7] = compt.diameter coords[8] = compt.diameter neuroMesh.coords = coords neuroMesh.preserveNumEntries = 1 createChemModel( neuroMesh ) dendCa = moose.element( '/n/chem/neuroMesh/Ca' ) assert dendCa.volume == compt.length * compt.diameter * compt.diameter dendKinaseEnzCplx = moose.element( '/n/chem/neuroMesh/Ca.kinase/enz/cplx' ) assert dendKinaseEnzCplx.volume == dendCa.volume # Make adaptors # Note that we can do this two ways: We can use an existing output # msg from the object, which will come whenever the object processes, # or the adapator can request the object for the field, which happens # whenever the adaptor processes. Here we illustrate both alternatives. adaptK = moose.Adaptor( '/n/chem/neuroMesh/adaptK' ) chemK = moose.element( '/n/chem/neuroMesh/kChan' ) elecK = moose.element( '/n/elec/compt/K' ) moose.connect( adaptK, 'requestOut', chemK, 'getConc', 'OneToAll' ) moose.connect( adaptK, 'output', elecK, 'setGbar', 'OneToAll' ) adaptK.scale = 0.3 # from mM to Siemens adaptCa = moose.Adaptor( '/n/chem/neuroMesh/adaptCa' ) chemCa = moose.element( '/n/chem/neuroMesh/Ca' ) elecCa = moose.element( '/n/elec/compt/ca' ) moose.connect( elecCa, 'concOut', adaptCa, 'input', 'OneToAll' ) moose.connect( adaptCa, 'output', chemCa, 'setConc', 'OneToAll' ) adaptCa.outputOffset = 0.0001 # 100 nM offset in chem conc adaptCa.scale = 0.05 # Empirical: 0.06 max to 0.003 mM def addPlot( objpath, field, plot ): assert moose.exists( objpath ) tab = moose.Table( '/graphs/' + plot ) obj = moose.element( objpath ) moose.connect( tab, 'requestOut', obj, field ) return tab def displayPlots(): for x in moose.wildcardFind( '/graphs/##[ISA=Table]' ): t = numpy.arange( 0, x.vector.size, 1 ) * x.dt pylab.plot( t, x.vector, label=x.name ) pylab.legend() pylab.show() def makeElecPlots(): graphs = moose.Neutral( '/graphs' ) elec = moose.Neutral( '/graphs/elec' ) addPlot( '/n/elec/compt', 'getVm', 'elec/Vm' ) addPlot( '/n/elec/compt/ca', 'getCa', 'elec/Ca' ) addPlot( '/n/elec/compt/K', 'getGk', 'elec/K_Gk' ) def makeChemPlots(): graphs = moose.Neutral( '/graphs' ) addPlot( '/n/chem/neuroMesh/Ca', 'getConc', 'chemCa' ) addPlot( '/n/chem/neuroMesh/kChan_p', 'getConc', 'chemkChan_p' ) addPlot( '/n/chem/neuroMesh/kChan', 'getConc', 'chemkChan' ) addPlot( '/n/chem/neuroMesh/Ca.kinase', 'getConc', 'activeKinase' ) def testCubeMultiscale( useSolver ): elecDt = 10e-6 chemDt = 1e-4 plotDt = 5e-4 plotName = 'cm.plot' if ( useSolver ): elecDt = 50e-6 chemDt = 2e-3 plotName = 'solve_cm.plot' makeModelInCubeMesh() makeChemPlots() makeElecPlots() ''' moose.setClock( 0, elecDt ) moose.setClock( 1, elecDt ) moose.setClock( 2, elecDt ) moose.setClock( 5, chemDt ) moose.setClock( 6, chemDt ) moose.setClock( 7, plotDt ) moose.setClock( 8, plotDt ) moose.useClock( 1, '/n/##[ISA=SpikeGen]', 'process' ) moose.useClock( 2, '/n/##[ISA=SynBase]','process') moose.useClock( 6, '/n/##[ISA=Adaptor]', 'process' ) moose.useClock( 7, '/graphs/#', 'process' ) moose.useClock( 8, '/graphs/elec/#', 'process' ) moose.useClock( 0, '/n/##[ISA=Compartment]', 'init' ) moose.useClock( 1, '/n/##[ISA=Compartment]', 'process' ) moose.useClock( 2, '/n/##[ISA=ChanBase],/n/##[ISA=SynBase],/n/##[ISA=CaConc]','process') moose.useClock( 5, '/n/##[ISA=PoolBase],/n/##[ISA=ReacBase],/n/##[ISA=EnzBase]', 'process' ) ''' if ( useSolver ): ksolve = moose.Ksolve( '/n/ksolve' ) stoich = moose.Stoich( '/n/stoich' ) stoich.compartment = moose.element( '/n/chem/neuroMesh' ) stoich.ksolve = ksolve stoich.path = '/n/##' ksolve.method = 'rk5' #moose.useClock( 5, '/n/ksolve', 'process' ) hsolve = moose.HSolve( '/n/hsolve' ) #moose.useClock( 1, '/n/hsolve', 'process' ) hsolve.dt = elecDt hsolve.target = '/n/compt' moose.reinit() moose.start( 1 ) displayPlots()
[docs]def main(): """ A toy compartmental neuronal + chemical model in just a cubic volume !""" testCubeMultiscale( 1 ) # change argument to 0 to run without solver.
if __name__ == '__main__': main() # cubeMeshSigNeur.py ends here