corticothalamic_spatial_integration/figs/figS5.m

clear % clear workspace
close all; % close all open figures
% load chR2V1_uInfo (Ntsr1 photoactivation)
load('/Users/busse/code/work/CTfeedbackSize/to_publish/data/chr2V1.mat')

%% plot parameters
mean_color = [218 165 35]/255; % color for mean
pinkEx = [255 20 147]/255; % color for example neurons
% percentage of layer thickness based on: Heumann D, Rabinowicz T (1977) Postnatal development of mouse cerebral neocortex. II Quantitative cytoarchitectonics of visual and auditory areas. J Hirnforsch 18: 483-500
%   Layers                  I    II/III    IV    V    VI    Total
%   Relative thickness (%)    13    26      13    24    24    100
% assumptions: thickness of V1 of 1 mm
% 0: base of layer 4
rel_layerDepth = [-520  -390  -130     0   240   480]; 
rel_layerCenters = [-455  -260   -65   120   360]; % add half of each layer thickness
%% S5b: example tuning curves


%% S5c: Ntsr1+ photoactivation across layers
rng(1);
layer_names = {'1', '2/3', '4', '5', '6'};
n_layers = numel(layer_names);

% compute means, CIs, SEMs across layers
layer_mean = nan(1, n_layers);
layer_sem = nan(1, n_layers);
layer_ci = nan(2, n_layers);
for ilayer = 1 : n_layers
    % find neurons in layer
    layer_log = strcmp({chr2V1_uInfo.layer}, layer_names(ilayer)); 
    % add 0.0001 to numerator and denominator to make log2 work for 0 firing rates
    layer_mean(ilayer) = mean(log2(([chr2V1_uInfo(layer_log).frL]+0.0001) ./ ...
        ([chr2V1_uInfo(layer_log).frC]+0.0001)));
    % compute confidence intervals via bootstrapping
    temp = bootstrp(100, @(x) {mean(log2(x))}, ...
        ([chr2V1_uInfo(layer_log).frL]+0.0001) ./([chr2V1_uInfo(layer_log).frC]+0.0001));    
    layer_ci(:, ilayer) = prctile([temp{:}], [2.5, 97.5]);
    layer_sem(ilayer) = var([temp{:}]);
    fprintf('Layer %s: mean %3.2f (%3.2f, %3.2f)\n', ...
        layer_names{ilayer}, layer_mean(ilayer), ...
        layer_ci(1, ilayer), layer_ci(2, ilayer));
end
fprintf('\n');

% make plot
xLim = 9;
layerY = [rel_layerDepth' rel_layerDepth']'/1000;
layerX = [repmat(-xLim,6,1) repmat(xLim,6,1)]';
modRats = log2(([chr2V1_uInfo.frL]+0.001) ./ ([chr2V1_uInfo.frC]+0.001));
modRats(modRats < -8) = -8.5;

figure('Name', 'Born et al. (2021), Fig. S5c'); hold on;
plot(modRats, [chr2V1_uInfo.relative_depth], 'ko', 'MarkerSize', 4);
plot(layer_mean, rel_layerCenters/1000, 'o', 'Color', mean_color, 'MarkerSize', 4);
line(layer_ci, [rel_layerCenters; rel_layerCenters]/1000, 'Color', mean_color)
% add layer boundaries
line(layerX, layerY, 'Color', 'k')
line([0 0], [-0.4 0.4], 'Color', 'k')
set(gca, 'YDir', 'reverse', 'XLim', [-xLim xLim], ...
    'YLim', [rel_layerDepth(1)/1000 rel_layerDepth(end)/1000])
xlabel('Fold change')
ylabel('Relative depth (mm)')

%% load uInfo for dLGN recordings
load('/Users/busse/code/work/CTfeedbackSize/to_publish/data/chr2LGN.mat');

%% 2e: example raster plots
%sponEx(1,1) = struct('mouse_id', 'Ntsr1-Cre_0039', 'mouse_counter', 91, 'series_num', 7, 'unit_id', 21);
%sponEx(1,2) = struct('mouse_id', 'Ntsr1-Cre_0039', 'mouse_counter', 91, 'series_num', 7, 'unit_id', 29);

exIdx = [43 48];
timelims{1} = [-0.5, 0.75];
timelims{2} = [-0.5, 1.325];

figure('Name', 'Born et al. (2021), Fig. S5e'); hold on;
% (1:31) trials for one particular experiment with 0.25 s light dur
plotRaster_pub(chr2LGN_uInfo(exIdx(1)).spikeTimes(1:31)', timelims{1});
line([0 0.25], [1 1], 'Color', 'k'); % light duration
set(gca, 'YColor', 'w')
ylabel('')
figure('Name', 'Born et al. (2021), Fig. S5f'); hold on;
% (32:100) trials for experiments with 0.5 s light dur
plotRaster_pub(chr2LGN_uInfo(exIdx(2)).spikeTimes(32:100)', timelims{2});
line([0 0.5], [1 1], 'Color', 'k');
set(gca, 'YColor', 'w')
ylabel('')

%% 2f: firing rate during spontaneous visual activity (corresponding to size 0)
% mean rate before onset of V1 Ntsr1+ stimulation and during V1 Ntsr1+ stimulation
chr2_meanBefore = arrayfun(@(u) nanmean(u.rateBefore), chr2LGN_uInfo);
chr2_meanDuring = arrayfun(@(u) nanmean(u.rateDuring), chr2LGN_uInfo);

% statistics: signrank treats nan as missing value (i.e. ignores it)
p_rates = signrank(chr2_meanBefore, chr2_meanDuring);
% change in rate, n
n_rates  = nnz(~isnan(chr2_meanBefore) & ~isnan(chr2_meanDuring));

fprintf('dLGN mean rate before %3.1f sp/s, during %3.1f sp/s V1 Ntsr1+ stimulation\n', ...
    mean(chr2_meanBefore), mean(chr2_meanDuring));
fprintf('\tp(signrank) = %0.2g (N = %d)\n', p_rates, n_rates);

% prepare for plotting
scale_out = 1.2; % outliers will be plotted at max+10%
ymin_r = 0;
ymaxVal_r = 20;
ymax_r = ymin_r + (ymaxVal_r-ymin_r)*scale_out;
chr2_meanBefore(chr2_meanBefore > ymaxVal_r) = ymax_r;
chr2_meanDuring(chr2_meanDuring > ymaxVal_r) = ymax_r;

% make the plot
figure('Name', 'Born et al. (2021), Fig. 3g'); hold on;
plot(chr2_meanBefore, chr2_meanDuring, 'ko', 'MarkerSize', 4);
hold on;
plot(chr2_meanBefore(exIdx), chr2_meanDuring(exIdx), 'ko', 'MarkerSize', 4, ...
    'MarkerEdgeColor', pinkEx, 'MarkerFaceColor', pinkEx);
plot(mean(chr2_meanBefore), mean(chr2_meanDuring), 'ko', 'MarkerSize', 4, ...
    'MarkerEdgeColor', mean_color, 'MarkerFaceColor', mean_color);
hline = refline(1,0); % line of unity slope
set(hline, 'color', 0.6*ones(1,3));
axis equal; axis square
set(gca, 'XLim', [0 ymax_r], 'YLim', [0 ymax_r], ...
    'XTick', [0:10:ymaxVal_r ymax_r], 'YTick', [0:10:ymaxVal_r ymax_r]);
xlabel('Control');
ylabel('L6 activation');
title('Firing rate (sp/s)');
%% 2g: burst ratios
% compute burst data

% compute burst tonic ratio before vs. during V1 Ntsr1+ stimulation
chr2_btBefore = arrayfun(@(u) nanmean(u.btRatioBefore), chr2LGN_uInfo);
chr2_btDuring = arrayfun(@(u) nanmean(u.btRatioDuring), chr2LGN_uInfo);
n_bts  = nnz(~isnan(chr2_btBefore) & ~isnan(chr2_btDuring));
pburst_pv = signrank(chr2_btBefore, chr2_btDuring); % treats nan as missign values

fprintf('\ndLGN burst ratio before and during V1 Ntsr1+ stimulation\n')
fprintf('\tmean before: %3.2f%%\n\tmean during: %3.2f%%\n\tsignrank p = %02g\n\tN = %d\n', ...
    nanmean(chr2_btBefore)*100, nanmean(chr2_btDuring)*100, ...
    pburst_pv, n_bts);

% compute burst lengths before vs. during Ntsr1+ photoactivation
chr2blenBefore = cell2mat(arrayfun(@(bt) cat(1, bt.blengthBefore{:}), ...
    chr2LGN_uInfo, 'uniformoutput', 0)');
chr2blenBefore(chr2blenBefore < 1) = []; % 0 burst length: tonic spikes
chr2blenDuring = cell2mat(arrayfun(@(bt) cat(1, bt.blengthDuring{:}), ...
    chr2LGN_uInfo, 'uniformoutput', 0)');
chr2blenDuring(chr2blenDuring < 1) = [];

[~,p_blen] = kstest2(chr2blenBefore, chr2blenDuring);
p_med_blen = ranksum(chr2blenBefore, chr2blenDuring);

fprintf('\ndLGN burst lengths before / during V1 Ntsr1+ stimulation\n')
fprintf('\tkstest: p = %0.2g\n\tranksum: p = %02g\n\tN = %d/%d\n', ...    
    p_blen, p_med_blen, numel(chr2blenBefore), numel(chr2blenDuring));

% prepare for plotting
btTicks = [10^-2 10^-1 10^0];
btTickLabels = [0.01 0.1 1];
% adjust for plotting
chr2_btBefore(chr2_btBefore == 0) = 10^-2.5;
chr2_btDuring(chr2_btDuring == 0) = 10^-2.5;

figure('Name', 'Born et al. (2021), Fig. S5h');
loglog(chr2_btBefore, chr2_btDuring, 'ko', 'MarkerSize', 4);
hold on;
loglog(chr2_btBefore(exIdx), chr2_btDuring(exIdx), 'ko', 'MarkerSize', 4, ...
    'MarkerEdgeColor', pinkEx, 'MarkerFaceColor', pinkEx);
loglog(nanmean(chr2_btBefore), nanmean(chr2_btDuring), 'ko', 'MarkerSize', 4, ...
    'MarkerEdgeColor', mean_color, 'MarkerFaceColor', mean_color);
hline = refline(1,0);
set(hline, 'color', 0.6*ones(1,3));
axis equal; axis square
set(gca, 'XLim', [10^-2.5 10^0], 'YLim', [10^-2.5 10^0], ...
    'XTick', btTicks, 'YTick', btTicks, ...
    'XTickLabel', btTickLabels, 'YTickLabel', btTickLabels, ...
    'XMinorTick','off', 'YMinorTick','off');
xlabel('Control');
ylabel('L6 activation');

% inset for distribution of burst lengths
figure('Name', 'Born et al. (2021), Fig. S5h inset'); hold on;
ecdf(chr2blenBefore)
ecdf(chr2blenDuring)
sh = get(gca, 'Children'); % to set the color of the individual ecdfs
set(sh(1), 'Color', [0 143 210]/255);
set(sh(2), 'Color', 'k');
ylabel('CDF')
set(gca, 'XLim', [2 inf], 'YTick', [0 1])
xlabel('sp/burst')