OpponencyModels/RunFigure6Models.m at master · ClarkLabCode/OpponencyModels · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
function RunFigure6Models()
%Run models from Figure 6 and S6 - JZV

%% Set parameters

% Add utils folder to MATLAB path
addpath('utils');

% Set parameters for models
[ params ] = SetFigure6ModelParameters();

% For TF sweep:
% tf = 2.^(-3:0.5:7)';

% For timetraces:
tf = 1;

% Set the contrasts to be tested
cont = [1/4, 1/2];

% Define stimulus names
numStim = 6;
inputNames = {'PD','ND','PD+ND','PD+OD', 'CISSumPD', 'CISSumND'};

% Define model names
modelNames = {'Rectified multiplier', 'Dynamic gain nonlinearity', 'Conductance nonlinearity', 'Sigmoidal LN', 'Half-quadratic LN'};

%% Run the models

% Get the number of TFs, the number of contrasts, and the number of models
numTf = length(tf);
numCont = length(cont);
numModel = length(modelNames);

% Allocate containers
meanResp = cell(numTf,numCont);
voltResp = cell(numTf,numCont);
modelResp = cell(numTf,numCont);
numResp = nan(numTf, numCont, numStim);
denResp = nan(numTf, numCont, numStim);

% Iterate over temporal frequencies
for tfInd = 1:numTf

    % Iterate over contrasts
    for contInd = 1:numCont
        tic;

        % Make filters
        [ filters ] = MakeFigure6Filters(params, tf(tfInd));

        % Run models
        [ meanResp{tfInd, contInd}, voltResp{tfInd, contInd}, modelResp{tfInd, contInd}, numResp(tfInd, contInd, :), denResp(tfInd, contInd, :) ] = ...
            ComputeFigure6ModelResponses(params, tf(tfInd), cont(contInd), filters);

        % Print a status update to the terminal
        fprintf('tf %d of %d, contrast %d of %d, %f\n', tfInd, numTf, contInd, numCont, toc);
    end
end

%% Visualize filters

% Make filters
params.numPeriods = 1;
[ filters ] = MakeFigure6Filters(params, 1);

% Visualize filters
PlotFigure6ModelFilters(filters);

%% Make bar plots

% Iterate over models
for ind = 1:numModel

    % Get the response for the selected model
    selResp = cellfun(@(x) x(:,ind), meanResp, 'uni', false);

    % Make bar plots
    if contains(modelNames{ind}, {'three', 'Three'})
        MakeModelBarPlots(selResp, tf, cont, modelNames{ind}, inputNames, numResp, denResp);
    else
        MakeModelBarPlots(selResp, tf, cont, modelNames{ind}, inputNames);
    end
end

%% Plot TF sweeps

% Only plot if multiple TFs were tested
if length(tf) > 1

    % Iterate over models
    for ind = 1:numModel

        % Get the response for the selected model
        selResp = cellfun(@(x) x(:,ind), meanResp, 'uni', false);

        % Plot TF sweeps
        if contains(modelNames{ind}, {'three', 'Three'})
            PlotModelTFSweep(selResp, tf, cont, modelNames{ind}, inputNames, numResp, denResp);
        else
            PlotModelTFSweep(selResp, tf, cont, modelNames{ind}, inputNames);
        end
    end
end

%% Plot time traces (averaged over spatial phase)

% Iterate over models
for ind = 1:numModel

    % Get the response for the selected model
    selResp = cellfun(@(x) x(:,:,:,ind), modelResp, 'uni', false);
    selMeanResp = cellfun(@(x) x(:,ind), meanResp, 'uni', false);

    % Make bar plots
    PlotModelTimeTracesAveragedOverSpatialPhase(params, selResp, selMeanResp, tf, cont, modelNames{ind}, inputNames);
end


%% Plot linearity analysis from Weinecke et al. 2018
ind = 3;

% Get the response for the selected model
selResp = cellfun(@(x) x(:,:,:,ind), voltResp, 'uni', false);

% Plot linearity analysis
PlotLinearityAnalysis(selResp, tf, cont, modelNames{ind});

end