% include figures_common
figures_common

% ****************************************
% FIGURE 2A: compare aCC in Cd, isopotential and ball-and-stick model
% firing activity.
% ****************************************

% PANEL: aCC recording in Cd
marley_Cd_20100204_cell2 = ...
    abf2voltage_clamp([ '../data/10204022.abf' ], ' Marley cell#2 - Cd^{2+}', ...
                      struct('scaleI', 1e3, 'iSteps', 1, 'ichan', 1:14));

% select only current steps 5, 25, and 45 pA consistent with Wei-Hsiang
% et al (2012); assuming bias current is -12 pA (actual value=-17 pA)
plotFigure(plot_abstract(setLevels(marley_Cd_20100204_cell2, [4 8 12]), ...
                         ' - 5, 25, 45 pA', ...
                         struct('fixedSize', [4 4], ...
                                'noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [100 650 NaN NaN])))

% PANEL: isopotential model
model_nap_spiking_fI_match_aCC_cc = ...
    XPP2current_clamp(...
      ['I-range_-12_42-by-13-gKs=50-gKf=24-gNa=100-gNaP=.8-c=4-gL=6.8.dat'], ...
      [10 510], vary_steps(-12, 42, 13)*1e-3, ...
      struct('threshold', 7, 'minInit2MaxAmp', 4, 'paramsStruct', ...
             struct('gL_nS', 6.8, 'gNa', 100, 'gNaP', .8, 'gKs_nS', 50, ...
                    'Cm_pA', 4)));

% with steps 5, 25, 45 pA as in Wei-Hsiang et al (2012):
plotFigure(plot_abstract(setLevels(model_nap_spiking_fI_match_aCC_cc, [4 8 12]), ...
                         ' - 5, 25, 45 pA', ...
                         struct('fixedSize', [4 4], ...
                                'noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 520 NaN NaN])))

% PANEL: ball and stick model
model_bs_spiking2_fI_match_aCC_3traces_cc = ...
    XPP2current_clamp(...
      ['I-range_-6_36-by-2-gaxon=1.2-gKf=1-gaKf=1000-gKs=1-gaKs=50-gNa=85-gNaP=.93-gleak=.3-richard-panel-B.dat'], ...
      [10 510], vary_steps(5, 45, 2)*1e-3, ...
      struct('threshold', 2.5, ...
             'downThreshold', -1, ...
             'minInit2MaxAmp', 2, ...
             'minMin2MaxAmp', 2));

plotFigure(plot_abstract(setLevels(model_bs_spiking2_fI_match_aCC_3traces_cc, [1:3]), ...
                         ' - 3 traces (0, 24, 44pA)', ...
                         struct('fixedSize', [6 4], ...
                                'noTitle', 1, ...
                                'relativeSizes', [1 3], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 520 NaN NaN])))
% Compare cell, iso model and bs model traces in one place for model paper
% traces
plotFigure(plot_stack({...
  plot_abstract(setLevels(marley_Cd_20100204_cell2, [4 8 12]), ...
                         ' - 5, 25, 45 pA', ...
                         struct('noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'ColorOrder', [0 0 0], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [100 650 -60 0])), ...
  plot_abstract(setLevels(model_nap_spiking_fI_match_aCC_cc, [4 8 12]), ...
                         ' - 5, 25, 45 pA', ...
                         struct('noTitle', 1, ...
                                'relativeSizes', [1 4], ...
                                'ColorOrder', [0 0 0], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 520 -60 0])), ...
  plot_abstract(setLevels(model_bs_spiking2_fI_match_aCC_3traces_cc, 1:3), ...
                         ' - 5, 25, 45 pA', ...
                         struct('noTitle', 1, ...
                                'ColorOrder', [0 0 0], ...
                                'relativeSizes', [1 4], ...
                                'curUnit', 'pA', ...
                                'xTicksPos', 'bottom', ...
                                'axisLimits', [0 520 -60 0]))}, [], 'x', ...
                      '', struct('yLabelsPos', 'left', 'yTicksPos', 'left', ...
                                 'fixedSize', [8 4], ...
                                 'ColorOrder', [0 0 0], ...
                                 'relativeSizes', [1 1 1])))

% ****************************************
% FIGURE 2B: compare f-I
% ****************************************

% load feb04_cell2_db 
load('../data/db_2011-02-04_cell2_aCC.mat');

% add bias to cip values (do this only one time after loading!)
feb04_cell2_db(:, 'cip_level_pA') = ...
    feb04_cell2_db(:, 'cip_level_pA') + 12;
% TODO: don't do additional correction below!

% find spikes for cell from database
cell_plot = ...
    set(setProp(plot_cip_stats(feb04_cell2_db(feb04_cell2_db(:, 'cadmium') == 1, :), 'aCC w/ Cd^{2+}', ...
                               'cip_level_pA', 'PulseSpikeRateISI'), ...
                'fixedSize', [2.5 2], 'axisLimits', [2 56 0 150], ...
                'legendLocation', 'SouthEast', 'ColorOrder', gray(1)), ...
        'axis_labels', {'current [pA]', 'firing rate [Hz]'});


% find spikes for iso model and make it into a Pandora database
% (shift by +12 pA like in recording)
model_nap_spiking_fI_match_aCC_db = ...
    params_tests_db(XPP2current_clamp(...
      ['I-range_-12_42-by-13-gKs=50-gKf=24-gNa=100-gNaP=.8-c=4-gL=6.8.dat'], ...
      [10 510], vary_steps(0, 52, 13)*1e-3, ...
      struct('threshold', 7, 'minInit2MaxAmp', 4, 'paramsStruct', ...
             struct('gL_nS', 6.8, 'gNa', 100, 'gNaP', .8, 'gKs_nS', 50, ...
                    'Cm_pA', 4))));


% compare fI curves of cell, iso, and bas models together:
% (shift models left 12 pA to simulate cell's negative offset)
plotFigure(plot_superpose({...
  cell_plot, ...
  plotScatter(model_nap_spiking_fI_match_aCC_db, ...
              'cip_level_pA', 'PulseSpikeRateISI', '', 'iso-model', ...
              struct('LineStyle', '--')), ...
  plot_XPP_fI_Irange(['I-range_-6.5_53.5_by-12-gaxon=1.3-gaKf=200-gaKs=700-gNa=180-gNaP=.01-gleak=.05-zi=2.1-eleak=-55-ealeak=-55-Ca=1.8-Cm=10-Ihold=-6.5.dat'], ...
                     vary_steps(-16.5+12, 43.5+12, 12)*1e-3, struct)}));

% ****************************************
% FIGURE 2C: compare delay-I
% ****************************************

plotFigure(plot_superpose({...
  set(setProp(plot_cip_stats(feb04_cell2_db(feb04_cell2_db(:, 'cadmium') == 1, :), 'aCC w/ Cd^{2+}', ...
                             'cip_level_pA', 'PulseFirstSpikeTime'), ...
              'fixedSize', [2.5 2], 'axisLimits', [-9 44 0 150], ...
              'legendLocation', 'NorthEast', 'ColorOrder', gray(1) ), ...
      'axis_labels', {'current [pA]', '1st spike delay [ms]'}), ...
  plotScatter(model_nap_spiking_fI_match_aCC_db, ...
              'cip_level_pA', 'PulseFirstSpikeTime', '', 'iso-model', ...
              struct('LineStyle', '--')), ...
  plot_XPP_delay_I_Irange(['I-range_-6.5_53.5_by-12-gaxon=1.3-gaKf=200-gaKs=700-gNa=180-gNaP=.01-gleak=.05-zi=2.1-eleak=-55-ealeak=-55-Ca=1.8-Cm=10-Ihold=-6.5.dat'], ...
                     vary_steps(-16.5, 43.5, 12)*1e-3, struct)}));


% ****************************************
% FIGURE 2D: compare f-V
% ****************************************

% for cell from database
cell_fV_plot = ...
    set(setProp(plot_cip_stats(feb04_cell2_db(feb04_cell2_db(:, 'cadmium') == 1, :), 'aCC w/ Cd^{2+}', ...
                               'PulsePotAvg', 'PulseSpikeRateISI'), ...
                'fixedSize', [2.5 2], 'axisLimits', [NaN NaN 0 150], ...
                'legendLocation', 'SouthEast', 'ColorOrder', gray(1)), ...
        'axis_labels', {'membrane voltage [mV]', 'firing rate [Hz]'});

plotFigure(plot_superpose({...
  cell_fV_plot, ...
  plotScatter(model_nap_spiking_fI_match_aCC_db, ...
              'PulsePotAvg', 'PulseSpikeRateISI', '', 'iso-model', ...
              struct('LineStyle', '--')), ...
  plot_XPP_rate_pulsepot_Irange(['I-range_-6.5_53.5_by-12-gaxon=1.3-gaKf=200-gaKs=700-gNa=180-gNaP=.01-gleak=.05-zi=2.1-eleak=-55-ealeak=-55-Ca=1.8-Cm=10-Ihold=-6.5.dat'], ...
                     vary_steps(-16.5, 43.5, 12)*1e-3, struct)}));


