Nov 9, 2015
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
Major cleanup and moderinization
|
|
|
61 |
/* Declaration of private functions */ |
|
62 |
/* Initialize the RNGs */ |
|
63 |
void set_RNG (void); |
|
64 |
|
|
65 |
/* Firing rates */ |
|
66 |
double get_Qt (int) const; |
|
67 |
double get_Qr (int) const; |
|
68 |
|
|
69 |
/* Synaptic currents */ |
|
70 |
double I_et (int) const; |
|
71 |
double I_gt (int) const; |
|
72 |
double I_er (int) const; |
|
73 |
double I_gr (int) const; |
|
74 |
|
|
75 |
/* Activation functions */ |
|
76 |
double m_inf_T_t (int) const; |
|
77 |
double m_inf_T_r (int) const; |
|
78 |
double m_inf_h (int) const; |
|
79 |
double tau_m_h (int) const; |
|
80 |
double P_h (int) const; |
|
81 |
double act_h (void)const; |
|
82 |
|
|
83 |
double m_inf_hs (int) const; |
|
84 |
double tau_m_hs (int) const; |
|
85 |
double tau_m_hf (int) const; |
|
86 |
|
|
87 |
/* Deactivation functions */ |
|
88 |
double h_inf_T_t (int) const; |
|
89 |
double h_inf_T_r (int) const; |
|
90 |
double tau_h_T_t (int) const; |
|
91 |
double tau_h_T_r (int) const; |
|
92 |
|
|
93 |
/* Currents */ |
|
94 |
double I_L_t (int) const; |
|
95 |
double I_L_r (int) const; |
|
96 |
double I_LK_t (int) const; |
|
97 |
double I_LK_r (int) const; |
|
98 |
double I_T_t (int) const; |
|
99 |
double I_T_r (int) const; |
|
100 |
double I_h (int) const; |
|
101 |
|
|
102 |
/* Noise functions */ |
|
103 |
double noise_xRK (int,int) const; |
|
104 |
double noise_aRK (int) const; |
|
105 |
|
|
106 |
/* Helper functions */ |
|
107 |
inline std::vector<double> init (double value) |
|
108 |
{return {value, 0.0, 0.0, 0.0, 0.0};} |
|
109 |
|
|
110 |
inline void add_RK (std::vector<double>& var) |
|
111 |
{var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6;} |
|
112 |
|
|
113 |
inline void add_RK_noise (std::vector<double>& var, unsigned noise) |
|
114 |
{var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6 + noise_aRK(noise);} |
|
115 |
|
|
116 |
/* Declaration and Initialization of parameters */ |
|
117 |
/* Membrane time in ms */ |
|
118 |
const double tau_t = 20; |
|
119 |
const double tau_r = 20; |
|
120 |
|
|
121 |
/* Maximum firing rate in ms^-1 */ |
|
122 |
const double Qt_max = 400.E-3; |
|
123 |
const double Qr_max = 400.E-3; |
|
124 |
|
|
125 |
/* Sigmoid threshold in mV */ |
|
126 |
const double theta_t = -58.5; |
|
127 |
const double theta_r = -58.5; |
|
128 |
|
|
129 |
/* Sigmoid gain in mV */ |
|
130 |
const double sigma_t = 6.; |
|
131 |
const double sigma_r = 6.; |
|
132 |
|
|
133 |
/* Scaling parameter for sigmoidal mapping (dimensionless) */ |
|
134 |
const double C1 = (M_PI/sqrt(3)); |
|
135 |
|
|
136 |
/* PSP rise time in ms^-1 */ |
|
137 |
const double gamma_e = 70E-3; |
|
138 |
const double gamma_g = 100E-3; |
|
139 |
|
|
140 |
/* Axonal flux time constant in ms^-1*/ |
|
141 |
const double nu = 120E-3; |
|
142 |
|
|
143 |
/* Membrane capacitance in muF/cm^2 */ |
|
144 |
const double C_m = 1.; |
|
145 |
|
|
146 |
/* Leak weight in aU */ |
|
147 |
const double g_L = 1.; |
|
148 |
|
|
149 |
/* Synaptic weights in ms */ |
|
150 |
const double g_AMPA = 1.; |
|
151 |
const double g_GABA = 1.; |
|
152 |
|
|
153 |
/* Conductivities */ |
|
154 |
/* Potassium leak current in mS/m^2 */ |
|
155 |
const double g_LK = 0.02; |
|
156 |
|
|
157 |
/* T current in mS/m^2 */ |
|
158 |
const double g_T_t = 3; |
|
159 |
const double g_T_r = 2.3; |
|
160 |
|
|
161 |
/* h current in mS/m^2 */ |
|
162 |
const double g_h = 0.06; |
|
163 |
|
|
164 |
/* Reversal potentials in mV */ |
|
165 |
/* Synaptic */ |
|
166 |
const double E_AMPA = 0; |
|
167 |
const double E_GABA = -70; |
|
168 |
|
|
169 |
/* Leak */ |
|
170 |
const double E_L_t = -70; |
|
171 |
const double E_L_r = -70; |
|
172 |
|
|
173 |
/* Potassium */ |
|
174 |
const double E_K = -100; |
|
175 |
|
|
176 |
/* I_T current */ |
|
177 |
const double E_Ca = 120; |
|
178 |
|
|
179 |
/* I_h current */ |
|
180 |
const double E_h = -40; |
|
181 |
|
|
182 |
/* Calcium parameters */ |
|
183 |
const double alpha_Ca = -51.8E-6; /* influx per spike in nmol */ |
|
184 |
const double tau_Ca = 10; /* calcium time constant in ms */ |
|
185 |
const double Ca_0 = 2.4E-4; /* resting concentration */ |
|
186 |
|
|
187 |
/* I_h activation parameters */ |
|
188 |
const double k1 = 2.5E7; |
|
189 |
const double k2 = 4E-4; |
|
190 |
const double k3 = 1E-1; |
|
191 |
const double k4 = 1E-3; |
|
192 |
const double n_P = 4; |
|
193 |
const double g_inc = 2; |
|
194 |
|
|
195 |
/* Noise parameters in ms^-1 */ |
|
196 |
const double mphi = 0E-3; |
|
197 |
const double dphi = 20E-3; |
|
198 |
double input = 0.0; |
|
199 |
|
|
200 |
/* Connectivities (dimensionless) */ |
|
201 |
const double N_rt = 3.; |
|
202 |
const double N_tr = 5.; |
|
203 |
const double N_rr = 25.; |
|
204 |
|
|
205 |
/* Connectivities from cortex (dimensionless) */ |
|
206 |
const double N_tp = 2.6; |
|
207 |
const double N_rp = 2.6; |
|
208 |
|
|
209 |
/* Pointer to cortical column */ |
|
210 |
Cortical_Column* Cortex; |
|
211 |
|
|
212 |
/* Parameters for SRK4 iteration */ |
|
213 |
const std::vector<double> A = {0.5, 0.5, 1.0, 1.0}; |
|
214 |
const std::vector<double> B = {0.75, 0.75, 0.0, 0.0}; |
|
215 |
|
|
216 |
/* Random number generators */ |
|
Sep 13, 2016
Major cleanup and moderinization
|
|
|
218 |
|
|
219 |
/* Container for noise */ |
|
220 |
std::vector<double> Rand_vars; |
|
221 |
|
|
222 |
/* Population variables */ |
|
223 |
std::vector<double> Vt = init(E_L_t), /* TC membrane voltage */ |
|
224 |
Vr = init(E_L_r), /* RE membrane voltage */ |
|
225 |
Ca = init(Ca_0), /* Calcium concentration of TC population */ |
|
226 |
s_et = init(0.0), /* PostSP from TC population to TC population */ |
|
227 |
s_er = init(0.0), /* PostSP from TC population to RE population */ |
|
228 |
s_gt = init(0.0), /* PostSP from RE population to TC population */ |
|
229 |
s_gr = init(0.0), /* PostSP from RE population to RE population */ |
|
230 |
y = init(0.0), /* axonal flux */ |
|
231 |
x_et = init(0.0), /* derivative of s_et */ |
|
232 |
x_er = init(0.0), /* derivative of s_er */ |
|
233 |
x_gt = init(0.0), /* derivative of s_gt */ |
|
234 |
x_gr = init(0.0), /* derivative of s_gr */ |
|
235 |
x = init(0.0), /* derivative of y */ |
|
236 |
h_T_t = init(0.0), /* inactivation of T channel */ |
|
237 |
h_T_r = init(0.0), /* inactivation of T channel */ |
|
238 |
m_h = init(0.0), /* activation of h channel */ |
|
239 |
m_h2 = init(0.0); /* activation of h channel bound with protein */ |
|
240 |
|
|
241 |
/* Data storage access */ |