%This script replicates Hasselmo et al., J. Neurosci. 1995.
%Plots neural activation "a" over time.  
%Plots interneuron activation "h" over time.
%Plots "a" vs. "h" in phase plot. 
%clear all;
%close all;

W=0.016;
H=0.06;
Wh=0.0042;
tha=8;
thh=8;
eta=0.01;
T=10000;
A=zeros(1,T);
a=zeros(1,T);
h=zeros(1,T);
for t=1:T;
    if t<5000; A(t)=0.1; else A(t)=0; end
    %aa=a(t); %No threshold on "a"
    if (a(t)-tha)>0; aa=a(t)-tha; else aa=0; end 
    %Comment out above line to remove threshold on "a".
    if (h(t)-thh)>0; hh=h(t)-thh; else hh=0; end
da=A(t)+W*aa-H*hh-eta*a(t);
dh=Wh*aa-eta*h(t);
a(t+1)=a(t)+da;
h(t+1)=h(t)+dh;
end

a(T)
h(T)

if max(a)>0
    maxfig=max(a);
else
    maxfig=50;
end

figure('Position',[40 40 300 300]);
axis([1 T 0 maxfig]); axis manual; hold on
figa=gcf; axa=gca; colormap(gray);
set(figa,'Name','a - activation','numbertitle','off');

figure('Position',[400 40 300 300]);
axis([1 T 0 maxfig]); axis manual; hold on
figb=gcf; axb=gca; colormap(gray);
set(figb,'Name','h - interneuron','numbertitle','off');

figure('Position',[760 40 300 300]);
axis([1 maxfig 0 maxfig]); axis manual; hold on
figc=gcf; axc=gca; colormap(gray);
set(figc,'Name','a vs. h phase plot','numbertitle','off');

axes(axa); plot(a); drawnow
axes(axb); plot(h); drawnow
axes(axc); plot(h,a); drawnow