unit_ids = nwb.units.id.data.load(); % array of unit ids represented within this
% Initialize trials & times Map containers indexed by unit_ids
unit_trials = containers.Map('KeyType',class(unit_ids),'ValueType','any');
unit_times = containers.Map('KeyType',class(unit_ids),'ValueType','any');
for i = 1:length(unit_ids)
row = nwb.units.getRow(unit_id, 'useId', true, 'columns', {'spike_times', 'trials'});
unit_trials(unit_id) = row.trials{1};
unit_times(unit_id) = row.spike_times{1};
sorted_ids = sort(unit_ids);
'ind', true,... % mask into xs and ys for this photostim
'duration', 0,... % in seconds
'ramp_offset', 0); % in seconds
% Initialize Map container of plotting data for each unit, stored as structure
Unit = containers.Map('KeyType',class(unit_ids),'ValueType','any');
'left_scatter', false,...
'right_scatter', false,...
'photostim', Photostim); % can have multiple photostim
trial = nwb.intervals_trials.getRow(unit_trial_id, 'useId', true,...
'columns', {'PoleInTime', 'PoleOutTime', 'CueTime', 'GoodTrials'});
unit_sample = trial.PoleInTime;
unit_delay = trial.PoleOutTime;
unit_response = trial.CueTime;
unit_good_trials = trial.GoodTrials;
response_threshold = 0.43;
expected_delay_offset = 1.3; % determined from figure 1a
expected_response_offset = 1.3;
expected_delay = unit_sample + expected_delay_offset;
expected_response = unit_delay + expected_response_offset;
good_delay = (unit_delay > expected_delay - delay_threshold) &...
(unit_delay < expected_delay + delay_threshold);
good_response = (unit_response > expected_response - response_threshold) &...
(unit_response < expected_response + response_threshold);
avg_sample = mean(unit_sample(good_delay & good_response));
avg_delay = mean(unit_delay(good_delay & good_response));
avg_response = mean(unit_response(good_delay & good_response));
xs = unit_spike_time - trial.start_time - avg_response;
curr_unit.left_scatter = logical(trial.HitL);
curr_unit.right_scatter = logical(trial.HitR);
curr_unit.sample = avg_sample - avg_response;
curr_unit.delay = avg_delay - avg_response;
curr_unit.photostim.ind = unit_no_stim;
SampleStim.ind = unit_sample_stim;
SampleStim.period = 'Sample';
SampleStim.duration = 0.5;
SampleStim.ramp_offset = 0.1;
curr_unit.photostim(end+1) = SampleStim;
early_stim_types = unique(unit_stim_type(unit_early_stim));
for i_early_types=1:length(early_stim_types)
early_type = early_stim_types(i_early_types);
EarlyStim.period = 'Early Delay';
EarlyStim.ind = early_type == unit_stim_type & unit_early_stim;
EarlyStim.duration = 0.5;
EarlyStim.ramp_offset = 0.1;
EarlyStim.duration = 0.8;
EarlyStim.ramp_offset = 0.2;
curr_unit.photostim(end+1) = EarlyStim;
MiddleStim.ind = unit_middle_stim;
MiddleStim.period = 'Middle Delay';
MiddleStim.duration = 0.5;
MiddleStim.ramp_offset = 0.1;
curr_unit.photostim(end+1) = MiddleStim;
Unit(unit_id) = curr_unit;
neuron_labels = [2, 3]; % neuron labels from Figure 1e
neuron_ids = [11, 2]; % neuron unit IDs corresponding to the Fig 1e labels
num_conditions = 4; % photostim conditions: nostim, sample, early, middle if applicable
num_neurons = length(neuron_ids);
% Inititalize data structures for each summary plot of categorized neural spike data at specified stimulus condition
ConditionPlot = struct(...
'right_scatter', RasterPlot,...
'left_scatter', RasterPlot,...
% Plot neural spike data for each neuron and stimulus condition in a subplot array: num_neurons (rows) x num_conditions (columns)
Neuron = Unit(neuron_ids(nn));
% Initialize structure with neural + stimulus condition data
CurrPlot = ConditionPlot;
CurrPlot.xlim = [min(Neuron.xs) max(Neuron.xs)];
CurrPlot.sample = Neuron.sample;
CurrPlot.delay = Neuron.delay;
CurrPlot.response = Neuron.response;
% Plot each neuron/condition
plot_row = (nn - 1) * num_conditions;
ax = subplot(num_neurons, num_conditions, plot_row + cc, 'Parent', fig);
Stim = Neuron.photostim(cc);
CurrPlot.stim_type = Stim.period;
if strcmp(Stim.period, 'none')
CurrPlot.label = sprintf('Neuron %d', neuron_labels(nn));
CurrPlot.psth_bin_window = 9;
CurrPlot.label = Stim.period;
CurrPlot.psth_bin_window = 2;
stim_left_scatter_ind = Stim.ind & Neuron.left_scatter;
stim_left_scatter_trials = Neuron.ys(stim_left_scatter_ind);
CurrPlot.left_scatter.xs = Neuron.xs(stim_left_scatter_ind);
[~,CurrPlot.left_scatter.ys] = ismember(stim_left_scatter_trials,unique(stim_left_scatter_trials));
stim_right_scatter_ind = Stim.ind & Neuron.right_scatter;
stim_right_scatter_trials = Neuron.ys(stim_right_scatter_ind);
CurrPlot.right_scatter.xs = Neuron.xs(stim_right_scatter_ind);
[~,CurrPlot.right_scatter.ys] = ismember(stim_right_scatter_trials,unique(stim_right_scatter_trials));
plot_condition(ax, CurrPlot);
function plot_condition(ax, ConditionPlot)
left_cdata = [1 0 0]; % red
right_cdata = [0 0 1]; % blue
% moving average over 200 ms as per figure 1e
hist_bin_width = 0.2 / ConditionPlot.psth_bin_window;
[left_psth_xs, left_psth_ys] =...
calculate_psth(ConditionPlot.left_scatter.xs, ConditionPlot.psth_bin_window, hist_bin_width);
[right_psth_xs, right_psth_ys] =...
calculate_psth(ConditionPlot.right_scatter.xs, ConditionPlot.psth_bin_window, hist_bin_width);
right_scatter_offset = min(ConditionPlot.right_scatter.ys);
right_scatter_height = max(ConditionPlot.right_scatter.ys) - right_scatter_offset;
left_scatter_offset = min(ConditionPlot.left_scatter.ys);
left_scatter_height = max(ConditionPlot.left_scatter.ys) - left_scatter_offset;
psth_height = max([left_psth_ys right_psth_ys]);
left_y_offset = hist_margin...
right_y_offset = scatter_margin...
subplot_height = right_y_offset...
+ right_scatter_offset...
plot(ax, left_psth_xs, left_psth_ys, 'Color', left_cdata);
plot(ax, right_psth_xs, right_psth_ys, 'Color', right_cdata);
ConditionPlot.left_scatter.xs,...
left_y_offset + ConditionPlot.left_scatter.ys,...
ConditionPlot.right_scatter.xs,...
right_y_offset + ConditionPlot.right_scatter.ys,...
% sample, delay, response lines
line(ax, repmat(ConditionPlot.sample, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
line(ax, repmat(ConditionPlot.delay, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
line(ax, repmat(ConditionPlot.response, 1, 2), [0 subplot_height],...
'Color', 'k', 'LineStyle', '--');
% blue bar for photoinhibition period
if ~strcmp(ConditionPlot.stim_type, 'none')
stim_height = subplot_height;
stim_width = 0.5; % seconds
% end time relative to 'go' cue as described in the paper.
switch ConditionPlot.stim_type
error('Invalid photostim period `%s`', ConditionPlot.stim_type);
stim_offset = ConditionPlot.response - stim_width - end_offset;
stim_offset, stim_height;...
stim_offset+stim_width, stim_height;...
stim_offset+stim_width, 0];
'Vertices', patch_vertices,...
'FaceColor', '#B3D3EC',... % light blue shading
title(ax, ConditionPlot.label);
xlabel(ax, 'Time (Seconds)');
ylabel(ax, 'Spikes s^{-1}')
yticks(ax, [0 max(10, round(psth_height, -1))]);
% legend(ax, [scatter_left_plot, scatter_right_plot], {'Left Lick', 'Right Lick'},...
% 'location', 'northwestoutside');
ax.XLim = ConditionPlot.xlim;
ax.YLim = [0 subplot_height];