Theta-Sequence-Network-Model/Simulate_Network.py at master · MattNolanLab/Theta-Sequence-Network-Model · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def Simulate_Network( Ni, Npc, v, PFcent, L, IntNoise):

	from brian import *

	# Note: Npc must be an integer multiple of Ni

	# Parameters of integrate and fire neuron

	tauI  = 40 * msecond      # membrane time constant, needs to be quite long to get theta spiking from dc input
	tauPC = 20 * msecond      # very variable estimates in the literature...
	Vt    = -50 * mvolt          # spike threshold
	Vr    = -70 * mvolt          # reset value (includes AHP)
	El    = -65 * mvolt          # resting potential (same as the reset)
	CmPC  = 155*pfarad          # place cell capacitance
	CmI   = 200*pfarad         # interneuron capacitance

	# Synapse time constants

	taue = 2 * msecond
	taui = 10 * msecond

	# Synapse reversal potentials

	Ee = 0 * mV
	Ei = -70 * mV

	# Other parameters

	f_th = 8*Hz   # septal theta frequency

	sigma = 40 * cm   # place field gaussian width (subthreshold field is larger than spiking field)

	Twidth = sigma / v  # place field width in time

	Trun = L/v # time to run along track


	# Place cell equations

	Idc_extPC =  (1.1  + 0.005 * v / (cm / second)) * (10 ** -4)  * uamp  # input current to E cells, Gaussian amplitude
	sigma_nE = (1.75  - 0.025* (v / (cm/second))) * mV


	eqsPC=('''
	    dVPC/dt = -(VPC - El)/tauPC - gi*(VPC - Ei)/CmPC + Idc/CmPC + sigma_nE * xi / tauPC**.5: volt
	    dgi/dt = -gi/taui        : uS
	    Idc = Idc_extPC* PF * second : amp

	    PF = exp(-(t-Tcent)**2/(2*Twidth**2)) * Hz : Hz
	    Tcent : second

	    ''')


	# Interneuron equations

	Idc_MS =  0.065 * v / (cm/second) * (10 ** -6) * uamp
	Idc0I  =  (7.93 + 0.003 * v / (cm/second)) * (10 ** -5) * uamp


	eqsI=Equations('''
			dVI/dt = -(VI - El )/tauI - ge*(VI - Ee)/CmI + Idc/CmI + sigma_n * xi / tauI**.5: volt
	      	        dge/dt = -ge/taue        : uS
			Idc = Idc0I - Idc_MS * cos(2*pi*f_th * t)  : amp

			sigma_n = IntNoise * volt : volt

		       ''')


	# Define cell groups


	PC = NeuronGroup(N=Npc, model = eqsPC,
		  	    threshold=Vt, reset=Vr )


	for i in range(Npc):

		PC.Tcent[i] = float(PFcent[i]) / v


	I = NeuronGroup(N=Ni, model = eqsI,
		    	  threshold=Vt, reset=Vr )


	# Define connections between cell groups

	we = 0.0005*uS
	wi = 0.025*uS  # set this so as to get realistic membrane oscillations outside place field


	CE = Connection(PC,I,'ge')
	CI = Connection(I,PC,'gi')

	# define E-I connections

	for j in range(Ni):

		CE[:,j] = 0*we/(Npc*10) # all to all component

		CE[(Npc/Ni*j):(Npc/Ni*(j+1)), j] = we # phase precessing component


	# define I-E connections

	for j in range(Ni):

		 CI[j, :] = 0*wi/(Ni*5) # all to all component

	   	 CI[j, (Npc/Ni*j):(Npc/Ni*(j+1))] = wi  # phase precessing component

	# Monitor state variables

	HPC = StateMonitor(PC, 'VPC', record=True)
	HI = StateMonitor(I, 'VI', record=True)
	SPC = SpikeCounter(PC)
	SI = SpikeCounter(I)
	MI = SpikeMonitor(I)
	MPC = SpikeMonitor(PC)


	# Simulate


	PC.VPC = Vr # intial membrane is at resting
	I.VI = Vr   # intial membrane is at resting

	run(Trun)


	return (HPC, HI, SPC, SI, MI, MPC)