From 3618469db8e8871965941299026a4dd530c44c35 Mon Sep 17 00:00:00 2001
From: Asia Jedrzejewska-Szmek <jjedrzej@gmu.edu>
Date: Thu, 26 Jan 2017 11:54:42 -0500
Subject: [PATCH] Script testing DifShells
Script creates a test model of a dendrite segment with one Ca
channel. DifShells are created. DifBuffers are created. DifShells are
connected to provide diffusion. DifBuffers are connected to provide
diffusion. DifBuffers are connected with DifShells to provide
buffer. MMPumps are created and connected to the outer DifShell. Ca
channel is connected to the outer Difshell to provide Ca influx and
outer DifShell is connected to Ca channel to provide ion
concentration. Segment's Vm, channel's Gk and Ik, DifShell's C and
DifBuffer's bBound are connected with output to save simulation
results.
---
moose-core/tests/python/chan_proto.py | 140 +++++++++++++++++
moose-core/tests/python/param_chan.py | 80 ++++++++++
moose-core/tests/python/params.py | 26 ++++
moose-core/tests/python/test_all.sh | 3 +-
moose-core/tests/python/test_difshells.py | 182 ++++++++++++++++++++++
5 files changed, 430 insertions(+), 1 deletion(-)
create mode 100644 moose-core/tests/python/chan_proto.py
create mode 100644 moose-core/tests/python/param_chan.py
create mode 100644 moose-core/tests/python/params.py
create mode 100644 moose-core/tests/python/test_difshells.py
diff --git a/moose-core/tests/python/chan_proto.py b/moose-core/tests/python/chan_proto.py
new file mode 100644
index 00000000..ae37024f
--- /dev/null
+++ b/moose-core/tests/python/chan_proto.py
@@ -0,0 +1,140 @@
+"""\
+Create general Channel Proto, pass in name, x and y power, and params
+
+Also, create the library of channels
+Might need a few other chan_proto types, such as
+ inf-tau channels
+ Ca dep channels
+chan_proto quires alpha and beta params for both activation and inactivation
+If no inactivation, just send in empty Yparam array.
+"""
+
+from __future__ import print_function, division
+import moose
+import numpy as np
+from util import NamedList
+
+SSTauChannelParams = NamedList('SSTauChannelParams', '''
+ Arate
+ A_B
+ A_C
+ Avhalf
+ Avslope
+ taumin
+ tauVdep
+ tauPow
+ tauVhalf
+ tauVslope''')
+
+AlphaBetaChannelParams = NamedList('AlphaBetaChannelParams', '''
+ A_rate
+ A_B
+ A_C
+ Avhalf
+ A_vslope
+ B_rate
+ B_B
+ B_C
+ Bvhalf
+ B_vslope''')
+
+ZChannelParams = NamedList('ZChannelParams', 'Kd power tau')
+
+
+ChannelSettings = NamedList('ChannelSettings', 'Xpow Ypow Zpow Erev name')
+
+def interpolate_values_in_table(tabA,V_0,l=40):
+ import param_chan
+ '''This function interpolates values in the table
+ around tabA[V_0]. '''
+ V = np.linspace(param_chan.VMIN, param_chan.VMAX, len(tabA))
+ idx = abs(V-V_0).argmin()
+ A_min = tabA[idx-l]
+ V_min = V[idx-l]
+ A_max = tabA[idx+l]
+ V_max = V[idx+l]
+ a = (A_max-A_min)/(V_max-V_min)
+ b = A_max - a*V_max
+ tabA[idx-l:idx+l] = V[idx-l:idx+l]*a+b
+ return tabA
+
+def fix_singularities(Params,Gate):
+ import param_chan
+
+ if Params.A_C < 0:
+
+ V_0 = Params.A_vslope*np.log(-Params.A_C)-Params.Avhalf
+
+ if V_0 > param_chan.VMIN and V_0 < param_chan.VMAX:
+ #change values in tableA and tableB, because tableB contains sum of alpha and beta
+ tabA = interpolate_values_in_table(Gate.tableA,V_0)
+ tabB = interpolate_values_in_table(Gate.tableB,V_0)
+ Gate.tableA = tabA
+ Gate.tableB = tabB
+
+ if Params.B_C < 0:
+
+ V_0 = Params.B_vslope*np.log(-Params.B_C)-Params.Bvhalf
+
+ if V_0 > param_chan.VMIN and V_0 < param_chan.VMAX:
+ #change values in tableB
+ tabB = interpolate_values_in_table(Gate.tableB,V_0)
+ Gate.tableB = tabB
+
+ return Gate
+
+#may need a CaV channel if X gate uses alpha,beta and Ygate uses inf tau
+#Or, have Y form an option - if in tau, do something like NaF
+def chan_proto(chanpath, params):
+ import param_chan
+
+ chan = moose.HHChannel(chanpath)
+ chan.Xpower = params.channel.Xpow
+ if params.channel.Xpow > 0:
+ xGate = moose.HHGate(chan.path + '/gateX')
+ xGate.setupAlpha(params.X + [param_chan.VDIVS, param_chan.VMIN, param_chan.VMAX])
+ xGate = fix_singularities(params.X, xGate)
+
+ chan.Ypower = params.channel.Ypow
+ if params.channel.Ypow > 0:
+ yGate = moose.HHGate(chan.path + '/gateY')
+ yGate.setupAlpha(params.Y + [param_chan.VDIVS, param_chan.VMIN, param_chan.VMAX])
+ yGate = fix_singularities(params.Y, yGate)
+ if params.channel.Zpow > 0:
+ chan.Zpower = params.channel.Zpow
+ zgate = moose.HHGate(chan.path + '/gateZ')
+ ca_array = np.linspace(param_chan.CAMIN, param_chan.CAMAX, param_chan.CADIVS)
+ zgate.min = param_chan.CAMIN
+ zgate.max = param_chan.CAMAX
+ caterm = (ca_array/params.Z.Kd) ** params.Z.power
+ inf_z = caterm / (1 + caterm)
+ tau_z = params.Z.tau * np.ones(len(ca_array))
+ zgate.tableA = inf_z / tau_z
+ zgate.tableB = 1 / tau_z
+ chan.useConcentration = True
+ chan.Ek = params.channel.Erev
+ return chan
+
+
+TypicalOneDalpha = NamedList('TypicalOneDalpha',
+ '''channel X Y Z=[] calciumPermeable=False calciumPermeable2=False''')
+
+_FUNCTIONS = {
+ TypicalOneDalpha: chan_proto,
+
+}
+
+#*params... passes the set of values not as a list but as individuals
+def make_channel(chanpath, params):
+ func = _FUNCTIONS[params.__class__]
+ return func(chanpath, params)
+
+def chanlib():
+ import param_chan
+ if not moose.exists('/library'):
+ moose.Neutral('/library')
+ #Adding all the channels to the library. *list removes list elements from the list,
+ #so they can be used as function arguments
+ chan = [make_channel('/library/'+key, value) for key, value in param_chan.ChanDict.items()]
+
+ return chan
diff --git a/moose-core/tests/python/param_chan.py b/moose-core/tests/python/param_chan.py
new file mode 100644
index 00000000..0db3f697
--- /dev/null
+++ b/moose-core/tests/python/param_chan.py
@@ -0,0 +1,80 @@
+# -*- coding: utf-8 -*-
+
+from util import NamedDict
+from chan_proto import (
+ SSTauChannelParams,
+ AlphaBetaChannelParams,
+ ZChannelParams,
+ ChannelSettings,
+ TypicalOneDalpha,)
+
+#chanDictSP.py
+#contains all gating parameters and reversal potentials
+# Gate equations have the form:
+#
+# y(x) = (rate + B * x) / (C + exp((x + vhalf) / vslope))
+#
+# OR
+# y(x) = tau_min + tau_vdep / (1 + exp((x + vhalf) / vslope))
+#
+# where x is membrane voltage and y is the rate constant
+#KDr params used by Sriram, RE paper1, Krp params used by RE paper 2
+#Parameters for Ca channels may need to be shifted - see Dorman model
+krev=-87e-3
+narev=50e-3
+carev=140e-3 #assumes CaExt=2 mM and CaIn=50e-3
+ZpowCDI=0
+
+VMIN = -100e-3
+VMAX = 50e-3
+VDIVS = 3000 #0.5 mV steps
+CAMIN=0.01e-3 #10 nM
+CAMAX=40e-3 #40 uM, might want to go up to 100 uM with spines
+CADIVS=3000 #10 nM steps
+
+#mtau: Ogata fig 5, no qfactor accounted in mtau, 1.2 will improve spike shape
+#activation minf fits Ogata 1990 figure 3C (which is cubed root)
+#inactivation hinf fits Ogata 1990 figure 6B
+#htau fits the main -50 through -10 slope of Ogata figure 9 (log tau), but a qfact of 2 is already taken into account.
+
+
+CaTparam = ChannelSettings(Xpow=3, Ypow=1, Zpow=1, Erev=carev, name='CaT')
+qfactCaT = 1
+
+#Original inactivation ws too slow compared to activation, made closder the alpha1G
+CaT_X_params = AlphaBetaChannelParams(A_rate = 5342.4*qfactCaT,
+ A_B = 2100.*qfactCaT,
+ A_C = 11.9,
+ Avhalf = 1e-3,
+ A_vslope = -12e-3,
+ B_rate = 289.7*qfactCaT,
+ B_B = 0.*qfactCaT,
+ B_C = 1.,
+ Bvhalf = 0.0969,
+ B_vslope = 0.0141)
+#CaT_Y_params = CaT_X_params
+CaT_Y_params = AlphaBetaChannelParams(A_rate = 0*qfactCaT,
+ A_B = -74.*qfactCaT,
+ A_C = 1.,
+ Avhalf = 9e-2,
+ A_vslope = 5e-3,
+ B_rate = 15*qfactCaT,
+ B_B = -1.5*qfactCaT,
+ B_C = 1.5,
+ Bvhalf = 0.05,
+ B_vslope = -15e-3)
+
+#These CDI params can be used with every channel, make ZpowCDI=2
+#If ZpowCDI=0 the CDI will not be used, power=-4 is to transform
+#(Ca/Kd)^pow/(1+(Ca/Kd)^pow) to 1/(1+(ca/Kd)^-pow)
+CDI_Z_params = ZChannelParams(Kd = 0.12e-3,
+ power = -4,
+ tau = 142e-3)
+
+#Dictionary of "standard" channels, to create channels using a loop
+#NaF doesn't fit since it uses different prototype form
+#will need separate dictionary for BK
+
+ChanDict = NamedDict(
+ 'ChannelParams',
+ CaT = TypicalOneDalpha(CaTparam, CaT_X_params, CaT_Y_params, CDI_Z_params, calciumPermeable=True))
diff --git a/moose-core/tests/python/params.py b/moose-core/tests/python/params.py
new file mode 100644
index 00000000..40696573
--- /dev/null
+++ b/moose-core/tests/python/params.py
@@ -0,0 +1,26 @@
+t_stop = 10
+dend_diameter = 2.2627398e-6
+dend_length = 1.131369936e-6
+Cm = 4.021231698e-12
+Rm = 1865100032
+Em = -0.07100000232
+Vm_0 = -0.0705
+dt = 50e-6
+spines_no = 0
+difshell_no = 2
+difshell_name = "Ca_shell"
+Ca_basal = 50e-6
+Ca_initial = Ca_basal*200
+dca = 200.0e-12
+difbuff_no = 1
+difbuff_name = "Buff"
+btotal = 80.0e-3
+kf = 0.028e6
+kb = 19.6
+d = 66e-12
+inject = 0.1e-9
+gbar = 1
+#MMpump
+km = 0.3e-3
+kcat = 85e-22
+pumps = 1
diff --git a/moose-core/tests/python/test_all.sh b/moose-core/tests/python/test_all.sh
index 9ec99b4e..338a20cb 100755
--- a/moose-core/tests/python/test_all.sh
+++ b/moose-core/tests/python/test_all.sh
@@ -6,7 +6,8 @@ test_scripts="./test_function.py
test_mumbl.py
test_pymoose.py
test_synchan.py
- test_vec.py"
+ test_vec.py
+ test_difshells.py"
for testFile in $test_scripts; do
echo "Executing $testFile"
diff --git a/moose-core/tests/python/test_difshells.py b/moose-core/tests/python/test_difshells.py
new file mode 100644
index 00000000..8e6da7c8
--- /dev/null
+++ b/moose-core/tests/python/test_difshells.py
@@ -0,0 +1,182 @@
+import moose
+import numpy as np
+import matplotlib.pyplot as plt
+import chan_proto
+
+from params import *
+
+def linoid(x,param):
+ #plt.figure()
+ #plt.plot(x,(param[2]+np.exp((V+param[3])/param[4])))
+ den = (param[2]+np.exp((V+param[3])/param[4]))
+ nom = (param[0]+param[1]*V)
+ return nom/den
+
+def add_difshell(comp,i,shell_thickness,shell_radius):
+ new_name = comp.path.split('[')[0]+'/'+comp.name+'/'+difshell_name+str(i)
+ dif = moose.DifShell(new_name)
+
+
+ #dif.C = Ca_initial
+ dif.Ceq = Ca_basal
+ dif.D = dca
+ dif.valence = 2
+ dif.leak = 0
+ dif.shapeMode = 0
+ dif.length = comp.length
+ dif.diameter = 2*shell_radius
+ dif.thickness = shell_thickness
+
+ return dif
+
+def add_difbuffer_to_dishell(comp,i,j,shell_thickness,shell_radius):
+ new_name = comp.path.split('[')[0]+'/'+comp.name+'/'+difshell_name+str(i)+'_'+difbuff_name+str(j)
+ buf = moose.DifBuffer(new_name)
+ buf.bTot = btotal
+ buf.shapeMode = 0
+ buf.length = comp.length
+ buf.diameter = 2*shell_radius
+ buf.thickness = shell_thickness
+ buf.kf = kf
+ buf.kb = kb
+ buf.D = d
+ return buf
+
+def add_difshells_and_buffers(comp):
+
+ if difshell_no < 1:
+ return [], []
+
+ difshell= []
+ shell_thickness = dend.diameter/difshell_no/2.
+ difbuffer = []
+ for i in range(difshell_no):
+ shell_radius = comp.diameter/2-i*shell_thickness
+ dif = add_difshell(comp,i,shell_thickness,shell_radius)
+ difshell.append(dif)
+
+ if i > 0 :
+ moose.connect(difshell[i-1],"outerDifSourceOut",difshell[i],"fluxFromOut")
+ moose.connect(difshell[i],"innerDifSourceOut",difshell[i-1],"fluxFromIn")
+
+ if difbuff_no > 0:
+ difbuffer.append([])
+ for j in range(difbuff_no):
+ buf = add_difbuffer_to_dishell(comp,i,j,shell_thickness,shell_radius)
+ difbuffer[i].append(buf)
+ moose.connect(difshell[i],"concentrationOut",buf,"concentration")
+ moose.connect(buf,"reactionOut",difshell[i],"reaction")
+ if i>0:
+ moose.connect(difbuffer[i-1][j],"outerDifSourceOut",difbuffer[i][j],"fluxFromOut")
+ moose.connect(difbuffer[i][j],"innerDifSourceOut",difbuffer[i-1][j],"fluxFromIn")
+
+ return difshell,difbuffer
+
+def addOneChan(chanpath,gbar,comp):
+
+ SA = np.pi * comp.length*comp.diameter
+ proto = moose.element('/library/'+chanpath)
+ chan = moose.copy(proto, comp, chanpath)
+ chan.Gbar = gbar * SA
+ #If we are using GHK AND it is a calcium channel, connect it to GHK
+ moose.connect(comp,'VmOut',chan,'Vm')
+ moose.connect(chan,"channelOut",comp,"handleChannel")
+ return chan
+
+
+if __name__ == '__main__':
+
+ for tick in range(0, 7):
+ moose.setClock(tick, dt)
+ moose.setClock(8,0.005) #set output clock
+ model = moose.Neutral('/model')
+ dend = moose.Compartment('/model/dend')
+ pulse = moose.PulseGen('/model/pulse')
+ data = moose.Neutral('/data')
+ vmtab = moose.Table('/data/dend_Vm')
+ gktab = moose.Table('/data/CaT_Gk')
+ iktab = moose.Table('/data/CaT_Ik')
+
+ dend.Cm = Cm
+ dend.Rm = Rm
+ dend.Em = Em
+ dend.initVm = Vm_0
+ dend.diameter = dend_diameter
+ dend.length = dend_length
+
+ pulse.delay[0] = 8.
+ pulse.width[0] = 500e-3
+ pulse.level[0] = inject
+ pulse.delay[1] = 1e9
+
+ m = moose.connect(pulse, 'output', dend, 'injectMsg')
+
+ moose.connect(vmtab, 'requestOut', dend, 'getVm')
+
+ chan_proto.chanlib()
+ chan = addOneChan('CaT',gbar,dend)
+
+
+ moose.connect(gktab, 'requestOut', chan, 'getGk')
+ moose.connect(iktab, 'requestOut', chan, 'getIk')
+ diftab = []
+ buftab = []
+ difs, difb = add_difshells_and_buffers(dend)
+ if pumps:
+ pump = moose.MMPump('/model/dend/pump')
+ pump.Vmax = kcat
+ pump.Kd = km
+ moose.connect(pump,"PumpOut",difs[0],"mmPump")
+ if difs:
+ moose.connect(chan,"IkOut",difs[0],"influx")
+ moose.connect(difs[0],'concentrationOut',chan,'concen')
+
+ for i,dif in enumerate(difs):
+ res_dif = moose.Table('/data/'+difshell_name+str(i))
+ diftab.append(res_dif)
+ moose.connect(diftab[i],'requestOut',dif,'getC')
+ if (difbuff_no):
+ buftab.append([])
+ for j,buf in enumerate(difb[i]):
+ res_buf = moose.Table('/data/'+difshell_name+str(i)+'_'+difbuff_name+str(j))
+ buftab[i].append(res_buf)
+ moose.connect(buftab[i][j],'requestOut',buf,'getBBound')
+
+
+ moose.reinit()
+ if not gbar:
+ for i,dif in enumerate(difs):
+ if i == 0:
+ dif.C = Ca_initial
+ else:
+ dif.C = 0
+ for j,dbuf in enumerate(difb[i]):
+ dbuf.bFree = dbuf.bTot
+
+ moose.start(t_stop)
+
+ t = np.linspace(0, t_stop, len(vmtab.vector))
+ fname = 'moose_results_difshell_no_'+str(difshell_no)+'_difbuffer_no_'+str(difbuff_no)+'_pump_'+str(pumps)+'_gbar_'+str(gbar)+'.txt'
+ print fname
+ header = 'time Vm Ik Gk'
+ number = 4+difshell_no*(difbuff_no+1)
+ res = np.zeros((len(t),number))
+ res[:,0] = t
+ res[:,1] = vmtab.vector
+ res[:,2] = iktab.vector
+ res[:,3] = gktab.vector
+
+ for i in range(difshell_no):
+ header += ' difshell_'+str(i)
+ res[:,4+i*(difbuff_no+1)] = diftab[i].vector
+ if (difbuff_no):
+ for j,buf in enumerate(buftab[i]):
+ res[:,4+i*(difbuff_no+1)+j+1] = buf.vector
+ header += ' difshell_'+str(i)+'_difbuff_'+str(j)
+ np.savetxt(fname,res,header=header, comments='')
+
+
+
+
+
+ # plt.show()
--
GitLab