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<ref_period) = ref_period;
intervals = round(isi/dt);
if num_isis > 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