%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Rate-response (at 1st order) of the Generalised EIF Model (GEM)
%
% This script runs in MATLAB (and hopefully OCTAVE) and will generate 
% panels B and C of figure 2 in the paper:
%
%  Dynamics of populations and networks of neurons with 
%  voltage-activated and calcium-activated currents  
%  Richardson MJE, Physical Review E  80: article-no 021928 (2009)
%
% Please reference this publication in any document that uses the code.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% COPYRIGHT: Magnus Richardson (2010).
% This code is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% See <http://www.gnu.org/licenses/> for a copy of  
% the GNU General Public License
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function GEM_fig_2BC
global C EL Ex Ee Ei gL gx ge0 ge1 gi0 sig0 DT VT Vre Vth V dV xi tx;

K=1000; % useful for converting rates from kHz to Hz

% parameters for the EIF model
C=1; gL=C/20; EL=-80;
DT=2; VT=-53; Vre=-60; Vth=0;

% synaptic parameters
gs0=2*gL;   Es0=-30; sig0=4;       
Ee=0;       Ei=-70;                
gi0=gs0*(Ee-Es0)/(Ee-Ei); 
ge0=gs0*(Es0-Ei)/(Ee-Ei); 
ge1=0.05*ge0;

% voltage-gated current parameters (most defined in function activ_var)
Vlb=-70; dV=0.01; V=[Vlb:dV:Vth]'; 
gx=2*gL;    Ex=-80;
[xi,tx]=activ_var(V);

lfth=0.1/1000; hfth=1000/1000; nfth=39;
fth=(exp(linspace(log(lfth),log(hfth),nfth))');

lfsim=1/1000; hfsim=100/1000; nfsim=5;
fsim=(exp(linspace(log(lfsim),log(hfsim),nfsim))');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The theory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[r0th,P0,x0,x0av,x0th]=findn;
[r1,re,rx,x1]=ge_mod(x0th,P0,fth);

Kar1=K*abs(r1); dr1=angle(r1)*360/(2*pi);
Kare=(ge1/gL)*K*abs(re); dre=angle(re)*360/(2*pi);
ax1=abs(x1); dx1=angle(x1)*360/(2*pi);

subplot(2,2,1); 
semilogx(fth*K,ax1,'-k');
xlabel('Freq. [Hz]'); 
ylabel('gating var. x, amplitude'); 
title(['2B']);
axis([min(K*fth) max(K*fth) 0 max(ax1)]);

subplot(2,2,3); 
semilogx(fth*K,dx1,'-k');
axis([min(K*fth) max(K*fth) -180 30]);
xlabel('Freq. [Hz]'); ylabel('gating var. x, phase [deg]'); 

subplot(2,2,2); 
semilogx(fth*K,Kar1,'-k',fth*K,Kare,'r:');
legend('theory','theory (static gating)',3);
xlabel(['Freq. [Hz]']); 
ylabel('rate r, amplitude [Hz]');
title(['2C']);
axis([min(K*fth) max(K*fth) 0 1.2*max(Kar1)]); 

subplot(2,2,4); 
semilogx(fth*K,dr1,'-k',fth*K,dre,'r:');
axis([min(K*fth) max(K*fth) -180 30]);
xlabel('Freq. [Hz]'); 
ylabel('rate r, phase [deg]'); 

drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optional simulations for comparison
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

yessim=1; 

if yessim
    
    T=4000; 
    dt1=0.05; 
    N=2000;
    r0sim=zeros(nfsim,1);
    arsim=zeros(nfsim,1);
    thrsim=zeros(nfsim,1);
    
    axsim=zeros(nfsim,1);
    thxsim=zeros(nfsim,1);
    
    for k=1:nfsim
        
        [r0sim(k),arsim(k),thrsim(k),axsim(k),thxsim(k)]=ge_mod_sim(fsim(k),T,dt1,N);
        
        subplot(2,2,1);
        semilogx(fth*K,ax1,'-k',fsim*K,axsim,'r*');
        xlabel('Freq. [Hz]');
        ylabel('gating var. x, amplitude');
        title(['2B']);
        axis([min(K*fth) max(K*fth) 0 max(ax1)]);
        
        subplot(2,2,3);
        semilogx(fth*K,dx1,'-k',fsim*K,thxsim*360/(2*pi),'r*');
        axis([min(K*fth) max(K*fth) -180 30]);
        xlabel('Freq. [Hz]'); ylabel('gating var. x, phase [deg]');
        
        subplot(2,2,2);
        semilogx(fth*K,Kar1,'-k',fth*K,Kare,'r:',fsim*K,arsim*K,'r*');
        legend('theory','theory (static gating)','simulations',3);
        xlabel(['Freq. [Hz]']);
        ylabel('rate r, amplitude [Hz]');
        title(['2C']);
        axis([min(K*fth) max(K*fth) 0 1.2*max(Kar1)]);
        
        subplot(2,2,4);
        semilogx(fth*K,dr1,'-k',fth*K,dre,'r:',fsim*K,thrsim*360/(2*pi),'r*');
        axis([min(K*fth) max(K*fth) -180 30]);
        xlabel('Freq. [Hz]');
        ylabel('rate r, phase [deg]');
        
        drawnow;
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  EIF steady-state using Threshold Integration of the Fokker-Planck Eq.
%  see Richardson, Phys. Rev. E 76: article-no 021919 (2007)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r0,P0,x0av]=steady_rate(x0)
global C EL Ex Ee Ei gL gx DT VT Vre xi tx V dV ge0 gi0 sig0;

n=length(V); kre=find(V==Vre);

G=(gL*(V-EL)+ge0*(V-Ee)+gi0*(V-Ei)+gx*x0*(V-Ex)-gL*DT*exp((V-VT)/DT))/(gL*sig0^2);
eG=exp(dV*G); F=(eG-1)./G; F(find(G==0))=dV; 

j0=zeros(n,1); j0(n)=1;
p0=zeros(n,1); p0(n)=0;
        
for k=n:-1:2

j0(k-1)=j0(k) - (k==kre+1);
p0(k-1)=p0(k)*eG(k) + j0(k)*(C/(gL*sig0^2))*F(k);

end

r0=1/(dV*sum(p0)); P0=r0*p0; 
x0av=sum(dV*P0.*xi./tx)/sum(dV*P0./tx);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Self-consistent steady-state solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r0,P0,x0,x0av,x0star]=findn

x0=[0:0.01:1]'; nx=length(x0);
x0av=zeros(nx,1);

for k=1:nx
    [r0,P0,x0av(k)]=steady_rate(x0(k));
end

[x0star,dummy]=intersects(x0,x0,x0av);
[r0,P0,dummy]=steady_rate(x0star);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  EIF 1st-order response from Fokker-Planck Eq. using Thresh. Int. Method
%  see Richardson, Phys. Rev. E 76: article-no 021919 (2007)
%  combined with the modulated intrinsic conductance gx. 
%  In this example the modulated afferent conductance is ge (excitation). 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r1,re,rx,x1]=ge_mod(x0,P0,freq)
global C EL Ex Ee Ei gL gx DT VT Vre xi tx V dV ge0 ge1 gi0 sig0;

dV=V(2)-V(1); n=length(V); 
kre=find(V==Vre);

w=2*pi*freq;    % NB this is in kHz 
nfreq=length(freq);

r1=zeros(nfreq,1); 
re=zeros(nfreq,1);
rx=zeros(nfreq,1);
x1=zeros(nfreq,1);

G=(gL*(V-EL)+ge0*(V-Ee)+gi0*(V-Ei)+gx*x0*(V-Ex)-gL*DT*exp((V-VT)/DT))/(gL*sig0^2);
eG=exp(dV*G); F=(eG-1)./G; F(find(G==0))=dV; 

for c=1:nfreq

    jr=zeros(n,1); jr(n)=1; pr=zeros(n,1); pr(n)=0;
    je=zeros(n,1); je(n)=0; pe=zeros(n,1); pe(n)=0;
    jx=zeros(n,1); jx(n)=0; px=zeros(n,1); px(n)=0;

    for k=n:-1:2

        jr(k-1)=jr(k) + dV*1i*w(c)*pr(k) - (k==kre+1);
        pr(k-1)=pr(k)*eG(k) + F(k)*(jr(k)*C)/(gL*sig0^2);

        je(k-1)=je(k) + dV*1i*w(c)*pe(k);
        pe(k-1)=pe(k)*eG(k) + F(k)*(je(k)*C/gL + (V(k)-Ee)*P0(k))/sig0^2;

        jx(k-1)=jx(k) + dV*1i*w(c)*px(k);
        px(k-1)=px(k)*eG(k) + F(k)*(jx(k)*C/gL + (V(k)-Ex)*P0(k))/sig0^2;

    end
    
    re(c)=-je(1)/jr(1); Pe=pe + re(c)*pr;
    rx(c)=-jx(1)/jr(1); Px=px + rx(c)*pr;
   
    xibyte=sum(dV*Pe.*xi./tx);
    xibytx=sum(dV*Px.*xi./tx);
    x0byte=x0*sum(dV*Pe./tx);
    x0bytx=x0*sum(dV*Px./tx);
    byt0=sum(dV*P0./tx);
       
    gammae1=ge1/gL;
    gammax=gx/gL;
    
    x1(c)=gammae1*(xibyte - x0byte)/(byt0 + i*w(c) - gammax*(xibytx-x0bytx));
    r1(c)=gammae1*re(c) + gammax*x1(c)*rx(c);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   calculates currents
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [I,IF]=currents(V,x,ge,gi0,dt)
global C EL Ex Ee Ei gL gx DT VT sig0;

IF=gL*sig0*sqrt(2*C/(gL*dt))*randn(1,length(V));    
I=gL*(V-EL)-gL*DT*exp((V-VT)/DT)+gx*x.*(V-Ex)-(ge*(Ee-V)+gi0*(Ei-V));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   the x_inf activation and timeconstant tau_x 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xiV,txV]=activ_var(V)

Dx=5; Vx=-50;                      
tx1=50; tx2=20; Vt=-50; Dt=30;     

xiV=1./(1+exp(-(V-Vx)/Dx));                          
txV=tx1 + tx2*exp(-(V-Vt).^2/(2*Dt));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   finds the intersections of two curves - only deals with one intersect
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xstar,astar]=intersects(x,y,z)

if y(1)<z(1)
    a=y; b=z;
else
    a=z; b=y;
end

s1=max(find(gt(b,a))); s2=s1+1;

if s2>length(x)
    s2=length(x);
    s1=s2-1;
end

x1=x(s1); x2=x(s2);
a1=a(s1); a2=a(s2);
b1=b(s1); b2=b(s2);

ma=(a2-a1)/(x2-x1);
mb=(b2-b1)/(x2-x1);

xstar=((b1-a1)+x1*(ma-mb))/(ma-mb);
astar=ma*(xstar-x1)+a1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    simulations 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r0,ar,thr,ax,thx]=ge_mod_sim(fkHz,T,dt1,N)
global C gL EL ge0 gi0 ge1 Vre Vth;
 
period=1/fkHz;                 % period in ms
nspp=ceil(period/dt1);         % number of steps per period
dt=period/nspp;                % the actual dt
np=ceil(T/period);             % the number of periods

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% oscillatory conductance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tp=dt*[0:nspp-1]';
ge1t=ge1*sin(2*pi*tp/period);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare the arrays
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r=zeros(nspp,1);
x=zeros(nspp,1);

nE=10;                  % number of examples
VE=zeros(nspp,nE);
xE=zeros(nspp,nE);

Vold=EL*ones(1,N); 
[xold,dummy]=activ_var(Vold);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a short equilibriation run without osc drive  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ns_therm=ceil(1000/dt);     

for k=1:ns_therm
    
    [I,IF]=currents(Vold,xold,ge0,gi0,dt);
    [xiV,txV]=activ_var(Vold);
                
    Vnew=Vold + dt*(IF-I)/C;
    xnew=xold + dt*(xiV-xold)./txV;
                
    spiked=find(Vnew>Vth);              
    Vnew(spiked)=Vre;
   
    Vold=Vnew;            % moves the system on
    xold=xnew;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loop over all the periods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for p=1:np

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % loop over a period
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for s=1:nspp

        [I,IF]=currents(Vold,xold,ge0+ge1t(s),gi0,dt);
        [xiV,txV]=activ_var(Vold);

        Vnew=Vold + dt*(IF-I)/C;
        xnew=xold + dt*(xiV-xold)./txV;

        spiked=find(Vnew>Vth);                                  
        r(s)=r(s)+length(spiked);
        Vold(spiked)=Vth;
        Vnew(spiked)=Vre;

        VE(s,:)=Vold(1:nE);
        xE(s,:)=xold(1:nE);

        Vold=Vnew;
        xold=xnew;

        x(s)=x(s)+sum(xnew);
    end

end

r=r/(N*np*dt);       
x=x/(N*np);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% extract the results for the fit for r
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r0=mean(r);
X=r-r0;

X1=(2*dt/period)*sum(X.*cos(2*pi*tp/period));
X2=(2*dt/period)*sum(X.*sin(2*pi*tp/period));

ar=sqrt(X1^2 + X2^2);

% put the angle between -180 and 180
shift=-((X2<0)&(X1>0) - (X2<0)&(X1<0))*pi;
thr=atan(X1/X2) + shift;

rfit=r0 + ar*sin(2*pi*tp/period + thr);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% extract the results for the fit for x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x0=mean(x);
X=x-x0;

X1=(2*dt/period)*sum(X.*cos(2*pi*tp/period));
X2=(2*dt/period)*sum(X.*sin(2*pi*tp/period));

ax=sqrt(X1^2 + X2^2);

% put the angle between -180 and 180
shift=-((X2<0)&(X1>0) - (X2<0)&(X1<0))*pi;
thx=atan(X1/X2) + shift;

xfit=x0 + ax*sin(2*pi*tp/period + thx);