%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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 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)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);