diff --git a/Cortical_Column.cpp b/Cortical_Column.cpp index 1c6e8c1..f1be3de 100644 --- a/Cortical_Column.cpp +++ b/Cortical_Column.cpp @@ -1,3 +1,25 @@ +/* +* Copyright (c) 2014 Michael Schellenberger Costa +* +* Permission is hereby granted, free of charge, to any person obtaining a copy +* of this software and associated documentation files (the "Software"), to deal +* in the Software without restriction, including without limitation the rights +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +* 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 +* 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 +* 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 +* THE SOFTWARE. +*/ + /****************************************************************************************************/ /* Functions of the cortical module */ /****************************************************************************************************/ @@ -7,15 +29,15 @@ /* Initialization of RNG */ /****************************************************************************************************/ void Cortical_Column::set_RNG(void) { - // the number of independent streams + /* Number of independent streams */ int N = 4; - // get the RNG + /* Create RNG for each stream */ for (int i=0; i #include #include -#include "macros.h" #include "Thalamic_Column.h" +#include "macros.h" using std::vector; - class Thalamic_Column; /****************************************************************************************************/ /* Typedefs for RNG */ /****************************************************************************************************/ -typedef boost::mt11213b ENG; // Mersenne Twister -typedef boost::normal_distribution DIST; // Normal Distribution -typedef boost::variate_generator GEN; // Variate generator +typedef boost::mt11213b ENG; /* Mersenne Twister */ +typedef boost::normal_distribution DIST; /* Normal Distribution */ +typedef boost::variate_generator GEN; /* Variate generator */ /****************************************************************************************************/ /* end */ /****************************************************************************************************/ @@ -29,30 +50,29 @@ typedef boost::variate_generator GEN; // Variate generator /****************************************************************************************************/ class Cortical_Column { public: - // Constructors + /* Constructors */ Cortical_Column(void) {set_RNG();} Cortical_Column(double* Par, double* Con) - :tau_e (Par[0]), theta_e (Par[1]), sigma_e (Par[2]), alpha_Na (Par[3]), - tau_Na (Par[4]), g_KNa (Par[5]), N_te (Con[2]), N_ti (Con[3]), - dphi (Par[6]) + :sigma_e (Par[0]), alpha_Na (Par[1]), tau_Na (Par[2]), dphi (Par[3]), + N_te (Con[2]), N_ti (Con[3]) {set_RNG();} - // get the pointer to the cortical module + /* Connect to the thalamic module */ void get_Thalamus(Thalamic_Column& T) {Thalamus = &T;} - // get axonal flux + /* Return axonal flux */ double get_phi (int N) const {_SWITCH((phi_e)); return var_phi_e;} - // Initialize the RNGs + /* Initialize the RNGs */ void set_RNG (void); - // Firing rates + /* Firing rates */ double get_Qe (int) const; double get_Qi (int) const; - // Currents + /* Currents */ double I_ee (int) const; double I_ei (int) const; double I_ie (int) const; @@ -61,112 +81,112 @@ class Cortical_Column { double I_L_i (int) const; double I_KNa (int) const; - // Potassium pump + /* Potassium pump */ double Na_pump (int) const; - // Noise function + /* Noise function */ double noise_xRK (int, int) const; - // ODE functions + /* ODE functions */ void set_RK (int); void add_RK (void); - // Data storage + /* Data storage access */ friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 2)); + /* Stimulation protocoll access */ + friend class Stim; + private: - // Population variables - vector 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_e = _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_e = _INIT(0.0); // derivative of phi_t - - // Random number generators + /* Population variables */ + vector 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_e = _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_e = _INIT(0.0); /* derivative of phi_t */ + + /* Random number generators */ vector MTRands; - // Container for noise + /* Container for noise */ vector Rand_vars; - // Declaration and Initialization of parameters - // Membrane time in ms + /* Declaration and Initialization of parameters */ + /* Membrane time in ms */ const double tau_e = 30; const double tau_i = 30; - // Maximum firing rate in ms^-1 + /* Maximum firing rate in ms^-1 */ const double Qe_max = 30.E-3; const double Qi_max = 60.E-3; - // Sigmoid threshold in mV + /* Sigmoid threshold in mV */ const double theta_e = -58.5; const double theta_i = -58.5; - // Sigmoid gain in mV + /* Sigmoid gain in mV */ const double sigma_e = 4; const double sigma_i = 6; - // Scaling parameter for sigmoidal mapping (dimensionless) + /* Scaling parameter for sigmoidal mapping (dimensionless) */ const double C1 = (3.14159265/sqrt(3)); - // parameters of the firing adaption - const double alpha_Na = 2; // Sodium influx per spike in mM ms - const double tau_Na = 1; // Sodium time constant in ms + /* parameters of the firing adaption */ + const double alpha_Na = 2; /* Sodium influx per spike in mM ms */ + const double tau_Na = 1; /* Sodium time constant in 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 + 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 + /* PSP rise time in ms^-1 */ const double gamma_e = 70E-3; const double gamma_i = 58.6E-3; - /* axonal flux time constant */ - const double nu = 60E-3; + /* Axonal flux time constant */ + const double nu = 120E-3; - - // Conductivities in mS/cm^-2 - // Leak current + /* Conductivities in mS/cm^-2 */ + /* Leak */ const double g_L = 1; - // KNa current + /* KNa */ const double g_KNa = 1.33; - // Connectivities (dimensionless) - const double N_ee = 120; - const double N_ei = 72; - const double N_ie = 90; - const double N_ii = 90; - const double N_te = 10; - const double N_ti = 10; - - // Reversal potentials in mV - // synaptic + /* Reversal potentials in mV */ + /* synaptic */ const double E_AMPA = 0; const double E_GABA = -70; - // Leak + /* Leak */ const double E_L_e = -66; const double E_L_i = -64; - // Potassium + /* Potassium */ const double E_K = -100; - // Noise parameters in ms^-1 + /* Noise parameters in ms^-1 */ const double mphi = 0E-3; - const double dphi = 30E-3;; + const double dphi = 60E-3;; double input = 0.0; + /* Connectivities (dimensionless) */ + const double N_ee = 120; + const double N_ei = 72; + const double N_ie = 90; + const double N_ii = 90; + const double N_te = 10; + const double N_ti = 10; + /* Pointer to thalamic column */ Thalamic_Column* Thalamus; - - friend class Stim; }; /****************************************************************************************************/ /* end */ diff --git a/Main.cpp b/Main.cpp index 39d1378..5a72229 100644 --- a/Main.cpp +++ b/Main.cpp @@ -1,3 +1,25 @@ +/* +* Copyright (c) 2014 Michael Schellenberger Costa +* +* Permission is hereby granted, free of charge, to any person obtaining a copy +* of this software and associated documentation files (the "Software"), to deal +* in the Software without restriction, including without limitation the rights +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +* 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 +* 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 +* 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 +* THE SOFTWARE. +*/ + /****************************************************************************************************/ /* Main file for compilation tests */ /* The Simulation requires the following boost libraries: Preprocessor */ @@ -9,13 +31,15 @@ #include "Thalamic_Column.h" #include "ODE.h" + /****************************************************************************************************/ /* Fixed simulation settings */ /****************************************************************************************************/ -extern const int T = 60; -extern const int res = 1E4; -extern const double dt = 1E3/res; -extern const double h = sqrt(dt); +extern const int T = 30; /* Simulation length s */ +extern const int res = 1E4; /* number of iteration steps per s */ +extern const int red = res/100; /* number of iterations that is saved */ +extern const double dt = 1E3/res; /* duration of a timestep in ms */ +extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */ /****************************************************************************************************/ /* end */ /****************************************************************************************************/ @@ -25,27 +49,27 @@ extern const double h = sqrt(dt); /* Main simulation routine */ /****************************************************************************************************/ int main(void) { - // Initializing the seeder. + /* Initializing the seeder */ srand(time(0)); - // Initializing the populations; + /* Initialize the populations */ Cortical_Column Cortex; Thalamic_Column Thalamus; - // Connect both modules + /* Connect both modules */ Cortex.get_Thalamus(Thalamus); Thalamus.get_Cortex(Cortex); - // takes the time of the simulation + /* Take the time of the simulation */ time_t start,end; time (&start); - // simulation + /* Simulation */ for (int t=0; t< T*res; ++t) { ODE (Cortex, Thalamus); } time (&end); - // time consumed by the simulation + /* Time consumed by the simulation */ double dif = difftime(end,start); std::cout << "simulation done!\n"; std::cout << "took " << dif << " seconds" << "\n"; diff --git a/ODE.h b/ODE.h index 2a61de2..44f647c 100644 --- a/ODE.h +++ b/ODE.h @@ -1,3 +1,25 @@ +/* +* Copyright (c) 2014 Michael Schellenberger Costa +* +* Permission is hereby granted, free of charge, to any person obtaining a copy +* of this software and associated documentation files (the "Software"), to deal +* in the Software without restriction, including without limitation the rights +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +* 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 +* 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 +* 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 +* THE SOFTWARE. +*/ + /****************************************************************************************************/ /* Implementation of the ODE solver */ /****************************************************************************************************/ @@ -9,12 +31,13 @@ /* Evaluation of SRK4 */ /****************************************************************************************************/ inline void ODE(Cortical_Column& Cortex, Thalamic_Column& Thalamus) { - // first calculating every ith RK moment - // has to be in order, 1th moment first + /* First calculate every ith RK moment. Has to be in order, 1th moment first */ for (int i=1; i<=4; ++i) { Cortex.set_RK(i); Thalamus.set_RK(i); } + + /* Add all moments */ Cortex.add_RK(); Thalamus.add_RK(); } diff --git a/Plots.m b/Plots.m index b807d4b..bca4beb 100644 --- a/Plots.m +++ b/Plots.m @@ -4,28 +4,21 @@ function Plots(T) if nargin == 0 - % fittet input - Input_N3 = [ 24; % tau_e - -64; % theta_e - 8.3; % sigma_e - 3.6; % alpha_Na - 2.7; % tau_Na - 0.63; % g_KNa - 50E-3]; % dphi + Input_N3 = [ 8.6; % sigma_e + 2.7; % alpha_Na + 3; % tau_Na + 60E-3]; % dphi - Input_N2 = [ 30; % tau_e - -58.5; % theta_e - 4.5; % sigma_e + Input_N2 = [ 4.6; % sigma_e 2; % alpha_Na - 1; % tau_Na - 1.33; % g_KNa - 30E-3]; % dphi + 1.2; % tau_Na + 60E-3]; % dphi - Con = [0; % N_et - 0; % N_er - 0; % N_te - 0]; % N_ti + Con = [0; % N_et + 0; % N_er + 0; % N_te + 0]; % N_ti % stimulation parameters % first number is the mode of stimulation @@ -40,7 +33,7 @@ function Plots(T) 8; % time between stimuli in s 550]; % time until stimuli after min in ms - T = 60; % duration of the simulation + T = 30; % duration of the simulation end [Ve, Vt] = TC(T, Input_N3, Con, var_stim); diff --git a/Stimulation.h b/Stimulation.h index 184f68f..325806f 100644 --- a/Stimulation.h +++ b/Stimulation.h @@ -1,3 +1,25 @@ +/* +* Copyright (c) 2014 Michael Schellenberger Costa +* +* Permission is hereby granted, free of charge, to any person obtaining a copy +* of this software and associated documentation files (the "Software"), to deal +* in the Software without restriction, including without limitation the rights +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +* 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 +* 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 +* 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 +* THE SOFTWARE. +*/ + /****************************************************************************************************/ /* Implementation of the stimulation protocol */ /****************************************************************************************************/ @@ -10,7 +32,7 @@ /****************************************************************************************************/ class Stim { public: - /* empty constructor for compiling */ + /* Empty constructor for compiling */ Stim(void); Stim(Cortical_Column& C, Thalamic_Column& T, double* var) @@ -18,25 +40,25 @@ class Stim { Thalamus = &T; setup(var);} - /* setup with respect to stimulation mode */ + /* Setup with respect to stimulation mode */ void setup (double* var_stim) { extern const int onset; extern const int res; extern const int dt; - /* mode of stimulation */ + /* Mode of stimulation */ mode = (int) var_stim[0]; - /* scale the stimulation strength from s^-1 to ms^-1 */ + /* Scale the stimulation strength from s^-1 to ms^-1 */ strength = var_stim[1] / 1000; - /* scale duration from ms to dt */ + /* Scale duration from ms to dt */ duration = (int) var_stim[2] * res / 1000; - /* scale the ISI from s to ms^-1 */ + /* Scale the ISI from s to ms^-1 */ ISI = (int) var_stim[3] * res; - /* scale time to stimulus from ms to dt */ + /* Scale time to stimulus from ms to dt */ time_to_stim= (int) var_stim[4] * res / 1000; if(mode==1) { @@ -49,22 +71,22 @@ class Stim { void check_stim (int time) { - /* check if stimulation should start */ + /* Check if stimulation should start */ switch (mode) { - /* no stimulation */ + /* No stimulation */ default: break; - /* periodic stimulation */ + /* Periodic stimulation */ case 1: - /* check if time is reached */ + /* Check if time is reached */ if(time == time_to_stim) { - /* switch stimulation on */ + /* Switch stimulation on */ stimulation_started = true; Thalamus->set_input(strength); - /* update the timer */ + /* Update the timer */ time_to_stim += ISI; @@ -74,9 +96,9 @@ class Stim { } break; - /* phase dependent up state stimulation */ + /* Phase dependent up state stimulation */ case 2: - /* search for threshold */ + /* Search for threshold */ if(!stimulation_started && !minimum_found && !threshold_crossed && time>correction) { if(Cortex->Ve[0]<=threshold) { threshold_crossed = true; @@ -84,7 +106,7 @@ class Stim { } } - /* search for minimum */ + /* Search for minimum */ if(threshold_crossed) { if(Cortex->Ve[0]>Ve_old) { threshold_crossed = false; @@ -96,12 +118,12 @@ class Stim { } } - /* wait until the stimulation should start */ + /* Wait until the stimulation should start */ if(minimum_found) { count_to_start++; - /* start stimulation after time_to_stim has passed */ + /* Start stimulation after time_to_stim has passed */ if(count_to_start==time_to_stim) { minimum_found = false; stimulation_started = true; @@ -112,9 +134,9 @@ class Stim { } break; - /* phase dependent down state stimulation */ + /* Phase dependent down state stimulation */ case 3: - /* search for threshold */ + /* Search for threshold */ if(!stimulation_started && !minimum_found && !threshold_crossed && time>correction) { if(Cortex->Ve[0]<=threshold) { threshold_crossed = true; @@ -122,7 +144,7 @@ class Stim { } } - /* search for minimum */ + /* Search for minimum */ if(threshold_crossed) { if(Cortex->Ve[0]>Ve_old) { threshold_crossed = false; @@ -134,7 +156,7 @@ class Stim { } } - /* start the stimulation */ + /* Start the stimulation */ if(minimum_found) { minimum_found = false; stimulation_started = true; @@ -144,11 +166,11 @@ class Stim { break; } - /* wait to switch the stimulation off */ + /* Wait to switch the stimulation off */ if(stimulation_started) { count_duration++; - /* switch stimulation off */ + /* Switch stimulation off */ if(count_duration==duration) { stimulation_started = false; count_duration = 0; @@ -157,6 +179,7 @@ class Stim { } } + /* Create MATLAB container for marker storage */ mxArray* get_marker(void) { mxArray* Marker = mxCreateDoubleMatrix(0, 0, mxREAL); mxSetM(Marker, 3); @@ -174,22 +197,22 @@ class Stim { private: /* Stimulation parameters */ - /* stimulation strength */ + /* Stimulation strength */ double strength = 0.0; - /* duration of the stimulation */ + /* Duration of the stimulation */ int duration = 500; - /* inter stimulus intervall in case of periodic stimulation */ + /* Inter stimulus intervall in case of periodic stimulation */ int ISI = 5E4; - /* threshold for phase dependent stimulation */ + /* Threshold for phase dependent stimulation */ double threshold = -80; - /* time until stimulus after minimum was found */ + /* Time until stimulus after minimum was found */ int time_to_stim = 5500; - /* mode of stimulation */ + /* Mode of stimulation */ /* 0 == none */ /* 1 == periodic */ /* 2 == phase dependent up state */ @@ -200,28 +223,28 @@ class Stim { /* Simulation on for TRUE and off for FALSE */ bool stimulation_started= false; - /* threshold has been reached */ + /* Threshold has been reached */ bool threshold_crossed = false; - /* minimum found */ + /* Minimum found */ bool minimum_found = false; - /* onset in timesteps to correct the given time of the markers */ + /* Onset in timesteps to correct the given time of the markers */ int correction = 10000; - /* counter for stimulation duration */ + /* Counter for stimulation duration */ int count_duration = 0; - /* counter after minimum */ + /* Counter after minimum */ int count_to_start = 0; - /* old voltage value for minimum detection */ + /* Old voltage value for minimum detection */ double Ve_old = 0; - /* pointer to cortical column */ + /* Pointer to cortical column */ Cortical_Column* Cortex; - /* pointer to thalamic column */ + /* Pointer to thalamic column */ Thalamic_Column* Thalamus; /* Data containers */ diff --git a/TC.cpp b/TC.cpp index aa61068..88254d8 100644 --- a/TC.cpp +++ b/TC.cpp @@ -1,3 +1,25 @@ +/* +* Copyright (c) 2014 Michael Schellenberger Costa +* +* Permission is hereby granted, free of charge, to any person obtaining a copy +* of this software and associated documentation files (the "Software"), to deal +* in the Software without restriction, including without limitation the rights +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +* 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 +* 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 +* 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 +* THE SOFTWARE. +*/ + /********************************************************************************************************/ /* Implementation of the simulation as MATLAB routine (mex compiler) */ /* mex command is given by: */ @@ -17,11 +39,11 @@ /****************************************************************************************************/ /* Fixed simulation settings */ /****************************************************************************************************/ -extern const int onset = 5; // time until data is stored in s -extern const int res = 1E4; // number of iteration steps per s -extern const int red = res/100; // number of iterations that is saved -extern const double dt = 1E3/res; // duration of a timestep in ms -extern const double h = sqrt(dt); // squareroot of dt for SRK iteration +extern const int onset = 5; /* time until data is stored in s */ +extern const int res = 1E4; /* number of iteration steps per s */ +extern const int red = res/100; /* number of iterations that is saved */ +extern const double dt = 1E3/res; /* duration of a timestep in ms */ +extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */ /****************************************************************************************************/ /* end */ /****************************************************************************************************/ @@ -31,35 +53,36 @@ extern const double h = sqrt(dt); // squareroot of dt for SRK iteration /* Simulation routine */ /****************************************************************************************************/ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - // Initializing the seeder. + /* Initialize the seeder */ srand(time(NULL)); - // Fetch inputs - const int T = (int) (mxGetScalar(prhs[0])); // duration of simulation in s - const int Time = (T+onset)*res; // total number of iteration steps - double* Param_Cortex = mxGetPr (prhs[1]); // parameters of cortical module - double* Connections = mxGetPr (prhs[2]); // parameters of thalamic module - double* var_stim = mxGetPr (prhs[3]); // parameters of stimulation protocol + /* Fetch inputs */ + const int T = (int) (mxGetScalar(prhs[0])); /* Duration of simulation in s */ + const int Time = (T+onset)*res; /* Total number of iteration steps */ + double* Param_Cortex = mxGetPr (prhs[1]); /* Parameters of cortical module */ + double* Connections = mxGetPr (prhs[2]); /* Parameters of thalamic module */ + double* var_stim = mxGetPr (prhs[3]); /* Parameters of stimulation protocol */ - // Initializing the populations; + /* Initialize the populations */ Cortical_Column Cortex(Param_Cortex, Connections); Thalamic_Column Thalamus(Connections); - // Linking both modules + /* Link both modules */ Cortex.get_Thalamus(Thalamus); Thalamus.get_Cortex(Cortex); - // Initialize the stimulation protocol + /* Initialize the stimulation protocol */ Stim Stimulation(Cortex, Thalamus, var_stim); - // setting up the data containers + /* Create data containers */ mxArray* Ve = SetMexArray(1, T*res/red); mxArray* Vt = SetMexArray(1, T*res/red); - // Pointer to the actual data block + /* Pointer to the actual data block */ double* Pr_Ve = mxGetPr(Ve); double* Pr_Vt = mxGetPr(Vt); + /* Simulation */ int count = 0; for (int t=0; t DIST; // Normal Distribution -typedef boost::variate_generator GEN; // Variate generator +typedef boost::mt11213b ENG; /* Mersenne Twister */ +typedef boost::normal_distribution DIST; /* Normal Distribution */ +typedef boost::variate_generator GEN; /* Variate generator */ /****************************************************************************************************/ /* end */ /****************************************************************************************************/ @@ -29,38 +51,38 @@ typedef boost::variate_generator GEN; // Variate generator /****************************************************************************************************/ class Thalamic_Column { public: - // Constructors + /* Constructors for compile check */ Thalamic_Column(void) {set_RNG();} - // Constructors - Thalamic_Column(double* Par) - : N_et (Par[0]), N_er (Par[1]) + /* Constructor for simulation */ + Thalamic_Column(double* Con) + : N_et (Con[0]), N_er (Con[1]) {set_RNG();} - // get the pointer to the cortical module + /* Get the pointer to the cortical module */ void get_Cortex (Cortical_Column& C) {Cortex = &C;} - // change the strength of input + /* Set strength of input */ void set_input (double I) {input = I;} - // get axonal flux + /* Get axonal flux */ double get_phi (int N) const {_SWITCH((phi_t)); return var_phi_t;} - // Initialize the RNGs + /* Initialize the RNGs */ void set_RNG (void); - // Firing rates + /* Firing rates */ double get_Qt (int) const; double get_Qr (int) const; - // Synaptic currents + /* Synaptic currents */ double I_et (int) const; double I_it (int) const; double I_er (int) const; double I_ir (int) const; - // Activation functions + /* Activation functions */ double m_inf_T_t (int) const; double m_inf_T_r (int) const; double tau_m_T_t (int) const; @@ -68,13 +90,13 @@ class Thalamic_Column { double m_inf_h (int) const; double tau_m_h (int) const; - // Deactivation functions + /* Deactivation functions */ double h_inf_T_t (int) const; double h_inf_T_r (int) const; double tau_h_T_t (int) const; double tau_h_T_r (int) const; - // Currents + /* Currents */ double I_L_t (int) const; double I_L_r (int) const; double I_LK_t (int) const; @@ -83,117 +105,117 @@ class Thalamic_Column { double I_T_r (int) const; double I_h (int) const; - // Noise functions + /* Noise functions */ double noise_xRK (int) const; - // ODE functions + /* ODE functions */ void set_RK (int); void add_RK (void); - // Data storage + /* Data storage access */ friend void get_data (int, Cortical_Column&, Thalamic_Column&, _REPEAT(double*, 2)); private: - // Population variables - vector Vt = _INIT(E_L_t), // TC membrane voltage - Vr = _INIT(E_L_r), // RE membrane voltage - Ca = _INIT(Ca_0), // Calcium concentration of TC population - Phi_tt = _INIT(0.0), // PostSP from TC population to TC population - Phi_tr = _INIT(0.0), // PostSP from TC population to RE population - Phi_rt = _INIT(0.0), // PostSP from RE population to TC population - Phi_rr = _INIT(0.0), // PostSP from RE population to RE population - phi_t = _INIT(0.0), // axonal flux - x_tt = _INIT(0.0), // derivative of Phi_tt - x_tr = _INIT(0.0), // derivative of Phi_tr - x_rt = _INIT(0.0), // derivative of Phi_rt - x_rr = _INIT(0.0), // derivative of Phi_rr - y_t = _INIT(0.0), // derivative of phi_t - h_T_t = _INIT(0.0), // inactivation of T channel - h_T_r = _INIT(0.0), // inactivation of T channel - m_T_t = _INIT(0.0), // activation of T channel - m_T_r = _INIT(0.0), // activation of T channel - m_h = _INIT(0.0), // activation of h channel - m_h2 = _INIT(0.0), // activation of h channel bound with protein - P_h = _INIT(0.0); // fraction of protein bound with calcium - - // Random number generators + /* Population variables */ + vector Vt = _INIT(E_L_t), /* TC membrane voltage */ + Vr = _INIT(E_L_r), /* RE membrane voltage */ + Ca = _INIT(Ca_0), /* Calcium concentration of TC population */ + Phi_tt = _INIT(0.0), /* PostSP from TC population to TC population */ + Phi_tr = _INIT(0.0), /* PostSP from TC population to RE population */ + Phi_rt = _INIT(0.0), /* PostSP from RE population to TC population */ + Phi_rr = _INIT(0.0), /* PostSP from RE population to RE population */ + phi_t = _INIT(0.0), /* axonal flux */ + x_tt = _INIT(0.0), /* derivative of Phi_tt */ + x_tr = _INIT(0.0), /* derivative of Phi_tr */ + x_rt = _INIT(0.0), /* derivative of Phi_rt */ + x_rr = _INIT(0.0), /* derivative of Phi_rr */ + y_t = _INIT(0.0), /* derivative of phi_t */ + h_T_t = _INIT(0.0), /* inactivation of T channel */ + h_T_r = _INIT(0.0), /* inactivation of T channel */ + m_T_t = _INIT(0.0), /* activation of T channel */ + m_T_r = _INIT(0.0), /* activation of T channel */ + m_h = _INIT(0.0), /* activation of h channel */ + m_h2 = _INIT(0.0), /* activation of h channel bound with protein */ + P_h = _INIT(0.0); /* fraction of protein bound with calcium */ + + /* Random number generators */ vector MTRands; - // Container for noise + /* Container for noise */ vector Rand_vars; - // Declaration and Initialization of parameters - // Membrane time in ms + /* Declaration and Initialization of parameters */ + /* Membrane time in ms */ const double tau_t = 30; const double tau_r = 30; - // Maximum firing rate in ms^-1 + /* Maximum firing rate in ms^-1 */ const double Qt_max = 400.E-3; const double Qr_max = 400.E-3; - // Sigmoid threshold in mV + /* Sigmoid threshold in mV */ const double theta_t = -45; const double theta_r = -45; - // Sigmoid gain in mV - const double sigma_t = 9; - const double sigma_r = 9; + /* Sigmoid gain in mV */ + const double sigma_t = 3; + const double sigma_r = 3; - // Scaling parameter for sigmoidal mapping (dimensionless) + /* Scaling parameter for sigmoidal mapping (dimensionless) */ const double C1 = (3.14159265/sqrt(3)); - // PSP rise time in ms^-1 + /* PSP rise time in ms^-1 */ const double gamma_e = 70E-3; const double gamma_i = 58.6E-3; /* axonal flux time constant */ - const double nu = 60E-3; + const double nu = 120E-3; - // Conductivities in mS/cm^-2 - // Leak current - const double g_L_t = 0.02; - const double g_L_r = 0.05; + /* Conductivities in mS/cm^-2 */ + /* Leak current */ + const double g_L_t = 0.6; + const double g_L_r = 1.5; - // Potassium leak current + /* Potassium leak current */ const double g_LK_t = 0.02; - const double g_LK_r = 0.01; + const double g_LK_r = 0.02; - // T current + /* T current */ const double g_T_t = 2.2; const double g_T_r = 2; /* h current */ const double g_h = 0.07; - // Connectivities (dimensionless) + /* Connectivities (dimensionless) */ const double N_tr = 2; const double N_rt = 5.5; const double N_rr = 5; const double N_et = 10; const double N_er = 10; - // Reversal potentials in mV - // synaptic + /* Reversal potentials in mV */ + /* Synaptic */ const double E_AMPA = 0; const double E_GABA = -70; - // Leak - const double E_L_t = -70; - const double E_L_r = -55; + /* Leak */ + const double E_L_t = -64; + const double E_L_r = -64; - // Potassium + /* Potassium */ const double E_K = -100; - // I_T current + /* I_T current */ const double E_Ca = 120; - // I_h current + /* I_h current */ const double E_h = -40; /* Calcium parameters */ - const double alpha_Ca = -50E-6; // influx per spike in nmol - const double tau_Ca = 5; // calcium time constant in ms - const double Ca_0 = 2E-4; // resting concentration + const double alpha_Ca = -50E-6; /* influx per spike in nmol */ + const double tau_Ca = 5; /* calcium time constant in ms */ + const double Ca_0 = 2E-4; /* resting concentration */ /* I_h activation parameters */ const double k1 = 2.5E7; @@ -203,7 +225,7 @@ class Thalamic_Column { const double n_P = 4; const double g_inc = 2; - // Noise parameters in ms^-1 + /* Noise parameters in ms^-1 */ const double mphi = 0E-3; const double dphi = 2E-3;; double input = 0.0; diff --git a/macros.h b/macros.h index 8fcf5b4..1baa7a1 100644 --- a/macros.h +++ b/macros.h @@ -1,3 +1,25 @@ +/* +* Copyright (c) 2014 Michael Schellenberger Costa +* +* Permission is hereby granted, free of charge, to any person obtaining a copy +* of this software and associated documentation files (the "Software"), to deal +* in the Software without restriction, including without limitation the rights +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +* 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 +* 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 +* 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 +* THE SOFTWARE. +*/ + /****************************************************************************************************/ /* Definition of all macros used */ /****************************************************************************************************/ @@ -9,6 +31,7 @@ #include #include + /****************************************************************************************************/ /* Macro for vector initialization */ /****************************************************************************************************/ diff --git a/saves.h b/saves.h index fc9f3c6..927a5ea 100644 --- a/saves.h +++ b/saves.h @@ -1,12 +1,36 @@ +/* +* Copyright (c) 2014 Michael Schellenberger Costa +* +* Permission is hereby granted, free of charge, to any person obtaining a copy +* of this software and associated documentation files (the "Software"), to deal +* in the Software without restriction, including without limitation the rights +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +* 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 +* 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 +* 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 +* THE SOFTWARE. +*/ + /****************************************************************************************************/ -/* functions for data storage */ +/* Functions for data storage */ /****************************************************************************************************/ #pragma once #include #include "Cortical_Column.h" #include "Thalamic_Column.h" -// saving the fluctuations of the populations +/****************************************************************************************************/ +/* Save data */ +/****************************************************************************************************/ inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thalamus, double* Ve, double* Vt) { Ve [counter] = Cortex.Ve [0]; Vt [counter] = Thalamus.Vt [0]; @@ -17,9 +41,8 @@ inline void get_data(int counter, Cortical_Column& Cortex, Thalamic_Column& Thal /****************************************************************************************************/ -/* function that creates a matlab data container */ +/* Create MATLAB data container */ /****************************************************************************************************/ -// function to create a MATLAB data container mxArray* SetMexArray(int N, int M) { mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL); mxSetM(Array, N);