%**********************************************************************************
%
%              CalC version 4.5.0, script "STF.decay.stage2.par"
%                     Victor Matveev, September 10, 2002
%
%      Supplemental material for the Letter to the Editor of the Biophys. J.
%           "New and Corrected Simulations of Synaptic Facilitation"
%                      V. Matveev, A. Sherman and R. Zucker
%__________________________________________________________________________________
%
%   This simulation script computes the decay time course of facilitation elicited
% by a 5-pulse train, as a function of the interval between the last two pulses
% (see Fig. 3 of the Letter), using the [Ca2+] files created by the "stage 1"
% simulation, "STF.decay.stage1.par". STF is measured as the ratio of the 5th
% reponse to the 1st response, minus one.
%
%   This script requires 2 command-line arguments, which are the same as the first
% two arguments of scripts "STF.growth.par" and "STF.decay.stage1.par" (they are
% only used to determine the names of [Ca2+] files created by stage 1 simulation):
%
% 1. The 1st argument specifies the parameter set used, and accepts values of
%     "A"  - for the quoted parameter set
%     "B"  - for the actual parameter set
%     "C"  - for the modified set
%     "D"  - for the parameter set with tortuosity
%
% 2. The 2nd argument specifies whether Fura-2 was included:
%     "with" or "w"      - Fura-2 is included
%     "without" or "wo"  - control simulation (no Fura-2)
%
% For example, to reproduce the control curve of Fig. 3D of the Letter, execute
% the following commands (see also the stage-1 script "STF.decay.stage1.par"):
%
%      calc STF.decay.stage1.par D without 10.0 800.0 25
%      calc STF.decay.stage2.par D without
%
% or, if using the xmgr(ace) data plotting application:
%
%      calc STF.decay.stage1.par D without 10.0 800.0 25
%      calc STF.decay.stage2.par D without | xmgr -pipe
%
%**********************************************************************************
%                 T H E   A L G O R I T H M   S U M M A R Y :
%                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% I. At 1st iteration (iteration = 0), perform these steps:
%
%   1. Find the steady-state of the system of ODEs describing the dynamics of the
%      Ca2+ binding model (equilibration).
%   2. Integrate the ODEs for the entire duration of the 5-pulse train, using the
%      [Ca2+] files generated by stage 1 script at 1st iteration. The interval
%      between the 4th and 5th APs corresponds to the 1st data point in the STF
%      decay time course curve.
%   3. Save the entire simulation state to the binary "dump" file at the end of the
%      inter-spike interval between the 4th and 5th pulses.
%   4. During the integration run, track the maximum of the release rate for the
%      1st and 5th pulses; this will yield the 1st data point in the STF decay time
%      course. Save the value of the 1st peak response to a file, to be used by
%      subsequent iterations.
%
% II. At each successive iteration, the sequence of operations mirrors the
%     corresponding steps in the "stage 1" simulation (see "STF.decay.stage2.par"):
%
%   1. Retreave the entire simulation state from the "dump" file, which was saved
%      there at the previous iteration.
%   2. Using the [Ca2+] files generated at the same iteration of stage 1 simulation,
%      integrate the ODEs for the duration of the interval between the 4th and 5th
%      pulses which, when combined with the corresponding interval in the previous
%      iteration, adds up to the full interval between the 4th and 5th pulses
%      corresponding to the current data point in the STF decay time course, plus
%      the duration of the final 5th pulse.
%   3. At the end of the last inter-spike interval (between 4th and 5th pulses),
%      save the entire simulation state to the binary "dump" file. In the course of
%      the simulation, this file will contain the state of the system at the end of
%      the progressively increasing interval between the 4th and 5th pulses.
%   4. Track the maximal release rate during time period corresponding to the last
%      action potential; using the value of 1st response saved at the 1st iteration,
%      calculate the next data point in the STF decay time course.
%
%==================================================================================
%                  B A S I C   C O N S T A N T S ,   E T C . :
%==================================================================================
% The following instruction tells the program that this script is only to solve the
% system of ordinary differential equations (ODEs) (the PDEs for the [Ca2+] have
% been already solved by the stage 1 simulation):
mode = ODE
% The next two lines simply define the constants corresponding to different user
% argument choices (we have to do this since boolean operations on strings are not
% yet implemented):
A = 1; B = 2; C = 3; D = 4
w = 1; with = 1; without = 2; wo = 2
% Define the path to the folder where the [Ca2+] data files were stored. This path
% should be the same as the path specified in the script "STF.decay.stage1.par"
% (you may also choose to set it using an additional command line argument, via
% statement "path = $4"):
path = ""
%   The two command-line arguments are only needed to determine the file prefix
% for the [Ca2+] data files produced by the stage 1 script:
if ($3 == with) then
  fString = "fura."      % a part of the filename prefix
else
  fString = "control."   % a part of the filename prefix
endif
% Now define the full prefix for all file names (e.g., "folder/actual.fura."):
if ($2 == A) then  prefix = path  "quoted."     fString;   endif
if ($2 == B) then  prefix = path  "actual."     fString;   endif
if ($2 == C) then  prefix = path  "modified."   fString;   endif
if ($2 == D) then  prefix = path  "tortuosity." fString;   endif
% For the parameter set with tortuosity, the affinity of the Y-site is 9 uM, not
% 3 uM as defined in the script "XY.ode.par" included below (it is important that
% Y.kon is redefined before the include statement: only the 1st definition of any
% constant is used; later definitions are discarded):
if ($2 == D) then  Y.kon = Y.koff / 9 ; endif
% Read the implementation of the Ca2+ binding part of the model from "XY.ode.par":
ode.script = path "XY.ode.par"
include ode.script
% The following file will contain the entire simulation data, to be stored at the
% end of each iteration and retreaved at the beginning of the next iteration:
dump.file = prefix "ode.dump"
% Next, the file that will hold the value of the response to the 1st pulse:
R1.file = prefix "R1"
% Next, the file that will hold the result - the STF decay time course:
Result.file = prefix "decay"
% The names of the [Ca2+] files will be formed by merging the file prefix, a suffix
% string "CaX."/"CaY.", and the value of the current interval between 4th and 5th
% action potentials ("dt.current" is the current interval between 4th and 5th pulses,
% which is varied in the main loop below):
CaX.file = prefix "CaX." dt.current "ms"
CaY.file = prefix "CaY." dt.current "ms"
% Now instruct the program to read variables CaX and CaY (appearing in the system
% of ODEs included from file "XY.ode.par") from these files :
CaX table CaX.file   % These instruction are performed at each iteration, which is
CaY table CaY.file   % why they appear here, outside of the "if" statement clauses
                     % following below
% Read the values of parameters "dt.first", "dt.last", and "dt.iterations", stored
% by the stage 1 script, "STF.decay.stage1.par":
log.file = prefix "log"
include log.file
dt.step = (dt.last - dt.first) / dt.iterations  % define the increment "dt.step"
%==================================================================================
%                     M A I N   S I M U L A T I O N   L O O P :
%==================================================================================
% Now we begin our loop over the increasing values of the time interval between the
% 4th and 5th pulses. The next statement is the main "for" loop command (notice that
% there is no loop termination statement, because the "for" statement can appear only
% once, and applies to this entire script; like almost all other statements, the
% "for" command can be put anywhere in the script file):
for iteration = 0 to dt.iterations step 1    % <---  THE MAIN LOOP STATEMENT
% The interval between the 4th and 5th pulses for the current iteration is given by
dt.current = dt.first + dt.step * iteration
if iteration == 0 then    %=========================================================
                          %  Instructions below are only performed at 1st iteration
                          %=========================================================
  % We must instruct the program to intialize the ODE variables to their steady-state
  % values, because at time zero the Ca2+ sensors should be in equilibrium with the
  % background [Ca2+]:
  equilibrate
  % At first iteration, the inter-pulse interval between 4th and 5th pulses ends at
  % time "T" given by
  T = 30 + dt.first
  % Total simulation time at the first iteration is T + duration of 5th pulse:
  Total = T + 1.2
  Run Total            % <-- Run the solver
  % Instruct the program to save the simulation variables at time T in the
  % "dump" file
  Export T dump.file
  % Define variables R1 and R2 that will track the maximal release rate during the
  % 1st and last pulses:
  R1 peak Release 0 1.2
  R5 peak Release T Total
  % Save the value of R1 to a file:
  print R1.file "R1 = " R1
else     % this else statement corresponds to "if iteration == 0"
  %===================================================================
  %  Instructions below are only performed at iteration > 0
  %===================================================================
  % 1. Retreave the simulation state from the "dump" file:
  Import dump.file
  % Read the amplitude of the response to the 1st pulse from a file created
  % at first iteration:
  include R1.file
  % Integrate the model ODEs for the duration given by the difference between
  % "dt.current" values at the current time step and the previous time step,
  % which is equal to dt.step, plus the duration of the 5th action potential:
  Total = dt.step + 1.2
  Run Total
  % Instruct the program to save the simulation variables at the end of the
  % silent period to the same simulation dump file:
  Export dt.step dump.file
  % Define the variable R5 that will track the maximal release rate during the
  % last pulse:
  R5 peak Release dt.step Total
endif      % this endif statement corresponds to "if iteration == 0"
%========================================================================
%    Instructions below are performed for all iterations:
%========================================================================
% Define the facilitation variable:
Facilitation := R5 / R1 - 1
% Save the result in file "Result.file", as a two-column ASCII file
% (arguments of "append" following the filename are concatenated):
append Result.file dt.current "  " Facilitation
% If the output is to be piped into the xmgr(ace) graphics application, uncomment
% the following lines:
% plot.method xmgr
% print stdout "@ g0.s0 POINT " dt.current ", " Facilitation
% print stdout "@ autoscale"
% print stdout "@ redraw"
%==================================================================================
%                                  T H E   E N D
%==================================================================================