Added data creation routines and plotting functions.

Figures/Create_Data.m

@@ -0,0 +1,47 @@
+function Create_Data()
+% This function creates the model data depicted in Figure 4-8 of
+%
+% 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).
+
+% Move to source folder(assuming it contains the Figures folder
+cd ..;
+
+% Check if the executable exists and compile if needed
+if(exist('Thalamus_mex.mesa64', 'file')==0)
+    mex CXXFLAGS="\$CXXFLAGS -std=c++11 -O3" TC_mex.cpp Cortical_Column.cpp Thalamic_Column.cpp;
+end
+
+% Add the path to the simulation routine
+addpath(pwd);
+
+% Go back into figures folder
+cd Figures;
+
+% NOTE: The data routines utilize various functions from the fieldtrip
+% toolbox http://www.fieldtriptoolbox.org/
+% You have to make sure, that the path to the fieldtrip routines is known
+% Edit the below statement to fit to your system.
+if(isempty(strfind(path, '/Tools/fieldtrip/preproc')))
+    addpath('/Tools/fieldtrip/preproc');
+end
+
+% Import the original data from Ngo et al 2013
+Import_Data(1);
+Import_Data(2);
+Import_Data(3);
+
+% Time series
+Data_Time_Series(1);
+Data_Time_Series(2);
+
+% Extract KC/SO averages
+Data_SO_Average(1);
+Data_SO_Average(2);
+ 
+% Data of rhytmic two-click stimulation
+Data_ERP_N3();
+end

Figures/Data_ERP_N3.m

@@ -0,0 +1,71 @@
+function Data_ERP_N3()
+% This function creates the model data depicted in Figure 7 of
+%
+% 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).
+
+% To ensure availability of the simulation routine and the utility functions 
+% from fieldtrip it should be called from Create_Data.m
+
+% Simulation parameters
+load('Data/Parameter_N3');   
+            
+% Simulation duration
+T = 3600;         % duration of the simulation
+
+% stimulation parameters
+% first number is the mode of stimulation
+% 0 == none
+% 1 == semi-periodic
+% 2 == phase dependend
+    
+var_stim    = [ 2;          % mode of stimulation
+                70;         % strength of the stimulus              in Hz (spikes per second)
+                80;       	% duration of the stimulus              in ms
+                5;          % time between stimulation events       in s  (ISI)
+                0;          % range of ISI                          in s  [ISI-range,ISI+range]  
+                2;          % Number of stimuli per event
+                1050;        % time between stimuli within a event   in ms         
+                450];       % time until stimuli after minimum      in ms
+ 
+
+% Model Output is given as: 
+% 1. Pyramidal membrane voltage         in mV
+% 2. Thalamic relay membrane voltage    in mV
+% 3. Thalamic calcium concentration     in mu mol
+% 4. Thalamic I_h activation            unitless
+% 5. Stimulation markers                in sampling rate
+[Vp, Vt, Ca, ah, Marker_Stim] = TC_mex(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim);
+
+% Power of spindle bands (BP filtered)
+Fs          = length(Vp)/T;
+Vp_FSP      = ft_preproc_hilbert(ft_preproc_bandpassfilter(Vp, Fs, [12, 15], 513, 'fir'), 'abs').^2;
+Vp_SSP      = ft_preproc_hilbert(ft_preproc_bandpassfilter(Vp, Fs, [09, 12], 513, 'fir'), 'abs').^2;
+
+% Check whether stimuli are too early
+Range_ERP   = [-1, 3];
+    
+Marker_Stim = Marker_Stim(Marker_Stim>-Range_ERP(1)*Fs);
+Marker_Stim = Marker_Stim(Marker_Stim< (T-Range_ERP(2))*Fs);
+
+% Define the matrices
+N_Stim      = length(Marker_Stim);
+time_event  = linspace(Range_ERP(1), Range_ERP(2), (Range_ERP(2)-Range_ERP(1))*Fs+1);
+Events      = zeros(length(time_event), N_Stim);
+Events_T    = zeros(length(time_event), N_Stim);
+Events_FSP  = zeros(length(time_event), N_Stim);
+Events_SSP  = zeros(length(time_event), N_Stim);
+
+for i=1:N_Stim
+    Events(:,i)     =  Vp((Marker_Stim(i)+Range_ERP(1)*Fs)+1:(Marker_Stim(i)+Range_ERP(2)*Fs+1));
+    Events_T(:,i)   =  Vt((Marker_Stim(i)+Range_ERP(1)*Fs)+1:(Marker_Stim(i)+Range_ERP(2)*Fs+1));
+    Events_FSP(:,i) =  Vp_FSP((Marker_Stim(i)+Range_ERP(1)*Fs)+1:(Marker_Stim(i)+Range_ERP(2)*Fs+1));
+    Events_SSP(:,i) =  Vp_SSP((Marker_Stim(i)+Range_ERP(1)*Fs)+1:(Marker_Stim(i)+Range_ERP(2)*Fs+1));
+end
+
+save('Data/ERP_Stim_Model', 'Events', 'Events_T', 'Events_FSP', 'Events_SSP', 'Marker_Stim');
+save('Data/Ca_Depletion', 'Vp', 'Vt', 'Ca', 'ah', 'Marker_Stim');
+end

Figures/Data_SO_Average.m

@@ -0,0 +1,73 @@
+function Data_SO_Average(type)
+% This function creates the model data depicted in Figure 6 of
+%
+% 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).
+
+% To ensure availability of the simulation routine and the utility functions 
+% from fieldtrip it should be called from Create_Data.m
+if nargin==0    
+    type = 2;
+end
+
+switch type
+    case 1;
+        load_model = 'Data/Time_Series_N2.mat';
+        load_data  = 'Data/KC_Average_data.mat';
+        savename   = 'Data/KC_Average.mat';
+    case 2;
+        load_model = 'Data/Time_Series_N3.mat';
+        load_data  = 'Data/SO_Average_data.mat';
+        savename   = 'Data/SO_Average.mat';
+end
+
+% Load data
+load(load_model, 'Vp');
+load(load_data);
+
+% Rename data
+mean_ERP_data = mean_ERP_sham;
+mean_FSP_data = mean_FSP_sham;
+sem_ERP_data  = sem_ERP_sham;
+sem_FSP_data  = sem_FSP_sham;
+
+% Process model Data
+T       = 3600;
+Fs      = length(Vp)/T;
+
+% Power of fast spindle band (BP filtered) and slow wave activity for peak
+% detection
+Vp_low      = ft_preproc_bandpassfilter(Vp, Fs, [0.25,4], 513,  'fir') + mean(Vp);
+Vp_FSP      = ft_preproc_hilbert(ft_preproc_bandpassfilter(Vp, Fs, [12, 15], 513, 'fir'), 'abs').^2;
+
+% Search for peaks
+[~, x_SO]   = findpeaks(-Vp_low, 'MINPEAKHEIGHT', 68, 'MINPEAKDISTANCE', 0.2*Fs);
+
+% Remove those events, that are too close to begin/end
+x_SO        = x_SO(x_SO<(T-2)*Fs); 
+x_SO        = x_SO(x_SO>    2*Fs); 
+
+% Set the variables 
+N_Stim      = length(x_SO);
+Range_ERP   = [-1.25, 1.25];
+Events      = zeros(length(time_events), N_Stim);
+Events_FSP  = zeros(length(time_events), N_Stim);
+
+% Segmentation
+for i=1:N_Stim
+    Events(:,i)     =  Vp    ((x_SO(i)+Range_ERP(1)*Fs)+1:(x_SO(i)+Range_ERP(2)*Fs+1));
+    Events_FSP(:,i) =  Vp_FSP((x_SO(i)+Range_ERP(1)*Fs)+1:(x_SO(i)+Range_ERP(2)*Fs+1));
+end
+
+% Averaging
+mean_ERP_model= mean(Events,    2); %#ok<*NASGU>
+mean_FSP_model= mean(Events_FSP,2);
+sd_ERP_model  = std (Events,    0, 2);
+sd_FSP_model  = std (Events_FSP,0, 2);
+
+% Save the data
+save(savename, 'N_Stim', 'time_events', 'mean_ERP_data', 'mean_FSP_data', 'sem_ERP_data', 'sem_FSP_data', 'mean_ERP_model', 'mean_FSP_model', 'sd_ERP_model', 'sd_FSP_model');
+end

Figures/Data_Time_Series.m

@@ -0,0 +1,40 @@
+function Data_Time_Series(type)
+% This function creates the model data depicted in Figure 4-8 of
+%
+% 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).
+
+% To ensure availability of the simulation routine and the utility functions 
+% from fieldtrip it should be called from Create_Data.m
+
+% Load the parameter settings
+if type == 1    
+    load('Data/Parameter_N2');           
+    Protocol_Name = 'N2';                              
+else    
+    load('Data/Parameter_N3');       
+    Protocol_Name = 'N3';                   
+end   
+ 
+% there is no stimulation so set var_stim to zero
+var_stim            = zeros(8,1);
+
+% Duration of the simulation
+T                   = 3600;         
+
+% Run the simulation (Can take some time)
+[Vp, Vt, Ca, ah, ~] = TC_mex(T, Param_Cortex, Param_Thalamus, Connectivity, var_stim);
+
+% Save the full data
+save(['Data/Time_Series_',Protocol_Name], 'Vp', 'Vt', 'Ca', 'ah');
+
+% Save a smaller snipplet for example time series plot
+Vp = Vp(1:3000);
+Vt = Vt(1:3000);
+Ca = Ca(1:3000);
+ah = ah(1:3000);
+save(['Data/Time_Series_Short_',Protocol_Name], 'Vp', 'Vt', 'Ca', 'ah');
+end

Figures/Import_Data.m

@@ -0,0 +1,44 @@
+function Import_Data(Type)
+% This function imports the experimental data from Ngo et al 2013 utilized
+% in
+%
+% 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).
+
+if nargin ==0
+    Type = 1;
+end
+
+load('Data/Orig/Experimental_Data');
+
+switch (Type)
+    case 1; 
+        Data = data.KC_Average; 
+        sn = 'Data/KC_Average_data';    
+    case 2; 
+        Data = data.SO_Average; 
+        sn = 'Data/SO_Average_data';   
+    case 3; 
+        Data = data.ERP_Average; 
+        sn = 'Data/ERP_Average_data';   
+end
+
+time_events     = Data(:,1); %#ok<*NASGU>
+
+mean_ERP        = Data(:,2);
+mean_ERP_sham   = Data(:,3);
+
+sem_ERP         = Data(:,4);
+sem_ERP_sham    = Data(:,5);
+
+mean_FSP        = Data(:,6);
+mean_FSP_sham   = Data(:,7);
+
+sem_FSP         = Data(:,8);
+sem_FSP_sham    = Data(:,9);
+
+save(sn, 'time_events', 'mean_ERP', 'mean_ERP_sham', 'sem_ERP', 'sem_ERP_sham', 'mean_FSP', 'mean_FSP_sham', 'sem_FSP', 'sem_FSP_sham');
+end

Figures/Pictogramm_TC.tex

@@ -0,0 +1,64 @@
+% command to run this on the terminal is lualatex -shell-escape Pictogram_Sleep_Regulation.tex
+% requires pgfplots
+
+\pdfpxdimen=1in
+\divide\pdfpxdimen by 600
+
+\documentclass[10pt]{standalone}
+\usepackage{tikz}
+%has to be compiled with ``lualatex -shell-escape Pictogram_TC.tex''
+\usetikzlibrary{external, positioning, backgrounds, arrows.meta,decorations.pathmorphing}
+
+% set up externalization
+%\tikzexternalize[prefix=tikz/]
+% WARNING the name of the Image CANNOT be the same as the tex file
+\tikzsetnextfilename{Pictogram_TC}
+\tikzset{external/system call={lualatex \tikzexternalcheckshellescape -halt-on-error
+-interaction=batchmode -jobname "\image" "\texsource";
+convert -density 600 -compress LZW "\image".pdf -flatten -alpha off -colorspace RGB -depth 8 "\image".tif;
+convert -density 600 "\image".tif  "\image".jpg;}}
+\tikzexternalize
+
+\begin{document}
+\begin{tikzpicture}[
+  s/.style  = {shorten >=1pt},
+  every loop/.style={shorten >=1pt},
+  pop/.style	 = {rectangle, draw, text width=10em, text centered, minimum height=3em, rounded corners=3mm},
+  noise_C/.style = {decorate, decoration={zigzag}, s,-{Circle[fill=blue!30]}},
+  noise_T/.style = {decorate, decoration={zigzag}, s,-{Circle[fill=green!30]}},
+  Ex/.style = {rounded corners, s, -{Circle[fill=blue!30]}},
+  Tc/.style = {rounded corners, s, -{Circle[fill=green!30]}},
+  In/.style = {rounded corners, s, -|},
+  ]
+
+  % Cortex populations
+  \node (ex)   at (0,0)	 [pop, fill=blue!30]	{Pyramidal (p)};
+  \node (in)   at (6,0)	 [pop, fill=yellow!30]	{Inhibitory (i)};
+  % Thalamus Populations
+  \node (tc)   at (6,-2) [pop, fill=green!30]	{Thalamocortical (t)};
+  \node (re)   at (0,-2) [pop, fill=yellow!30]	{Thalamic Reticular (r)};
+
+  % Connections within cortex
+  \draw
+  (		   ex)	    edge [Ex, loop above] node[below left, xshift = -0.4em]	{$N_{pp}$}	(		ex)
+  ([yshift= 0.2cm] ex.east) edge [Ex]		  node[above] 				{$N_{ip}$}	([yshift= 0.2cm]in.west)
+  ([yshift=-0.2cm] in.west) edge [In] 	    	  node[below] 				{$N_{pi}$}	([yshift=-0.2cm]ex.east)
+  (		   in)	    edge [In, loop above] node[below left, xshift = -0.4em]	{$N_{ii}$}	(		in)
+  (0.5,1.22)  		    edge [noise_C] 	  node[right] 				{$\phi_n$}	([xshift= 0.5cm]ex.north)
+  (6.5,1.22)  		    edge [noise_C] 	  node[right] 				{$\phi^{'}_n$}	([xshift= 0.5cm]in.north);
+
+  % Connections within thalamus
+  \draw
+  (6.5,-3.22)               edge [noise_T]	  node[right] 				{$\phi^{''}_{n}$}([xshift=0.5cm]  tc.south)
+  ([yshift= 0.15cm]tc.west) edge [Tc]	 	  node[above]				{$N_{rt}$}	([yshift= 0.15cm] re.east)
+  ([yshift=-0.15cm]re.east) edge [In]	 	  node[below]				{$N_{tr}$}	([yshift=-0.15cm] tc.west)
+  (                re)	    edge [In, loop below] node[above left, xshift = -0.4em]	{$N_{rr}$}      (                 re);
+
+  % TC connections
+  \path
+  (ex.south) 	  	    edge [Ex] 		  node[left]				{$N_{rp}$}	([yshift= 0.02cm] re.north)
+  (ex.south east) 	    edge [Ex, bend right=15, yshift =-0.2em] node[at start,below]{$N_{tp}$}	(tc.north west)
+  (tc.north) 	  	    edge [Tc]		  node[right]				{$N_{it}$}	([yshift=-0.02cm] in.south)
+  (tc.north west) 	    edge [Tc, bend right=15, yshift = 0.2em] node[at start,above]{$N_{pt}$}	(ex.south east);
+\end{tikzpicture}
+\end{document}