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
    • 401
    Answer it

    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()
            
     

 0 Answer(s)

Sign In
                           OR                           
                           OR                           
Register

Sign up using

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