Description of dsCalcCellProperties

0001 function stats = dsCalcCellProperties(data, varargin)
0002 %CALCCELLPROPERTIES - calculates the intrinsic electrophysiological properties of all cells in one or more populations
0003 %
0004 % This is designed to be used in conjunction with the experiment
0005 % dsProbeCellProperties which removes all connections from a model and produces
0006 % a data array of simulated data in response to a series of hyperpolarizing and
0007 % depolarizing pulses. This function is based on the DNSim experiment
0008 % "cell_pulses".
0009 %
0010 % Usage:
0011 %   stats = dsCalcCellProperties(data,'option1',option1,...)
0012 %
0013 % Inputs:
0014 %   - data: array of DynaSim data structures returned by dsProbeCellProperties
0015 %
0016 % Outputs:
0017 %   - stats.(pop).(measure) [1 x num_cells] (averaged over repetitions)
0018 %   - measures (intrinsic properties):
0019 %       RMP, V_thresh, R_in, tau_m, FI_slope, CV, AR24, AR_coefficient
0020 %       FR_min (threshrate), FR_min2 (steprate?), FR_max (steprate?)
0021 %     - AP morphology: AP_amp, AP_dur (spikewidth), AP_taur, AP_taud
0022 %                      Ih_relsag, Ih_abssag, hump, AHP_amp, AHP_dur
0023 %                      AHP_time2trough, ISI_median, AR23, AR13, ISI1,
0024 %                      min_ISI_median, ISI_step_median
0025 %
0026 % Algorithm walkthrough:
0027 %   From Iinj=0:
0028 %   RMP = (avg over 50-100% step | Iinj=0)
0029 %
0030 %   From largest hyperpolarizing step:
0031 %   Ih Sag = (Vend-Vmin)/|RMP-Vend|
0032 %       where Vmin=(min voltage during T sec hyperpolarizing step)
0033 %             Vend=(V at end of hyperpolarizing step)
0034 %   Ih abs sag = (Vend-Vmin)
0035 %
0036 %   From last subthreshold step:
0037 %   Hump = (Vmax-Vend)/|RMP-Vend|
0038 %       where Vmax=(max voltage during depolarizing step preceding spike)
0039 %             Vend=(V at end of step)
0040 %
0041 %   Across hyperpolarizing subthreshold steps:
0042 %   Rin = Input resistence (I/V slope) [4]
0043 %   taum = Membrane time constant: time for voltage relaxation to (1/e)Vmax
0044 %     where Vmax=max deflection from baseline on hyperpolarizing steps
0045 %     note: avg taum's over small drives that keep active currents silent (eg, 5-15pA)
0046 %
0047 %   Across suprathreshold steps with at least two spikes:
0048 %   FI slope [Hz/nA]: slope of f/I curve (firing freq vs injected current 0-140pA) (see [2])
0049 %
0050 %   Across suprathreshold steps >=60pA above first step with at least two spikes:
0051 %   AR coefficient: (slope of AR/I) where per step AR=ISI(1)/ISI(end)
0052 %   Note: AR = Adaptation ratio
0053 %   ISI_step_median = median ISI on step 60pA above first step w/ 2 spikes
0054 %   FR_step = mean FR on step 60pA above first step w/ 2 spikes
0055 %   ARif = max AR across steps >=60pA above first step w/ 2 spikes
0056 %
0057 %   From first suprathreshold T sec step:
0058 %   FRmin (threshrate) = (# spikes)/T
0059 %
0060 %   From first suprathreshold T sec step with at least two spikes:
0061 %   FRmin2 (steprate?) = (# spikes)/T
0062 %
0063 %   From first spike of first suprathreshold step: AP morphology
0064 %   Vthresh = V( crossing(dV/dt,20mV/ms) ) %10mV/ms) )
0065 %   AP_amp = (Vpeak-Vthresh)
0066 %   AP_taur = (time to rise from 10% to 90% between Vthresh and Vpeak)
0067 %   AP_taud = (time to decay from 10% to 90% between Vpeak and Vthresh)
0068 %   AP_dur (spikewidth) = (time between rising and falling (Vthresh+APamp/2) = (Vthresh+Vpeak)/2)
0069 %   AHP_amp = (Vbaseline-Vmin) where Vmin taken during repolarizing phase
0070 %   AHP_dur (AHP duration, half-width) = time between half-peak amplitude of AHP
0071 %                                     = (time between falling and rising (Vbaseline+AHPamp/2))
0072 %   AHP_time2trough = (time between falling Vbaseline and AHPamp)
0073 %   ADPamp? ADPdur?
0074 %
0075 %   From suprathreshold step 20pA above first step with at least two spikes:
0076 %   CV(ISIs over 30-100% step)
0077 %
0078 %   From last suprathreshold T sec step (or step at amp=max (eg, 140pA)):
0079 %   FR_max (steprate?): (# spikes)/T
0080 %   min_ISI_median
0081 %   max_ISI
0082 %   AR24 = ISI(2)/ISI(4)
0083 %
0084 % References for methods used:
0085 %   [1] Steffensen, Scott C., et al. "Electrophysiological characterization of
0086 %     GABAergic neurons in the ventral tegmental area." The Journal of
0087 %     neuroscience 18.19 (1998): 8003-8015.
0088 %   [2] Van Aerde, Karlijn I., et al. "Flexible spike timing of layer 5 neurons
0089 %     during dynamic beta oscillation shifts in rat prefrontal cortex." The
0090 %     Journal of physiology 587.21 (2009): 5177-5196.
0091 %   [3] Connors, BW, MJ Gutnick, DA Prince. "Electrophysiological properties of
0092 %     neocortical neurons in vitro." Journal of Neurophysiology 48.6 (1982):
0093 %     1302-1320.
0094 %   [4] Povysheva, Nadezhda V., et al. "Parvalbumin-positive basket
0095 %     interneurons in monkey and rat prefrontal cortex." Journal of
0096 %     neurophysiology 100.4 (2008): 2348-2360.
0097 %   [5] Gonz?lez-Burgos, Guillermo, et al. "Functional properties of fast
0098 %     spiking interneurons and their synaptic connections with pyramidal cells in
0099 %     primate dorsolateral prefrontal cortex." Journal of Neurophysiology 93.2
0100 %     (2005): 942-953.
0101 %   [6] Gorelova, Natalia, Jeremy K. Seamans, and Charles R. Yang. "Mechanisms
0102 %     of dopamine activation of fast-spiking interneurons that exert inhibition
0103 %     in rat prefrontal cortex." Journal of neurophysiology 88.6 (2002):
0104 %     3150-3166.
0105 %
0106 % Example:
0107 %   data = dsProbeCellProperties(model)
0108 %   stats = dsCalcCellProperties(data)
0109 %
0110 % See also: dsProbeCellProperties
0111 %
0112 % Author: Jason Sherfey, PhD <jssherfey@gmail.com>
0113 % Copyright (C) 2016 Jason Sherfey, Boston University, USA
0114 
0115 % Check inputs
0116 options=dsCheckOptions(varargin,{...
0117   'spike_threshold',1e-5,[],...
0118   'skip_time',10,[],... % time [ms] to exclude from detection
0119   'plot_flag',0,{0,1},...
0120   'equivalent_cells_flag',0,[],... % if true, only process one cell per pop
0121   'auto_gen_test_data_flag',0,{0,1},...
0122   },false);
0123 
0124 %% auto_gen_test_data_flag argin
0125 if options.auto_gen_test_data_flag
0126   varargs = varargin;
0127   varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
0128   varargs(end+1:end+2) = {'unit_test_flag',1};
0129   argin = [{data}, varargs]; % specific to this function
0130 end
0131 
0132 data=dsCheckData(data, varargin{:});
0133 model=data(1).model;
0134 time=data(1).time;
0135 num_times=length(time);
0136 
0137 % extract experiment parameters
0138 experiment_options=data(1).simulator_options.experiment_options;
0139 
0140 % stimulus timing
0141 onset=experiment_options.onset; %model.parameters.([pop_names{1} '_onset']);
0142 offset=experiment_options.offset; %model.parameters.([pop_names{1} '_offset']);
0143 tsel=find(time>=onset&time<=offset);
0144 
0145 % assume exactly one simulation per amplitude
0146 % if num_steps ~= length(data)
0147 %   error('there can only be one simulation per stimulus amplitude');
0148 % end
0149 
0150 % extract population info
0151 pop_names={model.specification.populations.name};
0152 num_pops=length(pop_names);
0153 
0154 % set num_cells=1 if equivalent_cells_flag
0155 if experiment_options.equivalent_cells_flag==1
0156   pop_sizes=ones(1,num_pops);
0157 else
0158   pop_sizes=[model.specification.populations.size];
0159 end
0160 
0161 % sort data by [data.(pop1)_TONIC] and sort amplitudes
0162 [amplitudes,I]=sort([data.([pop_names{1} '_TONIC'])]);
0163 data=data(I);
0164 num_steps=length(amplitudes);
0165 
0166 % get list of variables to analyze
0167 [vars,vars_pop]=dsSelectVariables(data(1).labels, varargin{:});
0168 
0169 % assume only one variable per population
0170 if length(unique(vars_pop))>num_pops
0171   [vars;vars_pop]
0172   error('there can only be one variable per population.');
0173 end
0174 
0175 % collect pulses
0176 % input=zeros(num_times,num_steps);
0177 % for s=1:num_steps
0178 %   input(:,s)=data(s).([pop_names{1} '_pulse'])(:,1);
0179 % end
0180 
0181 CF = (1e-6)/(1e-8);   % pA/um^2 => uA/cm^2. note: 1um2=1e-8cm2, 1pA=1e-6uA
0182 amplitudes_pA=amplitudes*experiment_options.membrane_area/CF;
0183 
0184 % initialize stats structure
0185 stats.experiment_options=experiment_options;
0186 stats.amplitudes=amplitudes;
0187 stats.time=time;
0188 stats.cell_results={}; % list of field names storing results for each cell
0189 stats.pop_results={};  % list of field names storing results for each population
0190 
0191 % analyze each population
0192 for p=1:num_pops
0193   % extract info for this population
0194   pop=pop_names{p};
0195   this_var=vars{strcmp(vars_pop,pop)};
0196   num_cells=pop_sizes(p);
0197   % preallocate intrinsic property measures (each will be an average over repetitions)
0198   stats.(pop).RMP     =nan(1,num_cells);
0199   stats.(pop).V_thresh=nan(1,num_cells);
0200   stats.(pop).R_in    =nan(1,num_cells);
0201   stats.(pop).tau_m   =nan(1,num_cells);
0202   stats.(pop).FI_slope=nan(1,num_cells);
0203   stats.(pop).CV      =nan(1,num_cells);
0204   stats.(pop).ISI_median=nan(1,num_cells);
0205   stats.(pop).ISI_step_median=nan(1,num_cells);
0206   stats.(pop).min_ISI_median=nan(1,num_cells);
0207   stats.(pop).max_ISI=nan(1,num_cells);
0208   stats.(pop).AR13    =nan(1,num_cells);
0209   stats.(pop).AR23    =nan(1,num_cells);
0210   stats.(pop).AR24    =nan(1,num_cells);
0211   stats.(pop).ARif    =nan(1,num_cells);
0212   stats.(pop).AR_coeff=nan(1,num_cells);
0213   stats.(pop).ISI1    =nan(1,num_cells);
0214   stats.(pop).FR_min  =nan(1,num_cells);
0215   stats.(pop).FR_min2 =nan(1,num_cells);
0216   stats.(pop).FR_step =nan(1,num_cells);
0217   stats.(pop).FR_max  =nan(1,num_cells);
0218   stats.(pop).AP_amp  =nan(1,num_cells);
0219   stats.(pop).AP_dur  =nan(1,num_cells);
0220   stats.(pop).AP_taur =nan(1,num_cells);
0221   stats.(pop).AP_taud =nan(1,num_cells);
0222   stats.(pop).Ih_relsag  =nan(1,num_cells);
0223   stats.(pop).Ih_abssag=nan(1,num_cells);
0224   stats.(pop).hump    =nan(1,num_cells);
0225   stats.(pop).AHP_amp =nan(1,num_cells);
0226   stats.(pop).AHP_dur =nan(1,num_cells);
0227   stats.(pop).AHP_time2trough=nan(1,num_cells);
0228   % determine whether cells are spiking or not on each step
0229   spike_times=cell(num_steps,num_cells);
0230   for c=1:num_cells
0231     for s=1:num_steps
0232       % select trace for this (pop,cell,step)
0233       X=data(s).(this_var)(tsel,c);
0234       % get spikes in this cell during step
0235       spike_inds=1+find((X(2:end)>=options.spike_threshold & X(1:end-1)<options.spike_threshold));
0236       spike_times{s,c}=time(tsel(spike_inds));
0237     end
0238     num_spikes=cellfun(@numel,spike_times(:,c));
0239     is_spiking=(num_spikes>0); % {0:subthreshold,1:suprathreshold}
0240 
0241     % Calculate measures of cell intrinsic properties
0242     % 1) calculate RMP (avg over repetitions) From steps with amp=0
0243     % RMP = (avg over 50-100% step | Iinj=0)
0244     step_sel=find(amplitudes==0); % select all simulations with amp=0
0245     time_sel=tsel(round(.5*numel(tsel)):end); % select 50-100% of step
0246     rmp=0;
0247     for s=1:length(step_sel)
0248       rmp=rmp+nanmean(data(step_sel(s)).(this_var)(time_sel,c));
0249     end
0250     stats.(pop).RMP(c)=rmp/length(step_sel);
0251 
0252     % 2) calculate Ih_relsag (avg over repetitions) From largest hyperpolarizing step
0253     % Ih Sag = (Vend-Vmin)/|RMP-Vend|
0254     %     where Vmin=(min voltage during T sec hyperpolarizing step)
0255     %           Vend=(V at end of hyperpolarizing step)
0256     step_sel=find(amplitudes==min(amplitudes) & amplitudes<0);
0257     sag=0; abssag=0;
0258     for s=1:length(step_sel)
0259       V=data(step_sel(s)).(this_var)(:,c);
0260       bl=V(tsel(1)-1);   % baseline value is point immediately before stim onset
0261       Vend=V(tsel(end)); % final point
0262       Vmin=min(V(tsel)); % minimum point of sag
0263       sag=sag+((Vend-Vmin)/abs(bl-Vend));
0264       abssag=abssag+(Vend-Vmin);
0265       %figure; plot(time,V); line(xlim,[Vmin Vmin]); line(xlim,[Vend Vend]); line(xlim,[bl bl]);
0266     end
0267     stats.(pop).Ih_relsag(c)=sag/length(step_sel);
0268     stats.(pop).Ih_abssag(c)=abssag/length(step_sel);
0269 
0270     % 3) calculate hump (avg over repetitions) From last subthreshold step
0271     % Hump = (Vmax-Vend)/|RMP-Vend|
0272     %     where Vmax=(max voltage during depolarizing step preceding spike)
0273     %           Vend=(V at end of step)
0274     step_sel=max(1,find(num_spikes>0,1,'first')-1);
0275     if ~isempty(step_sel)
0276       V=data(step_sel).(this_var)(:,c);
0277       bl=V(tsel(1)-1);   % baseline value is point immediately before stim onset
0278       Vend=V(tsel(end)); % final point
0279       Vmax=max(V(tsel)); % maximum point of hump
0280       stats.(pop).hump(c)=((Vmax-Vend)/abs(bl-Vend));
0281     end
0282 
0283     % 4) calculate (R_in,tau_m) (from avg over repetitions) Across hyperpolarizing subthreshold steps
0284     % Rin = Input resistence (I/V slope) [4]
0285     % taum = Membrane time constant: time for voltage relaxation to (1/e)Vmax
0286     %   where Vmax=max deflection from baseline on hyperpolarizing steps
0287     %   note: avg taum's over small drives that keep active currents silent (eg, 5-15pA)
0288     step_sel=find(num_spikes==0 & amplitudes'<0);
0289     amps=amplitudes(step_sel);
0290     uamps=unique(amps);
0291     nuamps=length(uamps);
0292     Veq=nan(1,nuamps);
0293     taum=nan(1,nuamps);
0294     time_sel=tsel(round(.5*numel(tsel)):end); % select 50-100% of step
0295     % loop over unique amplitudes
0296     for a=1:nuamps
0297       asel=find(amps==uamps(a)); % indices to subthreshold amps with this value
0298       % loop over repetitions of the same amplitude
0299       veq=0; tau=0;
0300       for l=1:length(asel)
0301         step=step_sel(asel(l)); % index to this simulation
0302         X=data(step).(this_var)(:,c);
0303         % store mean voltage over 50-100% step
0304         veq=veq+nanmean(X(time_sel));
0305         % calc membrane time constant based on this step (time to Vmax(1/e))
0306         bl=X(tsel(1)-1); % baseline value is point immediately before stim onset
0307         V=X-bl; % shift baseline to V=0
0308         if uamps(a)<0
0309           V=abs(V); % invert values so that deflection is positive
0310         end
0311         Vmax=V(tsel(end)); % max is final point during stim period before offset
0312         Vthresh=Vmax*(1/exp(1)); % voltage at 1/e
0313         Vdecay=V(tsel(end):end); % voltage decay after stim offset
0314         tthresh=time(tsel(end)+find(Vdecay<Vthresh,1,'first')-1); % time at which decay crosses 1/e
0315         tau=tau+(tthresh-offset); % time to drop to 1/e after stim offset
0316       end
0317       Veq(a)=veq/length(asel);
0318       if any(tau)
0319         taum(a)=tau/length(asel);
0320       end
0321     end
0322     % calc R_in (from slope of amp/Veq)
0323     if any(~isnan(Veq))
0324       sel=~isnan(Veq);
0325       P=polyfit(uamps(sel),Veq(sel),1);
0326       stats.(pop).R_in(c)=P(1);
0327     end
0328     if any(~isnan(taum))
0329       stats.(pop).tau_m(c)=taum(find(~isnan(taum),1,'first'));%1);%nanmedian(taum);
0330     end
0331 
0332     % 5) calculate FI_slope (from avg over repetitions) Across suprathreshold steps with at least two spikes
0333     % FI slope: slope of f/I curve (firing freq vs injected current 0-140pA) (see [2])
0334 %     step_sel=find(num_spikes>1 & amplitudes'>0);
0335 %     amps=amplitudes(step_sel);
0336 %     uamps=unique(amps);
0337 %     nuamps=length(uamps);
0338 %     FR=nan(1,nuamps);
0339 %     % loop over unique amplitudes
0340 %     for a=1:nuamps
0341 %       asel=find(amps==uamps(a)); % indices to subthreshold amps with this value
0342 %       % loop over repetitions of the same amplitude
0343 %       fr=0;
0344 %       for l=1:length(asel)
0345 %         step=step_sel(asel(l)); % index to this simulation
0346 %         fr=fr+(num_spikes(step)/((offset-onset)/1000));
0347 %       end
0348 %       FR(a)=fr/length(asel); % [Hz]
0349 %     end
0350 %     % calc FI_slope
0351 %     if any(~isnan(FR))
0352 %       sel=~isnan(FR);
0353 %       amps_nA=uamps(sel)*experiment_options.membrane_area/CF/1000;
0354 %       P=polyfit(amps_nA,FR(sel),1);
0355 %       stats.(pop).FI_slope(c)=P(1); % [Hz/nA]
0356 %       figure; plot(amps_nA,FR(sel)); xlabel('amps [nA]'); ylabel('FR [Hz]');
0357 %       [num_spikes(step_sel) amplitudes(step_sel)']
0358 %     end
0359 
0360     % 6) calculate ARs (avg over repetitions) & AR_coefficient Across
0361     %    suprathreshold steps >=60pA above first step with at least two spikes
0362     % AR coefficient: (slope of AR/I) where per step AR=ISI(1)/ISI(end)
0363     thresh_step=find(num_spikes>1,1,'first');
0364     if any(thresh_step)
0365       step_sel=find(amplitudes_pA>(amplitudes_pA(thresh_step)+60));
0366     else
0367       step_sel=[];
0368     end
0369     if any(step_sel)
0370       amps=amplitudes(step_sel);
0371       uamps=unique(amps);
0372       nuamps=length(uamps);
0373       AR=nan(1,nuamps);
0374       FR=nan(1,nuamps);
0375       % loop over unique amplitudes
0376       for a=1:nuamps
0377         asel=find(amps==uamps(a)); % indices to subthreshold amps with this value
0378         % loop over repetitions of the same amplitude
0379         ar=0; denom=0;
0380         fr=0;
0381         for l=1:length(asel)
0382           step=step_sel(asel(l)); % index to this simulation
0383           fr=fr+(num_spikes(step)/((offset-onset)/1000));
0384           ISI=diff(spike_times{step,c});
0385           if length(ISI)>1
0386             denom=denom+1;
0387   %           if length(ISI)>2
0388   %             ar=ar+(ISI(2)/ISI(end));
0389   %           else
0390               ar=ar+(ISI(1)/ISI(end));
0391   %           end
0392           end
0393         end
0394         AR(a)=ar/denom;
0395         FR(a)=fr/length(asel); % [Hz]
0396       end
0397       % calc ISI_median at 60pA above threshold
0398       spikes=spike_times{step_sel(1),c};
0399       ISI=diff(spikes);
0400       stats.(pop).ISI_step_median(c)=median(ISI);
0401       % store FR at 60pA above threshold
0402       stats.(pop).FR_step(c)=FR(1); % FR at 60pA above threshold
0403       if ~strcmp(reportUI,'matlab') && exist('nanmax') ~= 2 % 'nanmax is not in Octave's path
0404         try
0405           pkg load statistics; % trying to load octave forge 'statistics' package before using nanmax function
0406         catch
0407           error('nanmax function is needed, please install the statistics package from Octave Forge');
0408         end
0409       end
0410       stats.(pop).ARif(c)=nanmax(AR); % nanmean, nanmedian
0411       % calculate AR_coefficient
0412       if any(~isnan(AR))
0413         sel=~isnan(AR);
0414         P=polyfit(uamps(sel),AR(sel),1);
0415         stats.(pop).AR_coeff(c)=P(1);
0416         %figure; plot(uamps(sel),AR(sel)); xlabel('amps'); ylabel('AR');
0417       end
0418       % calc FI_slope
0419       if any(~isnan(FR))
0420         sel=~isnan(FR);
0421         amps_nA=uamps(sel)*experiment_options.membrane_area/CF/1000;
0422         P=polyfit(amps_nA,FR(sel),1);
0423         stats.(pop).FI_slope(c)=P(1); % [Hz/nA]
0424         %figure; plot(amps_nA,FR(sel)); xlabel('amps [nA]'); ylabel('FR [Hz]');
0425         %[num_spikes(step_sel) amplitudes(step_sel)']
0426       end
0427     end
0428 
0429     % 7) calculate FR_min From first suprathreshold T sec step
0430     % FRmin (threshrate) = (# spikes)/T
0431     step_sel=find(num_spikes>0,1,'first');
0432     if any(step_sel)
0433       stats.(pop).FR_min(c)=num_spikes(step_sel)/((offset-onset)/1000);
0434     end
0435 
0436     % 8) calculate FR_min2 From first suprathreshold T sec step with at least two spikes
0437     % FRmin2 (steprate?) = (# spikes)/T
0438     step_sel=find(num_spikes>1,1,'first');
0439     if any(step_sel)
0440       stats.(pop).FR_min2(c)=num_spikes(step_sel)/((offset-onset)/1000);
0441       stats.(pop).ISI1(c)=(spike_times{step_sel,c}(2)-spike_times{step_sel,c}(1));
0442     end
0443 
0444     % 9) calculate CV (avg over repetitions) From suprathreshold step 20pA above first step with at least two spikes
0445     % CV(ISIs over 30-100% step)
0446     tmin=time(tsel(round(.3*numel(tsel)))); % time at 30% through step
0447     thresh_step=find(num_spikes>1,1,'first');
0448     if any(thresh_step)
0449       step_sel=find(amplitudes_pA>(amplitudes_pA(thresh_step)+20));
0450     else
0451       step_sel=[];
0452     end
0453     if any(step_sel)
0454       spikes=spike_times{step_sel,c};
0455       spikes=spikes(spikes>tmin);
0456       ISI=diff(spikes);
0457       if any(ISI)
0458         stats.(pop).CV(c)=std(ISI)/mean(ISI);
0459       end
0460       spikes=spike_times{step_sel,c};
0461       ISI=diff(spikes);
0462       stats.(pop).ISI_median(c)=median(ISI);
0463       stats.(pop).max_ISI(c)=max(ISI);
0464     end
0465 
0466     % 10) calculate (FRmax,AR24) From last suprathreshold T sec step
0467     % FRmax (steprate?): (# spikes)/T
0468     % AR24 = ISI(2)/ISI(4)
0469     step_sel=find(num_spikes>4,1,'last');
0470     if any(step_sel)
0471       spikes=spike_times{step_sel,c};
0472       ISIs=diff(spikes)/1000;
0473       finst=1./(ISIs);
0474       stats.(pop).AR24(c)=finst(2)/finst(4);
0475 %       stats.(pop).AR23(c)=ISIs(2)/ISIs(3);
0476 %       stats.(pop).AR13(c)=ISIs(1)/ISIs(3);
0477       stats.(pop).FR_max(c)=num_spikes(step_sel)/((offset-onset)/1000);
0478       stats.(pop).min_ISI_median(c)=median(diff(spikes));
0479     end
0480     step_sel=find(num_spikes>3,1,'first');
0481     if any(step_sel)
0482       spikes=spike_times{step_sel,c};
0483       ISIs=diff(spikes)/1000;
0484       stats.(pop).AR23(c)=ISIs(2)/ISIs(3);
0485       stats.(pop).AR13(c)=ISIs(1)/ISIs(3);
0486     end
0487 
0488     % 11) calculate AP morphology From first spike of first suprathreshold step with at least two spikes
0489     step_sel=find(num_spikes>0,1,'first');
0490     if any(step_sel)
0491       X=double(data(step_sel).(this_var)(:,c));
0492       spks=spike_times{step_sel,c};
0493       spk_1=spks(1); % time of first spike
0494       if length(spks)>1
0495         spk_2=spks(2); % time of second spike
0496       else
0497         spk_2=inf;
0498       end
0499       pad=100; % ms, time before and after spike detection (make longer than expected AHP)
0500       spk_beg=max(spk_1-pad,onset);
0501       spk_end=min(spk_1+pad,spk_2); % make sure spike interval excludes the following spike
0502       spk_beg_i=nearest(time,spk_beg);
0503       spk_end_i=nearest(time,spk_end);
0504       spk_inds=spk_beg_i:spk_end_i; % indices for this spike
0505       t=time(spk_inds); % times around spike
0506       V=X(spk_inds);  % trace around spike
0507       % Vthresh = V( crossing(dV/dt,10mV/ms) )
0508       dVdt=diff(V)/(t(2)-t(1));
0509       dVdt=smooth(dVdt,round(1/(t(2)-t(1)))); % 1ms smoothing
0510       thresh_ind=find(dVdt>20,1,'first');%crossing(dVdt,t,10);
0511       Vthresh=V(thresh_ind(1));
0512       % APamp = (Vpeak-Vthresh)
0513       [pk,loc]=findpeaks(V,'NPeaks',1);
0514       Vpeak=V(loc);%max(V);
0515       Vpeak_i=loc;%find(V==Vpeak,1,'first');
0516       APamp=Vpeak-Vthresh;
0517       V10=Vthresh+.1*APamp; % 10% from Vthresh to Vpeak
0518       V90=Vthresh+.9*APamp; % 90% from Vthresh to Vpeak
0519       % depolarizing rising phase
0520       V10_i_r=1+find(V(2:end)>=V10 & V(1:end-1)<V10); % first
0521       V90_i_r=1+find(V(2:end)>=V90 & V(1:end-1)<V90); % second
0522       % APtaur = (time to rise from 10% to 90% between Vthresh and Vpeak)
0523       if any(V10_i_r) && any(V90_i_r)
0524         APtaur=t(V90_i_r(1))-t(V10_i_r(1));
0525       else
0526         APtaur=nan;
0527       end
0528       % repolarizing falling phase
0529       V10_i_d=1+find(V(1:end-1)>=V10 & V(2:end)<V10); % second
0530       V90_i_d=1+find(V(1:end-1)>=V90 & V(2:end)<V90); % first
0531       if numel(V10_i_d)>1 && numel(V10_i_d)>1
0532         V10_i_d=V10_i_d(end);
0533         V90_i_d=V90_i_d(end);
0534       end
0535       if isempty(V10_i_d)
0536         % set to the minimum point of the repolarizing phase
0537         V10=min(V(Vpeak_i:end));
0538         V10_i_d=find(V(Vpeak_i:end)==V10)+Vpeak_i-1;
0539       end
0540       % APtaud = (time to decay from 10% to 90% between Vpeak and Vthresh)
0541       if any(V10_i_d) && any(V90_i_d)
0542         APtaud=t(V10_i_d(1))-t(V90_i_d(1));
0543       else
0544         APtaud=nan;
0545       end
0546       % APdur (spikewidth) = (time between rising and falling (Vthresh+APamp/2) = (Vthresh+Vpeak)/2)
0547       V50=Vthresh+.5*APamp;
0548       V50_i_r=1+find(V(2:end)>=V50 & V(1:end-1)<V50);
0549       V50_i_d=1+find(V(1:end-1)>=V50 & V(2:end)<V50);
0550       if any(V50_i_d) && any(V50_i_r)
0551         APdur=t(V50_i_d(1))-t(V50_i_r(1));
0552       else
0553         APdur=nan;
0554       end
0555       % AHPamp = (Vbaseline-Vmin) where Vmin taken during repolarizing phase
0556       %bl=X(spk_inds(1)-1); % baseline voltage
0557       bl=mean(V);
0558       % find trough: look for first peak in the inverted post-spike trace
0559       [pk,loc]=findpeaks(-smooth(V(V10_i_d:end),100),'NPeaks',1);
0560       ahp_trough_i=V10_i_d+loc-1;
0561       Vmin=V(ahp_trough_i);
0562       %Vmin=-findpeaks(-V(V10_i_d:end),'NPeaks',1);
0563       %Vmin=max(findpeaks(-V(V10_i_d:end)));
0564       %Vmin=min(V(V10_i_d:end));
0565       AHPamp=abs(bl-Vmin);
0566       % AHPdur (AHP duration, half-width) = time between half-peak amplitude of AHP
0567       %                                   = (time between falling and rising (Vbaseline+AHPamp/2))
0568       ahp_V50=bl+.5*AHPamp;
0569       ahp_V50_i_d=1+find(V(1:end-1)>=ahp_V50 & V(2:end)<ahp_V50); % first
0570       ahp_V50_i_d(ahp_V50_i_d<V10_i_d)=[]; % remove crossings before repolarization
0571       ahp_V50_i_r=1+find(V(2:end)>=ahp_V50 & V(1:end-1)<ahp_V50); % second
0572       if any(ahp_V50_i_d)
0573         ahp_V50_i_r(ahp_V50_i_r<ahp_V50_i_d(1))=[]; % remove crossings before start of AHP
0574       end
0575       if any(ahp_V50_i_r) && any(ahp_V50_i_d)
0576         AHPdur=t(ahp_V50_i_r(1))-t(ahp_V50_i_d(1));
0577       else
0578         AHPdur=nan;
0579       end
0580       % AHP_time2trough = (time between falling Vbaseline and Vmin)
0581 %       ahp_bl_i_d=1+find(V(1:end-1)>=ahp_V50 & V(2:end)<ahp_V50); % first
0582 %       ahp_bl_i_d(ahp_bl_i_d<V10_i_d)=[]; % remove crossings before repolarization
0583       %ahp_trough_i=find(V==Vmin);
0584       %ahp_trough_i(ahp_trough_i<V10_i_d)=[]; % remove crossings before start
0585       if any(V10_i_d)
0586         AHPtime2trough=t(ahp_trough_i(1))-t(Vpeak_i);%t(V10_i_d(1));
0587         if options.plot_flag
0588           figure; plot(t,V);
0589           line(xlim,[V10 V10]); line(xlim,[Vmin Vmin]);
0590           line([t(ahp_trough_i(1)) t(ahp_trough_i(1))],ylim);
0591           line([t(Vpeak_i) t(Vpeak_i)],ylim,'color','r');
0592         end
0593       else
0594         AHPtime2trough=nan;
0595       end
0596       if options.plot_flag
0597         figure; plot(t,V);
0598         line(xlim,[bl bl]);
0599         line(xlim,[Vthresh Vthresh]);
0600         line(xlim,[Vpeak Vpeak]);
0601         line(xlim,[V10 V10]);
0602         line([t(V10_i_r(1)) t(V10_i_r(1))],ylim);
0603         line([t(V10_i_d(1)) t(V10_i_d(1))],ylim);
0604         line(xlim,[V50 V50]);
0605         line(xlim,[V90 V90]);
0606         line(xlim,[ahp_V50 ahp_V50],'color','r');
0607         line([t(ahp_V50_i_r(1)) t(ahp_V50_i_r(1))],ylim,'color','r');
0608         line([t(ahp_V50_i_d(1)) t(ahp_V50_i_d(1))],ylim,'color','r');
0609       end
0610       stats.(pop).V_thresh(c)=Vthresh;
0611       stats.(pop).AP_amp(c)=APamp;
0612       stats.(pop).AP_taur(c)=APtaur;
0613       stats.(pop).AP_taud(c)=APtaud;
0614       stats.(pop).AP_dur(c)=APdur;
0615       stats.(pop).AHP_amp(c)=AHPamp;
0616       stats.(pop).AHP_dur(c)=AHPdur;
0617       stats.(pop).AHP_time2trough(c)=AHPtime2trough;
0618     end
0619   end
0620   % collect data for this population across simulations
0621   % note: each simulation has a different input amplitude
0622   X=zeros(num_steps,num_times,num_cells);
0623   for s=1:num_steps
0624     X(s,:,:)=data(s).(this_var);
0625   end
0626   % calculate simple means
0627   if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path
0628     try
0629       pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function
0630     catch
0631       error('nanmean function is needed, please install the statistics package from Octave Forge');
0632     end
0633   end
0634   means=nanmean(X(:,tsel,:),2); % average over select times
0635   % store means in stats structure
0636   stats.([this_var '_mean_per_amp'])=squeeze(means);
0637   stats.([this_var '_pop_mean_per_amp'])=nanmean(means,3);
0638   stats.([this_var '_pop_mean_per_time'])=nanmean(X,3);
0639 end
0640 
0641 if options.plot_flag
0642   v=1;
0643   figure('position',[180 330 1450 530])
0644   subplot(2,2,1); % V(t)
0645   plot(time,stats.([vars{v} '_pop_mean_per_time']));
0646   xlabel('time [ms]'); ylabel(['<' strrep(vars{v},'_','\_') ', pop>']);
0647   legend(cellfun(@(x)['I=' num2str(x)],num2cell(amplitudes),'uni',0));
0648   subplot(2,2,3); % I(t)
0649   %plot(time,input); xlabel('time [ms]'); ylabel('I(t)');
0650   %legend(cellfun(@(x)['I=' num2str(x)],num2cell(amplitudes),'uni',0));
0651   subplot(2,2,2); % I/V
0652   plot(amplitudes,stats.([vars{v} '_mean_per_amp'])); hold on
0653   plot(amplitudes,stats.([vars{v} '_pop_mean_per_amp']),'k-','linewidth',5);
0654   xlabel('amplitude'); ylabel(['<' strrep(vars{v},'_','\_') ', time>']);
0655 end
0656 
0657 %% auto_gen_test_data_flag argout
0658 if options.auto_gen_test_data_flag
0659   argout = {stats}; % specific to this function
0660 
0661   dsUnitSaveAutoGenTestData(argin, argout);
0662 end
0663 
0664 end
0665 
0666 
0667 % check how Tallie calculated:
0668 % - threshrate (mean spike rate on tonic depol just above threshold)
0669 % - steprate (mean spike rate on final? step depol)
0670 % - AHP 5-20ms (Q: mean post-spike AHP amplitude over 5-20ms after trough?)
0671 % - SpikeAmp (same as [4])
0672 % - SpikeWidth (Spike width at half height) (same as [4])
0673 % - RMP
0674 
0675 % AHPtau (AHP duration, half-width) = time between half-peak amplitude of AHP
0676 
0677 % From Iinj=0:
0678 % RMP: resting membrane potential [mV] = avg V over step 150-500ms
0679 
0680 % From first spike:
0681 % Vthresh (spike threshold) [4]: level of voltage deflection exceeding 20mV/1ms %10mV/1ms
0682 % Peak amplitude [4]: (peak - threshold value)
0683 % Spike rise time [4]: time to rise from 10% to 90% of peak amplitude
0684 % Spike decay time [4]: time to decay from 10% to 90% of the amplitude b/w peak and spike threshold
0685 % Spike duration (half-width) [1,4]: "time between half-peak amplitude for the falling and rising edges"
0686 
0687 % Other measures from [4]:
0688 % Sag (Ih?) = (Vmin-Vend)/|RMP-Vend|
0689 %     where Vmin=(min voltage during 500ms hyperpolarizing step)
0690 %           Vend=(V at end of hyperpolarizing step)
0691 % Hump = (Vmax-Vend)/|RMP-Vend|
0692 %     where Vmax=(max voltage during depolarizing step preceding spike)
0693 %           Vend=(V at end of step)
0694 % Adaptation ratio AR = ISI1/ISIend = f(step amplitude)
0695 % AR coefficient = slope of (step amp, AR) for amp > 60pA above spike threshold
0696 % Coefficient of variation CV of ISIs was measured from spikes during
0697 %   150-500ms at depolarizing pulse 20pA above spike threshold
0698 
0699 % Measures from [2]:
0700 % FR: calculate from total spikes elicited by 1sec current injection at 140pA (see [2])
0701 % AR24: Adaptation ratio: ratio of inst. freq of 2nd and 4th spike intervals in train (see [2])
dsCalcCellProperties

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
