corticothalamic_spatial_integration/figs/fig2_BtoG.m

clear % clear workspace
close all; % close all open figures
% load pvV1_uInfo
load('/Users/busse/code/work/CTfeedbackSize/to_publish/data/pvV1.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
%% 2b: example tuning curves


%% 2c: V1 suppression 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
    layer_log = strcmp({pvV1_uInfo.layer}, layer_names(ilayer)) & ...
        ~[pvV1_uInfo.put_pv]; % find neurons in layer, but exclude the putative PVs
    % add 0.0001 to numerator and denominator to make log2 work for 0 firing rates
    layer_mean(ilayer) = mean(log2(([pvV1_uInfo(layer_log).frL]+0.0001) ./ ...
        ([pvV1_uInfo(layer_log).frC]+0.0001)));
    % compute confidence intervals via bootstrapping
    temp = bootstrp(100, @(x) {mean(log2(x))}, ...
        ([pvV1_uInfo(layer_log).frL]+0.0001) ./([pvV1_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(([pvV1_uInfo.frL]+0.001) ./ ([pvV1_uInfo.frC]+0.001));
modRats(modRats < -8) = -8.5;

figure('Name', 'Born et al. (2021), Fig. 3c'); hold on;
plot(modRats, [pvV1_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/pvdLGN.mat');

%% 2e: example raster plots

exIdx = [96 216];
timelims{1} = [-0.5, 0.75];
timelims{2} = [-0.5, 1.325];

figure('Name', 'Born et al. (2021), Fig. 3e_1'); hold on;
plotRaster_pub(pvLGN_uInfo(exIdx(1)).spikeTimes', timelims{1});
line([0 0.25], [1 1], 'Color', 'k');
set(gca, 'YColor', 'w')
ylabel('')
figure('Name', 'Born et al. (2021), Fig. 3e_2'); hold on;
plotRaster_pub(pvLGN_uInfo(exIdx(2)).spikeTimes', timelims{2});
line([0 1], [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 suppression and during V1 suppression
pv_meanBefore = arrayfun(@(u) nanmean(u.rateBefore), pvLGN_uInfo);
pv_meanDuring = arrayfun(@(u) nanmean(u.rateDuring), pvLGN_uInfo);

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

fprintf('dLGN mean rate before %3.1f sp/s, during %3.1f sp/s V1 suppression\n', ...
    mean(pv_meanBefore), mean(pv_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;
pv_meanBefore(pv_meanBefore > ymaxVal_r) = ymax_r;
pv_meanDuring(pv_meanDuring > ymaxVal_r) = ymax_r;

% make the plot
figure('Name', 'Born et al. (2021), Fig. 3f'); hold on;
plot(pv_meanBefore, pv_meanDuring, 'ko', 'MarkerSize', 4);
hold on;
plot(pv_meanBefore(exIdx), pv_meanDuring(exIdx), 'ko', 'MarkerSize', 4, ...
    'MarkerEdgeColor', pinkEx, 'MarkerFaceColor', pinkEx);
plot(mean(pv_meanBefore), mean(pv_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('V1 suppression');
title('Firing rate (sp/s)');
%% 2g: burst ratios
% compute burst data

% compute burst tonic ratio before vs. during V1 suppression
pv_btBefore = arrayfun(@(u) nanmean(u.btRatioBefore), pvLGN_uInfo);
pv_btDuring = arrayfun(@(u) nanmean(u.btRatioDuring), pvLGN_uInfo);
n_bts  = nnz(~isnan(pv_btBefore) & ~isnan(pv_btDuring));
pburst_pv = signrank(pv_btBefore, pv_btDuring); % treats nan as missign values

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

% compute burst lengths before vs. during V1 suppression
pvblenBefore = cell2mat(arrayfun(@(bt) cat(1, bt.blengthBefore{:}), ...
    pvLGN_uInfo, 'uniformoutput', 0)');
pvblenBefore(pvblenBefore < 1) = []; % 0 burst length: tonic spikes
pvblenDuring = cell2mat(arrayfun(@(bt) cat(1, bt.blengthDuring{:}), ...
    pvLGN_uInfo, 'uniformoutput', 0)');
pvblenDuring(pvblenDuring < 1) = [];

[~,p_blen] = kstest2(pvblenBefore, pvblenDuring);
p_med_blen = ranksum(pvblenBefore, pvblenDuring);

fprintf('\ndLGN burst lengths before / during V1 suppression\n')
fprintf('\tmedian before: %d\n\tmedian during: %d\n\tkstest: p = %0.2g\n\tranksum: p = %02g\n\tN = %d/%d\n', ...
    median(pvblenBefore), median(pvblenDuring), ...
    p_blen, p_med_blen, numel(pvblenBefore), numel(pvblenDuring));

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

figure('Name', 'Born et al. (2021), Fig. 3g');
loglog(pv_btBefore, pv_btDuring, 'ko', 'MarkerSize', 4);
hold on;
loglog(pv_btBefore(exIdx), pv_btDuring(exIdx), 'ko', 'MarkerSize', 4, ...
    'MarkerEdgeColor', pinkEx, 'MarkerFaceColor', pinkEx);
loglog(nanmean(pv_btBefore), nanmean(pv_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('V1 suppression');

% inset for distribution of burst lengths
figure('Name', 'Born et al. (2021), Fig. 3g inset'); hold on;
ecdf(pvblenBefore)
ecdf(pvblenDuring)
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')