仿真HH方程

发表于 |
字数统计: 734 | 阅读时长 ≈ 4

仿真恒流刺激

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
import numpy as np
import math
import matplotlib.pyplot as plt
plt.figure(figsize=(10,10))
result = []

DT = 0.001
t = np.linspace(0,50,int(50/DT)+1)   # 仿真时间   光 1200
Cm = 1   # 膜电容
I = 90   # 刺激电流 ?????电流不能太大
v = -112 # 静息电位
GK = 24 # GK=3.5;
GCl = 36 # GCl=16;
GCl_slow = 1
GKin = 20
GH = 0.28
n2 = 0
m2 = 0
nn2 = 0
s2 = 0
h2 = 1
k2 = 0
v2 = v

T0 = 25
Tinitial = 18
Tend = 4
Iin0 = 0.005 # Iin0=0.0193;
Iinmax = 3 # Iinmax=2.5;
K1 = 1
n_t1 = 3
Iexmax = 1
Km = 0.5
n_t2 = 2
Q = 10
P0 = 0.005
Kp = 0.5
Ca0 = 0.1   # 细胞质Ca离子浓度,uM
Ca2 = Ca0
KQ = math.log(Q)/10 
T1 = Tinitial
Rate = -0.8 # 摄氏度/秒
tc = 45 # 秒

for i in range(1,len(t)):
    #---------------------------------------------------------------------
    Ca1 = Ca2;
    T2 = T1+Rate*DT*math.exp(-t[i]/tc)
    temp = (-(Rate*math.exp(-t[i]/tc)-abs(Rate*math.exp(-t[i]/tc)))/2)**n_t1
    dCadt = Iin0+Iinmax*temp/(K1**n_t1+temp) - Iexmax*math.exp(KQ*(T1-T0))*(Ca1**n_t2)/(Km**n2+Ca1**n_t2)
    Ca2 = Ca1+DT*dCadt
    dIexmax = DT*P0*np.sign(Ca1-Ca0)*Ca1/(Kp+Ca1)
    Iexmax = Iexmax+dIexmax
    Itemperature = -0.6823*dCadt*450  # 单位没有转换 Itemperature=-0.0154*dCadt/0.00005;
    T1 = T2
    #---------------------------------------------------------------------
    n1 = n2
    m1 = m2
    h1 = h2 
    s1 = s2
    k1 = k2
    nn1 = nn2
    v1 = v2
    V = v1
    Ib = 2*(V+127.5) #   Ib=2*(V+130.5);
    am = 10.55*(V+60)/(1-math.exp(-(V+60)/7.029))
    bm = 44.32*math.exp(-V/98.2)
    ah = 38.35*math.exp(-V/41.39)
    bh = 7.249/(1+0.6061*math.exp(-V/26.58))
    an = (0.01812*V+2.598)/(1+0.5954*math.exp(-V/10.8))
    bn = 1.56*math.exp(-V/23.4) 
    AS = 0.02985*math.exp(V/144.6)   #慢Cl
    bs = 0.03542*math.exp(-V/91.67)
    ak = 0.01414*math.exp(-0.03175*V)/(1+math.exp(0.2434*V))#  内向K
    bk = 974.1*math.exp(0.04129*V)/(1+math.exp(0.7851*V))
    ann = 0.01*(V+60)/(1-math.exp(-(V+60)/2.8))
    bnn = 0.05*math.exp(-V/80.4)
    if i==1:
        n1 = an/(an+bn)
        m1 = 0 
        h1 = 1
        s1 = AS/(AS+bs)
        k1 = ak/(ak+bk) 
        nn1 = ann/(ann+bnn)
  
    N = n1
    M = m1
    H = h1
    S = s1  
    NN = nn1
    gK = GK*N*N
    gCl = GCl*M*M*H
    gCl_slow = GCl_slow*S; 
    # gKN=GK*NN*NN;% 复极化电导
    gKN=40*NN*NN # 最大电导40可调整，调大后，在较大的刺激下也可出现正常的波形
    gH = GH/(1+math.exp((65.53+V)/112.2))
    IH = gH*(V+231.8)

    Ik = gK*(V+53)
    ICl = gCl*(V-13.6)  # ICl=gCl*(V+23.6);
    ICl_slow = gCl_slow*(V-37.8)
    IKN = gKN*(V+115) #  IKN=gKN*(V+118);
    K = k1  # 内向K电流
    gKin =GKin*K 
    Ikin = gKin*(V+75) 
    kdot = (ak*(1-K)-bk*K)*DT
    k2=k1+kdot
    
    ndot = (an*(1-N)-bn*N)*DT
    nndot = (ann*(1-NN)-bnn*NN)*DT
    # IT=Ik+ICl+ICl_slow+Ikin+Ib+IH;
    IT=Ik+ICl+ICl_slow+IKN+Ikin+Ib+IH
    

#     if t[i]<800          %光刺激
#         Ilight=71*(1-math.exp(-t[i]/82))^3;
#     else
#         Ilight=71*math.exp(-(t[i]-800)/56);
#         #if rem(t[i],20)==0
#             #v1
#         #end
#     end
#     vdot=((Ilight-IT)/Cm)*DT;
   

    if t[i]<100:    #温度、恒流电刺激
        vdot = ((I-IT)/Cm)*DT
    #    vdot=(-(Itemperature+IT)/Cm)*DT;
    else:
        vdot = (-IT/Cm)*DT
    
       
    v2 = v1+vdot
    mdot = (am*(1-M)-bm*M)*DT
    n2 = n1+ndot
    m2 = m1+mdot
    sdot =(AS*(1-S)-bs*S)*DT
    s2 = s1+sdot
    nn2 = nn1+nndot
    if V<-20:
        h2=1
    else:
        hdot = (ah*(1-H)-bh*H)*DT
        h2=h1+hdot 
    
    '''
    if t[i]%0.01==0:
        result.append([t[i],v1])
    '''    
    
    result.append([t[i],v1])    
result = np.array(result)
plt.plot(result[:,0],result[:,1],c='b')
#plt.plot(result[:,0],result[:,2],'r-');  #仿真Ca离子浓度随时间的变化关系
plt.xlabel('时间(s)',fontproperties='SimHei')
plt.ylabel('膜电位(mV)',fontproperties='SimHei')