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