function  [sys, x0, str, ts]  = liandfrn_sim(t,x,u,flag,dt,thres,sigma,C,R,ref)
%liandfrn_sim: S-function that implements an integrate-and-fire neuron
%model with  leak and refractory period and normally distributed noise.
%The model uses discrete integration time steps. During a
%simulation the time course of the various variables is computed as follows:
%
%	The membrane voltage is updated using:
%
%	Vm_i+1 = Vm_i(1 - dt/RC) + (dt/C) * I_i if Vm_i < thres
%	Vm_i+1 = 0 if Vm_i >= thres or if reftime_i > 0
%
%	The occurence of a spike is determined from:
%
%	spk_i = 1  if Vm_i >=  thres
%	spk_i = 0  if Vm_i < thres
%
%	The state variable reftime is updated as follows:
%
%	reftime_i+1 = n if Vm_i >= thres
%	reftime_i+1 = reftime_i - 1 if reftime_i > 0
%	reftime_i+1 = 0 otherwise
%	where n = ceil(ref/dt) is the refractory time in units of the
%	          time step (rounded up)
%   The noise term is updated according to:
%   noise = sigma*randn
%
%	The input variable is:
%	I_i = input current at time step i
%
%
%	The output is a two dimensional vector [spk_i, Vm_i].
%	The first component spk_i is  equal to 0 if no spike occured
%	or 1 if a spike occured. The second component is the membrane
%	voltage Vm (in mV) at time i. 
%
%	The parameters are:
%	
%	dt  = sampling step (in msec)
%	thres = mean threshold value (in mV)
%	sigma = standard deviation of the noise term
%	C   = capacitance (in nF)
%	R   = resistance of the model cell (in MOhms)
%	ref = refractory period (in msec)
%
%

if (abs(flag) == 2)    	% Discrete state update
  if ( x(1) >= thres )	% threshold reached, next voltage value is reset
    sys(1) = 0;
    sys(2) = max(ceil(ref/dt)-1,0); %sets the refractory period
  elseif ( x(2) > 0 )	% in refractory period, wait
    sys(1) = 0;
    sys(2) = x(2) - 1;
  else	%integrates the voltage
    sys(1) = x(1)*(1-dt/(R*C)) + (dt/C) * u + sigma*randn;
    sys(2) = 0;
  end

elseif ( flag == 3 )	%output vector required
  if ( x(1) >= thres )	%fire a spike
    sys = [1, x(1)];
  else
    sys = [0, x(1)];
  end

elseif ( flag  == 0 )	% Initialization

  if ( (R*C) < dt )     % Check for stability
    disp(' ');
    disp(' ');
    disp('Error: the time constant RC computed from the resistance');
    disp('and the capacitance of the model is smaller than the time');
    disp('step. Execution aborted.');
    error('???');
  end    
    
  % Return system sizes
  sys(1) = 0;	% 0 continuous states
  sys(2) = 2;	% 2 discrete states: membrane voltage and refractory
		%   period indicator
  sys(3) = 2;	% 2 outputs: spike occurence and membrane voltage
  sys(4) = 1;	% 1 input: current
  sys(5) = 0;	% 0 roots
  sys(6) = 0;	% no direct feedthrough
  sys(7) = 1;	% 1 sample time: dt
  
  x0 = [0, 0];	% Initialization of the membrane voltage and
		% refractory time

  ts = [dt, 0];

else 

  sys = [];

end


str = [];