%% License Info
%A Model of Feedforward, Global, and Lateral Inhibition in the Locust Visual System.

%This model examines the architecture and function of inhibitory mechanisms in the 
%visual system of locusts, namely those involved in the processing of inputs to a
%key looming-sensitive neuron, the lobula giant movement detector (LGMD).

%Copyright (c) 2026, Erik Olson, Travis Wiens, and Jack Gray

%CITATION:
%When using the model code for scientific publications, cite the following work:
%Olson EGN, Wiens TK, Gray JR. A model of feedforward, global, and lateral inhibition 
%in the locust visual system predicts responses to looming stimuli. Biol Cybern. 2021 
%Jun;115(3):245-265. doi: 10.1007/s00422-021-00876-8. Epub 2021 May 16. PMID: 33997912.

%This program is free software: you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation, either version 3 of the License, or
%(at your option) any later version.

%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.

%You should have received a copy of the GNU General Public License
%along with this program.  If not, see <https://www.gnu.org/licenses/>.

%Contact: erik.olson@usask.ca

%% Code
%Creates Figure 11
%%
%Data Setup comes first.  Also, the data from Wang et al. (2018) needs to be loaded.

%Load data from Wang et al (2018) Fig. 4C
load('Wang_et_al_2018_FFI_Grating_Data.mat');

%Load data with no lateral inhibition
load('FFI_grating_no_lateral_Data.mat');
f_mean_no_lateral = f_mean;

%Load data for no global inhibition
load('FFI_grating_no_global_Data.mat');
f_mean_no_global = f_mean;

%Load data with no lateral inhibitory delay (and reduced lateral inhibitory gain)
load('FFI_grating_no_delay_Data.mat');
f_mean_no_delay = f_mean;

%Load data with no lateral inhibitory delay (and original gain value)
load('FFI_grating_no_delay_excess_gain_Data.mat');
f_mean_no_delay_excess_gain = f_mean;

%Load data for basic model response to grating
load('FFI_grating_Data.mat');
f_mean_base = f_mean;

%adjust model time values to match experiments (aligning t=0 to the grating
%appearance)
t_f = t_f-2;
%find the index where the grating appears in the experiment
grating_start_exp = find(t_f_exp>0);
grating_start_exp = grating_start_exp(1);
%find the experimental baseline (mean firing rate before grating
%appearance)
baseline_exp = mean(f_exp(1:grating_start_exp-1));

%find where the grating appears in the model
grating_start_model = find(t_f>0);
grating_start_model = grating_start_model(1);
%find the model baselines (mean firing rate before grating
%appearance)
baseline_model = mean(f_mean(1:grating_start_model-1));
baseline_no_lateral_model = mean(f_mean_no_lateral(1:grating_start_model-1));
baseline_no_global_model = mean(f_mean_no_global(1:grating_start_model-1));
baseline_no_delay_model = mean(f_mean_no_delay(1:grating_start_model-1));
baseline_no_delay_excess_gain_model = mean(f_mean_no_delay_excess_gain(1:grating_start_model-1));

%normalize to the baselines
f_exp_norm = f_exp-baseline_exp;

f_mean_norm = f_mean-baseline_model;
f_mean_no_lateral_norm = f_mean_no_lateral-baseline_no_lateral_model;
f_mean_no_global_norm = f_mean_no_global-baseline_no_global_model;
f_mean_no_delay_norm = f_mean_no_delay-baseline_no_delay_model;
f_mean_no_delay_excess_gain_norm = f_mean_no_delay_excess_gain-baseline_no_delay_excess_gain_model;

%%
%Plot Fig 11A
grating_plots = figure(1);
clf;
subplot(3,2,1:2);
hold on

y_bounds = [-20 180];

plot([0,0],y_bounds,'LineWidth',1,'LineStyle','--','Color',[0.8 0.8 0.8]);
plot([10,10],y_bounds,'LineWidth',1,'LineStyle','--','Color',[0.8 0.8 0.8]);
xlim([-2 20]);
ylim(y_bounds);
yticks(0:50:150);

experiment = plot(t_f_exp,f_exp_norm,'LineWidth',1,'Color',[0.8500 0.3250 0.0980],'DisplayName','Experiment');
model = plot(t_f,f_mean_norm,'LineWidth',1,'Color',[0 0.4470 0.7410],'DisplayName','Model');

text(0,160,'\leftarrow Grating appears')
text(10,160,{'\leftarrow Movement begins'})
legend([model,experiment],'Location','northeast');

ylabel({'Population Firing Rate for';'FFI Neurons Relative to';'Baseline [spikes/s]'});
xlabel('Time [s]');

title('A');

%%
%Plot Fig 11B
subplot(3,2,3);
hold on

y_bounds = [-20 170];

plot([0,0],y_bounds,'LineWidth',0.5,'LineStyle','--','Color',[0.8 0.8 0.8]);
plot([10,10],y_bounds,'LineWidth',0.5,'LineStyle','--','Color',[0.8 0.8 0.8]);
xlim([-2 20]);
ylim(y_bounds);
yticks(0:50:150);

experiment = plot(t_f_exp,f_exp_norm,'LineWidth',0.5,'Color',[0.8500 0.3250 0.0980]);
model = plot(t_f,f_mean_no_lateral_norm,'LineWidth',0.5,'Color',[0 0.4470 0.7410]);

ylabel({'Population Firing Rate for';'FFI Neurons Relative to';'Baseline [spikes/s]'});
xlabel('Time [s]');  

title('B');

%%
%Plot Fig 11C
subplot(3,2,4);
hold on

y_bounds = [-20 170];

plot([0,0],y_bounds,'LineWidth',0.5,'LineStyle','--','Color',[0.8 0.8 0.8]);
plot([10,10],y_bounds,'LineWidth',0.5,'LineStyle','--','Color',[0.8 0.8 0.8]);
xlim([-2 20]);
ylim(y_bounds);
yticks(0:50:250);

experiment = plot(t_f_exp,f_exp_norm,'LineWidth',0.5,'Color',[0.8500 0.3250 0.0980]);
model = plot(t_f,f_mean_no_global_norm,'LineWidth',0.5,'Color',[0 0.4470 0.7410]);

ylabel({'Population Firing Rate for';'FFI Neurons Relative to';'Baseline [spikes/s]'});
xlabel('Time [s]');  

title('C');

%%
%Plot Fig 11D
subplot(3,2,5);
hold on

y_bounds = [-20 100];

plot([0,0],y_bounds,'LineWidth',0.5,'LineStyle','--','Color',[0.8 0.8 0.8]);
plot([10,10],y_bounds,'LineWidth',0.5,'LineStyle','--','Color',[0.8 0.8 0.8]);
xlim([-2 20]);
ylim(y_bounds);
yticks(0:50:150);

experiment = plot(t_f_exp,f_exp_norm,'LineWidth',0.5,'Color',[0.8500 0.3250 0.0980]);
model = plot(t_f,f_mean_no_delay_excess_gain_norm,'LineWidth',0.5,'Color',[0 0.4470 0.7410]);

ylabel({'Population Firing Rate for';'FFI Neurons Relative to';'Baseline [spikes/s]'});
xlabel('Time [s]');  

title('D');

%%
%Plot Fig 11E
subplot(3,2,6);
hold on

y_bounds = [-20 100];

plot([0,0],y_bounds,'LineWidth',0.5,'LineStyle','--','Color',[0.8 0.8 0.8]);
plot([10,10],y_bounds,'LineWidth',0.5,'LineStyle','--','Color',[0.8 0.8 0.8]);
xlim([-2 20]);
ylim(y_bounds);
yticks(0:50:150);

experiment = plot(t_f_exp,f_exp_norm,'LineWidth',0.5,'Color',[0.8500 0.3250 0.0980]);
model = plot(t_f,f_mean_no_delay_norm,'LineWidth',0.5,'Color',[0 0.4470 0.7410]);

ylabel({'Population Firing Rate for';'FFI Neurons Relative to';'Baseline [spikes/s]'});
xlabel('Time [s]');  

title('E');

set(grating_plots,'Units','Inches');
set(grating_plots,'Position',[0 0 8 7]);
pos = get(grating_plots,'Position');
set(grating_plots,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);

savefig('grating_plots');
saveas(gcf,'Fig11.pdf');