diff --git a/moose-core/tests/python/chan_proto.py b/moose-core/tests/python/chan_proto.py
new file mode 100644
index 0000000000000000000000000000000000000000..ae37024f5de98402d8df87822af6258aaeef7b54
--- /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 0000000000000000000000000000000000000000..0db3f697226b9e30aa0a766d5e4b41ad8255cad8
--- /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 0000000000000000000000000000000000000000..406965738033096b2fd85fb6e4250c4decd88554
--- /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 9ec99b4e4ee2dda8faa1b4f5072ed0e3c2f99052..338a20cb746752c06064f4b6a3a4eb6d83a9f221 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 0000000000000000000000000000000000000000..8e6da7c8ccb8850ba98f529d8c690b8f21cad281
--- /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()