-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBackgroundStim.py
154 lines (124 loc) · 6.11 KB
/
BackgroundStim.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
import json
import numpy as np
# adapted from: https://github.com/BlueBrain/neurodamus/blob/2c7052096c22fb5183fdef120d080608748da48d/neurodamus/core/stimuli.py#L271
class addStim():
@staticmethod
def add_ornstein_uhlenbeck(tau, sigma, mean, duration, dt=0.025, seed=100000, plotFig=False):
from neuron import h
import numpy as np
"""
Adds an Ornstein-Uhlenbeck process with given correlation time,
standard deviation and mean value.
tau: correlation time [ms], white noise if zero
sigma: standard deviation [uS]
mean: mean value [uS]
duration: duration of signal [ms]
dt: timestep [ms]
"""
from math import sqrt, exp
"""Creates a default RNG, currently based on ACG"""
rng = h.Random(seed)
# rng = RNG() # Creates a default RNG
# if not self._rng:
# logging.warning("Using a default RNG for Ornstein-Uhlenbeck process")
# tvec = h.Vector()
# tvec.indgen(self._cur_t, self._cur_t + duration, dt) # time vector
tvec = h.Vector(np.linspace(0, duration, int(duration / dt)))
ntstep = len(tvec) # total number of timesteps
svec = h.Vector(ntstep, 0) # stim vector
noise = h.Vector(ntstep) # Gaussian noise
rng.normal(0.0, 1.0)
noise.setrand(rng) # generate Gaussian noise
if tau < 1e-9:
svec = noise.mul(sigma) # white noise
else:
mu = exp(-dt / tau) # auxiliar factor [unitless]
A = sigma * sqrt(1 - mu * mu) # amplitude [uS]
noise.mul(A) # scale noise by amplitude [uS]
# Exact update formula (independent of dt) from Gillespie 1996
for n in range(1, ntstep):
svec.x[n] = svec[n - 1] * mu + noise[n] # signal [uS]
svec.add(mean) # shift signal by mean value [uS]
svec = abs(svec)
for idx, val in enumerate(svec):
if val < 0:
print('original value is ' + str(val))
svec[idx] = 0
print('value is now ' + str(val))
if plotFig:
import matplotlib.pyplot as plt
plt.figure(figsize=(30, 5))
plt.plot(list(tvec), list(svec), 'k')
plt.savefig('test_fig_vec_OrnsteinUhlenbeck.png')
plt.figure(figsize=(30, 5))
plt.plot(list(tvec)[0:1000], list(svec)[0:1000], 'k')
plt.savefig('test_fig_vec_OrnsteinUhlenbeck_slice.png')
plt.figure()
plt.hist(list(svec), 100)
plt.savefig('test_fig_vec_OrnsteinUhlenbeck_hist.png')
return tvec, svec
def addNoiseGClamp(sim):
# from init import vecs_dict
import numpy as np
print('Creating Ornstein Uhlenbeck process to create noise conductance signal')
import math
vecs_dict = {}
for cell_ind, cell in enumerate(sim.net.cells):
vecs_dict.update({cell_ind: {'tvecs': {}, 'svecs': {}}})
cell_seed = sim.cfg.seeds['stim']+ cell.gid
for stim_ind, stim in enumerate(sim.net.cells[cell_ind].stims):
if 'NoiseSEClamp' in stim['label']:
mean = sim.net.params.NoiseConductanceParams[cell.tags['pop']]['g0']
sigma = sim.net.params.NoiseConductanceParams[cell.tags['pop']]['sigma']
tvec, svec = addStim.add_ornstein_uhlenbeck(
tau=10,
sigma=sigma,
mean=mean,
duration=stim['dur1'],
dt=0.05,
seed=cell_seed,
plotFig=False)
vecs_dict[cell_ind]['tvecs'].update({stim_ind: tvec})
vecs_dict[cell_ind]['svecs'].update({stim_ind: svec})
conductance_source = sim.net.cells[cell_ind].stims[stim_ind]['hObj']
# Play the vector into the rs variable of the ConductanceSource object
vecs_dict[cell_ind]['svecs'][stim_ind].play(conductance_source._ref_rs, vecs_dict[cell_ind]['tvecs'][stim_ind], True)
# vecs_dict[cell_ind]['svecs'][stim_ind].play(
# sim.net.cells[cell_ind].stims[stim_ind]['hObj']._ref_rs, vecs_dict[cell_ind]['tvecs'][stim_ind],
# True
# )
return sim, vecs_dict
def addNoiseIClamp(sim):
import numpy as np
print('\t>> Using Ornstein Uhlenbeck to add noise to the IClamp')
import math
# from CurrentStim import CurrentStim as CS
vecs_dict = {}
for cell_ind, cell in enumerate(sim.net.cells):
vecs_dict.update({cell_ind: {'tvecs': {}, 'svecs': {}}})
cell_seed = sim.cfg.seeds['stim']+ cell.gid
for stim_ind, stim in enumerate(sim.net.cells[cell_ind].stims):
if 'NoiseIClamp' in stim['label']:
# try:
mean = sim.net.params.NoiseIClampParams[cell.tags['pop']]['g0']
variance = sim.net.params.NoiseIClampParams[cell.tags['pop']]['sigma']
# print(cell.tags['pop'] + ' ' + str(variance) + ' ' + str(mean))
# print('mean noise: ', mean, ' nA')
# except:
# mean = 0
# print('except mean noise: ', mean, ' nA')
tvec, svec = BackgroundStim.add_ornstein_uhlenbeck(
tau= 10,
sigma=variance,
mean=mean,
duration=stim['dur'],
dt=0.025,
seed=cell_seed,
plotFig=False)
vecs_dict[cell_ind]['tvecs'].update({stim_ind: tvec})
vecs_dict[cell_ind]['svecs'].update({stim_ind: svec})
vecs_dict[cell_ind]['svecs'][stim_ind].play(
sim.net.cells[cell_ind].stims[stim_ind]['hObj']._ref_amp, vecs_dict[cell_ind]['tvecs'][stim_ind],
True
)
return sim, vecs_dict