Code cleanups.

Moved all noise generation into the module. Simplified most of the code.
miscco · 5a9762afd96d13a1e28a4daf31197f355bf6348e · miscco committed Mar 10, 2014 · Mar 10, 2014 · 5a9762a
1 parent 30534a9
commit 5a9762afd96d13a1e28a4daf31197f355bf6348e
Unified Split

Showing with 434 additions and 297 deletions.

+21 −13 Main.cpp

+26 −19 ODE.h

+11 −36 Plots.m

+95 −0 Stimulation.h

+92 −86 Thalamic_Column.cpp

+47 −39 Thalamic_Column.h

+44 −18 Thalamus.cpp

+22 −18 macros.h

+57 −57 parameters.h

+19 −11 saves.h
diff --git a/Main.cpp b/Main.cpp
@@ -1,34 +1,39 @@
-/*****************************************************************************************************/
-/**********************************		file for code verification	**********************************/
-/*****************************************************************************************************/
+/****************************************************************************************************/
+/*		Main file for compilation tests																*/
+/*		The Simulation requires the following boost libraries:	Preprocessor						*/
+/*																Random								*/
+/****************************************************************************************************/
 #include <iostream>
 #include <ctime>
-#include <vector>
 #include "Thalamic_Column.h"
 #include "ODE.h"
 
-using std::vector;
-
-extern const int T 		= 50;
+/****************************************************************************************************/
+/*										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);
+/****************************************************************************************************/
+/*										 		end			 										*/
+/****************************************************************************************************/
 
-// simulation of the thalamic model
-int main(void) {
-	// Initializing the seeder.
-	srand(time(0));
 
+/****************************************************************************************************/
+/*										Main simulation routine										*/
+/****************************************************************************************************/
+int main(void) {
 	// Initializing the populations;
-	Thalamic_Column Col;
+	Thalamic_Column Thalamus;
 
 	// takes the time of the simulation
 	time_t start,end;
 	time (&start);
 
 	// simulation
 	for (int t=0; t< T*res; ++t) {
-		ODE (Col);
+		ODE (Thalamus);
 	}
 	time (&end);
 
@@ -38,3 +43,6 @@ int main(void) {
 	std::cout << "took " << dif 	<< " seconds" << "\n";
 	std::cout << "end\n";
 }
+/****************************************************************************************************/
+/*										 		end			 										*/
+/****************************************************************************************************/
diff --git a/ODE.h b/ODE.h
@@ -1,13 +1,28 @@
-/*****************************************************************************************************/
-/***********************			ODE function of the system			******************************/
-/*****************************************************************************************************/
-#include <vector>
+/****************************************************************************************************/
+/*									Implementation of the ODE solver								*/
+/****************************************************************************************************/
+#pragma once
 #include "Thalamic_Column.h"
 
-/*****************************************************************************************************/
-/*****************************		parameters for integration of SRK4		**************************/
-/*****************************************************************************************************/
-// vectors needed for stochastic RK
+/****************************************************************************************************/
+/*										Evaluation of SRK4											*/
+/****************************************************************************************************/
+inline void ODE(Thalamic_Column& Thalamus) {
+	// first calculating every ith RK moment
+	// has to be in order, 1th moment first
+	for (int i=1; i<=4; ++i) {
+		Thalamus.set_RK(i);
+	}
+	Thalamus.add_RK();
+}
+/****************************************************************************************************/
+/*										 		end													*/
+/****************************************************************************************************/
+
+
+/****************************************************************************************************/
+/*									Parameters for SRK4	integration									*/
+/****************************************************************************************************/
 extern const vector<double> B1 = {0,
 								  0.626708569400000081728308032325,
 								  1.7296310295000001389098542858846,
@@ -16,14 +31,6 @@ extern const vector<double> B2 = {0,
 								  0.78000033203198970710445792065002,
 								  1.28727807507536762265942797967,
 								  0.44477273249350995909523476257164};
-/*****************************************************************************************************/
-/**********************************		 		end			 		**********************************/
-/*****************************************************************************************************/
-
-// function that evaluates ODE using stochastic Runge Kutta
-inline void ODE(Thalamic_Column& Col) {
-	for (auto i=1; i<=4; ++i) {
-		Col.set_RK(i);
-	}
-	Col.add_RK();
-}
+/****************************************************************************************************/
+/*										 		end													*/
+/****************************************************************************************************/
diff --git a/Plots.m b/Plots.m
@@ -4,14 +4,18 @@
 function Plots(T, onset)
 
 if nargin == 0
-    Con     = [10;             % N_tr
-                   20;             % N_rt
-                   40];            % N_rr    
-    T       = 30;      % duration of the simulation
-    onset   = 10;
+    Con     	= [10;		% N_tr
+               	   20;		% N_rt
+                   40];		% N_rr    
+
+    var_stim    = [ 0;          % strength of the stimulus in Hz (spikes per second)
+                    0;          % time between   stimuli in s    
+                    0;          % time until first stimuli in s
+                    0];		% duration of the stimulus in ms
+    T       	= 30;  		% duration of the simulation
 end
 
-[Vt, Vr, Cat, Car, Phi_tr, Phi_rt, I_T_t, I_T_r, I_KCa, I_CAN, I_h] = Thalamus(Con, T, onset);
+[Vt, Vr] = Thalamus(T, Con, var_stim);
 
 L        = max(size(Vt));
 timeaxis = linspace(0,T,L);
@@ -22,35 +26,6 @@ function Plots(T, onset)
 subplot(212), plot(timeaxis,Vr)
 title('Inhibitory membrane voltage'), xlabel('time in s'), ylabel('Vr in mV')
 
-% figure(1)
-% subplot(411), plot(timeaxis, Vt)
-% title('Thalamic relay membrane voltage'),  xlabel('time in s'), ylabel('V_{t} in mV')
-% subplot(412), plot(timeaxis, Vr)
-% title('Thalamic reticular membrane voltage'),  xlabel('time in s'), ylabel('V_{r} in mV')
-% subplot(413), plot(timeaxis, Cat)
-% title('TC Ca concentration'),                 xlabel('time in s'), ylabel('Ca in \muM')
-% subplot(414), plot(timeaxis, Car)
-% title('RE Ca concentration'),                 xlabel('time in s'), ylabel('Ca in \muM')
-
-%exportfig(gcf, 'TC_delta.png', 'Format', 'png', 'height', 11, 'Color', 'rgb', 'FontMode', 'fixed', 'FontSize', 22)
-% 
-% figure(2)
-% subplot(411), plot(timeaxis,I_T_r)
-% title('RE I_{T} current '),   xlabel('time in s'), ylabel('I_{T}   in \muA cm^{-2}')
-% subplot(412), plot(timeaxis,I_CAN)
-% title('RE I_{CAN} current '), xlabel('time in s'), ylabel('I_{CAN} in \muA cm^{-2}')
-% subplot(413), plot(timeaxis,I_KCa)
-% title('RE I_{KCa} current '), xlabel('time in s'), ylabel('I_{KCa} in \muA cm^{-2}')
-% subplot(414), plot(timeaxis,Phi_tr)
-% title('TC input to RE'),      xlabel('time in s'), ylabel('\Phi_{tr} in \muA cm^{-2}')
-% 
-% figure(4)
-% subplot(311), plot(timeaxis,I_T_t)
-% title('TC I_{T} current '), xlabel('time in s'), ylabel('I_{T} in \muA cm^{-2}')
-% subplot(312), plot(timeaxis,I_h)
-% title('TC I_{h} current '), xlabel('time in s'), ylabel('I_{h} in \muA cm^{-2}')
-% subplot(313), plot(timeaxis,Phi_rt)
-% title('RE input to TC'),    xlabel('time in s'), ylabel('\Phi_{rt} in \muA cm^{-2}')
 % 
 fs      = L/T;
 NFFT    = 2^nextpow2(L);
@@ -62,4 +37,4 @@ function Plots(T, onset)
 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')
+%save('Thalamus.mat','Vt','Vr')
diff --git a/Stimulation.h b/Stimulation.h
@@ -0,0 +1,95 @@
+/****************************************************************************************************/
+/*								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& C, double* var)
+	{Thalamus = &T;
+	 setup(var);}
+
+	// scaling from SI to simulation variables s -> ms
+	void setup		(double* var_stim) {
+		extern const int onset;
+		extern const int res;
+		extern const int dt;
+
+		// scale the stimulation strength with respect to ms^-1
+		strength = var_stim[0] / 1000;
+
+		// scale the stimulation variables with respect to simulation resolution
+		ISI 	 = (int) var_stim[1] * res;
+
+		// stimulation starts after the onset
+		start 	 = (int)(var_stim[2] + onset) *res;
+
+		// rescale duration with respect to dt
+		duration = (int) var_stim[3]/dt;
+	}
+
+	void check_stim	(int time) {
+		// check whether a stimulation should start
+		// the duration of the stimulation is ignored
+		if(time==(start + count_stim*ISI)){
+			// turn the stimulation on
+			mode = 1;
+			Thalamus->set_input(strength);
+
+		}
+
+		// check whether a stimulation should end
+		if(mode ==1 && count_dur ==duration) {
+			// turn off the stimulation
+			mode = 0;
+			Thalamus->set_input(0.0);
+
+			// add counter for stimulation occurrence
+			count_stim++ ;
+
+			// reset the stimulation counter
+			count_dur = 0;
+		}
+
+		// if stimulation is on track count its duration
+		if(mode==1){
+			count_dur++;
+		}
+	}
+
+private:
+
+	// stimulation strength
+	double strength = 0.0;
+
+	// inter stimulus intervall
+	int ISI = 0;
+
+	// onset until stimulation starts
+	int start = 0;
+
+	// duration of the stimulation
+	int duration = 0;
+
+	// counter for stimulation events
+	int count_stim = 0;
+
+	// counter for stimulation length
+	int count_dur  = 0;
+
+	// Simulation on for TRUE and off for FALSE
+	bool mode = 0;
+
+	// Pointer to thalamic module
+	Thalamic_Column* Thalamus;
+};
+/****************************************************************************************************/
+/*										 		end													*/
+/****************************************************************************************************/