function [Rtotal, Rseries, Rmem, Cmem, Ileak, current] = testpulse(Fs, scaling, gain, amplitude, pulsewidth, nrepeat)
%
% [Rtotal, Rseries, Rmem, Cmem, Ileak, current] =
% testpulse(Fs, scaling, gain, amplitude, pulsewidth, nrepeat)
% testpulse -- meant to be called from axontp
%
% Returned values:
% All resistances are returned in Ohms. Cmem is in Farads. Ileak
% and current are in pA.
%
% If nrepeat is > 1, then all returned values are averaged over
% nrepeat trials.
%
% Revision history:
% Implemented Rm, Cm calculation from Axon diagram 5/2003 (Neville),
% Cleaned up, added 2nd channel 11/2003 (Naveen),
% Fixed Ra mislabeling, fixed bugs with nrepeat 2/2004 (Neville)
% The STIM voltage array will be what we send through vstim.
% We make it here. We will send the same pulse to both the
% output channels.
pulse=amplitude*[zeros(pulsewidth,1); -ones(pulsewidth,1); zeros(pulsewidth,1)];
stim=[pulse pulse];
dV = amplitude * 10^-3;
avgCurrent = zeros(length(pulse),2);
Rtotal=zeros(1,2);
Rseries=zeros(1,2);
Rmem=zeros(1,2);
Cmem=zeros(1,2);
Ileak=zeros(1,2);
% In this loop, we send out the voltage pulse, and calculate
% the different parameters based on the output current.
for irepeat = 1:nrepeat
current=vstim(stim, scaling, gain, Fs, 1); % current in pA
% The bulk of the calculations below are based on figures and
% info in Axon Knowledge base article #652, found at:
% http://www.axon.com/mr_Axon_KB_Article.cfm?ArticleID=652
% As in Axon, I1 is steady state current during pulse and I2 is
% the steady state current before pulse
I1=mean(current(350:390,:));
I2=mean(current(80:160,:));
% We will calculate four variables: rt, rm, rs, and cm, for each
% of the two channels and then assign them to our returning variables
for ichannel = 1:2
Ileak(ichannel) = Ileak(ichannel) + mean(current(80:160,ichannel));
% We want to prevent "divide-by-zero" kind of errors, so we make sure
% the denominator isn't zero.
denom = mean(current(80:160,ichannel)) - mean(current(280:360, ichannel));
if denom == 0
rt = 0;
else
rt = amplitude / denom;
end
% Because amplitude is in mV, current is in pA, this returns
% Rtotal in gOhms. We will convert it to Ohms below.
% p will have a polynomial fit to the first little bit of the
% reponse to the downward step (which started at element 201
% of the current array)
p=polyfit(2:4,current(202:204,ichannel)',1); % linear extrapolation
p=polyfit(3:8,current(203:208,ichannel)',2); % quadratic extrapolation
% Convert measured Rt to Ohms
rt = rt * 10 ^9;
% Q1, Q2 and Qt are from the Axon diagram for Membrane Test
% Trying to find all points on the capacitive transient to
% calculate Q1
%
% Not sure if this is really the principled way to do this??
% We are trying to get all points where current > the steady
% steady state current during the pulse
% ASSUMPTION: Capacitive transient lasts until at least 202
% (true for all cells)
Q1_pts1 = find(current(202:400,ichannel) >= I1(ichannel));
numPointsAfterQ1 = length(Q1_pts1);
% By taking the second point, we allow some error... this ensures
% that we don't "stop including" points in Q1 too early.
% (e.g. some digitization mistake, etc.)
bestPointToEndTransient = 2;
% If there's no transient, this IF statement prevents errors
% (e.g. with a tight, tight seal in cell-attached mode)
if (numPointsAfterQ1 == 0)
Qt = 0;
else
if (numPointsAfterQ1 > bestPointToEndTransient)
Q1_lastpoint1 = Q1_pts1(bestPointToEndTransient);
else
Q1_lastpoint1 = Q1_pts1(numPointsAfterQ1);
end
% Calculate the area under the capacitive transient ABOVE I1
% This area is Q1. We are doing numerical integration here
% OLD NUMERICAL INTEGRATION (Rectangle method)
% Q1 = sum(abs(current((201:201+Q1_lastpoint1),ichannel)-I1(ichannel))/Fs);
% The formula above gives the area under the curve with rectangles
% at every timepoint. To convert this to trapezoids at every point,
% we have to subtract the areas of all the little traingles on top
% every rectangle. There are Q1_lastpoint1 such triangles, each of
% same width. Putting these together two at a time, we get Q1_lastpoint1/2
% rectangles, which cover (roughly) half the height of current(201).
% Because the curve stops at current(201+Q1_lastpoint1), we subtract the
% value of current(201+Q1_lastpoint1)/2.
rect_integral = sum(abs(current((201:201+Q1_lastpoint1),ichannel)-I1(ichannel)));
trap_subtraction = abs ( ((current(201,ichannel) - I1(ichannel)) / 2 ) - ...
((current(201+Q1_lastpoint1,ichannel) - I1(ichannel)) / 2 ) );
Q1 = (rect_integral - trap_subtraction) / Fs;
% Find the time constant, tau.
% At t= tau, the current is .368 of its final value.
% OLD, INCORRECT WAY: tauCurrent = .368 * p(end);
% (Some recordings were taken this way.)
%
% NEW: We take the value 63.2% of the way to the SS current
% from the peak of the transient spike as the current at tau.
tauCurrent = .632*(p(end) - I1(ichannel));
% Find indicies on the curve where current is greater than the
% current at tau (i.e. closer to the steady state current).
% Comparison sign is > because transient is a downward (-)
% deflection.
% +2 is because we started searching at 203 (so
% currentsAterTau=1 implies tau=(1+2)/Fs=3/Fs msec)
%
% Arbitrarily, we've chosen to take tau to be at the first
% recorded time point where the current *exceeded* the 63.2%
% current value. We could have just as well used the last time
% point where the current was *below* the 63.2% value.
currentsAfterTau = find(current(203:280,ichannel) > tauCurrent);
if length(currentsAfterTau > 0)
tau = (currentsAfterTau(1)+2)/Fs;
else
tau = 0;
end
% From Axon calculation
% abs measures the dI from the two steady-state currents (I2 & I1)
% Note: Including Q2 in Qt seems to make almost no difference
% in the final value of Cm.
Q2 = (abs(I2 - I1)) .* tau;
% Convert Q since I's are picoamps
Qt=(Q1+Q2(1)) *10^-12;
% Qt=(Q1) *10^-12;
end
% Derivations from Axon/Garcia-Diaz Membrane Test
% rs is the series resistance (pipette + access)
% rm is the membrane resistance
% cm is the membrane capacitance
%
% The Axon Rs seems less accurate than the transient-based rs
% calculation: rs2 = tau * dV / Qt;
rs = dV / ( (mean(current(160:199,ichannel)) - p(end)) * 10^-12);
rm = rt - rs;
cm = (Qt*rt)/(dV*rm);
% Another way to calculate rm. Gives similar results to rm above.
% rm2 = dV / (abs(p(end) - I1(1))*10^-12);
% An alternate way to estimate cm that seems less accurate
% on the model cell, true cap is 47 pF. cm=51pF while cm2=62pF
% cm2 = ((ra+rm)/(ra*rm))*tau;
Rtotal(ichannel) = Rtotal(ichannel) + rt;
Rseries(ichannel) = Rseries(ichannel) + rs;
Rmem(ichannel) = Rmem(ichannel) + rm;
Cmem(ichannel) = Cmem(ichannel) + cm;
% Old way of calculating Rseries (mislabeled as
% Raccess, which is really Rser - Rpip).
% Raccess(ichannel)=Raccess(ichannel) + amplitude / (mean(current(195:199,ichannel))-p(end));
end % for ichannel
avgCurrent = avgCurrent + current;
end % for irepeat
% Now we calculate the final values that we will send back
% Old calculation for Rseries
% Raccess = (10 ^9) * Raccess / nrepeat; % convert from GOhm to Ohms
Rtotal = Rtotal/nrepeat; % Rtotal is already in Ohms
Ileak = Ileak/nrepeat;
Rseries = Rseries/nrepeat;
Rmem = Rmem/nrepeat;
Cmem = Cmem/nrepeat;
current = avgCurrent/nrepeat;
%---------------------------------------------------------------------------