########################################################################## ########################################################################## ########################################################################## # # This script generates figure 3 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 function MakeFig3() # useful constants K=1000; r2d=360/(2*pi) ###################################################################### # The neuronal parameters ###################################################################### tau=20 vT=10 vre=5 vth=15 # factors of 3 dT=0.6 asA=0.2 asB=0.6 asC=1.8 ###################################################################### ###################################################################### # Get the data ###################################################################### ###################################################################### # some choices for the level of discretization #dt=0.01; nz=500; dv=0.001 # instantaneous dt=0.01; nz=1000; dv=0.001 # about 20secs #dt=0.002; nz=1000; dv=0.0001 # about 2mins #dt=0.002; nz=2000; dv=0.0001 # about 3mins #dt=0.002; nz=5000; dv=0.0001 # about 8mins # 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") # 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))) # exp spacing 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); # match the amplitudes so steady response the same gfA=0.025*ahA[1]/ahA[1] gfB=0.025*ahA[1]/ahB[1] gfC=0.025*ahA[1]/ahC[1] # get the stimulus: double sampled for the RK4 t2,RsAt2=thestimulus(dt/2,RsA,gfA) t2,RsBt2=thestimulus(dt/2,RsB,gfB) t2,RsCt2=thestimulus(dt/2,RsC,gfC); t,rAt,RsAt=numint(t2,RsAt2,asA,tau,dT,vT,vre,vth,nz) t,rBt,RsBt=numint(t2,RsBt2,asB,tau,dT,vT,vre,vth,nz) t,rCt,RsCt=numint(t2,RsCt2,asC,tau,dT,vT,vre,vth,nz) T=t[end]; ###################################################################### ###################################################################### # Make the figure ###################################################################### ###################################################################### # set up the figure fig1=figure(figsize=(10,4)); # figure parameters clf() lw=0.75 be=0.125 ms=3; labelfs=8 panelfs=13 tickfs=8 textfs=6 lw=0.75 be=0.125 ###################################################################### # Plot fig 3a ###################################################################### leA=0.075 pw=0.85 phA=0.2 ax1=PyPlot.axes([leA,0.725,pw,phA]) plot(t,RsAt,"b-",linewidth=lw) plot(t,RsBt,"g-",linewidth=lw) plot(t,RsCt,"r-",linewidth=lw) axis([0,550,0,2.5]) ylabel("Stimulus R(t) (kHz)",fontsize=labelfs) ax1[:text](20,1.7,"case (i) " *L"R_s=2.11"*"kHz",fontsize=textfs) ax1[:text](20,0.8,"case (ii) " *L"R_s=0.59"*"kHz",fontsize=textfs) ax1[:text](20,0.25,"case (iii) "*L"R_s=0.14"*"kHz",fontsize=textfs) ax1[:text](141,1.5,"Sharpening Gaussian impulses",fontsize=textfs) ax1[:text](295,1.5,"square pulse ",fontsize=textfs) ax1[:text](420,1.5,"100Hz high gamma",fontsize=textfs) ax1[:text](500,1.5,"200Hz ripple",fontsize=textfs) ax1[:spines]["top"][:set_visible](false) ax1[:spines]["right"][:set_visible](false) ax1[:tick_params](axis="both",labelsize=tickfs) fig1[:text](0.01,0.925, "A", fontsize=panelfs) ###################################################################### # Plot fig 3b ###################################################################### phB=0.2 ax2=PyPlot.axes([leA,0.45,pw,phB]) plot(t,K*rAt,"b-",linewidth=lw) plot(t,K*rBt,"g-",linewidth=lw) plot(t,K*rCt,"r-",linewidth=lw) axis([0,550,0,6.5]) plot([0,T],K*rA*[1,1],"b:",linewidth=lw) plot([0,T],K*rB*[1,1],"g:",linewidth=lw) plot([0,T],K*rC*[1,1],"r:",linewidth=lw); xlabel("Time (ms)",fontsize=labelfs); ylabel("Response r(t) (Hz)",fontsize=labelfs); plot([18,51.15],2.5*[1,1],"k:",linewidth=lw) ax2[:text](55,2.25,"33ms",fontsize=textfs) ax2[:text](48,1.0,"(i)",fontsize=textfs) ax2[:text](30,1.0,"(ii)",fontsize=textfs) ax2[:text](12,1.0,"(iii)",fontsize=textfs) ax2[:spines]["top"][:set_visible](false) ax2[:spines]["right"][:set_visible](false) ax2[:tick_params](axis="both",labelsize=tickfs) fig1[:text](0.01,0.65, "B", fontsize=panelfs) ###################################################################### # Plot fig 3c ###################################################################### pwC=0.175 Cgap=0.225 phC=0.2 lC1=leA lC2=lC1+Cgap lC3=lC2+Cgap lC4=lC3+Cgap ax3=PyPlot.axes([lC1,0.1,pwC,phC]) plot(t,K*rAt,"b-",linewidth=lw) plot(t,K*rBt,"g-",linewidth=lw) plot(t,K*rCt,"r-",linewidth=lw) plot(t .-40,K*rAt,"b-",linewidth=lw) plot(t .-40,K*rBt,"g-",linewidth=lw) plot(t .-40,K*rCt,"r-",linewidth=lw) plot(t .-80,K*rAt,"b-",linewidth=lw) plot(t .-80,K*rBt,"g-",linewidth=lw) plot(t .-80,K*rCt,"r-",linewidth=lw) axis([135,180,4.8,7]); xlabel("Time (ms)",fontsize=labelfs); ylabel("r(t) (Hz)",fontsize=labelfs); ax3[:text](141,6.9,"Gaussian impulse responses",fontsize=textfs) ax4=PyPlot.axes([lC2,0.1,pwC,phC]) plot(t,K*rAt,"b-",linewidth=lw) plot(t,K*rBt,"g-",linewidth=lw) plot(t,K*rCt,"r-",linewidth=lw); axis([285,310,4.8,6]); xlabel("Time (ms)",fontsize=labelfs); ylabel("r(t) (Hz)",fontsize=labelfs); ax4[:text](294,5.2,"square-pulse responses",fontsize=textfs) ax5=PyPlot.axes([lC3,0.1,pwC,phC]) plot(t,K*rAt,"b-",linewidth=lw) plot(t,K*rBt,"g-",linewidth=lw) plot(t,K*rCt,"r-",linewidth=lw); axis([415,475,3.5,6.5]); xlabel("Time (ms)",fontsize=labelfs); ylabel("r(t) (Hz)",fontsize=labelfs); ax5[:text](420,6.5,"100Hz high gamma responses",fontsize=textfs) ax6=PyPlot.axes([lC4,0.1,pwC,phC]) plot(t,K*rAt,"b-",linewidth=lw) plot(t,K*rBt,"g-",linewidth=lw) plot(t,K*rCt,"r-",linewidth=lw); axis([490,550,3.5,6.5]); xlabel("Time (ms)",fontsize=labelfs); ylabel("r(t) (Hz)",fontsize=labelfs); ax6[:text](500,6.5,"200Hz ripple responses",fontsize=textfs) for ax in (ax3,ax4,ax5,ax6) ax[:spines]["top"][:set_visible](false) ax[:spines]["right"][:set_visible](false) ax[:tick_params](axis="both",labelsize=tickfs) end fig1[:text](0.01,0.3, "C", fontsize=panelfs) ###################################################################### # Optional save ###################################################################### if false savefig("fig3.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=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 ######################################################################### ######################################################################### # 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 ######################################################################### ######################################################################### # Solves for an arbitrary waveform using an RK4 algorithm ######################################################################### ######################################################################### function numint(t2,Rst2,as,tau,dT,vT,vre,vth,nz) z,gz,vz,dvdz=ztransformB(tau,dT,vT,vth,nz) dz=z[2]-z[1] t=t2[1:2:end] Rst=Rst2[1:2:end] RstH=Rst2[2:2:end] # half step needed for the RK4 dt=t[2]-t[1] nt=length(t) dz=z[2]-z[1] r=zeros(nt) # initial conditions and initialisation muv=1 sigv=0.1 Pz=exp.(-(vz .-muv).^2/(2*sigv^2)).*dvdz Pz=Pz./sum(dz*Pz) # smooth the insertion at vre with a Gaussian (in voltage space) # this has little effect on result but improves convergence muvre=vre sigvre=0.1 hzre=exp.(-(vz .-muvre).^2/(2*sigvre^2)).*dvdz hzre=hzre/sum(dz*hzre) F1=exp.((vT .-vz)/as) G1=1.0 ./F1 # loop over time for k=1:nt-1 Jz,Kz,dJdz=getarray(z,F1,G1,gz,dz,Pz,Rst[k]) k1=-dt*dJdz + dt*hzre*Jz[end] Jz,Kz,dJdz=getarray(z,F1,G1,gz,dz,Pz+k1/2,RstH[k]) k2=-dt*dJdz + dt*hzre*Jz[end] Jz,Kz,dJdz=getarray(z,F1,G1,gz,dz,Pz+k2/2,RstH[k]) k3=-dt*dJdz + dt*hzre*Jz[end] Jz,Kz,dJdz=getarray(z,F1,G1,gz,dz,Pz+k3,Rst[k+1]) k4=-dt*dJdz + dt*hzre*Jz[end] Pz=Pz+k1/6+k2/3+k3/3+k4/6 r[k]=Jz[end] end r[end]=r[end-1] return t,r,Rst end ######################################################################### ######################################################################### # get arrays ######################################################################### ######################################################################### function getarray(z,F1,G1,gz,dz,Pz,Ret) Kz=dz*Ret*F1.*(cumsum(Pz.*G1)-Pz[1]*G1/2-Pz.*G1/2) Jz=Kz + gz.*Pz dJdz=[(Jz[2]-Jz[1]);(Jz[3:end]-Jz[1:end-2])/2;Jz[end]-Jz[end-1]]/dz; return Jz,Kz,dJdz end ######################################################################### ######################################################################### # the stimulus which needs to be multiplied by a gain factor ######################################################################### ######################################################################### function thestimulus(dt,Rs,gf) t=collect(0:dt:570) nt=length(t) x=zeros(nt) # number of extra spikes for the Gaussian ng=10 # Gaussian pulse 1 tgauss=140; FWHM=4; siggau=FWHM/sqrt(8*log(2)) x=x + ng*exp.(-(t .-tgauss).^2/(2*siggau^2))/sqrt(2*pi*siggau^2) # Gaussian pulse 2 tgauss=190; FWHM=2; siggau=FWHM/sqrt(8*log(2)) x=x + ng*exp.(-(t .-tgauss).^2/(2*siggau^2))/sqrt(2*pi*siggau^2) # Gaussian pulse 3 tgauss=240; FWHM=1; siggau=FWHM/sqrt(8*log(2)) x=x + ng*exp.(-(t .-tgauss).^2/(2*siggau^2))/sqrt(2*pi*siggau^2) # step x=x + 2*(t.>290).*(t.<350) # chirp 1 - fast gamma tch1=445; sch1=10; fch1=100/1000; ampch1=5; x=x+ampch1*cos.(2*pi*fch1*t).*exp.(-(t .-tch1).^2/(2*sch1^2)) # chirp 2 - ripples tch2=520; sch2=10; fch2=200/1000; ampch2=7; x=x+ampch2*cos.(2*pi*fch2*t).*exp.(-(t .-tch2).^2/(2*sch2^2)) # Multiply by prefactor and add steady-state Rst=Rs .+ gf*x Rst[findall(t.<5)] .=0 return t,Rst end ######################################################################### ######################################################################### # z transform to get v for a constant array over z ######################################################################### ######################################################################### function ztransformB(tau,dT,vT,vth,nz) # get the bounds for v vl,vu=MyfRoots.((0.0,vT+dT),vT,dT,10) vh=vth # add ve to v and sample at dv for the spline v=collect(range(vl,stop=vh,length=nz)) psi=exp.((v .-vT)/dT) zv=(psi .+v/vT)./(psi .+v/vT .+1) # z on the v grid z=collect(range(zv[1],stop=zv[end],length=nz)) dz=z[2]-z[1] # do one step using Halley method for finding roots vz=vT*ones(nz) for k=1:50 psi=exp.((vz .-vT)/dT) H=psi .+vz/vT .-z./(1 .-z) dHdv=psi/dT .+1/vT d2Hdv2=psi/dT^2 vz=vz-(2*H.*dHdv)./(2*dHdv.^2-H.*d2Hdv2) end if isnan(vz[end]) vz[end]=vth end if sum(isnan.(vz))>0 println("ERROR! there are some values of z that are NaN") end # now the other quantities we need on the z grid psi=exp.((vz .-vT)/dT) dzdv=(psi/dT .+1/vT)./(psi .+vz/vT .+1).^2 dvdz=1.0 ./dzdv gz=dzdv.*(dT*psi -vz)/tau return z,gz,vz,dvdz end