# 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 estimates: Cm, gL and EL for neurones when
# hyperpolarizing current steps are applied to them. In this version (2a),
# it is assumed that the leakage and the capacitive current completely
# determine the time evolution of the membrane pot. during the
# hyperpolarizing current steps.
# The estimation procedure is the lin. least square (LLSQ) method.
# It is assumed that a constant holding current Ihold is applied to keep
# the cell in steady-state at the membrane potential Eh. Then the
# equilibrium condition is
# gL*(Eh-EL)=Ihold. (1)
# During the hyperpolarizing step
# Cm*dV/dt+gL*(V-EL)=Ihold+Istep (2)
# applies. Subtracting (1) from (2), we get the equations for the LLSQ
# method, namely one eqn. at each sampling time point. dV/dt is computed
# (approximated) via the Chebyshev approximation of V(t).
#
# The neuronal input data is called burst, and its structure is given below.
# In this version (2a), a correction value for a possible d.c. shift of
# the input currents must be supplied, in order to obtain the unbiased value
# of Ihold. The input current can also be re-scaled by the scaling factor
# cifc (cf. Input below).
#
#
# Input:
# burst: data matrix whose 1st column is the sampled time; the subsequent
# K0 columns are the sampled voltage at different hyperpol. stimuli;
# and the last K0 columns are the sampled stimuli corresponding to
# the voltage traces, i.e. columns i and i+K0 correspond to each
# other (i>1). The no. of rows is the length of a data trace. The
# voltage values are given in mV, those of the current in nA;
# tsb: start time of the stimulus; all time data are given in ms;
# tse: end of the stimulus;
# tb: start time of the data segment used in the estimation, tb>tsb;
# te: end time of the data segment used in the estimation, te<tse;
# k0: the data column of a voltage trace selected (1<k0<K0+2);
# mch: no. of coeff.s in the Chebyshev approx.;
# dc_shift: d.c. shift of the current signal (nA); This is usually the mean
# value of a current trace with no stimulus and no holding current;
# cifc: correction factor for scaling the input current ci2.
#
# Output:
# Cm: membrane capacitance (pF);
# gL: leakage conductance (nS);
# cih: holding current applied (pA);
# Eh: holding potential (mV) at cih;
# EL: rev. pot. (mV) of the leakage current IL;
# Er: resting potential (mV), Er=EL;
# ci2: the constant stimulus current (pA);
# t: vector of the selected sampled time (ms) on [tb,te];
# v: vector of the selected sampled voltage data V corresponding to t (mV);
# vch: Chebyshev approx. vector of the voltage trace V (mV) on [tb,te];
# dv: Chebyshev approx. vector of dV/dt (mV/ms) on [tb,te];
# IL: leakage curr. (pA) on [tb,te].
# In addition sigma and the max. absolute error errmax=max(abs(Y-X*[Cm,gL]'))
# of the LLSQ estimation are displayed (cf. Octave's User's Guide).
#
# Internal variables:
# Y: right hand side of the lin. equation;
# X: coeff. matrix for Cm and gL.
#
#
function [Cm,gL,cih,Eh,EL,Er,ci2,t,v,vch,dv,IL]=in_estim2a(burst,\
tsb,tse,tb,te,k0,mch,dc_shift,cifc)
N0=rows(burst); # no. of data point in a trace
# First checking consistency of the data:
if (N0<50) error("In in_estim2a: too few data points!\n") endif
if ((tsb>=tse)||(tsb<burst(1,1))||(tse>burst(N0,1)))
printf("tsb=%f\t tse=%f\n",tsb,tse);
error("In in_estim2a: tsb or tse incorrect!\n")
endif
if ((tb<tsb)||(te>tse))
printf("tsb=%f\t tse=%f\t tb=%f\t te=%f\n",tsb,tse,tb,te);
error("In in_estim2a: tb or te is incorrect!\n")
endif
k1=floor((columns(burst)-1)/2); # no. of col.s of voltage data
k2=columns(burst)-2*k1-1;
if (k2>0) error("In in_estim2a: data are incomplete!\n") endif
if ((k0==1)||(k0>k1+1))
k0
error("In in_estim2a: k0 is incorrect!\n")
endif
if (mch<3)
mch
error("In in_estim2a: order of Chebyshev approx. is too low!\n")
endif
# Compute Eh (holding pot.) as average pre-stim. potential:
for n0=1:N0
if (tsb<=burst(n0,1)) break; endif # start time of stimulus
endfor
Eh=mean(burst(10:n0-20,k0)); # avoid transients
# Select data segment for the LLSQ, and find the end of the stimulus:
for n1=n0:N0
if (tb<=burst(n1,1)) break; endif # start time of selected data
endfor
for n2=n1:N0
if (te<=burst(n2,1)) break; endif # end time point of selected data
endfor
for n3=n2:N0
if (tse<=burst(n3,1)) break; endif # end of stimulus
endfor
N=n2-n1+1; # no. of data points in the selected segment
t=burst(n1:n2,1);
v=burst(n1:n2,k0);
# Compute the stim. curr. ci2 as average of the sampled values
# in [tsb,tse], and the pre-stimulus holding curr. cih
# avoiding transients:
ci2=mean(burst(n0+10:n3-10,k0+k1));
cih=mean(burst(10:n0-20,k0+k1));
ci2=1000*(ci2-cih); # in pA
# Correct the pre-stimulus holding current cih for dc shift:
cih=1000*(cih-dc_shift); # in pA
# Re-scale the input current if necessary:
if ((cifc !=1)&&(cifc !=0))
ci2=cifc*ci2;
cih=cifc*cih;
endif
# Compute the Chebyshev approx. of the voltage data and that of their
# time derivative dV/dt:
ck=cheb_linip(t,v,mch); # Chebyshev coeff.s
vch=chebev_vect(t(1),t(N),ck,t)';
dv=df_ch_vect(t,v,mch)';
# LLSQ procedure:
Y=ci2*ones(N,1);
X=zeros(N,2);
# The following expression contains the equlibrium condition at Er:
X(:,1)=dv;
X(:,2)=v-Eh;
[P,sig,R]=ols(Y,X);
sig
errmax=max(abs(R)) # max. abs. error of the LLSQ estim.
# Computing the output var.s:
Cm=P(1);
gL=P(2);
# Other output:
EL=Eh-cih/gL;
Er=EL;
IL=gL*(v-EL);
endfunction