function [in_n,in_t] = BG_GHS_inputs(dt,time_steps,n_inputs,switches,hz,n_neurons,neurons_per_channel,n_channels,ref_period,type) % BG_GHS_INPUTS spike train inputs % BG_GHS_INPUTS(DT,T,I,S,HZ,N,C,NC,R,FLAG) where DT is the simulation time-step (in seconds), T is the number of time-steps in the simulationm % I is the number of inputs, S is the array of switching times, HZ is the matrix of firing frequencies (length(S) x NC), N is the number % of simulated input neurons, C is the number of neurons per channel, NC is the number of channels, and FLAG is one of: % 'vivo' - Poisson process, irregular tonic firing % 'organo' - organotypic-like correlated tonic firing, specify a single HZ value or just use first HZ value in matrix % 'slow' - anaesthesia-like slow-wave firing % % % Mark Humphries & Rob Stewart 23/12/2004 switch type case 'vivo' old = 0; in_n = []; in_t = []; % in vivo normal tonic - Poisson processes for i = 1:length(switches) h = hz(i,:); %channel frequencies for this switch point time_seconds = (switches(i)-old)*dt; num_isis = ceil(time_seconds * max(h) * 2); isi = []; for j = 1:n_channels isi = [isi,random('exp',1/h(j),num_isis,neurons_per_channel)]; end isi(isi 1 spike_times = cumsum(intervals)+old; %only on a matrix else spike_times = intervals+old; end spike_times = spike_times'; [t1,t2] = find(spike_times<(switches(i)+2)); %shifted forward two (20/02/04) if ~isempty(t1) ind = sub2ind(size(spike_times),t1,t2); times = spike_times(ind); [in_t_temp,ord] = sort(times); %sort into time order in_n_temp = t1(ord); in_t = [in_t;in_t_temp]; in_n = [in_n;in_n_temp]; end old = switches(i); end in_n = uint32(in_n+n_neurons-1); in_t = uint32(in_t-2); %not quite sure why this is two... case 'organo' % organotypic - correlated tonic % NOTE: in_t must be in time-order isi = 1/hz(1)/dt; % ISI in time-steps regular_times = round(isi:isi:time_steps); % template train spike-times % create jittered versions for each input reg_mat = repmat(regular_times,n_inputs,1); jit_std = isi * 0.3; % 30% jitter - reduces with *increasing* Hz jitter = round(randn(n_inputs,length(regular_times)) * jit_std); % normally-distributed tonic firing - a la SNc neuron trains = reg_mat + jitter; % jittered trains for each input [in_t ord] = sort(trains(:)); % sort in to time order event_idxs = (1:n_inputs)+n_neurons-1; % input indices are after all other indices (shift back by 1 for MEX indexing) events = repmat(event_idxs',1,length(regular_times)); events = events(:); in_n = events(ord); % put input indices in same order as time of spike % remove all out-of-range events out_of_range = find(in_t < 0 | in_t > time_steps); in_t(out_of_range) = []; in_n(out_of_range) = []; in_n = uint32(in_n); in_t = uint32(in_t); case 'slow' % in vivo anaesthetic - slow-wave, like slow-wave sleep, at ~1Hz % up-period: correlated firing at ~12 Hz (Steriade et al 2001) [KG: % shouldn't it be 24Hz because thisis the *mean* rate over both states?] % down-period: complete silence % desynchronised: during ECoG recording, wave spontaneously desynchonises, suggesting that underlying cortical firing % is no longer correlated; however, Kasanetz et al. (2002) data shows that striatal neurons remain in their up-state during % this period. Therefore, if up/down state transition does follow cortical slow-wave exactly, as seems to be the case, then % this suggests that desynchronised ECoG corresponds to decorrelated cortical firing at a frequency at least equivalent to the % slow-wave up-period % To model: alternate periods of silence and correlated firing (a la organotypic) at ~ 1Hz, randomly adding % Poisson-generated section of desynchronisation if requested by user time_seconds = time_steps * dt; half_second = 0.5 * time_steps / time_seconds; in_t = []; in_n = []; slow_jitter = 2.5; % jitter - reduces with *increasing* Hz (try 2.6 with 24Hz input) %%%%%% KG mod: vary intensity of each coherent up-state in cortex - No %%%%%% spikes in each upstate is same for all inputs within each period slow_std_hz = 0.15; % * 100 % std with normally distributed rate for loop = 1:time_seconds*2 % silent on odd, correlated on even (so starts on silence) if ~mod(loop,2) % is even - do half second of correlated firing %%%%%% KG mod %%%%%%%%%%%% hz_eff = hz(1) .* (1 + slow_std_hz .* randn); isi = 1/hz_eff/dt; % ISI in time-steps %%%%%%% KG mod end %%%%%%% regular_times = round(isi:isi:half_second) + half_second.*(loop-1); % template train spike-times for this second % create jittered versions for each input reg_mat = repmat(regular_times,n_inputs,1); jit_std = isi * slow_jitter; % 30% jitter - reduces with *increasing* Hz jitter = round(randn(n_inputs,length(regular_times)) * jit_std); % normally-distributed tonic firing - a la SNc neuron trains = reg_mat + jitter; % jittered trains for each input [temp ord] = sort(trains(:)); % sort in to time order in_t = [in_t; temp]; event_idxs = (1:n_inputs)+n_neurons-1; % input indices are after all other indices (shift back by 1 for MEX indexing) events = repmat(event_idxs',1,length(regular_times)); events = events(:); temp = events(ord); % put input indices in same order as time of spike in_n = [in_n; temp]; end end % remove all out-of-range events out_of_range = find(in_t < 0 | in_t > time_steps); in_t(out_of_range) = []; in_n(out_of_range) = []; % sort ready for coincidence removal [ts tn] = sort(in_t); in_t = ts; in_n = in_n(tn); % remove duplicate events No_ins = length(in_t); if No_ins > 1 No_coincidents = 0; coincident_indices = []; input_t1 = in_t(1); input_n1 = in_n(1); for j=2:No_ins input_t2 = in_t(j); input_n2 = in_n(j); if input_t2 == input_t1 & input_n1 == input_n2 No_coincidents = No_coincidents + 1; coincident_indices = [coincident_indices j]; end input_t1 = input_t2; input_n1 = input_n2; end in_t(coincident_indices) = []; in_n(coincident_indices) = []; fprintf(1, 'There were %d coincident events removed from inputs\n', No_coincidents); end in_n = uint32(in_n); in_t = uint32(in_t); end % switch