-
Notifications
You must be signed in to change notification settings - Fork 2
/
run_sim.m
199 lines (175 loc) · 10.1 KB
/
run_sim.m
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
function run_sim(make_new_thalamic_input,make_new_thalamic_kernels,make_new_connectivity,includemodulationyn,includeSTDPyn, savename)
%% General settings
f = filesep;
addpath(genpath(['.']))
savefolder = ['Simulation_results' f savename '_data' f];
%% Specific settings thalamic input from Svoboda data
if make_new_thalamic_input
% What whisker data to use; here specific for Svoboda data
SvobodaStruct.loadfolder = ['Input data' f]; % folder where input data are stored
SvobodaStruct.animal = 'an171923'; % animal ID
SvobodaStruct.sessionvec = {'2012_06_04'}; % which sessions to load
SvobodaStruct.dataname = 'data_struct'; % addition on file names
SvobodaStruct.volume = 2; % trials for which recorded volume to use (see explanation Svoboda data)
SvobodaStruct.window.start = 'first touch'; % How to align trials ('first touch', 'first' or 'pole in reach')
SvobodaStruct.window.window = [-2000,4000]; % window (ms) around start time (above)
SvobodaStruct.trialvec = [8,8,9]; % (optional) which of the selected trials to use
savename_input = [savename '_' SvobodaStruct.animal '_' date];
SvobodaStruct.savename = savename;
% Thalamic spike trains (filter neurons responding to whiser data above)
SvobodaStruct.Nkernel_ba = 80; % # base angle kernels ('neurons')
SvobodaStruct.Nkernel_c = 80; % # curvature kernels ('neurons')
SvobodaStruct.Nkernel_m = 40; % # mixed kernels ('neurons')
else
dategen = input('Please give date input to load was generated: ', 's');
animalname = input('Please give animal identifier of input to load: ', 's');
savename_input = [savename '_' animalname '_' dategen];
end
if ~exist(savefolder, 'dir')
mkdir(savefolder)
end
%% Connectivity
if make_new_connectivity
disp('Generating connectivity')
% Make barrel-grid
Nbx = 3; % # of barrels 'x-direction'
Nby = 1; % # of barrels 'y-direction'
barrelstruct = cell(Nbx, Nby);
for nbx = 1:Nbx
for nby = 1:Nby
barrelstruct{nbx, nby}.xpos = (nbx-2)*300; % barrel position
barrelstruct{nbx, nby}.ypos = (nby-1)*300;
barrelstruct{nbx, nby}.Nthalamic = 200; % # filter neurons for this barrel
barrelstruct{nbx, nby}.mainbarrel = 3;
end
end
% Choose main and secondary barrels
barrelstruct{2,1}.mainbarrel = 1; % main
barrelstruct{1,1}.mainbarrel = 2; % secondary
barrelstruct{3,1}.mainbarrel = 2; % tertiary
% Generate connectivity
generate_connectivity(barrelstruct, savefolder, savename)
else
load([savefolder 'CMDMs_' savename], 'barrelstruct');
end
% reorganize the connectivity into a single list in order to speed up simulations.
ConMatfilename = ['CMDMs_' savename]; % connectivity matrix file
if exist([savefolder ConMatfilename '_ConData.mat'],'file') == 2
load([savefolder ConMatfilename '_ConData.mat'])
else
ConData = reorganize_conmat(savename ,savefolder, savefolder, ConMatfilename, includemodulationyn);
end
%% Make/load thalamic spike trains
if make_new_thalamic_input
% Make new spike trains (from Svoboda recordings)
SpikeTrainStruct = make_thalamic_spike_trains_svoboda_recordings(savefolder, savename_input, SvobodaStruct, barrelstruct, make_new_thalamic_kernels);
TSim = SvobodaStruct.window.window(2)-SvobodaStruct.window.window(1);
else
% Load existing spike trains
load([savefolder savename_input '_Thalamic_Spike_Trains']);
load([savefolder savename_input '_Thalamic_Kernels']);
TSim = length(SpikeTrainStruct{1}.PSTH{1,1})*(KernelStruct{1}.kerneltime(2)-KernelStruct{1}.kerneltime(1));
end
% Check
Input_spike_trains = check_reorganize_spike_trains(SpikeTrainStruct, barrelstruct);
clear SpikeTrainStruct KernelStruct SvobodaStruct
%% Simulation parameters
simdata.TSim = TSim; % (ms) Length of the simulation: can be single number (so same length for all trials) or an array of numbers
simdata.timestep = 0.1; % ms
simdata.Trials = 1:2; % This is the amount of times the simulation is run for the same initial conditions/parameters/trace etc
simdata.inputvec = [1,2]; % Which input spike traces (trials) from SpikeTrainStruct to use (optional, if empty all will be used).
if includeSTDPyn
simdata.STDP = true; % Turn on STDP
else
simdata.STDP = false; % Turn off STDP
end
% Parameters and variables to save for each trial:
simdata.whattosave = {'Neuron_Para', 'Tau_plas', 'cellinfo_all', 'cellinfo_input', ...
'modelsc', 'modelspt', 'inputsc', 'inputspikes', 'simLen', 'step', ...
'V0','U0','VT0','vr', 'V', 'seeds', 'timestepseed_input', ...
'timestepseed_model','simdata', 'final_connectivity', 'initial_connectivity'};
% Now left out (too large), but can be added:
% * initdata
% * MIn (whisker modulation for each cell, can easily be reconstructed)
% * U (recovery variable U for each cell)
% * VT (dynamic threshold trace for each cell)
%% Make whisker angles for direct modulation (optional)
if includemodulationyn
simdata.whiskermodulation = true; % Turn on direct whisker modulation
if make_new_thalamic_input
load([savefolder savename_input '_Thalamic_Spike_Trains'], 'WhiskerTrace')
WhPara = WhiskerPara_direct_modulation(WhiskerTrace, barrelstruct, simdata.timestep, savefolder, savename_input);
else
load([savefolder savename_input '_WhiskerModulation']);
end
else
simdata.whiskermodulation = false; % Turn off direct whisker modulation
WhPara = [];
end
clear WhiskerTrace barrelstruct
%% Initial conditions
% Initial parameters
initdata.Vrest = [-65]; % mean resting membrane potential simulated cells; can be list (multiple simulations)
initdata.stdVrest = [0]; % should have same size as Vrest
initdata.V0 = -65; % mean starting membrane potential simulated cells; can be list (multiple simulations)
initdata.stdV0 = 5; % should have same size as V0
initdata.u0 = 0; % mean starting u simulated cells; can be list (multiple simulations)
initdata.stdu0 = 0; % should have same size as u0
% Thresholds
initdata.Vthresdyn = 1; % whether to use a dynamic threshold (1) or not (0) for excitatory neurons
initdata.setVthres.type = 'pertype'; % how to set the threshold: 'distribution', 'pertype' or 'individual'
if strcmp(initdata.setVthres.type, 'distribution')
initdata.Vthres = -55; % mean threshold potential simulated cells; ; can be list (multiple simulations)
initdata.stdVthres = 0; % should have same size as Vthres
initdata.setVthres.nsim = length(initdata.Vthres);
elseif strcmp(initdata.setVthres.type, 'pertype')
% NB mean and std per celltype are set in ConData.Neuron_Vt; this is used
initdata.setVthres.nsim = 1; % multiple realizations can be done with Trials
initdata.Vthresvar = 0; % whether to use an initial distribution around the mean (1) or not (0)
elseif strcmp(initdata.setVthres.type, 'individual')
% set inidividual values, like shown in next section
end
% total # simulations
if isempty(simdata.inputvec)
inputvec = 1:length(Input_spike_trains);
else
inputvec = simdata.inputvec;
end
simdata.Nsim =length(inputvec)*length(initdata.Vrest)*initdata.setVthres.nsim*length(initdata.V0)*length(initdata.u0)*length(simdata.Trials);
% NB simulations are done in the order of for loops above (so Trials is innermost loop)
%% Set initial conditions for individual cells (optional)
% If you want to set parameters / initial conditions for cells specific, set this in
% setindcell, in a matrix of (Number of neurons NAll x total # simulations)
% Here for example: Vrest L4 is more depolarized
% Note: this overwrites anything that would be set with the mean and std of initdata (previous section)!
% NsimperVr = length(initdata.setVthres.nsim)*length(initdata.V0)*length(initdata.u0)*length(simdata.Trials);
% initdata.setindcell.Vrest = zeros(ConData.NAll, Nsim);
% % set all Vrest to initdata.Vrest
% for tt = 1:length(initdata.Vrest)
% initdata.setindcell.Vrest(:, (tt-1)*NsimperVr+1:tt*NsimperVr) = initdata.Vrest(tt)*ones(ConData.NAll,NsimperVr);
% end
% % make L4 cells more depolarized
% initdata.setindcell.Vrest(ConData.Cellinfo_All(:,4)<3,:) = initdata.setindcell.Vrest(ConData.Cellinfo_All(:,4)<3,:)+4.1*ones(sum(ConData.Cellinfo_All(:,4)<3),Nsim);
% initdata.setindcell.V0 = initdata.setindcell.Vrest;
%% External inputs (for instance random, optional)
initdata.ExternalInput = zeros(simdata.Nsim, ConData.NAll,simdata.TSim/simdata.timestep); % (optional) random inputs to each cell
%% Seeding (optional)
% NB should be as long as nr of trials, or non-existent
seeds.initseed = [3,3,3,3,3,3,2,2]; % for seeding initial values (if stdV>0 or stdu0>0)
seeds.runseed = [5,5,5,5,2,2,2,2]; % for seeding synaptic failures and amplitudes
%% Run simulation
run_trials(ConData, initdata, simdata, Input_spike_trains, seeds, WhPara)
%% Make 'luminescence' data: convolve and downsample, see Vogelstein et al 2009
disp('Calculating calcium data')
Para.frame_rate_c = 7/1000; % : sampling rate calcium calculations (kHz)
Para.tau_c = 500; % : calcium decay time constant (ms)
Para.Ca_b = 0.1; % : baseline calcium concentration (microM)
Para.A_c = 5; % : calcium concentration 'jump' for each spike
Para.sigma_c = 1; % : std noise for calcium signal
Para.alpha = 1; % : scale fluorescence signal
Para.beta =0; % : offset fluorescence signal
Para.sigma_f =1; % : std noise for luminescence
Para.tmax = TSim;
simulation_to_calcium( Para, savefolder,savename, initdata, simdata.Nsim);
%% Plot
plot_simulation(simdata.Nsim, initdata, ConData, savename_input)