 
% FHNNOISE integrates the FitzHugh Nagumo system 
% driven by sinusoidal forcing and colored noise. 
% The colored noise is the Ornstein-Uhlenbeck process eta(t). 
% Program averages results over navg realizations. 
% The program uses the integral Euler-Maruyama algorithm proposed in
% R.F. Fox et al., Physical Review A 38, 5938 (1988). 
% Solution can be seen using plot(t,v);
% An interspike interval histogram can be seen using plot(rhox,rhoy);
% The sequence of intervals can be plotted using plot(interv);
%
%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine
%
dt=0.0025;    % integration time step 
navg=5;    % number of sweeps (or ``realizations'') 
trans=1;  % number of time steps discarded as transients
% prmod=3;     % prmod: write to output solution vector u(1,tot) every 
% prmod time steps; this allows a plot of the solution that occupies less memory
tot=40000; %total number of integration time steps per realizations
intmax=100000; %maximum number of intervals generated 
a=0.5; % parameter a
b=0.15; % parameter b
d=1.0; % parameter d
ibias=0.04; %bias current
eps=0.005; % parameter epsilon 
amp=0.1; %amplitude of sine wave input
f=0.1; %frequency of sine wave input
dnz=0.0000001; %intensity of the Gaussian white noise process
tcor=0.01; %correlation time of the Ornstein-Uhlenbeck process
v0=0.0; %initial voltage
w0=0.0; %initial recovery state
thresh=0.5;  % spiking threshold
intervmin=0.4; %minimum time interval between two detected spikes (abs. refract. period)
t=ones(1,tot); %time vector
bins=200;  %number of bins to use for the interspike interval histogram 
rhoy=zeros(1,bins); %interspike interval histogram
transtime=trans*dt; %number of integration time steps discarded as transients
efac=exp(-dt/tcor); 
gausfac=sqrt(dnz*(1-efac*efac)/tcor);
v(1,1)=v0;
w(1,1)=w0;
nint=1;
interv=zeros(1,intmax);
for realiznumber=1:navg,
    realiznumber
    randn('state',sum(100*clock));
    eta=sqrt(dnz/tcor)*randn; %initial condition for Ornstein-Uhlenbeck process
    t1=trans*dt;
for i=1:tot,
    t(1,i)=(i-1)*dt;
end;    
v=v0*ones(1,tot);
w=w0*ones(1,tot);
for i=2:tot,
    gausx=randn*gausfac;
    etah=eta*efac+gausx;
    t(1,i)=(i-1)*dt;
    v(1,i)=v(1,i-1) +dt*(v(1,i-1)*(v(1,i-1)-a)*(1-v(1,i-1))-w(1,i-1)+ibias+amp*sin(2*pi*f*t(1,i))+etah)/eps;
    w(1,i)=w(1,i-1)+dt*(v(1,i-1)-d*w(1,i-1)-b);
    if (v(1,i)>thresh)&(v(1,i-1)<thresh), 
        int=t(1,i)-t1;
        if (int>intervmin)&(t(1,i)>transtime),
            interv(1,nint)=int;
            nint=nint+1;
            t1=t(1,i);
        end     
    end    
    eta=etah;
end;
end;
rhoy=hist(interv(1:nint),bins);
[rhodum,rhox]=hist(interv(1:nint),bins);
nint

