# 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 parameter estimates for INa, i.e.
# mNa_inf(V) as product of a steep Boltzmann curve and a polynomial
# of degree ng1<9 in V, and gamma_Na(V) or 1/gamma_Na(V) as
# a polynomial of degree ng0<9 in V whichever yields the better fit.
# Both gNa and gK are optimized by lin. lsq.
# In this version (df), IK is increased such that an "a-loop" appears
# with no self-intersection in the V-a plane. This is ensured by
# checking that a(t) values on the descending phase of the AP
# are greater than the corresponding ones on the ascending phase.
# Notation: tb=t(1)=ts1(length(ts1)), te=t(length(t)),
# tsb=ts1(1): stim. start.
#
#
# Input:
# t: time vector on [tb,te];
# v: V(t) on [tb,te];
# ci0: ci0=IK+INa on [tb,te];
# ts1: time from stim. start to AP: [tsb,tb];
# vs1: corresponding V vector;
# Er: resting pot.;
# ENa: Na rev.pot.;
# mNax0: (approx.) max. value of mNa_inf(V);
# ng0: degree of the polynom. for gamma_Na(V) or 1/gamma_Na(V);
# ng1: degree of the polynom. for maNa_inf(V);
# hNa: hNa(t) on [tb,te];
# hNas1: hNa(t) on [tsb,tb];
# gK: max. conductance of IK;
# IK: IK on [tb,te];
#
#
# Output:
# ivmx: index of the peak (max.) of v;
# inmx: index of t=tmx where mNa is maximal (has the largest positive peak);
# nne: index of the last point of the v-interval with v(nne) appr.=v(1);
# mNaa: 'original' mNa on [tb,te];
# vng: a subset of v(1:nne) where both mNa_inf(V) and gamma_Na(V)
# are computed and positive;
# tng: time vector corresponding to vng;
# mNasa: computed values of mNa_inf(V) on vng;
# gamNaa: computed values of gamma_Na(V) on vng;
# p_gamNa: coeff. vector of the polynom. estimate for gamma_Na(V)
# or 1/gamma_Na(V); the last entry of p_gamNa shows which case
# was accepted; +1: gamma_Na(V), -1: 1/gamma_Na(V);
# p_mNa: coeff. vector of the polynom. estimate for mNa_inf(V);
# q1: `slope' factor of the Boltzmann curve of mNa_inf(V);
# V0: Vhalf of the Boltzmann curve of mNa_inf(V);
# gNa: optimized value of gNa;
# da: da/dt where a=g0Na*mNa(t);
# mch2: order+1 of the Cheb. approx used to compute da;
# mNaes1: mNa(t) computed on [tsb,tb];
# mNass1: mNa_inf(V(t)) computed on [tsb,tb];
# gamNas1: gamma_Na(V(t)) computed on [tsb,tb];
# mNae: mNa(t) computed on [tb,te];
# mNas: mNa_inf(V(t)) computed on [tb,te];
# gamNa: gamma_Na(V(t)) computed on [tb,te];
# INas1: reconstructed INa(t) on [tsb,tb];
# INa: reconstructed INa(t) on [tb,te];
# gKe: new, optimized value of gK;
# IKe: new IK on [tb,te] computed with gKe;
# norm2_curr: quadratic error of ci0-IKe-INa;
# final_sig: min. sigma of the fit of gNa and gKe;
# final_err: maximal error of the fit;
#
#
# Internal variables:
# ci0mx1: limiting current value for IK;
# tvmx: t(ivmx) where V is maximal;
# bK: mK^4*(v-EK);
# dt1: sampling time interval, dt1=t(2)-t(1);
# thst: (varied) min. distance between t(inmx) and t(ivmx);
# thstmx: max. value of thst;
# cf: factor used to increase IK;
# dcf: (preset) increment of cf;
# ramn: relative min. distance of corresponding points in the "a-loop";
# atol: threshold for ramn;
# sig1: preceding error value for comparison with the actual one;
# g0Na: g0Na=gNa^(1/3);
# a: g0Na*mNa;
# mch2x: max. admissible order+1 of Cheb. approx.;
# t1: t(1:ivmx-1);
# v1: v(1:ivmx-1);
# t2: t(nne:-1:ivmx+1);
# v2: v(nne:-1:ivmx+1);
# a1,a2,da1,da2: sub-vectors of a and da corresponding to t1 and t2, resp.;
# Ha: coeff. matrix of the lin. eq.sys. for gamma(V) and g0Na*alpha(V)
# at each pair of identical V-values from v1 and v2;
# aNa0: g0Na*mNa_inf(V) as obtained from the solution to the lin. eq.sys.;
# G: coeff. matrix of the lsq. polyn. approx. for mNa_inf(V) and gamma_Na(V);
# G1: coeff. matrix of the lin. least square proc. for gNa and gK;
# sig, res: sigma and residual vector of the least square proc.s;
# inmx0,vng0,tng0 etc. : store the last `best' results;
#
#
# External procedures and functions:
# fm2df(): computes m(t), m_inf(V(t)), and gamma(V(t)) when m_inf(V) is
# a product of a Boltzmann curve and a polynomial in V, and
# gamma(V) or 1/gamma(V) is a polynomial in V.
# minf(): Boltzmnann curve;
function [ivmx,inmx,nne,mNaa,vng,tng,mNasa,gamNaa,p_gamNa,p_mNa,q1,V0,gNa,\
da,mch2,mNaes1,mNass1,gamNas1,mNae,mNas,gamNa,INas1,INa,gKe,\
IKe,norm2_curr,final_sig,final_err]=\
tc_estim7df(t,v,ci0,ts1,vs1,Er,ENa,mNax0,ng0,ng1,hNa,hNas1,gK,IK)
N0=length(t); # no. of data points on [tb,te]
# Check consistency of the data:
if ((length(v) !=N0)||(length(ci0)!=N0)||(length(IK)!=N0))
printf("data lengths: lt=%3d lv=%3d lci0=%3d lIK=%3d\n",\
N0,length(v),length(ci0),length(IK));
error("In tc_estim7df: Data are inconsistent\n");
endif
# Set and compute some test values for comparisons:
ci0mx1=2.0*max(abs(ci0)); # bound for IK
ivmx=max(find(1-sign(max(v)-v))); # index of Vmax
tvmx=t(ivmx); # where V is maximal
inmx=ivmx;
inmx0=0; # no useful results yet
bK=IK/gK; # mK^4*(v-EK)
atol=0.02; # minimal relative distance
dcf=0.02; # increment of cf (see below);
# Split the fct. V(t) into two monotonic (ascend. and descend.) parts:
V1=v(1);
nne=min(find(1+sign(V1-v(ivmx:N0))))+ivmx-1;
v1=v(1:ivmx-1);
t1=t(1:ivmx-1);
v2=v(nne:-1:ivmx+1);
t2=t(nne:-1:ivmx+1);
dt1=t(2)-t(1);
thst=2*dt1; # separation time between t_vmax and t_mNa_max
thst0=0;
thstmx=6*dt1+1e-6;
sig1=1e20;
final_sig=sig1;
final_err=sig1;
# Start a `big loop' which chooses the optimal thst=t(inmx)-t(ivmx)>0
while (thst<thstmx) # start of the `thst' loop
cf=1; # initial values
ramn=0;
# Test whether the peaks of V and mNa are well separated (Vmax must come first)
# and whether there is a proper "a-loop" (see above):
while ((t(inmx)-tvmx<thst)||(ramn<atol))
# Compute INa:
INa=ci0-cf*IK;
a=max(INa./(hNa.*(v-ENa)),0);
a=a.^(1/3);
# Find the largest positive peak of a:
mxa=max(a);
inmx=find(1-sign(mxa-a)); # index of mNa_max
# Exit loop if IK grows too large:
if (ci0mx1<max(cf*IK)) break; endif
cf=cf+dcf; # increase factor of IK
a1=a(1:ivmx-1);
a2=a(nne:-1:ivmx+1);
for ii=1:ivmx-1
vb=v1(ii);
a1a=a1(ii);
jj=0;
jj=min(find(1+sign(v2-vb)));
if (jj>1)
v2b=v2(jj);
a2b=a2(jj);
jj--;
v2a=v2(jj);
a2a=a2(jj);
else
ra(ii)=1e19;
continue;
endif
a1b=(vb-v2a)*(a2b-a2a)/(v2b-v2a)+a2a;
ra(ii)=(a1b-a1a)/mxa;
endfor
ramn=min(ra);
endwhile # end of ramn<atol loop
if (ci0mx1<max(cf*IK))
if (inmx0==0) # no useful results yet
printf("max V at %5.2f max mNa at %5.2f\n",t(ivmx),t(inmx));
printf("Inactivation of INa is conflicting with the data.\n");
printf("max_ci0=%13.5g max_IK=%13.5g min_INa=%13.5g\n",
max(ci0),max(IK),min(INa));
printf("In tc_estim7df: IK is too large\n");
endif
break; # exit `thst' loop, too!
endif
# Compute the actual diff. between the peaks of mNa and V:
thst=t(inmx)-tvmx;
# Find the corresponding interpolated values of a(t(V)) and da(t(V)):
mch2x=10;
for mch2=4:mch2x
da=df_ch_vect(t,a,mch2)'; # da/dt
if (da(inmx-1)*da(inmx+1)<0) break; endif
endfor
da1=da(1:ivmx-1); # da/dt on t1
da2=da(nne:-1:ivmx+1); # da/dt on t2
ng=0;
for ii=1:ivmx-1
vb=v1(ii);
a1a=a1(ii);
da1a=da1(ii);
jj=0;
jj=min(find(1+sign(v2-vb)));
if (jj>1)
v2b=v2(jj);
a2b=a2(jj);
da2b=da2(jj);
jj--;
v2a=v2(jj);
a2a=a2(jj);
da2a=da2(jj);
else
ii
continue;
endif
a1b=(vb-v2a)*(a2b-a2a)/(v2b-v2a)+a2a;
da1b=(vb-v2a)*(da2b-da2a)/(v2b-v2a)+da2a;
Ha=[1,-a1a;1,-a1b];
a_gam=Ha\[da1a;da1b];
gam=a_gam(2);
# Allow only positive solutions:
if ((gam>0)&&(a_gam(1)>0))
ng++;
vng(ng)=vb;
tng(ng)=t(ii);
aNa0(ng)=a_gam(1)/gam;
gamNaa(ng)=gam;
endif
endfor
for ii=ivmx+1:nne
vb=v(ii);
a2a=a(ii);
da2a=da(ii);
jj=0;
jj=min(find(1+sign(v1-vb)));
if (jj>1)
v1b=v1(jj);
a1b=a1(jj);
da1b=da1(jj);
jj--;
v1a=v1(jj);
a1a=a1(jj);
da1a=da1(jj);
else
ii
continue;
endif
a2b=(vb-v1a)*(a1b-a1a)/(v1b-v1a)+a1a;
da2b=(vb-v1a)*(da1b-da1a)/(v1b-v1a)+da1a;
Ha=[1,-a2a;1,-a2b];
a_gam=Ha\[da2a;da2b];
gam=a_gam(2);
# Allow only positive solutions:
if ((gam>0)&&(a_gam(1)>0))
ng++;
vng(ng)=vb;
tng(ng)=t(ii);
aNa0(ng)=a_gam(1)/gam;
gamNaa(ng)=gam;
endif
endfor
# Estimate gamma(V) or 1/gamma_Na(V) as a polynomial of degree ng0<9 in V:
# Usually ng0=8 or ng0=7;
clear G;
G(:,ng0+1)=ones(size(vng));
G(:,ng0)=vng;
for ng=ng0:-1:2
G(:,ng-1)=G(:,ng).*vng;
endfor
# gamma(V):
[p_gamNa1,sig,res]=ols(gamNaa,G);
p_gamNa1=[p_gamNa1;1];
# 1/gamma(V):
[p_gamNa2,sig,res]=ols(1./gamNaa,G);
p_gamNa2=[p_gamNa2;-1];
# Compute g0Na and mNa:
g0Na=max(aNa0)/mNax0; # max(mNa_inf(V)) is approx. mNax0
mNaa=a/g0Na; # `original' mNa(t)
mNasa=aNa0/g0Na; # computed mNa_inf(V) on vng
# Estimate par.s of the steep Boltzmann curve of mNa_inf(V):
V1=v(1);
V2=v(4);
q1=log((1/mNax0-1)/(1/mNasa(1)-1))/(V2-V1);
V0=log(1/mNasa(1)-1)/q1-V1;
# Estimate mNa_inf(V) as a polynomial of degree <9 in V:
# Usually ng1=8 or ng1=7;
clear G;
G(:,ng1+1)=ones(size(vng));
G(:,ng1)=vng;
for ng=ng1:-1:2
G(:,ng-1)=G(:,ng).*vng;
endfor
[p_mNa,sig,res]=ols(mNasa,G);
sig_mNa=sig
max_err=max(abs(res));
# Compute mNa(t), mNa_inf(V(t)), and gamma_Na(V(t)) on [tsb,tb] and [tb,te]
# with p_gamNa1:
bad_idx1=0;
[mNaes1,mNass1,gamNas1]=fm2df(ts1,vs1,0.,p_mNa,q1,V0,p_gamNa1);
m0=mNaes1(length(ts1));
[mNae,mNas,gamNa]=fm2df(t,v,m0,p_mNa,q1,V0,p_gamNa1);
# Estimate gNa and gK if possible:
if ((max(mNae)>=1.0)||(min(mNae)<0.)||\
(max(mNas)>=1.0)||(min(mNas)<0.)||(min(gamNa)<0.))
bad_idx1=1;
siga=1e18;
else
# Optimize gK and gNa:
clear G1;
G1(:,1)=bK;
G1(:,2)=mNae.^3.*hNa.*(v-ENa);
[gz1,siga,res]=ols(ci0,G1);
siga
max_erra=max(abs(res))
endif
# Compute mNa(t), mNa_inf(V(t)), and gamma_Na(V(t)) on [tsb,tb] and [tb,te]
# with p_gamNa2:
bad_idx2=0;
[mNaes1,mNass1,gamNas1]=fm2df(ts1,vs1,0.,p_mNa,q1,V0,p_gamNa2);
m0=mNaes1(length(ts1));
[mNae,mNas,gamNa]=fm2df(t,v,m0,p_mNa,q1,V0,p_gamNa2);
# Estimate gNa and gK if possible:
if ((max(mNae)>=1.0)||(min(mNae)<0.)||\
(max(mNas)>=1.0)||(min(mNas)<0.)||(min(gamNa)<0.))
bad_idx2=1;
sigb=1e18;
else
# Optimize gK and gNa:
clear G1;
G1(:,1)=bK;
G1(:,2)=mNae.^3.*hNa.*(v-ENa);
[gz2,sigb,res]=ols(ci0,G1);
sigb
max_errb=max(abs(res))
endif
# Ignore meaningless results or save good ones:
if ((bad_idx1)&&(bad_idx2))
thst=thst+dt1;
continue;
elseif (siga<sigb)
sig=siga;
max_err=max_erra;
gz=gz1;
p_gamNa=p_gamNa1;
else
sig=sigb;
max_err=max_errb;
gz=gz2;
p_gamNa=p_gamNa2;
endif
if (sig1>sig)
inmx0=inmx;
vng0=vng;
tng0=tng;
gamNaa0=gamNaa;
p_gamNa0=p_gamNa;
mNaa0=mNaa;
mNasa0=mNasa;
p_mNa0=p_mNa;
q10=q1;
V00=V0;
gz0=gz;
final_err=max_err;
final_sig=sig;
thst0=thst;
da0=da;
mch20=mch2;
sig1=sig;
endif
thst=thst+dt1;
endwhile # end of the `thst' loop
thst0
final_sig
final_err
# Set outputs to zero if no meaningful result:
if (thst0==0)
mNaa=0;vng=0;tng=0;mNasa=0;gamNaa=0;p_gamNa=0;p_mNa=0;
q1=0;V0=0;da=0;mch2=0;mNaes1=0;mNass1=0;gamNas1=0;
mNae=0;mNas=0;gamNa=0;
norm2_curr=final_err;
gKe=0; gNa=0; INa=0; IKe=0; INas1=0;
printf("No meaningful results!\n");
return;
endif
# Restore the `optimal' values:
gKe=gz0(1);
gNa=gz0(2);
inmx=inmx0;
vng=vng0;
tng=tng0;
gamNaa=gamNaa0;
p_gamNa=p_gamNa0;
mNaa=mNaa0;
mNasa=mNasa0;
p_mNa=p_mNa0;
q1=q10;
V0=V00;
da=da0;
mch2=mch20;
# Compute final mNa(t), mNa_inf(V(t)),gamma_Na(V(t)), INa(t) on [tsb,tb]
# and [tb,te] with best p_gamNa:
[mNaes1,mNass1,gamNas1]=fm2df(ts1,vs1,0.,p_mNa,q1,V0,p_gamNa);
m0=mNaes1(length(ts1));
[mNae,mNas,gamNa]=fm2df(t,v,m0,p_mNa,q1,V0,p_gamNa);
INa=gNa*mNae.^3.*hNa.*(v-ENa);
IKe=gKe*bK;
INas1=gNa*mNaes1.^3.*hNas1.*(vs1-ENa);
# Quadratic error:
norm2_curr=norm(ci0-IKe-INa)
endfunction