diff --git a/.cproject b/.cproject deleted file mode 100644 index 1fee099..0000000 --- a/.cproject +++ /dev/null @@ -1,177 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/.project b/.project deleted file mode 100644 index b0dd08a..0000000 --- a/.project +++ /dev/null @@ -1,83 +0,0 @@ - - - Steyn-Ross - - - - - - org.eclipse.cdt.managedbuilder.core.genmakebuilder - clean,full,incremental, - - - ?name? - - - - org.eclipse.cdt.make.core.append_environment - true - - - org.eclipse.cdt.make.core.autoBuildTarget - all - - - org.eclipse.cdt.make.core.buildArguments - - - - org.eclipse.cdt.make.core.buildCommand - make - - - org.eclipse.cdt.make.core.buildLocation - ${workspace_loc:/Steyn-Ross/Debug} - - - org.eclipse.cdt.make.core.cleanBuildTarget - clean - - - org.eclipse.cdt.make.core.contents - org.eclipse.cdt.make.core.activeConfigSettings - - - org.eclipse.cdt.make.core.enableAutoBuild - false - - - org.eclipse.cdt.make.core.enableCleanBuild - true - - - org.eclipse.cdt.make.core.enableFullBuild - true - - - org.eclipse.cdt.make.core.fullBuildTarget - all - - - org.eclipse.cdt.make.core.stopOnError - true - - - org.eclipse.cdt.make.core.useDefaultBuildCmd - true - - - - - org.eclipse.cdt.managedbuilder.core.ScannerConfigBuilder - full,incremental, - - - - - - org.eclipse.cdt.core.cnature - org.eclipse.cdt.core.ccnature - org.eclipse.cdt.managedbuilder.core.managedBuildNature - org.eclipse.cdt.managedbuilder.core.ScannerConfigNature - - diff --git a/Data_Storage.h b/Data_Storage.h new file mode 100644 index 0000000..6ebed5a --- /dev/null +++ b/Data_Storage.h @@ -0,0 +1,47 @@ +/* + * Copyright (c) 2015 University 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 + * 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. + * + * AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com + * + * Based on: A thalamocortical neural mass model of the EEG during NREM sleep and its response + * to auditory stimulation. + * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, + * JC Claussen. + * PLoS Computational Biology In Review (in review). + */ + +/****************************************************************************************************/ +/* Functions for data storage */ +/****************************************************************************************************/ +#pragma once +#include "Thalamic_Column.h" + +/****************************************************************************************************/ +/* Save data */ +/****************************************************************************************************/ +void get_data(int counter, Thalamic_Column& Col, double* Vt, double* Vr, double* ah) { + Vt [counter] = Col.Vt [0]; + Vr [counter] = Col.Vr [0]; + ah [counter] = Col.act_h (); +} +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ \ No newline at end of file diff --git a/Main.cpp b/Main.cpp deleted file mode 100644 index d38e8c3..0000000 --- a/Main.cpp +++ /dev/null @@ -1,86 +0,0 @@ -/* -* 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 */ -/* Random */ -/****************************************************************************************************/ -#include -#include -#include "Thalamic_Column.h" - -/****************************************************************************************************/ -/* Fixed simulation settings */ -/****************************************************************************************************/ -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 */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Constants for SRK4 iteration */ -/****************************************************************************************************/ -extern const vector B1 = {0, - 0.626708569400000081728308032325, - 1.7296310295000001389098542858846, - 1.2703689705000000831347506391467}; -extern const vector B2 = {0, - 0.78000033203198970710445792065002, - 1.28727807507536762265942797967, - 0.44477273249350995909523476257164}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Main simulation routine */ -/****************************************************************************************************/ -int main(void) { - /* Initialize the populations */ - Thalamic_Column Thalamus; - - /* Takes the time of the simulation */ - time_t start,end; - time (&start); - - /* Simulation */ - for (int t=0; t< T*res; ++t) { - Thalamus.iterate_ODE(); - } - time (&end); - - /* Time consumed by the simulation */ - double dif = difftime(end,start); - std::cout << "simulation done!\n"; - std::cout << "took " << dif << " seconds" << "\n"; - std::cout << "end\n"; -} -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/NM_Thalamus.pro b/NM_Thalamus.pro new file mode 100644 index 0000000..40a9042 --- /dev/null +++ b/NM_Thalamus.pro @@ -0,0 +1,19 @@ +TEMPLATE = app +CONFIG += console +CONFIG -= app_bundle +CONFIG -= qt + +TARGET = Thalamus.cpp + +SOURCES += Thalamic_Column.cpp \ + Thalamus.cpp \ + Thalamus_mex.cpp + +HEADERS += ODE.h \ + Data_Storage.h \ + Random_Stream.h \ + Thalamic_Column.h + +QMAKE_CXXFLAGS += -std=c++11 -O3 + +SOURCES -= Thalamus_mex.cpp diff --git a/Plots.m b/Plots.m index 7226504..2118f09 100644 --- a/Plots.m +++ b/Plots.m @@ -1,51 +1,24 @@ % mex command is given by: -% mex CXXFLAGS="\$CXXFLAGS -std=c++11" Thalamus.cpp Thalamic_Column.cpp +% mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" Thalamus_mex.cpp Thalamic_Column.cpp function Plots(T) if nargin == 0 - Con = [ 6; % sigma - 4; % N_tr - 4; % N_rt - 20]; % N_rr - + Con = [ 0.01; % g_h + 0.05]; % g_LK - % stimulation parameters - % first number is the mode of stimulation - % 0 == none - % 1 == periodic - % 2 == phase dependend up state - % 3 == phase dependend down state - - var_stim = [ 0; % mode of stimulation - 50; % strength of the stimulus in Hz (spikes per second) - 120; % duration of the stimulus in ms - 7; % time between stimuli in s - 1]; % time until stimuli after min in ms T = 30; % duration of the simulation end -[Vt, Vr, I_t, I_h] = Thalamus(T, Con, var_stim); +[Vt, Vr, ah] = Thalamus_mex(T, Con, var_stim); L = max(size(Vt)); timeaxis = linspace(0,T,L); -Fs = L/T; -[Pxx,f] = pwelch(Vt-mean(Vt), hamming(4*Fs), 2*Fs, [], Fs,'ConfidenceLevel', 0.9, 'power'); -n = find(f<=60, 1, 'last' ); - +% Plot the data figure(1) subplot(311), plot(timeaxis,Vt) -title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt in mV') +title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt [mV]') subplot(312), plot(timeaxis,Vr) -title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV') -subplot(313), plot(f(1:n),Pxx(1:n)) -title('Powerspectrum of Steyn-Ross model with pwelch'), xlabel('frequency in Hz'), ylabel('Power') -%save('Thalamus.mat','Vt','Vr') -% -% figure(2) -% subplot(311), plot(timeaxis,Vt) -% title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt in mV') -% subplot(312), plot(timeaxis,I_t) -% title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV') -% subplot(313), plot(timeaxis,I_h) -% title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV') \ No newline at end of file +title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr [mV]') +subplot(313), plot(timeaxis,ah) +title('Thalamic relay I_h activation'), xlabel('time in s'), ylabel('m_h') diff --git a/Random_Stream.h b/Random_Stream.h new file mode 100644 index 0000000..f5fb19a --- /dev/null +++ b/Random_Stream.h @@ -0,0 +1,80 @@ +/* + * Copyright (c) 2015 University 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 + * 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. + * + * AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com + * Stefanie Gareis: gareis@inb.uni-luebeck.de + */ + +/****************************************************************************************************/ +/* Random number streams */ +/****************************************************************************************************/ +#pragma once +#include + +/****************************************************************************************************/ +/* Struct for normal distribution */ +/****************************************************************************************************/ +struct random_stream_normal +{ + /* Random number engine: Mersenne-Twister */ + std::mt19937_64 mt; + /* Random number generator: Normal-distribution */ + std::normal_distribution norm_dist; + + /* Constructors */ + random_stream_normal(){} + random_stream_normal(double mean, double stddev) + : mt(rand()) , norm_dist(mean, stddev) + {} + + /* Overwrites the function-call operator "( )" */ + double operator( )(void) { + return norm_dist(mt); + } +}; +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ + +/****************************************************************************************************/ +/* Struct for uniform int distribution */ +/****************************************************************************************************/ +struct random_stream_uniform_int +{ + /* Random number engine: Mersenne-Twister */ + std::mt19937_64 mt; + /* Random number generator: Uniform integer-distribution */ + std::uniform_int_distribution<> uniform_dist; + + /* Constructors */ + random_stream_uniform_int(){} + random_stream_uniform_int(double lower_bound, double upper_bound) + : mt(rand()) , uniform_dist(lower_bound, upper_bound) + {} + + /* Overwrites the function-call operator "( )" */ + double operator( )(void) { + return uniform_dist(mt); + } +}; +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/Stimulation.h b/Stimulation.h deleted file mode 100644 index b00fa70..0000000 --- a/Stimulation.h +++ /dev/null @@ -1,169 +0,0 @@ -/* -* 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 */ -/****************************************************************************************************/ -#pragma once -#include "Thalamic_Column.h" - -/****************************************************************************************************/ -/* Stimulation object */ -/****************************************************************************************************/ -class Stim { -public: - /* Empty constructor for compiling */ - Stim(void); - - Stim(Thalamic_Column& T, double* var) - { Thalamus = &T; - setup(var);} - - /* 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 = (int) var_stim[0]; - - /* Scale the stimulation strength from s^-1 to ms^-1 */ - strength = var_stim[1] / 1000; - - /* Scale duration from ms to dt */ - duration = (int) var_stim[2] * res / 1000; - - /* Scale the ISI from s to ms^-1 */ - ISI = (int) var_stim[3] * res; - - /* Scale time to stimulus from ms to dt */ - time_to_stim= (int) var_stim[4] * res / 1000; - - if(mode==1) { - time_to_stim = (onset+1) * res; - } - - correction = onset * res; - - } - - void check_stim (int time) { - - /* Check if stimulation should start */ - switch (mode) { - - /* No stimulation */ - default: - break; - - /* Periodic stimulation */ - case 1: - /* Check if time is reached */ - if(time == time_to_stim) { - /* Switch stimulation on */ - stimulation_started = true; - Thalamus->set_input(strength); - - /* Update the timer */ - time_to_stim += ISI; - - /* Update the marker */ - marker_stimulation.push_back(time - correction); - } - break; - } - - /* Wait to switch the stimulation off */ - if(stimulation_started) { - count_duration++; - - /* Switch stimulation off */ - if(count_duration==duration) { - stimulation_started = false; - count_duration = 0; - Thalamus->set_input(0.0); - } - } - } - - /* Create MATLAB container for marker storage */ - mxArray* get_marker(void) { - mxArray* Marker = mxCreateDoubleMatrix(0, 0, mxREAL); - mxSetM(Marker, 1); - mxSetN(Marker, marker_stimulation.size()); - mxSetData(Marker, mxMalloc(sizeof(double)*marker_stimulation.size())); - double* Pr_Marker = mxGetPr(Marker); - for(unsigned i=0; i marker_stimulation; -}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index f25e461..82c9525 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -1,44 +1,56 @@ /* -* 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. -*/ + * Copyright (c) 2015 University 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 + * 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. + * + * AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com + * + * Based on: A thalamocortical neural mass model of the EEG during NREM sleep and its response + * to auditory stimulation. + * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, + * JC Claussen. + * PLoS Computational Biology In Review (in review). + */ /****************************************************************************************************/ /* Functions of the thalamic module */ /****************************************************************************************************/ #include "Thalamic_Column.h" - /****************************************************************************************************/ /* Initialization of RNG */ /****************************************************************************************************/ void Thalamic_Column::set_RNG(void) { - /* Number of independent streams */ - int N = 2; + extern const double dt; + /* Number of independent random variables */ + int N = 1; /* Create RNG for each stream */ for (int i=0; i B1, B2; - double n = 1 / h * (B1[N-1] * Rand_vars[0] + B2[N-1] * Rand_vars[1]); - return n; +double Thalamic_Column::noise_xRK(int N, int M) const{ + return gamma_e * gamma_e * (Rand_vars[2*M] + Rand_vars[2*M+1]/std::sqrt(3))*B[N]; +} + +double Thalamic_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 */ @@ -247,27 +251,21 @@ double Thalamic_Column::noise_xRK(int N) const{ /****************************************************************************************************/ void Thalamic_Column::set_RK (int N) { extern const double dt; - _SWITCH((Ca) - (Phi_tt)(Phi_tr)(Phi_rt)(Phi_rr) - (x_tt) (x_tr) (x_rt) (x_rr) - (h_T_t) (h_T_r) - (m_h) (m_h2) (P_h)) - Vt [N] = dt*(-(I_L_t(N) + I_et(N) + I_it(N))/tau_t - (I_LK_t(N) + I_T_t(N) + I_h(N))); - Vr [N] = dt*(-(I_L_r(N) + I_er(N) + I_ir(N))/tau_r - (I_LK_r(N) + I_T_r(N))); - Ca [N] = dt*(alpha_Ca * I_T_t(N) - (var_Ca - Ca_0)/tau_Ca); - h_T_t [N] = dt*(h_inf_T_t(N) - var_h_T_t)/tau_h_T_t(N); - h_T_r [N] = dt*(h_inf_T_r(N) - var_h_T_r)/tau_h_T_r(N); - m_h [N] = dt*((m_inf_h(N) * (1 - var_m_h2) - var_m_h)/tau_m_h(N) - k3 * var_P_h * var_m_h + k4 * var_m_h2); - m_h2 [N] = dt*(k3 * var_P_h * var_m_h - k4 * var_m_h2); - P_h [N] = dt*(k1 * pow(var_Ca, n_P) * (1 - var_P_h) - k2 * var_P_h); - Phi_tt [N] = dt*(var_x_tt); - Phi_tr [N] = dt*(var_x_tr); - Phi_rt [N] = dt*(var_x_rt); - Phi_rr [N] = dt*(var_x_rr); - x_tt [N] = dt*(pow(gamma_e, 2) * (noise_xRK(N) - var_Phi_tt) - 2 * gamma_e * var_x_tt); - x_tr [N] = dt*(pow(gamma_e, 2) * (N_tr * get_Qt(N) - var_Phi_tr) - 2 * gamma_e * var_x_tr); - x_rt [N] = dt*(pow(gamma_i, 2) * (N_rt * get_Qr(N) - var_Phi_rt) - 2 * gamma_i * var_x_rt); - x_rr [N] = dt*(pow(gamma_i, 2) * (N_rr * get_Qr(N) - var_Phi_rr) - 2 * gamma_i * var_x_rr); + Vt [N+1] = Vt [0] + A[N]*dt*(-(I_L_t(N) + I_et(N) + I_gt(N))/tau_t - (I_LK_t(N) + I_T_t(N) + I_h(N))); + Vr [N+1] = Vr [0] + A[N]*dt*(-(I_L_r(N) + I_er(N) + I_gr(N))/tau_r - (I_LK_r(N) + I_T_r(N))); + Ca [N+1] = Ca [0] + A[N]*dt*(alpha_Ca * I_T_t(N) - (Ca[N] - Ca_0)/tau_Ca); + h_T_t [N+1] = h_T_t[0] + A[N]*dt*(h_inf_T_t(N) - h_T_t[N])/tau_h_T_t(N); + h_T_r [N+1] = h_T_r[0] + A[N]*dt*(h_inf_T_r(N) - h_T_r[N])/tau_h_T_r(N); + m_h [N+1] = m_h [0] + A[N]*dt*((m_inf_h(N) * (1 - m_h2[N]) - m_h[N])/tau_m_h(N) - k3 * P_h(N) * m_h[N] + k4 * m_h2[N]); + m_h2 [N+1] = m_h2 [0] + A[N]*dt*(k3 * P_h(N) * m_h[N] - k4 * m_h2[N]); + s_tt [N+1] = s_tt [0] + A[N]*dt*(x_tt[N]); + s_tr [N+1] = s_tr [0] + A[N]*dt*(x_tr[N]); + s_rt [N+1] = s_rt [0] + A[N]*dt*(x_rt[N]); + s_rr [N+1] = s_rr [0] + A[N]*dt*(x_rr[N]); + x_tt [N+1] = x_tt [0] + A[N]*dt*(pow(gamma_e, 2) * ( - s_tt[N]) - 2 * gamma_e * x_tt[N]); + x_tr [N+1] = x_tr [0] + A[N]*dt*(pow(gamma_e, 2) * (N_tr * get_Qt(N) - s_tr[N]) - 2 * gamma_e * x_tr[N]); + x_rt [N+1] = x_rt [0] + A[N]*dt*(pow(gamma_g, 2) * (N_rt * get_Qr(N) - s_rt[N]) - 2 * gamma_g * x_rt[N]); + x_rr [N+1] = x_rr [0] + A[N]*dt*(pow(gamma_g, 2) * (N_rr * get_Qr(N) - s_rr[N]) - 2 * gamma_g * x_rr[N]); } /****************************************************************************************************/ /* end */ @@ -278,23 +276,21 @@ void Thalamic_Column::set_RK (int N) { /* Function that adds all SRK terms */ /****************************************************************************************************/ void Thalamic_Column::add_RK(void) { - extern const double h; - Vt [0] += (Vt [1] + Vt [2] * 2 + Vt [3] * 2 + Vt [4])/6; - Vr [0] += (Vr [1] + Vr [2] * 2 + Vr [3] * 2 + Vr [4])/6; - Ca [0] += (Ca [1] + Ca [2] * 2 + Ca [3] * 2 + Ca [4])/6; - Phi_tt [0] += (Phi_tt [1] + Phi_tt [2] * 2 + Phi_tt [3] * 2 + Phi_tt [4])/6; - Phi_tr [0] += (Phi_tr [1] + Phi_tr [2] * 2 + Phi_tr [3] * 2 + Phi_tr [4])/6; - Phi_rt [0] += (Phi_rt [1] + Phi_rt [2] * 2 + Phi_rt [3] * 2 + Phi_rt [4])/6; - Phi_rr [0] += (Phi_rr [1] + Phi_rr [2] * 2 + Phi_rr [3] * 2 + Phi_rr [4])/6; - x_tt [0] += (x_tt [1] + x_tt [2] * 2 + x_tt [3] * 2 + x_tt [4])/6 + pow(gamma_e, 2) * h * Rand_vars[0]; - x_tr [0] += (x_tr [1] + x_tr [2] * 2 + x_tr [3] * 2 + x_tr [4])/6; - x_rt [0] += (x_rt [1] + x_rt [2] * 2 + x_rt [3] * 2 + x_rt [4])/6; - x_rr [0] += (x_rr [1] + x_rr [2] * 2 + x_rr [3] * 2 + x_rr [4])/6; - h_T_t [0] += (h_T_t [1] + h_T_t [2] * 2 + h_T_t [3] * 2 + h_T_t [4])/6; - h_T_r [0] += (h_T_r [1] + h_T_r [2] * 2 + h_T_r [3] * 2 + h_T_r [4])/6; - m_h [0] += (m_h [1] + m_h [2] * 2 + m_h [3] * 2 + m_h [4])/6; - m_h2 [0] += (m_h2 [1] + m_h2 [2] * 2 + m_h2 [3] * 2 + m_h2 [4])/6; - P_h [0] += (P_h [1] + P_h [2] * 2 + P_h [3] * 2 + P_h [4])/6; + Vt [0] =(-3*Vt [0] + 2*Vt [1] + 4*Vt [2] + 2* Vt [3] + Vt [4])/6; + Vr [0] =(-3*Vr [0] + 2*Vr [1] + 4*Vr [2] + 2* Vr [3] + Vr [4])/6; + Ca [0] =(-3*Ca [0] + 2*Ca [1] + 4*Ca [2] + 2* Ca [3] + Ca [4])/6; + s_tt [0] =(-3*s_tt [0] + 2*s_tt [1] + 4*s_tt [2] + 2* s_tt [3] + s_tt [4])/6; + s_tr [0] =(-3*s_tr [0] + 2*s_tr [1] + 4*s_tr [2] + 2* s_tr [3] + s_tr [4])/6; + s_rt [0] =(-3*s_rt [0] + 2*s_rt [1] + 4*s_rt [2] + 2* s_rt [3] + s_rt [4])/6; + s_rr [0] =(-3*s_rr [0] + 2*s_rr [1] + 4*s_rr [2] + 2* s_rr [3] + s_rr [4])/6; + x_tt [0] =(-3*x_tt [0] + 2*x_tt [1] + 4*x_tt [2] + 2* x_tt [3] + x_tt [4])/6; + x_tr [0] =(-3*x_tr [0] + 2*x_tr [1] + 4*x_tr [2] + 2* x_tr [3] + x_tr [4])/6; + x_rt [0] =(-3*x_rt [0] + 2*x_rt [1] + 4*x_rt [2] + 2* x_rt [3] + x_rt [4])/6; + x_rr [0] =(-3*x_rr [0] + 2*x_rr [1] + 4*x_rr [2] + 2* x_rr [3] + x_rr [4])/6; + h_T_t [0] =(-3*h_T_t[0] + 2*h_T_t[1] + 4*h_T_t[2] + 2* h_T_t [3] + h_T_t [4])/6; + h_T_r [0] =(-3*h_T_r[0] + 2*h_T_r[1] + 4*h_T_r[2] + 2* h_T_r [3] + h_T_r [4])/6; + m_h [0] =(-3*m_h [0] + 2*m_h [1] + 4*m_h [2] + 2* m_h [3] + m_h [4])/6; + m_h2 [0] =(-3*m_h2 [0] + 2*m_h2 [1] + 4*m_h2 [2] + 2* m_h2 [3] + m_h2 [4])/6; /* Generate noise for the next iteration */ for (unsigned i=0; i #include -#include -#include -#include -#include "macros.h" +#include "Random_Stream.h" using std::vector; /****************************************************************************************************/ -/* Typedefs for RNG */ +/* Macro for vector initialization */ /****************************************************************************************************/ -typedef boost::mt11213b ENG; /* Mersenne Twister */ -typedef boost::normal_distribution DIST; /* Normal Distribution */ -typedef boost::variate_generator GEN; /* Variate generator */ +#ifndef _INIT +#define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0} +#endif /****************************************************************************************************/ /* end */ /****************************************************************************************************/ @@ -53,14 +58,19 @@ class Thalamic_Column { {set_RNG();} /* Constructor for simulation */ - Thalamic_Column(double* Con) - : sigma_t (Con[0]), sigma_r (Con[0]), - N_tr (Con[1]), N_rt (Con[2]), N_rr (Con[3]) + Thalamic_Column(double* Par) + : g_LK (Par[1]), g_h (Par[0]), + N_tr (Par[2]), N_rt (Par[3]), N_rr (Par[4]) {set_RNG();} - /* Set strength of input */ - void set_input (double I) {input = I;} + /* Iterate one time step through SRK4 */ + void iterate_ODE (void); + + /* Data storage access */ + friend void get_data (int, Thalamic_Column&, double*, double*, double*); +private: + /* Declaration of private functions */ /* Initialize the RNGs */ void set_RNG (void); @@ -70,15 +80,17 @@ class Thalamic_Column { /* Synaptic currents */ double I_et (int) const; - double I_it (int) const; + double I_gt (int) const; double I_er (int) const; - double I_ir (int) const; + double I_gr (int) const; /* Activation functions */ double m_inf_T_t (int) const; double m_inf_T_r (int) const; double m_inf_h (int) const; double tau_m_h (int) const; + double P_h (int) const; + double act_h (void)const; /* Deactivation functions */ double h_inf_T_t (int) const; @@ -96,40 +108,12 @@ class Thalamic_Column { double I_h (int) const; /* Noise functions */ - double noise_xRK (int) const; + double noise_xRK (int,int) const; + double noise_aRK (int) const; /* ODE functions */ void set_RK (int); void add_RK (void); - void iterate_ODE (void); - - /* Data storage access */ - friend void get_data (int, Thalamic_Column&, _REPEAT(double*, 4)); - -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 */ - 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 */ - h_T_t = _INIT(0.0), /* inactivation of T channel */ - h_T_r = _INIT(0.0), /* inactivation 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 */ - vector Rand_vars; /* Declaration and Initialization of parameters */ /* Membrane time in ms */ @@ -149,27 +133,29 @@ class Thalamic_Column { const double sigma_r = 6; /* Scaling parameter for sigmoidal mapping (dimensionless) */ - const double C1 = (3.14159265/sqrt(3)); + const double C1 = (M_PI/sqrt(3)); /* PSP rise time in ms^-1 */ const double gamma_e = 70E-3; - const double gamma_i = 100E-3; + const double gamma_g = 100E-3; + + /* Conductivities */ + /* Leak in aU */ + const double g_L = 1.; - /* Conductivities in mS/cm^-2 */ - /* Leak current */ - const double g_L_t = 1; - const double g_L_r = 1; + /* Synaptic conductivity in ms */ + const double g_AMPA = 1.; + const double g_GABA = 1.; - /* Potassium leak current */ - const double g_LK_t = 0.02; - const double g_LK_r = 0.02; + /* Potassium leak current in mS/m^2 */ + const double g_LK = 0.02; - /* T current */ + /* T current in mS/m^2 */ const double g_T_t = 3; const double g_T_r = 2.3; - /* h current */ - const double g_h = 0.05; + /* h current in mS/m^2 */ + const double g_h = 0.051; /* Reversal potentials in mV */ /* Synaptic */ @@ -190,9 +176,9 @@ class Thalamic_Column { const double E_h = -40; /* Calcium parameters */ - const double alpha_Ca = -52E-6; /* influx per spike in nmol */ + const double alpha_Ca = -51.8E-6; /* influx per spike in nmol */ const double tau_Ca = 10; /* calcium time constant in ms */ - const double Ca_0 = 2.4E-4; /* resting concentration */ + const double Ca_0 = 2.4E-4; /* resting concentration */ /* I_h activation parameters */ const double k1 = 2.5E7; @@ -204,15 +190,40 @@ class Thalamic_Column { /* Noise parameters in ms^-1 */ const double mphi = 0E-3; - const double dphi = 10E-3; + const double dphi = 0E-3; double input = 0.0; /* Connectivities (dimensionless) */ - const double N_tr = 4; - const double N_rt = 4; - const double N_rr = 20; + const double N_tr = 3; + const double N_rt = 5; + const double N_rr = 30; + + /* Parameters for SRK4 iteration */ + const vector A = {0.5, 0.5, 1.0, 1.0}; + const vector B = {0.75, 0.75, 0.0, 0.0}; + + /* Random number generators */ + vector MTRands; + + /* Container for noise */ + vector Rand_vars; - friend class Stim; + /* 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 */ + s_tt = _INIT(0.0), /* PostSP from TC population to TC population */ + s_tr = _INIT(0.0), /* PostSP from TC population to RE population */ + s_rt = _INIT(0.0), /* PostSP from RE population to TC population */ + s_rr = _INIT(0.0), /* PostSP from RE population to RE population */ + x_tt = _INIT(0.0), /* derivative of s_tt */ + x_tr = _INIT(0.0), /* derivative of s_tr */ + x_rt = _INIT(0.0), /* derivative of s_rt */ + x_rr = _INIT(0.0), /* derivative of s_rr */ + h_T_t = _INIT(0.0), /* inactivation of T channel */ + h_T_r = _INIT(0.0), /* inactivation of T channel */ + m_h = _INIT(0.0), /* activation of h channel */ + m_h2 = _INIT(0.0); /* activation of h channel bound with protein */ }; /****************************************************************************************************/ /* end */ diff --git a/Thalamus.cpp b/Thalamus.cpp index f355a25..e2156d8 100644 --- a/Thalamus.cpp +++ b/Thalamus.cpp @@ -1,44 +1,46 @@ /* -* 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. -*/ + * Copyright (c) 2015 University 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 + * 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. + * + * AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com + * + * Based on: A thalamocortical neural mass model of the EEG during NREM sleep and its response + * to auditory stimulation. + * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, + * JC Claussen. + * PLoS Computational Biology In Review (in review). + */ /****************************************************************************************************/ -/* Implementation of the simulation as MATLAB routine (mex compiler) */ -/* mex command is given by: */ -/* mex CXXFLAGS="\$CXXFLAGS -std=gnu++0x -fpermissive" Thalamus.cpp Thalamic_Column.cpp */ -/* The Simulation requires the following boost libraries: Preprocessor */ -/* Random */ +/* Main file for compilation tests */ /****************************************************************************************************/ -#include "mex.h" -#include "matrix.h" +#include +#include #include "Thalamic_Column.h" -#include "Stimulation.h" -#include "saves.h" /****************************************************************************************************/ /* Fixed simulation settings */ /****************************************************************************************************/ -extern const int onset = 15; /* time until data is stored in s */ +typedef std::chrono::high_resolution_clock::time_point timer; +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 */ /****************************************************************************************************/ @@ -47,72 +49,28 @@ extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */ /****************************************************************************************************/ -/* Constants for SRK4 iteration */ +/* Main simulation routine */ /****************************************************************************************************/ -extern const vector B1 = {0, - 0.626708569400000081728308032325, - 1.7296310295000001389098542858846, - 1.2703689705000000831347506391467}; -extern const vector B2 = {0, - 0.78000033203198970710445792065002, - 1.28727807507536762265942797967, - 0.44477273249350995909523476257164}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Simulation routine */ -/* lhs defines outputs */ -/* rhs defines inputs */ -/****************************************************************************************************/ -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - /* 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_Thalamus = mxGetPr (prhs[1]); /* Parameters of cortical module */ - double* var_stim = mxGetPr (prhs[2]); /* Parameters of stimulation protocol */ - +int main(void) { /* Initialize the populations */ - Thalamic_Column Thalamus(Param_Thalamus); - - /* Initialize the stimulation protocol */ - Stim Stimulation(Thalamus, var_stim); + Thalamic_Column Thalamus = Thalamic_Column(); - /* Create data containers */ - mxArray* Vt = SetMexArray(1, T*res/red); - mxArray* Vr = SetMexArray(1, T*res/red); - mxArray* I_T = SetMexArray(1, T*res/red); - mxArray* I_h = SetMexArray(1, T*res/red); - - /* Pointer to the actual data block */ - double* Pr_Vt = mxGetPr(Vt); - double* Pr_Vr = mxGetPr(Vr); - double* Pr_I_T = mxGetPr(I_T); - double* Pr_I_h = mxGetPr(I_h); + /* Take the time of the simulation */ + timer start,end; /* Simulation */ - int count = 0; - for (int t=0; t=onset*res && t%red==0){ - get_data(count, Thalamus, Pr_Vt, Pr_Vr, Pr_I_T, Pr_I_h); - ++count; - } } + end = std::chrono::high_resolution_clock::now(); - /* Output of the simulation */ - plhs[0] = Vt; - plhs[1] = Vr; - plhs[2] = I_T; - plhs[3] = I_h; -return; + /* Time consumed by the simulation */ + double dif = 1E-3*std::chrono::duration_cast( end - start ).count(); + std::cout << "simulation done!\n"; + std::cout << "took " << dif << " seconds" << "\n"; + std::cout << "end\n"; } /****************************************************************************************************/ -/* end */ +/* end */ /****************************************************************************************************/ diff --git a/Thalamus_mex.cpp b/Thalamus_mex.cpp new file mode 100644 index 0000000..ef858fb --- /dev/null +++ b/Thalamus_mex.cpp @@ -0,0 +1,114 @@ +/* + * Copyright (c) 2015 University 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 + * 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. + * + * AUTHORS: Michael Schellenberger Costa: mschellenbergercosta@gmail.com + * + * Based on: A thalamocortical neural mass model of the EEG during NREM sleep and its response + * to auditory stimulation. + * M Schellenberger Costa, A Weigenand, H-VV Ngo, L Marshall, J Born, T Martinetz, + * JC Claussen. + * PLoS Computational Biology In Review (in review). + */ + +/****************************************************************************************************/ +/* Implementation of the simulation as MATLAB routine (mex compiler) */ +/* mex command is given by: */ +/* mex CXXFLAGS="\$CXXFLAGS -std=c++11" Thalamus_mex.cpp Thalamic_Column.cpp */ +/* The Simulation requires the following boost libraries: Random */ +/****************************************************************************************************/ +#include "mex.h" +#include "matrix.h" +#include "Data_Storage.h" +mxArray* SetMexArray(int N, int M); + +/****************************************************************************************************/ +/* Fixed simulation settings */ +/****************************************************************************************************/ +extern const int onset = 15; /* time until data is stored in s */ +extern const int res = 1E3; /* 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 */ +/****************************************************************************************************/ + +/****************************************************************************************************/ +/* Simulation routine */ +/* lhs defines outputs */ +/* rhs defines inputs */ +/****************************************************************************************************/ +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { + /* 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_Thalamus = mxGetPr (prhs[1]); /* Parameters of cortical module */ + + /* Initialize the populations */ + Thalamic_Column Thalamus(Param_Thalamus); + + /* Create data containers */ + mxArray* Vt = SetMexArray(1, T*res/red); + mxArray* Vr = SetMexArray(1, T*res/red); + mxArray* ah = SetMexArray(1, T*res/red); + + /* Pointer to the actual data block */ + double* Pr_Vt = mxGetPr(Vt); + double* Pr_Vr = mxGetPr(Vr); + double* Pr_ah = mxGetPr(ah); + + /* Simulation */ + int count = 0; + for (int t=0; t=onset*res && t%red==0){ + get_data(count, Thalamus, Pr_Vt, Pr_Vr, Pr_ah); + ++count; + } + } + + /* Output of the simulation */ + plhs[0] = Vt; + plhs[1] = Vr; + plhs[2] = ah; +return; +} +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ + + +/****************************************************************************************************/ +/* Create MATLAB data container */ +/****************************************************************************************************/ +mxArray* SetMexArray(int N, int M) { + mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL); + mxSetM(Array, N); + mxSetN(Array, M); + mxSetData(Array, mxMalloc(sizeof(double)*M*N)); + return Array; +} +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ diff --git a/macros.h b/macros.h deleted file mode 100644 index 1baa7a1..0000000 --- a/macros.h +++ /dev/null @@ -1,81 +0,0 @@ -/* -* 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 */ -/****************************************************************************************************/ -#pragma once -#include -#include -#include -#include -#include -#include - - -/****************************************************************************************************/ -/* Macro for vector initialization */ -/****************************************************************************************************/ -#define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0} -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Macros for calculation of nth RK term */ -/****************************************************************************************************/ -#define _GET_VARNAMES(r, data, i, elem) BOOST_PP_COMMA_IF( i ) BOOST_PP_CAT(data, elem) -#define _SET_RK1(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0]); -#define _SET_RK2(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[1] * 0.5); -#define _SET_RK3(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[2] * 0.5); -#define _SET_RK4(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[3]); - -#define _SWITCH(...) \ - double BOOST_PP_SEQ_FOR_EACH_I(_GET_VARNAMES, var_, __VA_ARGS__); \ - switch(N) { \ - default: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK1, var_, __VA_ARGS__) \ - break; \ - case 2: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK2, var_, __VA_ARGS__) \ - break; \ - case 3: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK3, var_, __VA_ARGS__) \ - break; \ - case 4: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK4, var_, __VA_ARGS__) \ - break; \ - } -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Macros for repeated strings */ -/****************************************************************************************************/ -#define _REPEAT1(z, num, elem) BOOST_PP_COMMA_IF(num) elem -#define _REPEAT(elem, num) BOOST_PP_REPEAT(num , _REPEAT1 , elem) -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/saves.h b/saves.h deleted file mode 100644 index d089033..0000000 --- a/saves.h +++ /dev/null @@ -1,56 +0,0 @@ -/* -* 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 */ -/****************************************************************************************************/ -#pragma once -#include "Thalamic_Column.h" - -/****************************************************************************************************/ -/* Save data */ -/****************************************************************************************************/ -inline void get_data(int counter, Thalamic_Column& Col, double* Vt, double* Vr, double* I_T, - double* I_h) { - Vt [counter] = Col.Vt [0]; - Vr [counter] = Col.Vr [0]; - I_T [counter] = -Col.I_T_t (0); - I_h [counter] = -Col.I_h (0); -} -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Create MATLAB data container */ -/****************************************************************************************************/ -mxArray* SetMexArray(int N, int M) { - mxArray* Array = mxCreateDoubleMatrix(0, 0, mxREAL); - mxSetM(Array, N); - mxSetN(Array, M); - mxSetData(Array, mxMalloc(sizeof(double)*M*N)); - return Array; -} -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/