TITLE Essig's random channel model
:  National Biomedical Simulation Resource, 12-17-90

DEFINE ARRAYSIZE 101

PARAMETER
{
  max_time = 1.	: With your value of m and k, the "channel" decays in less
		: than 1 time unit (msec?), so I changed max_time and the
		: equivalent in INDEPENDENT from 10 to 1.
  gate0 = 10
  cutoff = 5
  m = 10
  k = 80
  d = 0
  r = 500	: I put in a real number here to see the effect of the
		: random coefficient.
}

STEPPED
{
  repetitions = 5,200
}

INDEPENDENT
{
  time FROM 0 TO 1 WITH 100	: Upper limit of 10 changed to 1
}

PLOT current VS time

STATE
{
  gate FROM -20 TO 20
}

ASSIGNED
{
  delta_time
  RANDOM FROM 0 TO 1
  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)
  FROM count = 1 TO repetitions
  {
    RANDOM = normrand(0.5, 0.2)		: You have to have a statement like
					: this to put random numbers in the
					: variable RANDOM. Putting the
					: statement here changes RANDOM each
					: individual channel simulation.
    gate = gate0
    gate' = 0.				: Since you have a 2nd order
					: differential equation, you also
					: have to set the value of the first
					: derivative of gate for each channel
					: simulation as well as gate itself.
    time = 0.
    WHILE (time <= max_time)
    {
      index = time * (ARRAYSIZE-1) / max_time + 0.5
      SOLVE rates METHOD runge
      IF (gate > cutoff)
      {
        accum[index] = accum[index] + 1
        time = time + delta_time
      }
      ELSE
      {
	time = max_time + 1
      }
    }
  }
}

DERIVATIVE rates
{
  gate'' = (-1/m) * (k * gate + d * gate' + r * RANDOM)
}

BREAKPOINT
{
  LOCAL index
  index = time * (ARRAYSIZE-1) / max_time + 0.5
  current = accum[index] / repetitions
}