-
Notifications
You must be signed in to change notification settings - Fork 3
/
iv_model.py
336 lines (265 loc) · 9.88 KB
/
iv_model.py
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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
"""
Single neuron circuit model
Circuit consists of a parallel interconnection of an arbitrary number of either
'Current' or 'Conductance' elements
@author: Luka
"""
from numpy import tanh, exp
import numpy as np
from scipy.integrate import solve_ivp, BDF
def sigmoid(x, k = 1):
return 1 / (1 + exp(-k * (x)))
class EulerSolver():
"""
ODE solver using the basic Euler step
args:
odesys: f(t,y) -> dy/dt = f(t,y)
t0: initial time
y0: initial state
dt: time step
methods:
step: iterate a single simulation step
"""
def __init__(self, odesys, t0, y0, dt):
self.odesys = odesys
self.t = t0
self.y = y0
self.dt = dt
def step(self):
self.y += self.odesys(self.t,self.y)*self.dt
self.t += self.dt
# Return error message for compatibility with scipy solvers
errorMessage = False
return errorMessage
class System():
"""
Parent class implementing basic simulation methods
methods:
sys: f(i_app, y) -> dV/dt = f(i_app, y)
set_solver: set the ODE solver and simulation parameters
step: iterate a single simulation step and return next (t,y)
"""
def __init__(self):
self.y0 = []
def sys(self, i_app, y):
pass
def set_solver(self, solver, i_app, t0, sstep, dt = 1):
def odesys(t, y):
return self.sys(i_app(t), y)
if (solver == "Euler"):
self.solver = EulerSolver(odesys, t0, self.y0, dt)
elif (solver == "BDF"):
self.solver = BDF(odesys, t0, self.y0, np.inf, max_step = sstep)
else:
raise ValueError("Undefined solver")
def step(self):
msg = self.solver.step()
t = self.solver.t
y = self.solver.y
if msg:
raise ValueError('Solver terminated with message: %s ' % msg)
return t,y
def simulate(self, trange, i_app, method = "Default", max_step = np.inf):
def odesys(t, y):
return self.sys(i_app(t), y)
if (method == "Default"):
sol = solve_ivp(odesys, trange, self.y0, max_step = max_step)
return sol
class SingleTimescaleElement():
"""
Parent class for elements depending on a single filtered voltage Vx
args:
neuron: pointer to the neuron containing the element
timescale: time constant of the first-order filter
v0: initial condition
v_index: index of Vx, assigned when interconnected in a circuit
methods:
out: Iout for input V
outx: Iout using the appropriate Vx (when interconnected)
IV: IV curve of the element in timescale tau
-> return out(V) if tau <= timescale
-> return out(Vrest) if tau > timescale
"""
vx0 = None # Default initial conditions for first-order filters
def __init__(self, neuron, timescale, v0):
self.neuron = neuron
self.timescale = timescale
self.v0 = v0
self.v_index = None
if (timescale == 0) and (v0 is not None):
raise ValueError("Initial condition of an instantaneous element "
"cannot be set")
self._add_timescale(neuron.timescales, neuron.y0)
def _add_timescale(self, timescales, y0):
"""
Add a first-order filter for Vx, or associate the element with an
existing Vx
"""
if (self.timescale in timescales) and (self.v0 is None):
self.v_index = timescales.index(self.timescale)
else:
timescales.append(self.timescale)
self.v_index = len(timescales) - 1
if (self.v0):
y0.append(self.v0)
else:
y0.append(self.vx0)
def out(self, V):
"""
Steady-state output depending on the input V
"""
return V
def outx(self, y):
"""
Output depending on the appropriate Vx
"""
return self.out(y[self.v_index])
def IV(self, V, tau, Vrest = 0):
"""
I-V curve of the element in timescale tau
"""
if (self.timescale <= tau):
return self.out(V)
else:
return self.out(Vrest)
class Neuron(System):
"""
Parallel interconnection of current or conductance elements
C dV/dT = - sum(I_x) + Iapp
where I_x is the output current of each current/conductance element
kwargs:
C: membrane capacitance
v0: initial voltage
vx0: default initial state for current/conductance elements
methods:
add_current: add current element
add_conductance: add conductance element
IV: IV curve in timescale tau and resting voltage Vrest
IV_ss: steady-state IV curve
get_init_conditions: return y0
i_sum: sum(Ix) for all conductance/circuit elements
"""
# Membrane capacitor value + init conditions
_stdPar = {'C': 1, 'v0': -1.9, 'vx0': -1.8}
@property
def stdPar(self):
return type(self)._stdPar
def __init__(self, **kwargs):
self.__dict__.update(self.stdPar) # Default circuit parameters
self.__dict__.update(kwargs) # Modify circuit parameters
SingleTimescaleElement.vx0 = self.vx0 # Default initial conditions
self.timescales = [0] # Timescales of membrane voltage + all filters
self.y0 = [self.v0] # Initial conditions
self.elements = [] # List containing all circuit elements
def add_current(self, a, voff, timescale, v0 = None):
I = self.CurrentElement(self, a, voff, timescale, v0)
self.elements.append(I)
return I
def add_conductance(self, g_max, E_rev = 0):
I = self.ConductanceElement(self, g_max, E_rev)
self.elements.append(I)
return I
def IV(self, V, tau, Vrest = 0):
I = 0
for el in self.elements:
I += el.IV(V, tau, Vrest)
return I
def IV_ss(self, V):
I = 0
for el in self.elements:
I += el.out(V)
return I
def get_init_conditions(self):
return np.array(self.y0)
def i_sum(self, y):
"""
Returns total internal current
"""
s = 0
for el in self.elements:
s += el.outx(y)
return s
def sys(self, i_app, y):
"""
Returns the state vector update
y[0] = membrane voltage
y[1],y[2],... = Element first-order filters, in order of definition
"""
dy = []
dvmem = (i_app - self.i_sum(y)) / self.C
dy.append(dvmem)
# First-order filters
for index, tau in enumerate(self.timescales[1:]):
dy.append((y[0] - y[index+1]) / tau)
return np.array(dy)
class CurrentElement(SingleTimescaleElement):
"""
Current element of the form:
Iout = a * tanh(Vx - voff).
taux * dVx/dt = Vmem - Vx
"""
def __init__(self, neuron, a, voff, timescale, v0):
super().__init__(neuron, timescale, v0)
self.a = a
self.voff = voff
def out(self, V):
return (self.a * tanh(V - self.voff))
def update_a(self, a):
self.a = a
def update_voff(self, voff):
self.voff = voff
class ConductanceElement:
"""
Single conductance element consisting of multiple gates:
Iout = g_max * x1 * x2 * ... * xn * (V - E_rev)
*args: [x1,x2,...,xn] = gates
"""
def __init__(self, neuron, g_max, E_rev = 0):
self.neuron = neuron
self.g_max = g_max
self.E_rev = E_rev
self.gates = []
class Gate(SingleTimescaleElement):
"""
Single gating variable with sigmoidal activation/inactivation:
x = S(k*(Vx - voff))
taux * dVx/dt = Vmem - Vx
"""
def __init__(self, neuron, k, voff, timescale, v0):
super().__init__(neuron, timescale, v0)
self.k = k
self.voff = voff
def out(self, V):
return sigmoid(V - self.voff, self.k)
def update_voff(self, voff):
self.voff = voff
def update_k(self, k):
self.k = k
# Add a gating variable to the conductance element
def add_gate(self, k, voff, timescale, v0 = None):
x = self.Gate(self.neuron, k, voff, timescale, v0)
self.gates.append(x)
return x
def out(self, V):
iout = self.g_max * (V - self.E_rev)
for x in self.gates:
iout *= x.out(V)
return iout
def outx(self, y):
iout = self.g_max * (y[0] - self.E_rev)
for x in self.gates:
iout *= x.outx(y)
return iout
def update_g_max(self, g_max):
self.g_max = g_max
def update_E_rev(self, E_rev):
self.E_rev = E_rev
def IV(self, V, tau, Vrest = 0):
I = self.g_max * (V - self.E_rev)
gates_fast = [g for g in self.gates if g.timescale <= tau]
gates_slow = [g for g in self.gates if g.timescale > tau]
for g in gates_fast:
I *= g.out(V)
for g in gates_slow:
I *= g.out(Vrest)
return I