Permalink
Newer
Older
100644 332 lines (271 sloc) 11.2 KB
Apr 24, 2014 @miscco Cleanup. Made comments line save
1 /*
Nov 9, 2015 @miscco Updated the code with correct license.
2 * Copyright (c) 2015 University of Lübeck
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 *
22 * AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com
23 *
24 * Based on: A thalamocortical neural mass model of the EEG during NREM sleep and its response
25 * to auditory stimulation.
26 * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz,
27 * JC Claussen.
Sep 13, 2016 @miscco Major cleanup and moderinization
28 * PLoS Computational Biology http://dx.doi.org/10.1371/journal.pcbi.1005022
Nov 9, 2015 @miscco Updated the code with correct license.
29 */
Apr 24, 2014 @miscco Cleanup. Made comments line save
30
Sep 13, 2016 @miscco Major cleanup and moderinization
31 /******************************************************************************/
32 /* Implementation of the stimulation protocol */
33 /******************************************************************************/
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
34 #pragma once
Sep 13, 2016 @miscco Major cleanup and moderinization
35 #include <vector>
36
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
37 #include "Cortical_Column.h"
Sep 13, 2016 @miscco Major cleanup and moderinization
38 #include "Random_Stream.h"
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
39 #include "Thalamic_Column.h"
Sep 13, 2016 @miscco Major cleanup and moderinization
40
41 /******************************************************************************/
42 /* Stimulation object */
43 /******************************************************************************/
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
44 class Stim {
45 public:
Sep 13, 2016 @miscco Major cleanup and moderinization
46 /* Constructor with references and stimulation variables */
47 Stim(Cortical_Column& C, Thalamic_Column& T, double* var)
48 { Cortex = &C;
49 Thalamus = &T;
50 setup(var);}
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
51
Sep 13, 2016 @miscco Major cleanup and moderinization
52 /* Initialize stimulation class with respect to stimulation mode */
53 void setup (double* var_stim);
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
54
Sep 13, 2016 @miscco Major cleanup and moderinization
55 /* Check whether stimulation should be started/stopped */
56 void check_stim (int time);
Oct 1, 2014 @miscco Updated implementation of the Thalamo-cortical model
57 private:
Sep 13, 2016 @miscco Major cleanup and moderinization
58 /* Mode of stimulation */
59 /* 0 == none */
60 /* 1 == semi-periodic */
61 /* 2 == phase dependent */
62 int mode = 0;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
63
Sep 13, 2016 @miscco Major cleanup and moderinization
64 /* Default values already in dt: E1==ms, E4==s */
65 /* Stimulation strength */
66 double strength = 0.0;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
67
Sep 13, 2016 @miscco Major cleanup and moderinization
68 /* Duration of the stimulation */
69 int duration = 120E1;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
70
Sep 13, 2016 @miscco Major cleanup and moderinization
71 /* Interval between different Stimulus events */
72 int ISI = 5E4;
Aug 21, 2014 @miscco Tuned the parameters of both the thalammic and the cortical column.
73
Sep 13, 2016 @miscco Major cleanup and moderinization
74 /* Range of Inter Stimulus Interval */
75 int ISI_range = 1E4;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
76
Sep 13, 2016 @miscco Major cleanup and moderinization
77 /* Number of stimuli in case of multiple stimuli per event */
78 int number_of_stimuli = 1;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
79
Sep 13, 2016 @miscco Major cleanup and moderinization
80 /* Time until next stimulus */
81 /* Function varies between different stimulation modes */
82 int time_to_stimuli = 350E1;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
83
Sep 13, 2016 @miscco Major cleanup and moderinization
84 /* Time between stimuli in case of multiple stimuli per event */
85 int time_between_stimuli = 1050E1;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
86
Sep 13, 2016 @miscco Major cleanup and moderinization
87 /* Threshold for phase dependent stimulation */
88 double threshold = -68;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
89
Sep 13, 2016 @miscco Major cleanup and moderinization
90 /* Internal variables */
91 /* Simulation on for TRUE and off for FALSE */
92 bool stimulation_started = false;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
93
Sep 13, 2016 @miscco Major cleanup and moderinization
94 /* If threshold has been reached */
95 bool threshold_crossed = false;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
96
Sep 13, 2016 @miscco Major cleanup and moderinization
97 /* If minimum was found */
98 bool minimum_found = false;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
99
Sep 13, 2016 @miscco Major cleanup and moderinization
100 /* If a stimulation event has occurred and there is a minimal time (pause) until the next one */
101 bool stimulation_paused = false;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
102
Sep 13, 2016 @miscco Major cleanup and moderinization
103 /* If burst mode is enabled */
104 bool burst_enabled = false;
Jan 9, 2015 @miscco Fixed some of the parameters to get a better SO_Average for N2. However,
105
Sep 13, 2016 @miscco Major cleanup and moderinization
106 /* In case of bursted stimulation start the bursts */
107 bool burst_started = true;
Dec 4, 2014 @miscco Imported the project into QTcreator,
108
Sep 13, 2016 @miscco Major cleanup and moderinization
109 /* Length of a burst stimulus */
110 int burst_length = 10;
Dec 4, 2014 @miscco Imported the project into QTcreator,
111
Sep 13, 2016 @miscco Major cleanup and moderinization
112 /* Length of silence between burst stimuli */
113 int burst_ISI = 10;
Dec 4, 2014 @miscco Imported the project into QTcreator,
114
Sep 13, 2016 @miscco Major cleanup and moderinization
115 /* Onset in time steps where no data is recorded, so stimulation has to be delayed */
116 int onset_correction = 10E4;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
117
Sep 13, 2016 @miscco Major cleanup and moderinization
118 /* Counter for number of stimuli that occurred within a stimulation event */
119 int count_stimuli = 1;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
120
Sep 13, 2016 @miscco Major cleanup and moderinization
121 /* Counter for number of bursts per stimuli */
122 int count_bursts = 0;
Dec 4, 2014 @miscco Imported the project into QTcreator,
123
Sep 13, 2016 @miscco Major cleanup and moderinization
124 /* Counter for stimulation duration since started*/
125 int count_duration = 0;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
126
Sep 13, 2016 @miscco Major cleanup and moderinization
127 /* Counter after minimum was found */
128 int count_to_start = 0;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
129
Sep 13, 2016 @miscco Major cleanup and moderinization
130 /* Counter for time between two stimulation events (with multiple tones) */
131 int count_pause = 0;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
132
Sep 13, 2016 @miscco Major cleanup and moderinization
133 /* Old voltage value for minimum detection */
134 double Vp_old = 0.0;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
135
Sep 13, 2016 @miscco Major cleanup and moderinization
136 /* Pointer to columns */
137 Cortical_Column* Cortex;
138 Thalamic_Column* Thalamus;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
139
Sep 13, 2016 @miscco Major cleanup and moderinization
140 /* Data containers */
141 std::vector<int> marker_stimulation;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
142
Sep 13, 2016 @miscco Major cleanup and moderinization
143 /* Random number generator in case of semi-periodic stimulation */
Sep 26, 2016 @miscco Fixed a minor casting bug in random stream implementation and renamed…
144 randomStreamUniformInt Uniform_Distribution = randomStreamUniformInt(0, 0);
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
145
Sep 13, 2016 @miscco Major cleanup and moderinization
146 /* Create MATLAB container for marker storage */
147 friend mxArray* get_marker(Stim &stim);
148 };
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
149
Sep 13, 2016 @miscco Major cleanup and moderinization
150 /******************************************************************************/
151 /* Function definitions */
152 /******************************************************************************/
Oct 1, 2014 @miscco Updated implementation of the Thalamo-cortical model
153 void Stim::setup (double* var_stim) {
Sep 13, 2016 @miscco Major cleanup and moderinization
154 extern const int onset;
155 extern const int res;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
156
Sep 13, 2016 @miscco Major cleanup and moderinization
157 /* Set the onset onset_correction for the marker */
158 onset_correction = onset * res;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
159
Sep 13, 2016 @miscco Major cleanup and moderinization
160 /* Mode of stimulation */
161 mode = (int) var_stim[0];
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
162
Sep 13, 2016 @miscco Major cleanup and moderinization
163 /* Scale the stimulation strength from s^-1 (Hz) to ms^-1 */
164 strength = var_stim[1] / 1000;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
165
Sep 13, 2016 @miscco Major cleanup and moderinization
166 /* Scale duration from ms to dt */
167 duration = (int) var_stim[2] * res / 1000;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
168
Sep 13, 2016 @miscco Major cleanup and moderinization
169 /* Scale the inter stimulus event interval from s to dt */
170 ISI = (int) var_stim[3] * res;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
171
Sep 13, 2016 @miscco Major cleanup and moderinization
172 /* Scale inter stimulus event interval range from s to dt */
173 ISI_range = (int) var_stim[4] * res;
Mar 1, 2014 @miscco New implementation of the thalamo-cortical model
174
Sep 13, 2016 @miscco Major cleanup and moderinization
175 /* Number of stimuli per Stimulus event */
176 number_of_stimuli = (int) var_stim[5];
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
177
Sep 13, 2016 @miscco Major cleanup and moderinization
178 /* Scale time_between_stimuli from ms to dt */
179 time_between_stimuli = (int) var_stim[6] * res / 1000;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
180
Sep 13, 2016 @miscco Major cleanup and moderinization
181 /* Scale the length of burst_length and burst_ISI from ms to dt*/
182 burst_length = (int) 2 * res / 1000;
183 burst_ISI = (int) 28 * res / 1000;
Dec 4, 2014 @miscco Imported the project into QTcreator,
184
Sep 13, 2016 @miscco Major cleanup and moderinization
185 /* If ISI is fixed do not create RNG */
186 if (mode == 1) {
187 /* Set first time_to_stimuli to 1 sec after onset */
188 time_to_stimuli = (int) (onset+1) * res;
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
189
Sep 13, 2016 @miscco Major cleanup and moderinization
190 /* If ISI is random create RNG */
191 if (ISI_range != 0){
192 /* Generate uniform distribution */
Sep 26, 2016 @miscco Fixed a minor casting bug in random stream implementation and renamed…
193 Uniform_Distribution = randomStreamUniformInt(ISI-ISI_range, ISI+ISI_range);
Sep 13, 2016 @miscco Major cleanup and moderinization
194 }
195 } else {
196 /* In case of phase dependent stimulation, time_to_stim is the time from minimum detection to start of stimulation */
197 /* Scale time_to_stimuli from ms to dt */
198 time_to_stimuli = (int) var_stim[7] * res / 1000;
199 }
Dec 4, 2014 @miscco Imported the project into QTcreator,
200 }
Apr 2, 2014 @miscco Updated the TC model to latest code changes.
201
Oct 1, 2014 @miscco Updated implementation of the Thalamo-cortical model
202 void Stim::check_stim (int time) {
Sep 13, 2016 @miscco Major cleanup and moderinization
203 /* Check if stimulation should start */
204 switch (mode) {
205
206 /* No stimulation */
207 default:
208 break;
209
210 /* Semi-periodic stimulation */
211 case 1:
212 /* Check if stimulation time is reached */
213 if(time == time_to_stimuli) {
214 /* Switch stimulation on */
215 stimulation_started = true;
216 Thalamus->set_input(strength);
217
218 /* Add marker for the first stimuli in the event */
219 if(count_stimuli == 1) {
220 marker_stimulation.push_back(time - onset_correction);
221 }
222
223 /* Check if multiple stimuli should be applied */
224 if (count_stimuli < number_of_stimuli) {
225 /* Update the timer with respect to time between stimuli */
226 time_to_stimuli += time_between_stimuli;
227 count_stimuli++;
228 }
229 /* After last stimulus in event update the timer with respect to (random) ISI*/
230 else {
231 time_to_stimuli += (ISI_range==0)? ISI : Uniform_Distribution();
232
233 /* Reset the stimulus counter for next stimulation event */
234 count_stimuli = 1;
235 }
236 }
237 break;
238
239 /* Phase dependent stimulation */
240 case 2:
241 /* Search for threshold */
242 if(!stimulation_started && !minimum_found && !threshold_crossed && time>onset_correction && !stimulation_paused) {
243 if(Cortex->Vp[0]<=threshold) {
244 threshold_crossed = true;
245 }
246 }
247
248 /* Search for minimum */
249 if(threshold_crossed) {
250 if(Cortex->Vp[0]>Vp_old) {
251 threshold_crossed = false;
252 minimum_found = true;
253 Vp_old = 0;
254 } else {
255 Vp_old = Cortex->Vp[0];
256 }
257 }
258
259 /* Wait until the stimulation should start */
260 if(minimum_found) {
261 /* Start stimulation after time_to_stimuli has passed */
262 if(count_to_start==time_to_stimuli + (count_stimuli-1) * time_between_stimuli) {
263 stimulation_started = true;
264 Thalamus->set_input(strength);
265
266 /* Add marker for the first stimuli in the event */
267 if(count_stimuli == 1) {
268 marker_stimulation.push_back(time - onset_correction);
269 }
270
271 /* Check if multiple stimuli should be applied */
272 if (count_stimuli < number_of_stimuli) {
273 /* Update the number of stimuli */
274 count_stimuli++;
275 } else {
276 /* After last stimulus in event pause the stimulation */
277 minimum_found = false;
278 stimulation_paused = true;
279 count_to_start = 0;
280
281 /* Reset the stimulus counter for next stimulation event */
282 count_stimuli = 1;
283 }
284 }
285 /* Update counter */
286 count_to_start++;
287 }
288 break;
289 }
290
291 /* Actual stimulation protocols */
292 if(stimulation_started) {
293 /* Wait to switch the stimulation off */
294 if(count_duration==duration) {
295 stimulation_started = false;
296 burst_started = true;
297 count_duration = 0;
298 count_bursts = 0;
299 Thalamus->set_input(0.0);
300 }
301
302 count_duration++;
303 count_bursts++;
304
305 /* Switch stimulation on and off wrt burst parameters */
306 if(burst_enabled) {
307 if(burst_started) {
308 if(count_bursts%burst_length==0) {
309 count_bursts = 0;
310 burst_started = false;
311 Thalamus->set_input(0.0);
312 }
313 } else {
314 if(count_bursts%burst_ISI==0) {
315 count_bursts = 0;
316 burst_started = true;
317 Thalamus->set_input(strength);
318 }
319 }
320 }
321 }
322
323 /* Wait if there is a pause between stimulation events */
324 if(stimulation_paused) {
325 if(count_pause == ISI) {
326 stimulation_paused = false;
327 count_pause = 0;
328 }
329 count_pause++;
330 }
Dec 4, 2014 @miscco Imported the project into QTcreator,
331 }