%  This script is used to generate simulation data using the HH
%propagated action potential software.  The extracellular potassium 
%concentration is varied and the approximate refractory period is determined
% by varying stimulus interpulse times. All other parameters remain unchanged.

%  With each setting of the concentration, the simulation is run and the
%resulting data is checked to see if a second action potential occurred. 
%This is determined by checking whether the propagation of the potential down
%the length of the axon is decrement-free, as would be expected in an action
%potential. When the refractory period is determined, the resulting data is 
%saved to data files which can then be used to display the results in the 
%Comparisons figure.

%  Before using the script, make sure that the Hodgkin-Huxley
%biophysics software has been started.
%If it hasn't been started, start the software using SOFTCELL.M
%      softcell('pap');

%APscr('set') will list the parameters that can be changed

%  Since a lot of files will be created, it is best to store them in a
%separate folder.  Examples of foldername are
%foldername = 'C:\data\temp'; %Windows
%foldername = '/user/tsb/data/temp'; %UNIX
%foldername = 'MacHD:matlab:data:temp'; %Mac
foldername = pwd;

% First reset all parameters to their defaults:
APstim('default');
APparams('default');
APnum('default');

% Make sure the duration of the simulation is long enough to get in two action
% potentials:
APscr('set','dur',10);

%  Set up a loop that samples through a concentration range.
%In this case, we use the range c_K_o = [0.11 10.11 20.11 30.11] mmol/L
for c_K_o = 0.11:10:30.11

  %set the model's temperature parameter
  newconc = APscr('set','c_K_o',c_K_o);

  %give yourself a note: indicate in the console what the settings are
  disp(['Concentration is set to ' num2str(newconc)]);

  % Set up a loop to figure out the refractory period by sampling through a
  %range of interpulse periods (the parameter that we are changing here
  %is the start time of pulse2) from 0 until we find the refractory period
				
  no_second_AP = 1;
  t = -0.1;

  while (no_second_AP == 1),
    t = t + 0.1;
    APscr('set','stimulus','pulse2','amplitude',0.05);
    APscr('set','stimulus','pulse2','duration',0.5);
    newInterval = APscr('set','stimulus','pulse2','start',t);
    disp(['Interpulse interval is set to ' num2str(newInterval)]);

    %perform simulation
    APscr('start');

    Vm1 = APscr('get','V_{ m}',1.5,':');
    % only look at values of Vm _after_ the first AP (from ~2ms on)
    peak1 = max(Vm1(200:end));
    Vm2 = APscr('get','V_{ m}',2.5,':');
    peak2 = max(Vm2(200:end));

    % Compare the peaks of Vm at z=1.5cm and z=2.5cm to see if the propagation
    % is decrement-free, and therefore an AP
    if ((peak1 - peak2) < 10) & (min([peak1;peak2]) > 10),
      disp(['Action potential detected.']);
      no_second_AP = 0;
    else
      disp(['No action potential detected.']);
    end
    
  end 
  
  %create a unique filename that contains the simulation data
  filename = ['/c_K_o' num2str(c_K_o) 't' num2str(t) '.mat']; 

  %save the simulation data to the file
  APscr('save','simulation',foldername,filename);
end

%That's it -- all the data has been saved.

%Now, put up all the data simultaneously by following these steps:
%1. call up the "Graphics" figure,
%2. using "Setup", select whichever variable(s) to plot,
%3. still inside "Setup HH plot", select the source of "Get values from"
%   to be the "Files with parameters".
% 3a. using "Select files", add the files to the list "Use these files".
% 3b. select "Overlay" to view the selected variable(s) in all the datafiles.