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