Thank you for your help. I appreciate it. I tried to implement what you suggested, and it works. In the following code the callback is a wrapper python function and at each time step, z1 is calculated and transferred to the .mod file.
The only problem is that at each time step I want my python code to use the updated alphan, betan, and n for z1 calculation. I thought that
recording_callback = (wrapper, (seg.alphan_Kv1 , seg.betan_Kv1, seg.n_Kv1 , 0 ))
h.CVode().extra_scatter_gather(0, recording_callback)
would update alphan, betan and n per time step. But it seems like it does not. I am wondering how the updated alphan, betan and n from the mod file can be transferred to the python code per time step.
Thanks again.
Code: Select all
from neuron import h, gui
h.load_file("mosinit.hoc")
KHat = zeros((5,5)) # Initial values set to 0
def findnoiseterm(func):
def inner(alphan,betan,n,noise):
global KHat
KNoise = randn(5,5)
AK = array(( (-4*alphan, betan, 0, 0, 0),
(4*alphan, -3*alphan, 0, 0, 0),
(0, 3*alphan, -2*alphan, 3*betan, 0),
(0, 0, 2*alphan, -alphan-3*betan, 4*betan),
(0, 0, 0, alphan, -4*betan) ))
X = array(((1-n)**4, 4*n*(1-n)**3, 6*n**2*(1-n)**2, 4*n**3*(1-n), n**4))
DK= 1/20 * array(( (4*alphan*X[0] + betan*X[1],-(4*alphan*X[0]) , 0, 0, 0),
(-(4*alphan*X[0]+betan*X[1]),(4*alphan*X[1]+(3*alphan+ betan)*X[1]), -(2*betan*X[2] + 3*alphan*X[1]), 0,0),
(0, -(2*betan*X[2]), (3*alphan*X[1]+(2*alphan)*X[2]+3*betan*X[3]),-(3*betan*X[3]), 0),
(0, 0, -(3*betan*X[3]), (2*alphan*X[2]+(3*betan)*X[3]+4*betan*X[4]),-(4*betan*X[4])),
(0, 0, 0, -(4*betan*X[4]), (alphan*X[3]) ) ))
u,s,v = linalg.svd(DK)
SK = u*sqrt(s)*v
KHat += h.dt*dot(AK,KHat) + sqrt(h.dt)*dot(SK,KNoise[:,:])
noise = KHat[4,4]
return noise, func(alphan,betan,n, noise)
return inner
@findnoiseterm
def wrapper(alphan,betan,n, noise):
seg.z1_Kv1=noise
for sec in h.allsec():
if h.ismembrane('Kv1',sec = sec):
for seg in sec:
recording_callback = (wrapper, (seg.alphan_Kv1 , seg.betan_Kv1, seg.n_Kv1 , 0 ))
h.CVode().extra_scatter_gather(0, recording_callback)
h.finitialize(-65)
h.fadvance()