%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Steady state of the Generalised EIF Model (GEM)
%
% This script runs in MATLAB (and hopefully OCTAVE) and will generate 
% panels B,C,D and E of figure 1 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_1BCDE
global C EL Ex Ee Ei gL gx DT VT Vre Vth V xi tx;

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

% parameters for the EIF model
C=1; gL=C/20; EL=-80;           % leak current
DT=2; VT=-53; Vre=-60; Vth=0;   % spike currents, threshold and reset      

% general synaptic parameters
gs0=2*gL;  Ee=0; Ei=-70; sig0=4;                

% voltage-gated current parameters (most defined in function activ_var)
Vlb=-100; dV=0.05; V=[Vlb:dV:Vth]';
gx=2*gL;    Ex=-80;
[xi,tx]=activ_var(V);   % get the x_inf(V) and tau_x(V)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a plots showing the non-linear xi and tx (part of Fig 1B)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(4,2,1);  
plot(V,tx,'k-'); 
xlabel('V [mV]'); 
ylabel('tau_x(V) [ms]');
title(['Panel 1B']);
axis([-80 -20 0 80]); 

subplot(4,2,3); 
plot(V,xi,'k-');
xlabel('V [mV]'); 
ylabel('x_{inf}(V)');
axis([-80 -20 0 1.1]); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the example (used for Fig 1C)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fix the mean reversal to get ge0 and gi0
Es0EX=-30;     
ge0=gs0*(Es0EX-Ei)/(Ee-Ei);
gi0=gs0*(Ee-Es0EX)/(Ee-Ei);

[r0thEX,P0,x0,x0av,x0thEX]=findn(ge0,gi0,sig0);

subplot(2,2,3); 
plot(x0,x0,x0thEX,x0thEX,'o',x0,x0av,'b:');
title(['Panel 1C']);
text(0.5,0.35,['x_0=',num2str(x0thEX,2)]); 
text(0.1,0.9,['example E_{s0}=',num2str(Es0EX),'mV']);
axis([0 1 0 1]); 
xlabel('gating variable x_0^{(in)}'); 
ylabel('gating variable x_0^{(out)}');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% run over tonic synaptic drive to give Figs 1D and 1E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Es0=[-50:2.5:-20]'; 
nEs0=length(Es0);

r0th=zeros(nEs0,1);  
x0th=zeros(nEs0,1);

r0sim=zeros(nEs0,1); 
x0sim=zeros(nEs0,1);

yessim=1;  % optional simulations

for k=1:nEs0
    
    ge0=gs0*(Es0(k)-Ei)/(Ee-Ei);
    gi0=gs0*(Ee-Es0(k))/(Ee-Ei);
    
    [r0th(k),P0,x0,x0av,x0th(k)]=findn(ge0,gi0,sig0);
    
    if yessim
        [r0sim(k),x0sim(k)]=sim(ge0,gi0,sig0);
    end
    
    subplot(2,2,2);
    plot(Es0,K*r0th,Es0,K*r0sim,'r+',Es0EX,K*r0thEX,'o');
    legend('theory','simulation','example',2);
    xlabel('Es0 [mV]'); ylabel('rate [Hz]'); title('Panel 1D')
    
    subplot(2,2,4);
    plot(Es0,x0th,Es0,x0sim,'r+',Es0EX,x0thEX,'o');
    xlabel('Es0 [mV]'); ylabel('x0th'); title('Panel 1E');
    
    drawnow;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  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(ge0,gi0,x0,sig0)
global C EL Ex Ee Ei gL gx DT VT Vre xi tx V;

dV=V(2)-V(1); 
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);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   this gets the self-consistent solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r0,P0,x0,x0av,x0star]=findn(ge0,gi0,sig0)

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

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

[x0star,~]=intersects(x0,x0,x0av);
[r0,P0,~]=steady_rate(ge0,gi0,x0star,sig0);

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

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)-(ge0*(Ee-V)+gi0*(Ei-V));

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

Dx=5; Vx=-50;                           % voltage-dependent gating
xiV=1./(1+exp(-(V-Vx)/Dx));

tx1=50; tx2=20; Vt=-50; Dt=30;          % voltage-dependent time constant
txV=tx1 + tx2*exp(-(V-Vt).^2/(2*Dt));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  finds the intersection 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,x0]=sim(ge0,gi0,sig0)
global C Vre Vth EL;

T=7000; 
dt=0.05; 
ns=ceil(T/dt); 
t=dt*[1:ns]';

V=zeros(ns,1);
x=zeros(ns,1);
r=zeros(ns,1);

V(1)=EL;    %initial conditions
x(1)=0;

for k=1:ns-1
 
    [I,IF]=currents(V(k),x(k),ge0,gi0,sig0,dt);
    [xiV,txV]=activ_var(V(k));
    
    V(k+1)=V(k) + dt*(IF-I)/C;  
    x(k+1)=x(k) + dt*(xiV-x(k))./txV;
        
    if V(k+1)>Vth
        r(k+1)=1;
        V(k+1)=Vre;
    end

end

trelax=2000;             % time to allow transients to decay
srelax=ceil(trelax/dt);

x0=mean(x(srelax:ns));   % calculate means with transients ignored
r0=sum(r(srelax:ns))/(t(ns)-trelax);