心臓刺激伝導系の数理モデルのPythonでのシミュレーション

import math
import numpy as np
import matplotlib.pyplot as plt 
from scipy.integrate import odeint

class NobleModel():
    C_m = 12.0 #membrane capacitance, in uF/cm^2
    g_Na = 400.0 #Sodium (Na) maximum conductances, in mS/cm^2
    g_K = 1.2 #Postassium (K) maximum conductances, in mS/cm^2
    g_Cl = 0.075 #Leak maximum conductances, in mS/cm^2
    E_Na = 40.0 #Sodium (Na) Nernst reversal potentials, in mV
    E_K = -100.0 #Postassium (K) Nernst reversal potentials, in mV
    E_Cl = -60 #Leak Nernst reversal potentials, in mV
    
    dt = 0.1 #ms
    T = 1000 #ms
    t = np.arange(0.0, T, dt) #The time to integrate over

    def alpha_m(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 0.1*(V + 48) / (1 - math.exp(-(V + 48)/15))

    def beta_m(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 0.12*(V + 8) / (math.exp((V + 8)/5) - 1)

    def alpha_h(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 0.17*math.exp(-(V + 90)/20)

    def beta_h(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 1 / (1 + math.exp(-(V + 42)/10))

    def alpha_n(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 1e-4*(V + 50) / (1 - math.exp(-(V + 50)/10))

    def beta_n(self, V):
        """Channel gating kinetics. Functions of membrane voltage"""
        return 2e-3*math.exp(-(V + 90)/80)

    def I_Na(self, V, m, h):
        return (self.g_Na * (m**3) * h + 0.14)* (V - self.E_Na)

    def I_K(self, V, n):
        gK1 = self.g_K * np.exp(-(V + 90)/50) + 0.015*np.exp((V + 90)/60)
        gK2 = self.g_K * (n**4)
        return (gK1 + gK2) * (V - self.E_K)

    def I_Cl(self, V):
        return self.g_Cl * (V - self.E_Cl)

    @staticmethod
    def dALLdt(X, t, self):
        V, m, h, n = X

        dVdt = - (self.I_Na(V, m, h) + self.I_K(V, n) + self.I_Cl(V)) / self.C_m
        dmdt = self.alpha_m(V)*(1.0-m) - self.beta_m(V)*m
        dhdt = self.alpha_h(V)*(1.0-h) - self.beta_h(V)*h
        dndt = self.alpha_n(V)*(1.0-n) - self.beta_n(V)*n
        return dVdt, dmdt, dhdt, dndt

    def simulation(self):
        initX = [-75, 0.06, 0.7, 0.3] # V, m, h, n
        X = odeint(self.dALLdt, initX, self.t, args=(self,))
        V = X[:,0]
        m = X[:,1]
        h = X[:,2]
        n = X[:,3]
        ina = self.I_Na(V, m, h)
        ik = self.I_K(V, n)
        icl = self.I_Cl(V)

        plt.figure(figsize=(6,6))
        plt.subplot(3,1,1)
        plt.title('Noble model')
        plt.plot(self.t*1e-3, V, 'k')
        plt.ylabel('V (mV)')

        plt.subplot(3,1,2)
        plt.plot(self.t*1e-3, ina, 'c', label='$I_{Na}$')
        plt.plot(self.t*1e-3, ik, 'y', label='$I_{K}$')
        plt.plot(self.t*1e-3, icl, 'm', label='$I_{L}$')
        plt.ylabel('Current ($\mu$A/cm$^2$)')
        plt.legend()

        plt.subplot(3,1,3)
        plt.plot(self.t*1e-3, m, 'r', label='m')
        plt.plot(self.t*1e-3, h, 'g', label='h')
        plt.plot(self.t*1e-3, n, 'b', label='n')
        plt.ylabel('Gating Value')
        plt.xlabel('t (sec)')
        plt.legend()
        
        plt.tight_layout()
        plt.savefig("Noble_model.png")
        #plt.show()

if __name__ == '__main__':
    model = NobleModel()
    model.simulation()