-
Notifications
You must be signed in to change notification settings - Fork 1
/
in_estim5c.m
289 lines (264 loc) · 8.62 KB
/
in_estim5c.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
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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
# THIS SOFTWARE COMES WITH NO WARRANTY WHATSOEVER, EXPRESSED OR IMPLIED.
# USE IT AT YOUR OWN RISK!
#
# By T.I. Toth, Cardiff University, U.K.; 1996-2002
#
#
#
# This procedure computes the inital estimates of gNa, and the par.s of
# alpha_Na(V) and beta_Na(V) in the form, shown below.
# It uses the same order of Chebyshev approx. as during the computaion
# of dV/dt in in_estim4bc().
# Notation: tb=t(1)>tsb and te=t(length(t))<tse
# (tsb: stim. start, tse: stim. end)
#
#
# Input:
# t: time vector on [tb,te];
# v: sampled V(t) on [tb,te];
# ci0: ci0=IK+INa on [tb,te];
# mch: no. of polynom. terms in the Chebyshev approx. of V;
# ts1: sampled time on [tsb,tb];
# vs1: sampled V(t) on [tsb,tb];
# Er: resting pot.;
# hNa: hNa(t) on [tb,te];
# hNas1: hNa(t) on [tsb,tb];
# ENa: Na rev.pot.;
# gammn: pre-defined minimal value of gamma_Na(V)=1/taum_Na(V) (activ.rate);
#
# Output:
# ka: minimal index in vector t where ci0<0, index of tka;
# icmn: index of tcmn where ci0 is minimal (has a negative peak);
# gNa: estimated max. conductance of INa;
# alpNa0: vector of non-negative alpha_Na values on [tka,tcmn];
# valp: corresponding V vector;
# betNa0: vector of non-negative beta_Na values on [tka,tcmn];
# vbet: corresponding V vector;
# lambda: parameter in the solution to alpha_Na and beta_Na;
# alpNae: alpha_Na(V) estimated as ap*V1/(exp(V1)-1) where V1=alp1*V+alp2
# on V=v;
# p_alpNa: param. vector [ap,alp1,alp2];
# betNae: beta_Na(V) estimated as bp*V1/(exp(V1)-1) where V1=bet1*V+bet2
# on V=v;
# p_betNa: param. vector [bp,bet1,bet2];
# mNaes1: mNa(t) on [tsb,tb];
# mNass1: mNa_inf(V(t)) on [tsb,tb];
# gamNas1: gamma_Na(V(t)) on [tsb,tb];
# mNae: mNa(t) on [tb,te];
# mNas: mNa_inf(V(t)) on [tb,te];
# gamNa: gamma_Na(V(t)) on [tb,te];
# INas1: reconstructed INa(t) on [tsb,tb];
# INa: reconstructed INa(t) on [tb,te];
#
#
# Internal variables:
# ci0mn: minimum of ci0(t);
# tcmn: tcmn=t(icmn), end point of the estimation interval [tka,tcmn];
# tka: tka=t(ka) starting point of the estimation interval [tka,tcmn];
# N1: length of the selected sub-interval [tka,tcmn], N1=icmn-ka+1;
# mch2: the same as mch;
# a: a(t)=[INa/(hNa*(v-ENa))]^(1/3)=gNa^(1/3)*mNa(t) on [tka,tcmn];
# c_a: vector of Cheb. coeffs. of a;
# ach: Chebyshev approximate of a(t) on [tka,tcmn];
# av_err_a_cheb: `average' error of the Cheb approx. of a(t);
# Ha: coeff. matrix using c_a;
# D: derivative matrix of the Cheb. approx.;
# g0Na: g0Na=gNa^(1/3);
# c_alpNa0: vector of Cheb. coeffs. of alpha_Na0;
# c_betNa0: vector of Cheb. coeffs. of beta_Na0;
# G: coeff. matrix of polynom. lin. lsq.;
# p_alp: coeff vector [alp1,alp2];
# p_bet: coeff vector [bet1,bet2];
# sig, res: sigma val. and residual vector of the lin. least square proc.;
# m0: steady-state value of mNa at Er;
# m00: dummy var.;
#
#
# External procedures and functions:
# cheb_linip(): computes the Cheb. coeffs.;
# chebev_vect(): computes fct values from the Cheb. coeffs.;
# hgv3(): creates coeff. mtx from the Cheb coeffs.;
# d_cch(): derivative matrix of the Cheb. coeffs.;
# f_eq1(): nonlin.eqn. given as 1+b*x=exp(x) used in fsolve();
# fm1(): computes the actual values of the (in)activ. over a time vector
# t when alpha(V) and beta(V) are a*x/(exp(x)-1), where x=b1*V+b2,
# as well as the time course of m_inf(V(t)) and gamma(V(t))=1/tau;
function [ka,icmn,gNa,alpNa0,valp,betNa0,vbet,lambda,alpNae,p_alpNa,\
betNae,p_betNa,mNaes1,mNass1,gamNas1,mNae,mNas,gamNa,INas1,INa]=\
in_estim5c(t,v,ci0,mch,ts1,vs1,Er,hNa,hNas1,ENa,gammn)
N0=length(t); # no. of data points on [tb,te]
# Check consistency of the data:
if ((length(v) !=N0)||(length(ci0)!=N0)||(length(hNa)!=N0))
printf("data lengths: lt=%3d lv=%3d lci0=%3d lhNa=%3d\n",\
N0,length(v),length(ci0),length(hNa));
error("In in_estim5c: Data are inconsistent\n");
endif
# Find the negative peak of ci0:
ci0mn=min(ci0);
icmn=find(1+sign(ci0mn-ci0));
tcmn=t(icmn); # end point of the estim. interval
for k=icmn:-1:1
if (ci0(k)>=0) break; endif
endfor
ka=k+1;
tka=t(ka); # starting point of the estim. interval
N1=icmn-ka+1; # no. of points in the estim. interv.
mch2=mch; # the same as used for V;
a=ci0(ka:icmn)./(hNa(ka:icmn).*(v(ka:icmn)-ENa)); # gNa*mNa^3
a=a.^(1/3);
c_a=cheb_linip(t(ka:icmn),a, mch2); # Cheb.coeff.vect of a
ach=chebev_vect(tka,tcmn,c_a,t(ka:icmn))';
av_err_a_cheb=max(abs(a-ach))*length(a)/sum(a) # average error
Ha=hgv3(mch2-1,c_a);
Ha=Ha(1:mch2,:);
D=d_cch(tka,tcmn,mch2); # derivative matrix
dca=D*c_a;
g0Na=1.15*max(eig(Ha)); # 1st estimate of gNa^(1/3)
c_alpNa0=(g0Na*eye(mch2)-Ha)\dca;
c_betNa0=Ha\dca;
alpNa0=chebev_vect(tka,tcmn,c_alpNa0,t(ka:icmn))';
betNa0=chebev_vect(tka,tcmn,c_betNa0,t(ka:icmn))';
# Use only positive values for alpNa0 and betNa0 (select the corresponding
# values of v accordingly):
ka2=0; kb2=0;
for k=1:N1
al0=alpNa0(k);
if (al0>0)
ka2++;
alpNa0(ka2)=al0;
valp(ka2)=v(ka+ka2-1);
endif
be0=betNa0(k);
if (be0>0)
kb2++;
betNa0(kb2)=be0;
vbet(kb2)=v(ka+kb2-1);
endif
endfor
alpNa0=alpNa0(1:ka2);
betNa0=betNa0(1:kb2);
al0=0;
be0=0;
# Estimate par. values of alpha_Na and beta_Na;
# Coeff. matrix of the lin.sq. est. for alpha_Na:
G(:,2)=ones(size(valp));
G(:,1)=valp;
# First solve a nonlin.eqn. at each value of alpNa0:
global ap ri;
sq_err0=1e10;
dcf0=0.5*abs(min(alpNa0(2:ka2)-alpNa0(1:ka2-1)));
cf0mx=2*(max(alpNa0)+dcf0-al0)+1e-10;
cf0=0.7*(min(alpNa0)-al0);
while(cf0<cf0mx)
ap=cf0;
for ii=1:ka2
ri=alpNa0(ii)-al0;
if (ap==ri)
y0=0;
info=1;
elseif (ap>ri)
[y0,info]=fsolve("f_eq1",[10]);
else
[y0,info]=fsolve("f_eq1",[-ri/ap]);
endif
if (info!=1)
perror("fsolve",info)
printf("Error during estim. of alpNa0 at cf0=%5.2f\n",cf0);
break;
endif
ya(ii)=y0;
endfor
[p_alp,sig,res]=ols(ya,G);
# Compute the square error of the estimation
va1=polyval(p_alp,valp);
alpNa1=ap*va1./(exp(va1)-1);
sq_err1=sumsq(alpNa0-al0-alpNa1);
if (sq_err1<sq_err0)
p_alp0=p_alp;
sig_alp=sig;
res0=max(abs(res));
ap0=ap;
sq_err0=sq_err1;
cf0a=cf0;
endif
cf0=cf0+dcf0;
endwhile
printf("al0=%6.3f ap0=%6.3f alp1=%5.2f alp2=%5.2f\n",al0,ap0,p_alp0);
printf("sig_alp=%13.6g res=%13.6g\n",sig_alp,res0);
# Now solve a nonlin. eqn. at each value of betNa0:
# Re-define coeff. matrix of the lin.sq. est. for beta_Na:
clear G ya;
G(:,2)=ones(size(vbet));
G(:,1)=vbet;
sq_err0=1e10;
dcf0=0.5*abs(min(betNa0(2:kb2)-betNa0(1:kb2-1)));
cf0mx=2*(max(betNa0)+dcf0-be0)+1e-10;
cf0=0.7*(min(betNa0)-be0);
while(cf0<cf0mx)
ap=cf0;
for ii=1:kb2
ri=betNa0(ii)-be0;
if (ap==ri)
y0=0;
info=1;
elseif (ap>ri)
[y0,info]=fsolve("f_eq1",[10]);
else
[y0,info]=fsolve("f_eq1",[-ri/ap]);
endif
if (info!=1)
perror("fsolve",info)
printf("Error during estim. of betNa0 at cf0=%5.2f\n",cf0);
break;
endif
ya(ii)=y0;
endfor
[p_bet,sig,res]=ols(ya,G);
# Compute the square error of the estimation
va1=polyval(p_bet,vbet);
betNa1=ap*va1./(exp(va1)-1);
sq_err1=sumsq(betNa0-be0-betNa1);
if (sq_err1<sq_err0)
p_bet0=p_bet;
sig_bet=sig;
res0=max(abs(res));
bp0=ap;
sq_err0=sq_err1;
cf0b=cf0;
endif
cf0=cf0+dcf0;
endwhile
printf("be0=%6.3f bp0=%6.3f bet1=%5.2f bet2=%5.2f\n",be0,bp0,p_bet0);
printf("sig_bet=%13.6g res=%13.6g\n",sig_bet,res0);
# Compute the estimated alpha_Na0(V) and beta_Na0(V) for V=v:
va2=polyval(p_alp0,v);
alpNa0e=ap0*va2./(exp(va2)-1)+al0;
va2=polyval(p_bet0,v);
betNa0e=bp0*va2./(exp(va2)-1)+be0;
# Set lambda=10:
lambda=10;
alpNae=(1+lambda)*alpNa0e;
betNae=lambda*betNa0e;
# Compute new lambda such that gam0=gammn:
gam0=min(alpNae+betNae);
igam0=min(find(1+sign(gam0-(alpNae+betNae))));
# Re-define alphaNae and betaNae with the new lambda:
lambda=lambda+(gammn-gam0)/(alpNa0e(igam0)+betNa0e(igam0));
alpNae=(1+lambda)*alpNa0e;
betNae=lambda*betNa0e;
# Define par. vectors of alpNa_(V), beta_Na(V) for computing mNa(t):
p_alp2=[(1+lambda)*ap0;p_alp0;(1+lambda)*al0];
p_bet2=[lambda*bp0;p_bet0;lambda*be0];
p_alpNa=p_alp2(1:3); # usual par. set of alpha_Na(V)
p_betNa=p_bet2(1:3); # usual par. set of beta_Na(V)
# Compute mNa(t), mNa_inf(V(t)), and gamma_Na(V(t)) during ts1 and t:
[m00,m0,gamNa0]=fm1(0,Er,0,p_alp2,p_bet2);
[mNaes1,mNass1,gamNas1]=fm1(ts1,vs1,m0,p_alp2,p_bet2);
m0=mNaes1(length(ts1));
[mNae,mNas,gamNa]=fm1(t,v,m0,p_alp2,p_bet2);
# Compute the estimated INa:
m00=mNae.^3.*hNa.*(v-ENa);
gNa=ci0mn/min(m00);
INa=gNa*m00;
INas1=gNa*mNaes1.^3.*hNas1.*(vs1-ENa);
endfunction