Join the social network of Tech Nerds, increase skill rank, get work, manage projects...

• # How to Combine Stochastic Hodgkin Huxley Model Code with Voltage Dependent Calcium Channel

• 0
• 0
• 0
• 0
• 0
• 0
• 0
• 190

I have written a python script for Hodgkin huxley model as well as for a voltage dependent calcium channel. I want to combine both the code together and want that the voltage given by hodgkin huxley model should be the input for voltage dependent calcium channels.

here's a look at my code:

import numpy as np
import matplotlib.pyplot as plt

#defining inputs
Cm = 1.0
gK = 36.0
VK = -12.0
gl = 0.3
Vl = 10.6
gNa = 120.0
VNa = 120.0

#rate equation for pottasium channels
def alpha_n(V):
return (0.01*(10.0 - V))/(np.exp(1.0-(0.1*V)) - 1.0)
def beta_n(V):
return 0.125 * np.exp(-V/80.0)

#rate equations for sodium channels
def alpha_m(V):
return 0.1 * (25-V) / (np.exp((25-V)/10)-1)
def beta_m(V):
return 4 * np.exp(-V/18)
def alpha_h(V):
return 0.07 * np.exp(-V/20)
def beta_h(V):
return 1 / (np.exp((30-V)/10)+1)
def alpha1(V):
return 4.04*np.exp(V/49.14)
def alpha2(V):
return 6.70*np.exp(V/42.08)
def alpha3(V):
return 4.39*np.exp(V/55.31)
def alpha4(V):
return 17.33*np.exp(V/26.55)
def beta1(V):
return 2.88*np.exp(-V/49.14)
def beta2(V):
return 6.30*np.exp(-V/42.08)
def beta3(V):
return 8.16*np.exp(-V/55.31)
def beta4(V):
return 1.84*np.exp(-V/26.55)

#gillepsie algorithm
def condition(V,Kin,Nain,t,dt):

tau = t
Ks = Kin
Nas = Nain

while(tau<(t+dt)):

r = np.zeros((28))
r[0] = 3*alpha_m(V)*Nas[0,0]
r[1] = r[0] + 2*alpha_m(V)*Nas[1,0]
r[2] = r[1] + 1*alpha_m(V)*Nas[2,0]
r[3] = r[2] + 3*beta_m(V)*Nas[3,0]
r[4] = r[3] + 2*beta_m(V)*Nas[2,0]
r[5] = r[4] + 1*beta_m(V)*Nas[1,0]
r[6] = r[5] + 1*alpha_h(V)*Nas[0,0]
r[7] = r[6] + 1*alpha_h(V)*Nas[1,0]
r[8] = r[7] + 1*alpha_h(V)*Nas[2,0]
r[9] = r[8] + 1*alpha_h(V)*Nas[3,0]
r[10] = r[9] + 1*beta_h(V)*Nas[0,1]
r[11] = r[10] + 1*beta_h(V)*Nas[1,1]
r[12] = r[11] + 1*beta_h(V)*Nas[2,1]
r[13] = r[12] + 1*beta_h(V)*Nas[3,1]
r[14] = r[13] + 3*alpha_m(V)*Nas[0,1]
r[15] = r[14] + 2*alpha_m(V)*Nas[1,1]
r[16] = r[15] + 1*alpha_m(V)*Nas[2,1]
r[17] = r[16] + 3*beta_m(V)*Nas[3,1]
r[18] = r[17] + 2*beta_m(V)*Nas[2,1]
r[19] = r[18] + 1*beta_m(V)*Nas[1,1]
r[20] = r[19] + 4*alpha_n(V)*Ks[0]
r[21] = r[20] + 3*alpha_n(V)*Ks[1]
r[22] = r[21] + 2*alpha_n(V)*Ks[2]
r[23] = r[22] + 1*alpha_n(V)*Ks[3]
r[24] = r[23] + 4*beta_n(V)*Ks[4]
r[25] = r[24] + 3*beta_n(V)*Ks[3]
r[26] = r[25] + 2*beta_n(V)*Ks[2]
r[27] = r[26] + 1*beta_n(V)*Ks[1]

R = r[27]
dt = -np.log(np.random.uniform(0,1))/R
tau +=dt

while(tau<(t+dt)):

r = np.random.randn()*R
if(r<r[0]):
Nas[0,0] -= 1
Nas[1,0] += 1
elif(r<r[1]):
Nas[1,0] -= 1
Nas[2,0] += 1
elif(r<r[2]):
Nas[2,0] -= 1
Nas[3,0] += 1
elif(r<r[3]):
Nas[3,0] -= 1
Nas[2,0] += 1
elif(r<r[4]):
Nas[2,0] -= 1
Nas[1,0] += 1
elif(r<r[5]):
Nas[1,0] -= 1
Nas[0,0] += 1
elif(r<r[6]):
Nas[0,0] -= 1
Nas[0,1] += 1
elif(r<r[7]):
Nas[1,0] -= 1
Nas[1,1] += 1
elif(r<r[8]):
Nas[2,0] -= 1
Nas[2,1] += 1
elif(r<r[9]):
Nas[3,0] -= 1
Nas[3,1] += 1
elif(r<r[10]):
Nas[0,1] -= 1
Nas[0,0] += 1
elif(r<r[11]):
Nas[1,1] -= 1
Nas[1,0] += 1
elif(r<r[12]):
Nas[2,1] -= 1
Nas[2,0] += 1
elif(r<r[13]):
Nas[3,1] -= 1
Nas[3,0] += 1
elif(r<r[14]):
Nas[0,1] -= 1
Nas[1,1] += 1
elif(r<r[15]):
Nas[1,1] -= 1
Nas[2,1] += 1
elif(r<r[16]):
Nas[2,1] -= 1
Nas[3,1] += 1
elif(r<r[17]):
Nas[3,1] -= 1
Nas[2,1] += 1
elif(r<r[18]):
Nas[2,1] -= 1
Nas[1,1] += 1
elif(r<r[19]):
Nas[1,1] -= 1
Nas[0,1] += 1
elif(r<r[20]):
Ks[0] -= 1
Ks[1] += 1
elif(r<r[21]):
Ks[1] -= 1
Ks[2] += 1
elif(r<r[22]):
Ks[2] -= 1
Ks[3] += 1
elif(r<r[23]):
Ks[3] -= 1
Ks[4] += 1
elif(r<r[24]):
Ks[4] -= 1
Ks[3] += 1
elif(r<r[25]):
Ks[3] -= 1
Ks[2] += 1
elif(r<r[26]):
Ks[2] -= 1
Ks[1] += 1
elif(r<r[27]):
Ks[1] -= 1
Ks[0] += 1

Ko = Ks
Nao = Nas
return Ko,Nao

def condition1(V,Cain,t,dt):
Cas = Cain
tau = t

while(tau<(t+dt)):
r = np.zeros((8))
r[0] = 1*alpha1(V)*Cas[0]
r[1] = r[0] + 1*alpha2(V)*Cas[1]
r[2] = r[1] + 1*alpha3(V)*Cas[2]
r[3] = r[2] + 1*alpha4(V)*Cas[3]
r[4] = r[3] + 1*beta4(V)*Cas[4]
r[5] = r[4] + 1*beta3(V)*Cas[3]
r[6] = r[5] + 1*beta2(V)*Cas[2]
r[7] = r[6] + 1*beta1(V)*Cas[1]

R = r[7]
dt = -np.log(np.random.uniform(0,1))/R
tau += dt

while(tau<(t+dt)):
r = np.random.randn()*R
if(r<r[0]):
Cas[0] -= 1
Cas[1] += 1
elif(r<r[1]):
Cas[1] -= 1
Cas[2] += 1
elif(r<r[2]):
Cas[2] -= 1
Cas[3] += 1
elif(r<r[3]):
Cas[3] -= 1
Cas[4] += 1
elif(r<r[4]):
Cas[4] -= 1
Cas[3] += 1
elif(r<r[5]):
Cas[3] -= 1
Cas[2] += 1
elif(r<r[6]):
Cas[2] -= 1
Cas[1] += 1
else:
Cas[1] -= 1
Cas[0] += 1

Cao = Cas
return Cao

#defining a class object
class Model(object):
def __init__(self):
ta = np.arange(0.0,100.01,0.001) #time array
Cin = np.zeros(len(ta)) #input current
for i,t in enumerate(ta):
if(5<=t<=15):
Cin[i] = 10.0
s = 'ODE' # input model

self.ta = ta
self.Cin = Cin
self.s = s

def __call__(self):
self.Final_model()
self.plot()
self.subunit_plot()

def Final_model(self):
t = self.ta
Cin = self.Cin
s = self.s

dt = t[1] - t[0] #step size
n_steps = len(t) #iterate over
solve = n_steps - 1 # at which to solve

v_array = np.zeros((n_steps)) #voltage array
n_array_4 = np.zeros((n_steps)) #n subunit array
m_array_4 = np.zeros((n_steps)) #m subunit array
h_array_4 = np.zeros((n_steps)) #h subunit array
c0_array = np.zeros((n_steps))
c1_array = np.zeros((n_steps))
c2_array = np.zeros((n_steps))
c3_array = np.zeros((n_steps))
o_array = np.zeros((n_steps))

t0 = t[0] #initial time
V = v_array[0] #initial potential
n = alpha_n(V)/(alpha_n(V) + beta_n(V)) #initial n
m = alpha_m(V)/(alpha_m(V) + beta_m(V)) #initial m
h = alpha_h(V)/(alpha_h(V) + beta_h(V)) #initial h
C0 = c0_array[0]
C1 = c1_array[0]
C2 = c2_array[0]
C3 = c3_array[0]
O = o_array[0]

NK = 150 #number of pottasium channels
NNa = 150 #number of sodium channels

noise_n = lambda V,n: np.sqrt((alpha_n(V)*(1-n) + beta_n(V)*n)/NK)*np.random.randn() #n subunit noise
noise_m = lambda V,m: np.sqrt((alpha_m(V)*(1-m) + beta_m(V)*m)/NNa)*np.random.randn() #m subunit noise
noise_h = lambda V,h: np.sqrt((alpha_h(V)*(1-h) + beta_h(V)*h)/NNa)*np.random.randn() #h subunit noise

#Markov chain
if(s=="MC"):
Na_mc = np.zeros((4,2))
Na_mc[0,0] = np.floor(NNa*(1-m)**3*(1-h))
Na_mc[1,0] = np.floor(NNa*3*m*(1-m)**2*(1-h))
Na_mc[2,0] = np.floor(NNa*3*m**2*(1-m)*(1-h))
Na_mc[3,0] = np.floor(NNa*m**3*(1-h))
Na_mc[0,1] = np.floor(NNa*(1-m)**3*h)
Na_mc[1,1] = np.floor(NNa*3*m*(1-m)**2*h)
Na_mc[2,1] = np.floor(NNa*3*m**2*(1-m)*h)
Na_mc[3,1] = NNa - sum(sum(Na_mc))

K_mc = np.zeros(5)
K_mc[0] = np.floor(NK*(1-n)**4)
K_mc[1] = np.floor(NK*4*n*(1-n)**3)
K_mc[2] = np.floor(NK*6*n**2*(1-n)**2)
K_mc[3] = np.floor(NK*4*n**3*(1-n))
K_mc[4] = np.floor(NK*n**4)

for i in range(1,n_steps):
I = Cin[i-1] #input current
n += dt*(alpha_n(V)*(1-n) - beta_n(V)*n) + noise_n(V,n)*np.sqrt(dt) #updated n subunit with noise
m += dt*(alpha_m(V)*(1-m) - beta_m(V)*m) + noise_m(V,m)*np.sqrt(dt) #updated m subunit with noise
h += dt*(alpha_h(V)*(1-h) - beta_h(V)*h) + noise_h(V,h)*np.sqrt(dt) #updated h subunit with noise

if(s=='MC'):
Na_mc, K_mc = condition(V,Na_mc,K_mc,t0,dt)

C0_ = (beta1(V)*C1) - (alpha1(V)*C0)
C1_ = (alpha1(V)*C0) + (beta2(V)*C2) - (alpha2(V)+beta1(V))*C1
C2_ = (alpha2(V)*C1) + (beta3(V)*C3) - (alpha3(V)+beta2(V))*C2
C3_ = (alpha3(V)*C2) + (beta4(V)*O) - (alpha4(V)+beta3(V))*C3
O_ = (alpha4(V)*C3) - (beta4(V)*O)
C0 += dt*C0_
C1 += dt*C1_
C2 +=  dt*C2_
C3 += dt*C3_
O += dt*O_

V_ = (I/Cm) - ((gNa*(V-VNa)*(np.power(m,3)*h))/Cm) - ((gK*(V-VK)*(np.power(n,4)))/Cm) - (gl*(V-Vl))/Cm
V += dt*V_ + np.sqrt(dt)/Cm #updated potential
v_array[i] = V
n_array_4[i] = n
m_array_4[i] = m
h_array_4[i] = h
c0_array[i] = C0
c1_array[i] = C1
c2_array[i] = C2
c3_array[i] = C3
o_array[i] = O

self.v_array = v_array
self.n_array_4 = n_array_4
self.m_array_4 = m_array_4
self.h_array_4 = h_array_4
self.c0_array = c0_array
self.c1_array = c1_array
self.c2_array = c2_array
self.c3_array = c3_array
self.o_array = o_array

def plot(self):
t = self.ta
V = self.v_array
plt.plot(t,V)
plt.xlabel("Time(ms)")
plt.ylabel("Spikes")
plt.show()

def subunit_plot(self):
t = self.ta
n_array_4 = self.n_array_4
m_array_4 = self.m_array_4
h_array_4 = self.h_array_4
plt.plot(t,n_array_4,'b',label="n subunit")
plt.plot(t,m_array_4,'y',label="m sububit")
plt.plot(t,h_array_4,'k',label="h subunit")
plt.xlabel("Time(ms)")
plt.ylabel("Subunits")
plt.xlim(0,10)
plt.legend()
plt.show()

c0_array = self.c0_array
c1_array = self.c1_array
c2_array = self.c2_array
c3_array = self.c3_array
o_array = self.o_array
plt.plot(t,c0_array)
plt.plot(t,c1_array)
plt.plot(t,c2_array)
plt.plot(t,c3_array)
plt.plot(t,o_array)
plt.show()

OR
OR
Register

OR
Fill out the form below and instructions to reset your password will be emailed to you:

• Hire
• Post Projects

### Post Projects

• All at 0 Cost ....
• Post Tech Job
• Select Best Bidder
• Track the Project
• Approve Work and Pay safely
• Browse Nerds
• Work
• Find Projects Find Projects
• UI Design and UX
• Software Engineering
View more...
View less...
• Marketing
• General
View more...
View less...
• Manage
• Company Company

### Manage Company

• All at 0 Cost ....
• Manage Company and Employee Profiles
• Company wide Employee Productivity Reports
• Knowledge Sharing and Collaboration Tools
• Get Sales Lead and Bid for Tech Projects
• Send Invoices and Receive Payment Safely
• Learn
• Nerd Digest Nerd Digest
• UI Design and UX
• Software Engineering
View more...
View less...
• Marketing
• General
View more...
View less...
• Tech Q & A Tech Q & A
• UI Design and UX
• Software Engineering
View more...
View less...
• Marketing
• General
View more...
View less...