-
Notifications
You must be signed in to change notification settings - Fork 24
/
vclamp.m
119 lines (103 loc) · 2.82 KB
/
vclamp.m
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
%Program for simulating full model with coupling
%importing data
data = importdata('data/input.csv'); %Not sure if these are the best paramaters...
% Model paramaters
Cmem = data(1);
gKS = data(2)*Cmem;
gKF = data(3)*Cmem;
gCa = data(4)*Cmem;
gL = data(5)*Cmem;
VKS = data(6);
VKF = data(7);
VCa = data(8);
VL = 10e-3;%data(9);
Vhalf_n = data(10);
Vhalf_p = data(11);
Vhalf_q = data(12);
Vhalf_e = data(13);
Vhalf_f = data(14);
Cahalf_h = data(15)*1e-9;
k_n = data(16);
k_p = data(17);
k_q = data(18);
k_e = data(19);
k_f = data(20);
k_h = data(21)*1e-9;
T_n = data(22);
T_p = data(23);
T_q = data(24);
T_e = data(25);
T_f = data(26);
T_Ca = data(27);
alphaCa = data(28);
thiCa = 6.1e-6/(T_Ca*gCa);
%Simulation Paramaters
deltat = 0.01e-3;
duration = 0.03; %********************* Duration Duration ***************
numpoints = round(duration/deltat);
numtests = 11;
xaxis = (deltat:deltat:duration)*1e3;
%Input paramaters
onset = round(0.002/deltat);
offset = round(0.022/deltat);
%Variable Declaration
V = zeros(numtests,numpoints);
I_j = zeros(numtests+1,1);
I_mem = zeros(numtests,numpoints);
Ca = zeros(numtests,numpoints);
n = zeros(numtests,1);
p = zeros(numtests,1);
q = zeros(numtests,1);
e = zeros(numtests,1);
f = zeros(numtests,1);
h = zeros(numtests,1);
%input initialization
for i = 1:numtests
for j = 1:numpoints
V(i,j) = -70e-3;
end
end
Vstim = 40e-3;
for i = 1:numtests
for j = onset:offset
V(i,j) = Vstim;
end
Vstim = Vstim - 10e-3;
end
%variable initialization
for j = 1:numtests
Ca(j,1) = 0;
n(j,1) = x_inf(V(j,1),Vhalf_n,k_n);
p(j,1) = x_inf(V(j,1),Vhalf_p,k_p);
q(j,1) = x_inf(V(j,1),Vhalf_q,k_q);
e(j,1) = x_inf(V(j,1),Vhalf_e,k_e);
f(j,1) = x_inf(V(j,1),Vhalf_f,k_f);
end
%start of simulation
for j = 1:numtests
for i = 2:numpoints
dn = (x_inf(V(j,i-1),Vhalf_n,k_n) - n(j,1))/T_n;
n(j,1) = n(j,1) + dn*deltat;
dp = (x_inf(V(j,i-1),Vhalf_p,k_p) - p(j,1))/T_p;
p(j,1) = p(j,1) + dp*deltat;
dq = (x_inf(V(j,i-1),Vhalf_q,k_q) - q(j,1))/T_q;
q(j,1) = q(j,1) + dq*deltat;
de = (x_inf(V(j,i-1),Vhalf_e,k_e) - e(j,1))/T_e;
e(j,1) = e(j,1) + de*deltat;
df = (x_inf(V(j,i-1),Vhalf_f,k_f) - f(j,1))/T_f;
f(j,1) = f(j,1) + df*deltat;
h(j,1) = x_inf(Ca(j,i-1),Cahalf_h,k_h);
H_rec(1,i) = h(j,1);
IKS = gKS*n(j,1)*(V(j,i-1) - VKS);
IKF = gKF*p(j,1)^4*q(j,1)*(V(j,i-1) - VKF);
ICa = gCa*e(j,1)^2*f(j,1)*(1 + (h(j,1) - 1)*alphaCa)*(V(j,i-1) - VCa);
IL = gL*(V(j,i-1) - VL);
dCa = -(Ca(j,i-1)/T_Ca + thiCa*ICa);
Ca(j,i) = Ca(j,i-1) + dCa*deltat;
I_mem(j,i) = (IKS + IKF + ICa);
end
plot(xaxis,I_mem(j,:).*1e9,'k')
hold on
end
ylabel('I_{mem} (nA)')
xlabel('Time (ms)')