##########################################################################
##########################################################################
##########################################################################
#
# This script generates figure 2 of the paper:
#
# "Spike shape and synaptic-amplitude distribution interact to set ...
# the high-frequency firing-rate response of neuronal populations"
# MJE Richardson (2018) to appear in Physical Review E.
#
# The code is in Julia version 1.0 (https://julialang.org/)
#
##########################################################################
#
# COPYRIGHT: Magnus Richardson (2018).
# 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 version 3 of the License.
#
# 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.
#
#
##########################################################################
##########################################################################
##########################################################################
using PyPlot
using Random
using SpecialFunctions
function MakeFig2()
Random.seed!(1)
K=1000.0
######################################################################
# The neuronal parameters
######################################################################
tau=20
vT=10
vre=5
vth=30
dT=0.6
asA=0.2
asB=0.6
asC=1.8
dv=0.0001;
# finds the synaptic rate for a given output rate for cases A, B and C
r0aim=5.0/1000
rA,RsA=Rintorout(r0aim,tau,vT,dT,vre,vth,asA,dv)
rB,RsB=Rintorout(r0aim,tau,vT,dT,vre,vth,asB,dv)
rC,RsC=Rintorout(r0aim,tau,vT,dT,vre,vth,asC,dv)
println("rA=$(round(rA*K;digits=4))Hz and RsA=$(round(RsA;digits=4))kHz")
println("rB=$(round(rB*K;digits=4))Hz and RsB=$(round(RsB;digits=4))kHz")
println("rC=$(round(rC*K;digits=4))Hz and RsC=$(round(RsC;digits=4))kHz")
######################################################################
######################################################################
# Generate the data for each panel
######################################################################
######################################################################
######################################################################
# Generate data for figure 2A
######################################################################
# methods fig for the paper
asE=asC
RsE=RsC
tE,RstE,vE,rtE,tEb,rtEsimb,fE,PE,TE=nicemethodsfig(dv,vth,vre,RsE,asE,tau,vT,dT);
######################################################################
# Generate data for figure 2B
######################################################################
# ThIn method for steady-state and modulation
nfth=33; lfth=0.1/1000; hfth=1000/1000;
fth=exp.(collect(range(log(lfth),stop=log(hfth),length=nfth)))
w=2*pi*fth
rhA=zeros(Complex{Float64},nfth)
rhB=zeros(Complex{Float64},nfth)
rhC=zeros(Complex{Float64},nfth)
for k=1:nfth
~,rhA[k],~,~,~=EshotEIFThIn(dv,vth,vre,RsA,asA,tau,vT,dT,fth[k],1)
~,rhB[k],~,~,~=EshotEIFThIn(dv,vth,vre,RsB,asB,tau,vT,dT,fth[k],1)
~,rhC[k],~,~,~=EshotEIFThIn(dv,vth,vre,RsC,asC,tau,vT,dT,fth[k],1)
end
ahA=abs.(rhA); ahB=abs.(rhB); ahC=abs.(rhC)
phA=angle.(rhA); phB=angle.(rhB); phC=angle.(rhC)
# calculate the frequency domain asymptotics
# case A as
dT
~,vC,PC,~=EshotEIFThIn0(dv,vth,vre,RsC,asC,tau,vT,dT);
Is=sum(dv*PC.*exp.((vC .-vT)/asC))/(rC*tau);
r1thasymC=rC[1]*tau*Is*gamma(1+dT/asC)./(1im*w*tau).^(dT/asC);
asymA=abs.(r1thasymA);
thasymA=-ones(length(w))*pi/2;
asymB=abs.(r1thasymB); thasymB=-pi/2 .+atan.(pi./(2*log.(c*w*tau)));
asymC=abs.(r1thasymC); thasymC=-ones(length(w))*pi*dT/(2*asC);
######################################################################
######################################################################
# Now plot the figures
######################################################################
######################################################################
# set up the figure
fig2=figure(figsize=(10,3.5));
# figure parameters
# useful constants
K=1000;
r2d=360/(2*pi)
lw=0.75
be=0.125
ms=3;
labelfs=8
panelfs=13
tickfs=8
textfs=6
######################################################################
# Plot Figure 2A
######################################################################
pwA=0.25
phAA=0.125
leA=0.075
ax1=PyPlot.axes([leA,0.75,pwA,phAA])
plot(tE,RstE,"k-",linewidth=lw)
plot([tE[1],tE[end]],RsC*[1,1],"k:",linewidth=lw)
ax1[:set_yticks]([0,0.1,0.2])
axis([TE-2*PE,TE,0,0.2])
ylabel(L"R_s(t)~"*" (kHz)",fontsize=labelfs)
ax2=PyPlot.axes([leA,0.45,pwA,0.175])
plot(tE,vE,linewidth=lw)
ylabel("Voltage (mV)",fontsize=labelfs)
ax2[:set_yticks]([0,10,20,30])
axis([TE-2*PE,TE,-5,30])
ax3=PyPlot.axes([leA,be,pwA,0.225])
plot(tE,K*rtE,linewidth=lw,"k")
plot([tE[1],tE[end]],5*[1,1],"k:",linewidth=lw)
fill_between(tEb,0,K*rtEsimb,linewidth=lw,color="lightgray")
plot(tEb,K*rtEsimb,linewidth=lw,color="gray")
axis([TE-2*PE,TE,0,7.5]);
ax3[:set_yticks]([0,2.5,5,7.5])
ylabel("Rate "*L"r(t)"*" (Hz)",fontsize=labelfs)
xlabel("Time (ms)",fontsize=labelfs)
for ax in (ax1,ax2,ax3)
ax[:spines]["top"][:set_visible](false)
ax[:spines]["right"][:set_visible](false)
ax[:tick_params](axis="both",labelsize=tickfs)
end
fig2[:text](0.02,0.925, "A", fontsize=panelfs)
######################################################################
# Plot Figure 2B
######################################################################
pwB=0.225
phBB=0.325
leB=0.45
ax4=PyPlot.axes([leB,0.575,pwB,phBB])
loglog(K*fth,abs.(rhA/rhA[1]),"b-",linewidth=lw)
loglog(K*fth,abs.(rhB/rhB[1]),"g-",linewidth=lw)
loglog(K*fth,abs.(rhC/rhC[1]),"r-",linewidth=lw)
loglog(K*fth,asymA/ahA[1],"b:",linewidth=lw)
loglog(K*fth,asymB/ahB[1],"g:",linewidth=lw)
loglog(K*fth,asymC/ahC[1],"r:",linewidth=lw);
axis([K*lfth,K*hfth,0.02,1.2]);
ylabel("Normalised amplitude",fontsize=labelfs)
ax4[:text](100,0.035,"case (i)",fontsize=textfs)
ax4[:text](100,0.025,L"a_s<\delta_\mathrm{T}",fontsize=textfs)
ax4[:text](500,0.12,"case (ii)",fontsize=textfs)
ax4[:text](800,0.08,L"a_s=\delta_\mathrm{T}",fontsize=textfs)
ax4[:text](300,0.5,"case (iii)",fontsize=textfs)
ax4[:text](300,0.35,L"a_s>\delta_\mathrm{T}",fontsize=textfs)
ax5=PyPlot.axes([leB+0.05,0.68,0.06,0.15])
loglog(K*fth,K*abs.(rhA),"b-",linewidth=lw)
loglog(K*fth,K*abs.(rhB),"g-",linewidth=lw)
loglog(K*fth,K*abs.(rhC),"r-",linewidth=lw)
loglog(K*fth,K*asymA,"b:",linewidth=lw)
loglog(K*fth,K*asymB,"g:",linewidth=lw)
loglog(K*fth,K*asymC,"r:",linewidth=lw)
ylabel("Amp. (Hz)",fontsize=7)
ax5[:set_yticks]([1,10,100])
xlabel("Freq. (Hz)",fontsize=7,labelpad=0)
ax5[:set_xticks]([1,1000])
#label_params(axis="x", pad=-1)
yticks(rotation=90,fontsize=6)
xticks(fontsize=6)
axis([1,K*hfth,1.0,100]);
#xlabel("f (Hz)",fontsize=labelfs)
ss=findall(K*fth.>100);
ax6=PyPlot.axes([leB,be,pwB,phBB])
semilogx(K*fth,phA*r2d,"b-",linewidth=lw)
semilogx(K*fth,phB*r2d,"g-",linewidth=lw)
semilogx(K*fth,phC*r2d,"r-",linewidth=lw)
semilogx(K*fth,thasymA*r2d,"b:",linewidth=lw)
semilogx(K*fth[ss],thasymB[ss]*r2d,"g:",linewidth=lw)
semilogx(K*fth,thasymC*r2d,"r:",linewidth=lw)
axis([K*lfth,K*hfth,-120,0]);
xlabel("Modulation frequency (Hz)",fontsize=labelfs)
ylabel("Phase (deg)",fontsize=labelfs)
ax6[:set_yticks]([0,-30,-60,-90,-120])
ax6[:text](100,-70,"(i)",fontsize=textfs)
ax6[:text](200,-52,"(ii)",fontsize=textfs)
ax6[:text](400,-25,"(iii)",fontsize=textfs)
for ax in (ax4,ax6)
ax[:spines]["top"][:set_visible](false)
ax[:spines]["right"][:set_visible](false)
ax[:tick_params](axis="both",labelsize=tickfs)
end
fig2[:text](0.375,0.925, "B", fontsize=panelfs)
######################################################################
# Plot Figure 2C
######################################################################
pwC=0.15
phCC=0.325
leC=0.8
xx=collect(0:0.01:10)
yy=1*(xx.<1.0)+(xx.>=1.0)./xx
zz=-90*(xx.<1.0)-(xx.>=1.0)*90 ./xx
ax7=PyPlot.axes([leC,0.575,pwC,phCC])
plot(xx,yy,"k-",linewidth=lw)
plot(asA/dT,1,"bo",markersize=ms)
plot(asB/dT,1,"go",markersize=ms)
plot(asC/dT,dT/asC,"ro",markersize=ms)
fill_between([0,1],[0,0],[1.5,1.5],color="lightgray")
xlabel(L"a_s / \delta_\mathrm{T}",fontsize=labelfs)
ylabel("Exponent "*L"\beta",fontsize=labelfs)
axis([0,4,0,1.5])
ax7[:text](0.1,0.4, "range of ", fontsize=textfs)
ax7[:text](0.1,0.25, "diffusion ", fontsize=textfs)
ax7[:text](0.1,0.1, "approximation", fontsize=textfs)
ax7[:text](0.3,1.1,"(i)",fontsize=textfs)
ax7[:text](1.1,1.1,"(ii)",fontsize=textfs)
ax7[:text](3,0.45,"(iii)",fontsize=textfs)
ax8=PyPlot.axes([leC,be,pwC,phCC])
plot(xx,zz,"k-",linewidth=lw)
plot(asA/dT,-90,"bo",markersize=ms)
plot(asB/dT,-90,"go",markersize=ms)
plot(asC/dT,-90*dT/asC,"ro",markersize=ms)
fill_between([0,1],[0,0],[-120,-120],color="lightgray")
xlabel(L"a_s / \delta_\mathrm{T}",fontsize=labelfs)
ylabel("phase (deg)",fontsize=labelfs)
axis([0,4,-120,0])
ax8[:set_yticks]([0,-30,-60,-90,-120])
ax8[:text](0.3,-105,"(i)",fontsize=textfs)
ax8[:text](1.1,-105,"(ii)",fontsize=textfs)
ax8[:text](3,-45,"(iii)",fontsize=textfs)
for ax in (ax7,ax8)
ax[:spines]["top"][:set_visible](false)
ax[:spines]["right"][:set_visible](false)
ax[:tick_params](axis="both",labelsize=tickfs)
end
fig2[:text](0.725,0.925, "C", fontsize=panelfs)
######################################################################
# Optional save
######################################################################
if true
savefig("fig2.pdf")
end
end
#########################################################################
#########################################################################
#########################################################################
# Subfunctions called by the main function
#########################################################################
#########################################################################
#########################################################################
#########################################################################
#########################################################################
# finds the steady-state input rates for a particular output rate
#########################################################################
#########################################################################
function Rintorout(r0aim,tau,vT,dT,vre,vth,as,dv)
# first run
Rsl=0.01; Rsh=5; Rsn=30;
RsL=collect(range(Rsl,stop=Rsh,length=Rsn))
rL=zeros(Rsn)
for k=1:Rsn
rL[k],~,~,~=EshotEIFThIn0(dv,vth,vre,RsL[k],as,tau,vT,dT)
end
s1=maximum(findall(rL.r0aim)); R2=RsL[s2]
# second run
Rsl=R1; Rsh=R2; Rsn=30;
RsL=collect(range(Rsl,stop=Rsh,length=Rsn))
rL=zeros(Rsn)
for k=1:Rsn
rL[k],~,~,~=EshotEIFThIn0(dv,vth,vre,RsL[k],as,tau,vT,dT)
end
s1=maximum(findall(rL.r0aim)); R2=RsL[s2]; r2=rL[s2]
Rstar=R1+(r0aim-r1)*(R2-R1)/(r2-r1);
rstar,~,~,~=EshotEIFThIn0(dv,vth,vre,Rstar,as,tau,vT,dT)
return rstar,Rstar
end
#########################################################################
#########################################################################
# ThIn for steady-state EIF with one shot-noise process
#########################################################################
#########################################################################
function EshotEIFThIn0(dv,vth,vre,Re,ae,tau,vT,dT)
# get the fixed points
vl,vu=MyfRoots.((0.0,vT+dT),vT,dT,10)
# set up the voltage
v=collect(vl+dv/2:dv:vth-dv/2)
F=(dT*exp.((v .-vT)/dT)-v)/tau
n=length(v)
# above and below the reset
krep=minimum(findall(v.>vre))
krem=krep-1
# above and below upper fixed point
nup=minimum(findall(v.>vu))
num=nup-1
# above the lower fixed point (NB should be entry 1)
nlp=minimum(findall(v.>vl))
p=zeros(n)
je=zeros(n)
# useful quantities at the upper fixed point
dFdvu=(exp((vu .-vT)/dT)-1)/tau
jeu=1
pu=(jeu/ae)/(Re+dFdvu)
#improve the initial conditions around vu - good for stability
je[num]=jeu-(vu-v[num])*(Re*pu-jeu/ae);
je[nup]=jeu+(v[nup]-vu)*(Re*pu-jeu/ae);
#########################################################################
# integrate vs <-----<< vu ..... vth
#########################################################################
for k=num:-1:nlp+1
p[k]=((v[k]>vre)-je[k])/F[k]
je[k-1]=je[k]-dv*(Re*p[k]-je[k]/ae)
end
p[nlp]=((v[nlp]>vre)-je[nlp])/F[nlp];
#########################################################################
# integrate vs ....... vu >>--------> vth
#########################################################################
for k=nup:n-1
p[k]=((v[k]>vre)-je[k])/F[k]
je[k+1]=je[k]+dv*(Re*p[k]-je[k]/ae)
end
p[n]=((v[n]>vre)-je[n])/F[n]
r=1.0/sum(p*dv)
P=r*p
Je=r*je
return r,v,P,Je
end
#########################################################################
#########################################################################
# Halley's method for lower and upper fixed points
#########################################################################
#########################################################################
function MyfRoots(v1,vT,dT,n)
v=zeros(n)
v[1]=v1
for k=1:n-1
expk=exp((v[k]-vT)/dT)
f0=dT*expk-v[k]
f1=expk-1
f2=(1/dT)*expk
v[k+1]=v[k]-2*f0*f1/(2*f1^2-f0*f2)
end
return v[end]
end
#########################################################################
#########################################################################
# binning function
#########################################################################
#########################################################################
function binnit(dt,t,tbin,x)
sbin=Int(ceil(tbin/dt))
sbin2=Int(sbin/2)
tb=t[sbin2:sbin:end-sbin2]
Cx=dt*cumsum(x)
Cxb=Cx[sbin:sbin:end]
xb=[Cxb[1];diff(Cxb)]/tbin
return tb,xb
end
#########################################################################
#########################################################################
# nice methods fig
#########################################################################
#########################################################################
function nicemethodsfig(dv,vth,vre,RsE,asE,tau,vT,dT)
dtE=0.05
TE=500
tE=collect(0:dtE:TE)
ntE=length(tE)
fE=10.0/1000
PE=1.0/fE
gfE=0.02
RstE=RsE .+ gfE*sin.(2*pi*fE*tE)
# get the analytical result
rE,rhE,~,~,~=EshotEIFThIn(dv,vth,vre,RsE,asE,tau,vT,dT,fE,1)
ahE=abs(rhE)
phE=angle(rhE)
rtE=rE .+ gfE*ahE*sin.(2*pi*fE*tE .+phE)
# do some sims
N=2000
tbinE=2
nrunsE=10
~,vE,rtEsim=EIFshotRtsim(vth,vre,tau,vT,dT,tE,RstE,asE,N)
for k=2:nrunsE
~,vEx,x=EIFshotRtsim(vth,vre,tau,vT,dT,tE,RstE,asE,N)
rtEsim=rtEsim+x
vE=[vE vEx]
end
tEb,rtEsimb=binnit(dtE,tE,tbinE,rtEsim)
rtEsimb=rtEsimb/nrunsE
return tE,RstE,vE,rtE,tEb,rtEsimb,fE,PE,TE
end
#########################################################################
#########################################################################
# ThIn for modulated shot EIF
#########################################################################
#########################################################################
function EshotEIFThIn(dv1,vth,vre,Re,ae,tau,vT,dT,fkHz,r1thG)
# f to w
w=2*pi*fkHz
# get the fixed points
vl,vu=MyfRoots.((0.0,vT+dT),vT,dT,10)
# match grid to fixed points
su=Int(ceil((vu-vl)/dv1)+1) # this is the unstable fixed point grid number
dv=(vu-vl)/(su-1)
v=collect(vl:dv:vth)
n=length(v)
sre=Int(ceil((vre-vl)/dv)+1) # first point above vre
# the forcing term and useful functions of it (on half sites)
vh=v .+dv/2
fh=(dT*exp.((vh .-vT)/dT)-vh)/tau
X=dv./fh
Xw=exp.(-1im*w*X)
Xe=exp.(-Re*X)
Xe[su]=0
vh=v .-dv/2
fh=(dT*exp.((vh .-vT)/dT)-vh)/tau
X=dv./fh
Yw=exp.(1im*w*X)
Ye=exp.(Re*X)
Ye[su]=0
#########################################################################
# steady-state first
#########################################################################
# the value an vu is 1
je=ones(n)
# integrate vs ....... vu >>--------> vth
for k=su:n-1
je[k+1]=exp(-dv/ae)*(je[k]*Xe[k]+(1-Xe[k]))
end
# integrate vs <-----<< vu ..... vth
for k=su:-1:2
je[k-1]=exp(dv/ae)*(je[k]*Ye[k]+(k>sre)*(1-Ye[k]));
end
# calculation of p using djedv+je/ae=ReP
djedv=[diff(je);0]/dv
p=(djedv+je/ae)/Re
# normalization
r=ae*Re/sum(je*dv)
P=r*p
Je=r*je
#########################################################################
# Now the modulation
#########################################################################
Aj,Aje=ones(Complex{Float64},n),ones(Complex{Float64},n)
Ej,Eje=zeros(Complex{Float64},n),zeros(Complex{Float64},n)
#########################################################################
# integrate from vlb to vs
# vs ...... vu ------> vth
#
# for the A and E trial solns at vu we have
# Aj=1 Aje=1
# Ej=0 Eje=0
#########################################################################
for k=su:n-1
Aj[k+1]=Aj[k]*Xw[k]+Aje[k]*(1-Xw[k])
Aje[k+1]=exp(-dv/ae)*(Aje[k]*Xe[k] + Aj[k]*(1-Xe[k]))
Ej[k+1]=Ej[k]*Xw[k]+Eje[k]*(1-Xw[k])
Eje[k+1]=exp(-dv/ae)*(Eje[k]*Xe[k] + Ej[k]*(1-Xe[k])) + dv*P[k]
end
# at threshold require rj=a*Aj=1 and ej=aa*Aj+Ej=0
a=1/Aj[n]; rj=a*Aj; rje=a*Aje;
aa=-Ej[n]/Aj[n]; ej=aa*Aj+Ej; eje=aa*Aje+Eje;
#########################################################################
# now integrate from vu to vs
# vs <----- vu ..... vth
#########################################################################
r1thGA=r1thG*0.999
r1thGB=r1thG*1.111
jA,jB=zeros(Complex{Float64},n),zeros(Complex{Float64},n)
jeA,jeB=zeros(Complex{Float64},n),zeros(Complex{Float64},n)
jA[su]=r1thGA*rj[su]+ej[su]
jB[su]=r1thGB*rj[su]+ej[su]
jeA[su]=r1thGA*rje[su]+eje[su]
jeB[su]=r1thGB*rje[su]+eje[su]
for k=su:-1:2
jA[k-1]=jA[k]*Yw[k] + jeA[k]*(1-Yw[k]) - r1thGA*(k==sre)
jeA[k-1]=exp(dv/ae)*jeA[k]*Ye[k] + (jA[k])*(1-Ye[k]) - dv*P[k]
jB[k-1]=jB[k]*Yw[k] + jeB[k]*(1-Yw[k]) - r1thGB*(k==sre)
jeB[k-1]=exp(dv/ae)*jeB[k]*Ye[k] + (jB[k])*(1-Ye[k]) - dv*P[k]
rj[k-1]=rj[k]*Yw[k]+rje[k]*(1-Yw[k])-(k==sre)
rje[k-1]=exp(dv/ae)*rje[k]*Ye[k]+(rj[k])*(1-Ye[k])
ej[k-1]=ej[k]*Yw[k]+eje[k]*(1-Yw[k])
eje[k-1]=exp(dv/ae)*eje[k]*Ye[k]+ej[k]*(1-Ye[k])-dv*P[k]
end
re=(r1thGA*jB[1]-r1thGB*jA[1])/(jB[1]-jA[1])
j=zeros(Complex{Float64},n)
je=zeros(Complex{Float64},n)
# check this below
j[1:su]=(jA[1:su]*jB[1]-jB[1:su]*jA[1])/(jB[1]-jA[1])
je[1:su]=(jeA[1:su]*jB[1]-jeB[1:su]*jA[1])/(jB[1]-jA[1])
j[su+1:end]=rj[su+1:end]*re+ej[su+1:end]
je[su+1:end]=rje[su+1:end]*re+eje[su+1:end]
return r,re,v,j,je
end
#########################################################################
#########################################################################
# response to a time course t, R(t)
# uses a second-order RKG
#########################################################################
#########################################################################
function EIFshotRtsim(vth,vre,tau,vT,dT,t,Rt,ae,N)
dt=t[2]-t[1]
ns=length(t)
Xt=Rt*dt
S=zeros(ns) # spike time series
v=zeros(ns) # voltage holder
for n=1:N
# create the noise vector
se=ae*(rand(ns).vth
v[k]=50; # decorative spike
S[k+1]=S[k+1]+1
v[k+1]=vre
end
end
end
r=S/(dt*N)
return t,v,r
end