'''
PYR network model
@author: Ferguson et al. (2015) J. Comput. Neurosci.
Figure 4, Network Size 1K
'''
import brian_no_units
from brian import *
defaultclock.dt = 0.02*ms
seed(3)
C = 115 * pF
vr = -61.8 * mV
vpeak = 22.6 * mV
c = -65.8 * mV
N = 1000
klow = 0.1 * nS/mV
khigh = 3.3 * nS/mV
a = 0.0012 /ms
d = 10 * pA
vt = -57.0 * mV
b = 3 * nS
Ishift= 0 * pA
g = 1.425 * nS
alpha = 2000 /second
Beta = 333.33 /second
s_inf =alpha/(alpha+Beta)
tau_s =1/(alpha+Beta)
Ee=-15 *mV
Iapp = 80*pA
Iapp_std = 15*pA
eqs_pyr = """
Iext : amp
Isyn = g*(v-Ee)*s_sum : amp
s_sum : 1
k=(v<vt)*klow+(v>=vt)*khigh : (siemens/volt)
du/dt = a*(b*(v-vr)-u) : amp
dv/dt = (k*(v-vr)*(v-vt)+Ishift+Iext-Isyn-u)/C : volt
"""
PYR = NeuronGroup(N, model=eqs_pyr, reset ="v = c; u += d" , threshold="v>=vpeak")
PYR.Iext = randn(len(PYR))*Iapp_std + Iapp
S=Synapses(PYR,PYR,
model="""ds0/dt=-(s0-s_inf)/(tau_s) :1
ds1/dt=-s1*Beta :1
s_tot=clip(s0,0,s1) :1""",
pre="""s0=s1; s1=(s_inf*(1-exp(-0.001*second/tau_s))+s1*exp(-0.001*second/tau_s))*exp(0.001*Beta*second)""")
PYR.s_sum=S.s_tot
S.s0=0
S.s1=0
PYR.v = rand(len(PYR))*0.01 -0.065
S[:,:]='(rand()<0.01)'
PYR_v = StateMonitor(PYR,'v',record=[0])
PYR_spktimes = SpikeMonitor(PYR, record=True)
duration = 2 * second
net =Network(PYR,S,PYR_v,PYR_spktimes)
net.run(duration)
figure(1)
plot(PYR_v.times/ms,PYR_v[0]/mV)
xlabel("Time (ms)")
ylabel("PYR Membrane Potential (mV)")
title('PYR network with gpyr=%.3f nS, Iapp=%d pA, IappSD=%d pA'%(g/nS,Iapp/pA,Iapp_std/pA))
figure(2)
raster_plot(PYR_spktimes)
xlabel("Time (ms)")
ylabel("PYR Cell Number")
title('gpyr=%.3f nS, Iapp=%d pA, IappSD=%d pA'%(g/nS,Iapp/pA,Iapp_std/pA))
show()