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