/NM_TC

Permalink
Browse files

Update of RK iteration 

Updated to new scheme after Rößler2010
No more macro shenanigans
Replaced pow(*, N) function with N times product
  • Loading branch information...
1 parent 89fc8c1 commit 2a0c42c0878d21dcb644836252bae3cd751d42c1 @miscco committed Aug 7, 2015
Unified Split
Showing with 405 additions and 533 deletions.
  1. +55 −62 Cortical_Column.cpp
  2. +62 −49 Cortical_Column.h
  3. +1 −17 ODE.h
  4. +13 −21 Plots.m
  5. +16 −17 TC.cpp
  6. +1 −1 TC_model.pro
  7. +83 −93 Test_Stimulation.m
  8. +74 −102 Thalamic_Column.cpp
  9. +88 −76 Thalamic_Column.h
  10. +0 −83 macros.h
  11. +12 −12 saves.h
View
117 Cortical_Column.cpp
@@ -29,16 +29,21 @@
/* Initialization of RNG */
/****************************************************************************************************/
void Cortical_Column::set_RNG(void) {
- /* Number of independent streams */
- int N = 4;
+ extern const double dt;
+ /* Number of independent random variables */
+ int N = 2;
/* Create RNG for each stream */
for (int i=0; i<N; ++i){
- /* Add the RNG */
- MTRands.push_back({ENG(rand()), DIST (mphi, sqrt(dphi))});
+ /* Add the RNG for I_{l}*/
+ MTRands.push_back({ENG(rand()), DIST (0.0, dphi*dt)});
+
+ /* Add the RNG for I_{l,0} */
+ MTRands.push_back({ENG(rand()), DIST (0.0, dt)});
/* Get the random number for the first iteration */
- Rand_vars.push_back(MTRands[i]());
+ Rand_vars.push_back(MTRands[2*i]());
+ Rand_vars.push_back(MTRands[2*i+1]());
}
}
/****************************************************************************************************/
@@ -51,15 +56,13 @@ void Cortical_Column::set_RNG(void) {
/****************************************************************************************************/
/* Pyramidal firing rate */
double Cortical_Column::get_Qe (int N) const{
- _SWITCH((Ve))
- double q = Qe_max / (1 + exp(-C1 * (var_Ve - theta_e) / sigma_e));
+ double q = Qe_max / (1 + exp(-C1 * (Ve[N] - theta_e) / sigma_e));
return q;
}
/* Inhibitory firing rate */
double Cortical_Column::get_Qi (int N) const{
- _SWITCH((Vi))
- double q = Qi_max / (1 + exp(-C1 * (var_Vi - theta_i) / sigma_i));
+ double q = Qi_max / (1 + exp(-C1 * (Vi[N] - theta_i) / sigma_i));
return q;
}
/****************************************************************************************************/
@@ -72,28 +75,24 @@ double Cortical_Column::get_Qi (int N) const{
/****************************************************************************************************/
/* Excitatory input to pyramidal population */
double Cortical_Column::I_ee (int N) const{
- _SWITCH((Ve)(Phi_ee))
- double I = var_Phi_ee * (var_Ve - E_AMPA);
+ double I = y_ee[N] * (Ve[N] - E_AMPA);
return I;
}
/* Inhibitory input to pyramidal population */
double Cortical_Column::I_ie (int N) const{
- _SWITCH((Ve)(Phi_ie))
- double I = var_Phi_ie * (var_Ve - E_GABA);
+ double I = y_ie[N] * (Ve[N] - E_GABA);
return I;
}
/* Excitatory input to inhibitory population */
double Cortical_Column::I_ei (int N) const{
- _SWITCH((Vi)(Phi_ei))
- double I = var_Phi_ei * (var_Vi - E_AMPA);
+ double I = y_ei[N] * (Vi[N] - E_AMPA);
return I;
}
/* Inhibitory input to inhibitory population */
double Cortical_Column::I_ii (int N) const{
- _SWITCH((Vi)(Phi_ii))
- double I = var_Phi_ii * (var_Vi - E_GABA);
+ double I = y_ii[N] * (Vi[N] - E_GABA);
return I;
}
/****************************************************************************************************/
@@ -106,23 +105,20 @@ double Cortical_Column::I_ii (int N) const{
/****************************************************************************************************/
/* Leak current of pyramidal population */
double Cortical_Column::I_L_e (int N) const{
- _SWITCH((Ve))
- double I = g_L * (var_Ve - E_L_e);
+ double I = g_L * (Ve[N] - E_L_e);
return I;
}
/* Leak current of inhibitory population */
double Cortical_Column::I_L_i (int N) const{
- _SWITCH((Vi))
- double I = g_L * (var_Vi - E_L_i);
+ double I = g_L * (Vi[N] - E_L_i);
return I;
}
/* Sodium dependent potassium current */
double Cortical_Column::I_KNa (int N) const{
- _SWITCH((Ve)(Na))
- double w_KNa = 0.37/(1+pow(38.7/var_Na, 3.5));
- double I_KNa = g_KNa * w_KNa * (var_Ve - E_K);
+ double w_KNa = 0.37/(1+pow(38.7/Na[N], 3.5));
+ double I_KNa = g_KNa * w_KNa * (Ve[N] - E_K);
return I_KNa;
}
/****************************************************************************************************/
@@ -134,8 +130,7 @@ double Cortical_Column::I_KNa (int N) const{
/* Potassium pump */
/****************************************************************************************************/
double Cortical_Column::Na_pump (int N) const{
- _SWITCH((Na))
- double Na_pump = R_pump*( pow(var_Na, 3)/(pow(var_Na, 3)+3375) - pow(Na_eq, 3)/(pow(Na_eq, 3)+3375));
+ double Na_pump = R_pump*(Na[N]*Na[N]*Na[N]/(Na[N]*Na[N]*Na[N]+3375) - Na_eq*Na_eq*Na_eq/(Na_eq*Na_eq*Na_eq+3375));
return Na_pump;
}
/****************************************************************************************************/
@@ -147,10 +142,11 @@ double Cortical_Column::Na_pump (int N) const{
/* RK noise scaling */
/****************************************************************************************************/
double Cortical_Column::noise_xRK(int N, int M) const{
- extern const double h;
- extern const vector<double> B1, B2;
- double n = 1 / h * (B1[N-1] * Rand_vars[2*M] + B2[N-1] * Rand_vars[2*M+1]);
- return n;
+ return gamma_e * gamma_e * (Rand_vars[2*M] + Rand_vars[2*M+1]/std::sqrt(3))*B[N];
+}
+
+double Cortical_Column::noise_aRK(int M) const{
+ return gamma_e * gamma_e * (Rand_vars[2*M] - Rand_vars[2*M+1]*std::sqrt(3))/4;
}
/****************************************************************************************************/
/* end */
@@ -162,21 +158,19 @@ double Cortical_Column::noise_xRK(int N, int M) const{
/****************************************************************************************************/
void Cortical_Column::set_RK (int N) {
extern const double dt;
- _SWITCH((Phi_ee)(Phi_ei)(Phi_ie)(Phi_ii)(phi)
- (x_ee) (x_ei) (x_ie) (x_ii) (y))
- Ve [N] = dt*(-(I_L_e(N) + I_ee(N) + I_ie(N))/tau_e - I_KNa(N));
- Vi [N] = dt*(-(I_L_i(N) + I_ei(N) + I_ii(N))/tau_i);
- Na [N] = dt*(alpha_Na * get_Qe(N) - Na_pump(N))/tau_Na;
- Phi_ee [N] = dt*(var_x_ee);
- Phi_ei [N] = dt*(var_x_ei);
- Phi_ie [N] = dt*(var_x_ie);
- Phi_ii [N] = dt*(var_x_ii);
- phi [N] = dt*(var_y);
- x_ee [N] = dt*(pow(gamma_e, 2) * (N_ee * get_Qe(N) + noise_xRK(N, 0) + N_te * Thalamus->get_phi(N) - var_Phi_ee) - 2 * gamma_e * var_x_ee);
- x_ei [N] = dt*(pow(gamma_e, 2) * (N_ei * get_Qe(N) + noise_xRK(N, 1) + N_ti * Thalamus->get_phi(N) - var_Phi_ei) - 2 * gamma_e * var_x_ei);
- x_ie [N] = dt*(pow(gamma_i, 2) * (N_ie * get_Qi(N) - var_Phi_ie) - 2 * gamma_i * var_x_ie);
- x_ii [N] = dt*(pow(gamma_i, 2) * (N_ii * get_Qi(N) - var_Phi_ii) - 2 * gamma_i * var_x_ii);
- y [N] = dt*(pow(nu, 2) * ( get_Qe(N) - var_phi) - 2 * nu * var_y);
+ 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]);
+ 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]);
}
/****************************************************************************************************/
/* end */
@@ -187,24 +181,23 @@ void Cortical_Column::set_RK (int N) {
/* Function that adds all SRK terms */
/****************************************************************************************************/
void Cortical_Column::add_RK(void) {
- extern const double h;
- Ve [0] += (Ve [1] + Ve [2] * 2 + Ve [3] * 2 + Ve [4])/6;
- Vi [0] += (Vi [1] + Vi [2] * 2 + Vi [3] * 2 + Vi [4])/6;
- Na [0] += (Na [1] + Na [2] * 2 + Na [3] * 2 + Na [4])/6;
- Phi_ee [0] += (Phi_ee [1] + Phi_ee[2] * 2 + Phi_ee[3] * 2 + Phi_ee[4])/6;
- Phi_ei [0] += (Phi_ei [1] + Phi_ei[2] * 2 + Phi_ei[3] * 2 + Phi_ei[4])/6;
- Phi_ie [0] += (Phi_ie [1] + Phi_ie[2] * 2 + Phi_ie[3] * 2 + Phi_ie[4])/6;
- Phi_ii [0] += (Phi_ii [1] + Phi_ii[2] * 2 + Phi_ii[3] * 2 + Phi_ii[4])/6;
- phi [0] += (phi [1] + phi [2] * 2 + phi [3] * 2 + phi [4])/6;
- x_ee [0] += (x_ee [1] + x_ee [2] * 2 + x_ee [3] * 2 + x_ee [4])/6 + pow(gamma_e, 2) * h * Rand_vars[0];
- x_ei [0] += (x_ei [1] + x_ei [2] * 2 + x_ei [3] * 2 + x_ei [4])/6 + pow(gamma_e, 2) * h * Rand_vars[2];
- x_ie [0] += (x_ie [1] + x_ie [2] * 2 + x_ie [3] * 2 + x_ie [4])/6;
- x_ii [0] += (x_ii [1] + x_ii [2] * 2 + x_ii [3] * 2 + x_ii [4])/6;
- y [0] += (y [1] + y [2] * 2 + y [3] * 2 + y [4])/6;
-
- /* Generat noise for the next iteration */
+ Ve [0] = (-3*Ve [0] + 2*Ve [1] + 4*Ve [2] + 2*Ve [3] + Ve [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;
+ 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_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 [0] = (-3*x [0] + 2*x [1] + 4*x [2] + 2*x [3] + x [4])/6;
+
+ /* Generate noise for the next iteration */
for (unsigned i=0; i<Rand_vars.size(); ++i) {
- Rand_vars[i] = MTRands[i]();
+ Rand_vars[i] = MTRands[i]() + input;
}
}
/****************************************************************************************************/
View
111 Cortical_Column.h
@@ -30,7 +30,6 @@
#include <boost/random/normal_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include "Thalamic_Column.h"
-#include "macros.h"
using std::vector;
class Thalamic_Column;
@@ -46,6 +45,17 @@ typedef boost::variate_generator<ENG,DIST> GEN; /* Variate generator */
/****************************************************************************************************/
+/* Macro for vector initialization */
+/****************************************************************************************************/
+#ifndef _INIT
+#define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0}
+#endif
+/****************************************************************************************************/
+/* end */
+/****************************************************************************************************/
+
+
+/****************************************************************************************************/
/* Implementation of the cortical module */
/****************************************************************************************************/
class Cortical_Column {
@@ -62,56 +72,55 @@ class Cortical_Column {
/* Connect to the thalamic module */
void get_Thalamus(Thalamic_Column& T) {Thalamus = &T;}
- /* Return axonal flux */
- double get_phi (int N) const {_SWITCH((phi)); return var_phi;}
-
/* Data storage access */
- friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 4));
+ friend void get_data (int, Cortical_Column&, Thalamic_Column&, vector<double*>);
- /* ODE functions */
- void set_RK (int);
- void add_RK (void);
+ /* ODE functions */
+ void set_RK (int);
+ void add_RK (void);
/* Stimulation protocol access */
- friend class Stim;
+ friend class Stim;
+ friend class Thalamic_Column;
private:
- /* Initialize the RNGs */
- void set_RNG (void);
-
- /* Firing rates */
- double get_Qe (int) const;
- double get_Qi (int) const;
-
- /* Currents */
- double I_ee (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_KNa (int) const;
-
- /* Potassium pump */
- double Na_pump (int) const;
-
- /* Noise function */
- double noise_xRK (int, int) const;
-
- /* Population variables */
+ /* Initialize the RNGs */
+ void set_RNG (void);
+
+ /* Firing rates */
+ double get_Qe (int) const;
+ double get_Qi (int) const;
+
+ /* Currents */
+ double I_ee (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_KNa (int) const;
+
+ /* Potassium pump */
+ double Na_pump (int) const;
+
+ /* Noise functions */
+ 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 */
- Phi_ee = _INIT(0.0), /* PostSP from excitatory to excitatory population */
- Phi_ei = _INIT(0.0), /* PostSP from excitatory to inhibitory population */
- Phi_ie = _INIT(0.0), /* PostSP from inhibitory to excitatory population */
- Phi_ii = _INIT(0.0), /* PostSP from inhibitory to inhibitory population */
- phi = _INIT(0.0), /* axonal flux */
- x_ee = _INIT(0.0), /* derivative of Phi_ee */
- x_ei = _INIT(0.0), /* derivative of Phi_ei */
- x_ie = _INIT(0.0), /* derivative of Phi_ie */
- x_ii = _INIT(0.0), /* derivative of Phi_ii */
- y = _INIT(0.0); /* derivative of phi */
+ 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<GEN> MTRands;
@@ -141,9 +150,9 @@ class Cortical_Column {
/* parameters of the firing adaption */
const double alpha_Na = 2; /* Sodium influx per spike in mM ms */
- const double tau_Na = 1.7; /* Sodium time constant in ms */
+ const double tau_Na = 1.7; /* Sodium time constant in ms */
- const double R_pump = 0.09; /* Na-K pump constant in mM/ms */
+ const double R_pump = 0.09; /* Na-K pump constant in mM/ms */
const double Na_eq = 9.5; /* Na-eq concentration in mM */
/* PSP rise time in ms^-1 */
@@ -166,27 +175,31 @@ class Cortical_Column {
const double E_GABA = -70;
/* Leak */
- const double E_L_e = -64;
+ const double E_L_e = -64;
const double E_L_i = -64;
/* Potassium */
const double E_K = -100;
/* Noise parameters in ms^-1 */
const double mphi = 0E-3;
- const double dphi = 60E-3;
+ const double dphi = 20E-1;
double input = 0.0;
/* Connectivities (dimensionless) */
- const double N_ee = 115;
+ const double N_ee = 115;
const double N_ei = 72;
- const double N_ie = 90;
- const double N_ii = 90;
+ const double N_ie = 90;
+ const double N_ii = 90;
const double N_te = 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};
+ const vector<double> B = {0.75, 0.75, 0.0, 0.0};
};
/****************************************************************************************************/
/* end */
Oops, something went wrong.

0 comments on commit 2a0c42c

Please sign in to comment.