/NM_TC

Permalink
Browse files

Switched to the random_stream header for RNGs and updated the naming … 

…convention.
  • Loading branch information...
1 parent 1cba4bc commit 60a768b1e3a4401750b60d31dbddcda8eedc3274 @miscco committed Feb 1, 2016
Unified Split
Showing with 189 additions and 170 deletions.
  1. +39 −39 Cortical_Column.cpp
  2. +51 −46 Cortical_Column.h
  3. +1 −1 Data_Storage.h
  4. +8 −8 NM_TC.pro
  5. +1 −1 ODE.h
  6. +0 −1 Random_Stream.h
  7. +7 −5 Stimulation.h
  8. +3 −0 TC.cpp
  9. +1 −0 TC_mex.cpp
  10. +31 −29 Thalamic_Column.cpp
  11. +47 −40 Thalamic_Column.h
View
78 Cortical_Column.cpp
@@ -1,5 +1,5 @@
/*
- * Copyright (c) 2015 University of Lübeck
+ * Copyright (c) 2015 UniVprsity of Lübeck
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
@@ -8,12 +8,12 @@
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
- * The above copyright notice and this permission notice shall be included in
+ * The aboVp copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
- * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVpNT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
@@ -44,10 +44,10 @@ void Cortical_Column::set_RNG(void) {
/* Create RNG for each stream */
for (int i=0; i<N; ++i){
/* Add the RNG for I_{l}*/
- MTRands.push_back(random_stream_normal(0.0, dphi*dt));
+ MTRands.push_back(random_stream_normal(0.0, dphi*dt));
/* Add the RNG for I_{l,0} */
- MTRands.push_back(random_stream_normal(0.0, dt));
+ MTRands.push_back(random_stream_normal(0.0, dt));
/* Get the random number for the first iteration */
Rand_vars.push_back(MTRands[2*i]());
@@ -63,8 +63,8 @@ void Cortical_Column::set_RNG(void) {
/* Firing Rate functions */
/****************************************************************************************************/
/* Pyramidal firing rate */
-double Cortical_Column::get_Qe (int N) const{
- double q = Qe_max / (1 + exp(-C1 * (Ve[N] - theta_e) / sigma_e));
+double Cortical_Column::get_Qp (int N) const{
+ double q = Qp_max / (1 + exp(-C1 * (Vp[N] - theta_p) / sigma_p));
return q;
}
@@ -82,25 +82,25 @@ double Cortical_Column::get_Qi (int N) const{
/* Synaptic currents */
/****************************************************************************************************/
/* Excitatory input to pyramidal population */
-double Cortical_Column::I_ee (int N) const{
- double I = y_ee[N] * (Ve[N] - E_AMPA);
+double Cortical_Column::I_ep (int N) const{
+ double I = g_AMPA * s_ep[N] * (Vp[N] - E_AMPA);
return I;
}
/* Inhibitory input to pyramidal population */
-double Cortical_Column::I_ie (int N) const{
- double I = y_ie[N] * (Ve[N] - E_GABA);
+double Cortical_Column::I_gp (int N) const{
+ double I = g_GABA * s_gp[N] * (Vp[N] - E_GABA);
return I;
}
/* Excitatory input to inhibitory population */
double Cortical_Column::I_ei (int N) const{
- double I = y_ei[N] * (Vi[N] - E_AMPA);
+ double I = g_AMPA * s_ei[N] * (Vi[N] - E_AMPA);
return I;
}
/* Inhibitory input to inhibitory population */
-double Cortical_Column::I_ii (int N) const{
- double I = y_ii[N] * (Vi[N] - E_GABA);
+double Cortical_Column::I_gi (int N) const{
+ double I = g_GABA * s_gi[N] * (Vi[N] - E_GABA);
return I;
}
/****************************************************************************************************/
@@ -112,21 +112,21 @@ double Cortical_Column::I_ii (int N) const{
/* Current functions */
/****************************************************************************************************/
/* Leak current of pyramidal population */
-double Cortical_Column::I_L_e (int N) const{
- double I = g_L * (Ve[N] - E_L_e);
+double Cortical_Column::I_L_p (int N) const{
+ double I = g_L * (Vp[N] - E_L_p);
return I;
}
/* Leak current of inhibitory population */
-double Cortical_Column::I_L_i (int N) const{
- double I = g_L * (Vi[N] - E_L_i);
+double Cortical_Column::I_L_g (int N) const{
+ double I = g_L * (Vi[N] - E_L_g);
return I;
}
/* Sodium dependent potassium current */
double Cortical_Column::I_KNa (int N) const{
double w_KNa = 0.37/(1+pow(38.7/Na[N], 3.5));
- double I_KNa = g_KNa * w_KNa * (Ve[N] - E_K);
+ double I_KNa = g_KNa * w_KNa * (Vp[N] - E_K);
return I_KNa;
}
/****************************************************************************************************/
@@ -166,19 +166,19 @@ double Cortical_Column::noise_aRK(int M) const{
/****************************************************************************************************/
void Cortical_Column::set_RK (int N) {
extern const double dt;
- Ve [N+1] = Ve [0] + A[N] * dt*(-(I_L_e(N) + I_ee(N) + I_ie(N))/tau_e - I_KNa(N));
- Vi [N+1] = Vi [0] + A[N] * dt*(-(I_L_i(N) + I_ei(N) + I_ii(N))/tau_i);
- Na [N+1] = Na [0] + A[N] * dt*(alpha_Na * get_Qe(N) - Na_pump(N))/tau_Na;
- y_ee[N+1] = y_ee[0] + A[N] * dt*(x_ee[N]);
- y_ei[N+1] = y_ei[0] + A[N] * dt*(x_ei[N]);
- y_ie[N+1] = y_ie[0] + A[N] * dt*(x_ie[N]);
- y_ii[N+1] = y_ii[0] + A[N] * dt*(x_ii[N]);
+ Vp [N+1] = Vp [0] + A[N] * dt*(-(I_L_p(N) + I_ep(N) + I_gp(N))/tau_p - I_KNa(N));
+ Vi [N+1] = Vi [0] + A[N] * dt*(-(I_L_g(N) + I_ei(N) + I_gi(N))/tau_i);
+ Na [N+1] = Na [0] + A[N] * dt*(alpha_Na * get_Qp(N) - Na_pump(N))/tau_Na;
+ s_ep[N+1] = s_ep[0] + A[N] * dt*(x_ep[N]);
+ s_ei[N+1] = s_ei[0] + A[N] * dt*(x_ei[N]);
+ s_gp[N+1] = s_gp[0] + A[N] * dt*(x_gp[N]);
+ s_gi[N+1] = s_gi[0] + A[N] * dt*(x_gi[N]);
y [N+1] = y [0] + A[N] * dt*(x [N]);
- x_ee[N+1] = x_ee[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ee * get_Qe(N) + N_te * Thalamus->y[N] - y_ee[N]) - 2 * gamma_e * x_ee[N]) + noise_xRK(N, 0);
- x_ei[N+1] = x_ei[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ei * get_Qe(N) + N_ti * Thalamus->y[N] - y_ei[N]) - 2 * gamma_e * x_ei[N]) + noise_xRK(N, 1) ;
- x_ie[N+1] = x_ie[0] + A[N] * dt*(pow(gamma_i, 2) * (N_ie * get_Qi(N) - y_ie[N]) - 2 * gamma_i * x_ie[N]);
- x_ii[N+1] = x_ii[0] + A[N] * dt*(pow(gamma_i, 2) * (N_ii * get_Qi(N) - y_ii[N]) - 2 * gamma_i * x_ii[N]);
- x [N+1] = x [0] + A[N] * dt*(pow(nu, 2) * ( get_Qe(N) - y [N]) - 2 * nu * x [N]);
+ x_ep[N+1] = x_ep[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ep * get_Qp(N) + N_tp * Thalamus->y[N] - s_ep[N]) - 2 * gamma_e * x_ep[N]) + noise_xRK(N, 0);
+ x_ei[N+1] = x_ei[0] + A[N] * dt*(pow(gamma_e, 2) * (N_ei * get_Qp(N) + N_ti * Thalamus->y[N] - s_ei[N]) - 2 * gamma_e * x_ei[N]) + noise_xRK(N, 1) ;
+ x_gp[N+1] = x_gp[0] + A[N] * dt*(pow(gamma_g, 2) * (N_ip * get_Qi(N) - s_gp[N]) - 2 * gamma_g * x_gp[N]);
+ x_gi[N+1] = x_gi[0] + A[N] * dt*(pow(gamma_g, 2) * (N_ii * get_Qi(N) - s_gi[N]) - 2 * gamma_g * x_gi[N]);
+ x [N+1] = x [0] + A[N] * dt*(pow(nu, 2) * ( get_Qp(N) - y [N]) - 2 * nu * x [N]);
}
/****************************************************************************************************/
/* end */
@@ -189,18 +189,18 @@ void Cortical_Column::set_RK (int N) {
/* Function that adds all SRK terms */
/****************************************************************************************************/
void Cortical_Column::add_RK(void) {
- Ve [0] = (-3*Ve [0] + 2*Ve [1] + 4*Ve [2] + 2*Ve [3] + Ve [4])/6;
+ Vp [0] = (-3*Vp [0] + 2*Vp [1] + 4*Vp [2] + 2*Vp [3] + Vp [4])/6;
Vi [0] = (-3*Vi [0] + 2*Vi [1] + 4*Vi [2] + 2*Vi [3] + Vi [4])/6;
Na [0] = (-3*Na [0] + 2*Na [1] + 4*Na [2] + 2*Na [3] + Na [4])/6;
- y_ee[0] = (-3*y_ee[0] + 2*y_ee[1] + 4*y_ee[2] + 2*y_ee[3] + y_ee[4])/6;
- y_ei[0] = (-3*y_ei[0] + 2*y_ei[1] + 4*y_ei[2] + 2*y_ei[3] + y_ei[4])/6;
- y_ie[0] = (-3*y_ie[0] + 2*y_ie[1] + 4*y_ie[2] + 2*y_ie[3] + y_ie[4])/6;
- y_ii[0] = (-3*y_ii[0] + 2*y_ii[1] + 4*y_ii[2] + 2*y_ii[3] + y_ii[4])/6;
+ s_ep[0] = (-3*s_ep[0] + 2*s_ep[1] + 4*s_ep[2] + 2*s_ep[3] + s_ep[4])/6;
+ s_ei[0] = (-3*s_ei[0] + 2*s_ei[1] + 4*s_ei[2] + 2*s_ei[3] + s_ei[4])/6;
+ s_gp[0] = (-3*s_gp[0] + 2*s_gp[1] + 4*s_gp[2] + 2*s_gp[3] + s_gp[4])/6;
+ s_gi[0] = (-3*s_gi[0] + 2*s_gi[1] + 4*s_gi[2] + 2*s_gi[3] + s_gi[4])/6;
y [0] = (-3*y [0] + 2*y [1] + 4*y [2] + 2*y [3] + y [4])/6;
- x_ee[0] = (-3*x_ee[0] + 2*x_ee[1] + 4*x_ee[2] + 2*x_ee[3] + x_ee[4])/6 + noise_aRK(0);
+ x_ep[0] = (-3*x_ep[0] + 2*x_ep[1] + 4*x_ep[2] + 2*x_ep[3] + x_ep[4])/6 + noise_aRK(0);
x_ei[0] = (-3*x_ei[0] + 2*x_ei[1] + 4*x_ei[2] + 2*x_ei[3] + x_ei[4])/6 + noise_aRK(1);
- x_ie[0] = (-3*x_ie[0] + 2*x_ie[1] + 4*x_ie[2] + 2*x_ie[3] + x_ie[4])/6;
- x_ii[0] = (-3*x_ii[0] + 2*x_ii[1] + 4*x_ii[2] + 2*x_ii[3] + x_ii[4])/6;
+ x_gp[0] = (-3*x_gp[0] + 2*x_gp[1] + 4*x_gp[2] + 2*x_gp[3] + x_gp[4])/6;
+ x_gi[0] = (-3*x_gi[0] + 2*x_gi[1] + 4*x_gi[2] + 2*x_gi[3] + x_gi[4])/6;
x [0] = (-3*x [0] + 2*x [1] + 4*x [2] + 2*x [3] + x [4])/6;
/* Generate noise for the next iteration */
View
97 Cortical_Column.h
@@ -60,8 +60,8 @@ class Cortical_Column {
{set_RNG();}
Cortical_Column(double* Param, double* Con)
- :sigma_e (Param[0]), g_KNa (Param[1]), dphi (Param[2]),
- N_te (Con[2]), N_ti (Con[3])
+ :sigma_p (Param[0]), g_KNa (Param[1]), dphi (Param[2]),
+ N_tp (Con[2]), N_ti (Con[3])
{set_RNG();}
/* Connect to the thalamic module */
@@ -79,20 +79,21 @@ class Cortical_Column {
friend class Thalamic_Column;
private:
+ /* Declaration of private functions */
/* Initialize the RNGs */
void set_RNG (void);
/* Firing rates */
- double get_Qe (int) const;
+ double get_Qp (int) const;
double get_Qi (int) const;
/* Currents */
- double I_ee (int) const;
+ double I_ep (int) const;
double I_ei (int) const;
- double I_ie (int) const;
- double I_ii (int) const;
- double I_L_e (int) const;
- double I_L_i (int) const;
+ double I_gp (int) const;
+ double I_gi (int) const;
+ double I_L_p (int) const;
+ double I_L_g (int) const;
double I_KNa (int) const;
/* Potassium pump */
@@ -102,46 +103,25 @@ class Cortical_Column {
double noise_xRK (int,int) const;
double noise_aRK (int) const;
- /* Population variables */
- vector<double> Ve = _INIT(E_L_e), /* excitatory membrane voltage */
- Vi = _INIT(E_L_i), /* inhibitory membrane voltage */
- Na = _INIT(Na_eq), /* Na concentration */
- y_ee = _INIT(0.0), /* PostSP from excitatory to excitatory population */
- y_ei = _INIT(0.0), /* PostSP from excitatory to inhibitory population */
- y_ie = _INIT(0.0), /* PostSP from inhibitory to excitatory population */
- y_ii = _INIT(0.0), /* PostSP from inhibitory to inhibitory population */
- y = _INIT(0.0), /* axonal flux */
- x_ee = _INIT(0.0), /* derivative of y_ee */
- x_ei = _INIT(0.0), /* derivative of y_ei */
- x_ie = _INIT(0.0), /* derivative of y_ie */
- x_ii = _INIT(0.0), /* derivative of y_ii */
- x = _INIT(0.0); /* derivative of y */
-
- /* Random number generators */
- vector<random_stream_normal> MTRands;
-
- /* Container for noise */
- vector<double> Rand_vars;
-
/* Declaration and Initialization of parameters */
/* Membrane time in ms */
- const double tau_e = 30;
+ const double tau_p = 30;
const double tau_i = 30;
/* Maximum firing rate in ms^-1 */
- const double Qe_max = 30.E-3;
+ const double Qp_max = 30.E-3;
const double Qi_max = 60.E-3;
/* Sigmoid threshold in mV */
- const double theta_e = -58.5;
+ const double theta_p = -58.5;
const double theta_i = -58.5;
/* Sigmoid gain in mV */
- const double sigma_e = 4;
+ const double sigma_p = 4;
const double sigma_i = 6;
/* Scaling parameter for sigmoidal mapping (dimensionless) */
- const double C1 = (3.14159265/sqrt(3));
+ const double C1 = (M_PI/sqrt(3));
/* parameters of the firing adaption */
const double alpha_Na = 2; /* Sodium influx per spike in mM ms */
@@ -152,16 +132,20 @@ class Cortical_Column {
/* PSP rise time in ms^-1 */
const double gamma_e = 70E-3;
- const double gamma_i = 58.6E-3;
+ const double gamma_g = 58.6E-3;
/* Axonal flux time constant */
const double nu = 120E-3;
- /* Conductivities in mS/cm^-2 */
- /* Leak */
- const double g_L = 1;
+ /* Conductivities */
+ /* Leak in aU*/
+ const double g_L = 1.;
- /* KNa */
+ /* Synaptic conductivities in ms */
+ const double g_AMPA = 1.;
+ const double g_GABA = 1.;
+
+ /* KNa in mS/cm^2 */
const double g_KNa = 1.33;
/* Reversal potentials in mV */
@@ -170,8 +154,8 @@ class Cortical_Column {
const double E_GABA = -70;
/* Leak */
- const double E_L_e = -64;
- const double E_L_i = -64;
+ const double E_L_p = -64;
+ const double E_L_g = -64;
/* Potassium */
const double E_K = -100;
@@ -182,19 +166,40 @@ class Cortical_Column {
double input = 0.0;
/* Connectivities (dimensionless) */
- const double N_ee = 115;
+ const double N_ep = 115;
const double N_ei = 72;
- const double N_ie = 90;
+ const double N_ip = 90;
const double N_ii = 90;
- const double N_te = 2.5;
+ const double N_tp = 2.5;
const double N_ti = 2.5;
/* Pointer to thalamic column */
Thalamic_Column* Thalamus;
- /* Interation matrix for SRK4 */
- const vector<double> A = {0.5, 0.5, 1.0, 1.0};
+ /* Parameters for SRK4 iteration */
+ const vector<double> A = {0.5, 0.5, 1.0, 1.0};
const vector<double> B = {0.75, 0.75, 0.0, 0.0};
+
+ /* Random number generators */
+ vector<random_stream_normal> MTRands;
+
+ /* Container for noise */
+ vector<double> Rand_vars;
+
+ /* Population variables */
+ vector<double> Vp = _INIT(E_L_p), /* excitatory membrane voltage */
+ Vi = _INIT(E_L_g), /* inhibitory membrane voltage */
+ Na = _INIT(Na_eq), /* Na concentration */
+ s_ep = _INIT(0.0), /* PostSP from excitatory to excitatory population */
+ s_ei = _INIT(0.0), /* PostSP from excitatory to inhibitory population */
+ s_gp = _INIT(0.0), /* PostSP from inhibitory to excitatory population */
+ s_gi = _INIT(0.0), /* PostSP from inhibitory to inhibitory population */
+ y = _INIT(0.0), /* axonal flux */
+ x_ep = _INIT(0.0), /* derivative of s_ep */
+ x_ei = _INIT(0.0), /* derivative of s_ei */
+ x_gp = _INIT(0.0), /* derivative of s_gp */
+ x_gi = _INIT(0.0), /* derivative of s_gi */
+ x = _INIT(0.0); /* derivative of y */
};
/****************************************************************************************************/
/* end */
View
2 Data_Storage.h
@@ -38,7 +38,7 @@
/****************************************************************************************************/
/* Save data */
/****************************************************************************************************/
-inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, vector<double*> Data) {
+void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, vector<double*> Data) {
Data[0][counter] = Cortex.Ve [0];
Data[1][counter] = Thalamus.Vt [0];
Data[2][counter] = Thalamus.Ca [0];
View
16 NM_TC.pro
@@ -6,16 +6,16 @@ CONFIG -= qt
TARGET = TC.cpp
SOURCES += Cortical_Column.cpp \
- TC.cpp \
- TC_mex.cpp \
+ TC.cpp \
+ TC_mex.cpp \
Thalamic_Column.cpp
-HEADERS += Cortical_Column.h \
- Data_Storage.h \
- ODE.h \
- Random_Stream.h \
- Stimulation.h \
- Thalamic_Column.h
+HEADERS += Cortical_Column.h \
+ Data_Storage.h \
+ ODE.h \
+ Random_Stream.h \
+ Stimulation.h \
+ Thalamic_Column.h
QMAKE_CXXFLAGS += -std=c++11 -O3
View
2 ODE.h
@@ -38,7 +38,7 @@
/****************************************************************************************************/
/* Evaluation of SRK4 */
/****************************************************************************************************/
-inline void ODE(Cortical_Column& Cortex, Thalamic_Column& Thalamus) {
+void ODE(Cortical_Column& Cortex, Thalamic_Column& Thalamus) {
/* First calculate every ith RK moment. Has to be in order, 1th moment first */
for (int i=0; i<4; ++i) {
Cortex.set_RK(i);
View
1 Random_Stream.h
@@ -21,7 +21,6 @@
*
* AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com
* Stefanie Gareis: gareis@inb.uni-luebeck.de
- *
*/
/****************************************************************************************************/
Oops, something went wrong.

0 comments on commit 60a768b

Please sign in to comment.