%% 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
navgs=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.1;     % bias current
eps=0.005;      % parameter epsilon
amp=0.001; % 0.1;        % amplitude of sine wave input
f=0.1;          % frequency of sine wave input
dnz=0.0000001; % 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;      % ode initial condition v
w(1,1)=w0;      % ode initial condition w
nint=1;
interv=zeros(1,intmax);
for realiznumber=1:navgs
    realiznumber
    rng(sum(100*clock)); % initialize the random number generator
    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