TITLE Assembly of random channels
:  National Biomedical Simulation Resource, 12-7-90
:  A simple channel model with a random rate coefficient is solved repeatedly
:    and the average behavior shown at the conclusion. The "channel" model is a
:    simple decay reaction, with the gate variable being set to 1 initially and
:    the channel closing when gate drops below a threshold. The rate coefficient
:    for the decay is chosen from a set of random numbers with a normal
:    distribution, i.e. a known mean and standard deviation. The output
:    current is the average of the currents for all the individual channels,
:    an exponential decay from 1 toward zero.  This file uses a stepped
:    variable to do two complete simulations, one with a small number of
:    channels (showing "steps" of current) and one with a large number of
:    repetitions showing a smooth current curve. Note that this SCoP program
:    does not plot individual current curves, so there is no output until
:    all the individual solutions are completed and the average current is
:    calculated.

DEFINE ARRAYSIZE 101    : ARRAYSIZE determines the time resolution of the
			: current curve for output and cannot be changed at
			: run time.

PARAMETER
{
  max_time = 10.        : To change the upper limit of time for a simulation,
			: max_time must be set equal to the upper limit for
			: time in the INDEPENDENT block or the upper limit
			: for the time variable at run time.
  gate0 = 1
  cutoff = 0.3
}

STEPPED
{
  repetitions = 5,200
}

INDEPENDENT
{
  time FROM 0 TO 10 WITH 100    : These numbers govern only the final display
				: of average current.
}

PLOT current VS time

STATE
{
  gate FROM 0 TO 1.5
}

ASSIGNED
{
  delta_time
  decay_rate
  current FROM 0 TO 1.5
  accum[ARRAYSIZE]
}

INITIAL
{
  LOCAL index

  FROM count = 1 TO ARRAYSIZE
  {
    accum[count - 1] = 0
  }
  delta_time = max_time / (ARRAYSIZE-1) : The time step must fit with the
					: desired final time and the number
					: of storage locations in the accum
					: array.  If a smaller time step is
					: necessary, recompile with a larger
					: value for ARRAYSIZE.
  FROM count = 1 TO repetitions         : This outer loop solves the channel
					: model repeatedly.
  {
    decay_rate = normrand(0.5, 0.2)
    gate = gate0
    time = 0.                           : Because this simulation requires
					: repeated single channel solutions
					: for one average curve, we have to     
					: write our own solution loop for
					: each channel.
    WHILE (time <= max_time)
    {
      index = time * (ARRAYSIZE-1) / max_time + 0.5
					: index is converted to an integer
					: below when it is used to access
					: an element of the vector accum.
					: The 0.5 is for rounding.
      SOLVE rates
      IF (gate > cutoff)
      {
	accum[index] = accum[index] + 1
      }
      time = time + delta_time
    }
  }
}

KINETIC rates
{
~ gate ->  (decay_rate)                 : A more sophisticated channel model
					: can easily be substituted here.
}

BREAKPOINT
{
  LOCAL index                           : This block serves only to display
					: the average current at the end of
					: the simulation.
  index = time * (ARRAYSIZE-1) / max_time + 0.5
  current = accum[index] / repetitions

}