%% 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
%Generates Figures 8, 9 and 10
%%

%Data Preparation.  Data from Wang et al. 2018 Figure S7 also needs to be
%loaded.

%(Model data is already processed into firing rate via histogram and
%Gaussian convolution.)

%Load model data for looming with l/|v| = 20 ms
load('FFI_l_over_v_20ms_Data.mat');
f_20ms = f_output;
f_singles_20ms = f_all;
sigma_f_20ms = std(f_all,0,2);
theta_20ms = theta;
V_trace_20ms = V_trace;
t_trace = t_raw;

%Load model data for looming with l/|v| = 40 ms
load('FFI_l_over_v_40ms_Data.mat')
f_40ms = f_output;
sigma_f_40ms = std(f_all,0,2);
theta_40ms = theta;

%Load model data for looming with l/|v| = 80 ms
load('FFI_l_over_v_80ms_Data.mat')
f_80ms = f_output;
sigma_f_80ms = std(f_all,0,2);
theta_80ms = theta;

%Load model data for looming with l/|v| = 20 ms and no global inhibition
load('FFI_l_over_v_20ms_no_global_Data.mat')
f_n_g_20ms = f_output;
sigma_f_n_g_20ms = std(f_all,0,2);
theta_n_g_20ms = theta;

%Load model data for looming with l/|v| = 20 ms and no lateral inhibition
load('FFI_l_over_v_20ms_no_lateral_Data.mat')
f_n_l_20ms = f_output;
sigma_f_n_l_20ms = std(f_all,0,2);
theta_n_l_20ms = theta;

%Clean up unnecessary variables
clear V_trace f_all f_output l_over_v raw t_raw theta

%Load data from Wang et al. 2018 Figure S7
load('Wang_et_al_2018_FFI_Looming_Data.mat');

close all

%% Linear Fit
%As data from Wang et al is from multiple animals, first normalize each
%animal's firing data to its peak firing rate.
%Create blank arrays to hold the normalized data:
f_exp_norm_20ms = NaN(size(f_exp_20ms));
f_exp_norm_40ms = NaN(size(f_exp_40ms));
f_exp_norm_80ms = NaN(size(f_exp_80ms));

%For each animal
for i = 1:6
    %Find the animal's peak rate across all three l/v values
    max_animal = max([max(f_exp_20ms(:,i)),max(f_exp_40ms(:,i)),max(f_exp_80ms(:,i))]);
    
    %Then normalize
    f_exp_norm_20ms(:,i) = f_exp_20ms(:,i)/max_animal;
    f_exp_norm_40ms(:,i) = f_exp_40ms(:,i)/max_animal;
    f_exp_norm_80ms(:,i) = f_exp_80ms(:,i)/max_animal;
end

%Next, resample experimental values so that they are evenly distributed 
%with respect to subtense angle, to avoid over-weighting of early parts of
%the looming stimuli.  Model values have already been sampled in this way.
theta_q_20ms = linspace(theta_20ms(1),theta_20ms(end),1001)';
f_exp_resampled_20ms = interp1(theta_exp_20ms,f_exp_norm_20ms,theta_q_20ms);

theta_q_40ms = linspace(theta_40ms(1),theta_40ms(end),1001)';
f_exp_resampled_40ms = interp1(theta_exp_40ms,f_exp_norm_40ms,theta_q_40ms);

theta_q_80ms = linspace(theta_80ms(1),theta_80ms(end),1001)';
f_exp_resampled_80ms = interp1(theta_exp_80ms,f_exp_norm_80ms,theta_q_80ms);

%concatenate experimental theta values, and add a column of ones for 
%fitting of a linear function
theta_exp_overall = [theta_q_20ms;theta_q_40ms;theta_q_80ms];
theta_exp_overall = [theta_exp_overall,ones(size(theta_exp_overall))];

%concatenate experimental firing rate data, and take the mean across the
%six locusts
f_exp_overall = mean([f_exp_resampled_20ms;f_exp_resampled_40ms;f_exp_resampled_80ms],2);

%Fit a linear function to experimental data
beta_exp = (theta_exp_overall'*theta_exp_overall)\(theta_exp_overall'*f_exp_overall);

%Find R-squared for the experimental data linear fit
SS_tot_exp = (f_exp_overall-mean(f_exp_overall))'*(f_exp_overall-mean(f_exp_overall));
SS_res_exp = (f_exp_overall-theta_exp_overall*beta_exp)'*(f_exp_overall-theta_exp_overall*beta_exp);
R_sq_exp = 1-SS_res_exp/SS_tot_exp;

%concatenate model theta values, and add a column of ones for 
%fitting of a linear function
theta_model_overall = [theta_20ms;theta_40ms;theta_80ms];
theta_model_overall = [theta_model_overall,ones(size(theta_model_overall))];

%concatenate model firing rate data
f_model_overall = [f_20ms;f_40ms;f_80ms];

%Fit a linear function to experimental data
beta_model = (theta_model_overall'*theta_model_overall)\(theta_model_overall'*f_model_overall);

%Find R-squared for the experimental data linear fit
SS_tot_model = (f_model_overall-mean(f_model_overall))'*(f_model_overall-mean(f_model_overall));
SS_res_model = (f_model_overall-theta_model_overall*beta_model)'*(f_model_overall-theta_model_overall*beta_model);
R_sq_model = 1-SS_res_model/SS_tot_model;


%% Plot Figure 8
l_over_v_plots = figure('Position', [0 0 800 700]);
%Fig.8 A, l/|v| = 20 ms
plot_firing(theta_20ms,f_20ms,sigma_f_20ms,theta_exp_20ms,f_exp_20ms,{'A1';'l/|v| = 20 ms'},l_over_v_plots,1,beta_model);
%Fig.8 B, l/|v| = 40 ms
plot_firing(theta_40ms,f_40ms,sigma_f_40ms,theta_exp_40ms,f_exp_40ms,{'A2';'l/|v| = 40 ms'},l_over_v_plots,2,beta_model);
%Fig.8 C, l/|v| = 80 ms
plot_firing(theta_80ms,f_80ms,sigma_f_80ms,theta_exp_80ms,f_exp_80ms,{'A3';'l/|v| = 80 ms'},l_over_v_plots,3,beta_model);

%8 D, showing a single trace at l/|v| = 20 ms
plot_single(theta_20ms,f_singles_20ms(:,1),theta_exp_20ms,f_exp_20ms,{'B';'l/|v| = 20 ms'},l_over_v_plots,4,beta_model);

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

savefig('l_over_v_plots');
saveas(gcf,'Fig8.pdf');

%% Figure 9
% Generate figure for voltage traces from 9 neurons across the eye centerline
% during a l/|v| = 20 ms loom
voltage_plots = figure('Position', [0 0 800 400]);
voltage_traces = plot(t_trace-1,V_trace_20ms(:,27:35),'LineWidth',0.5);
ylabel({'Membrane Voltages for';'Feedforward Inhibitory Neurons [V]'});
xlabel('Time Relative to Collision [s]');

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

savefig('voltage_plots');
saveas(gcf,'Fig9.pdf');

%% Figure 10
%Generate figures for model output during a l/|v| = 20 ms loom with
%multiple inhibition schemes
inhibition_plots = figure('Position', [0 0 800 400]);

%First plot values unnormalized (Fig 10A)
subplot(1,2,1);
hold on
complete_model = plot(theta_20ms,f_20ms,'LineWidth',2,'DisplayName','Complete Model');
no_lateral = plot(theta_n_l_20ms,f_n_l_20ms,'LineWidth',2,'DisplayName','No Lateral Inhibition');
no_global = plot(theta_n_g_20ms,f_n_g_20ms,'LineWidth',2,'DisplayName','No Global Inhibition');
ylabel({'Population Firing Rate for';'Feedforward Inhibitory Neurons [spikes/s]'});
xlabel('Stimulus Angular Size, theta [deg]');        
legend([complete_model,no_lateral,no_global],'Location','northwest');
title('A');

%Then plot values normalized (Fig 10B); note division of data by their maximum
%values
subplot(1,2,2);
hold on
f_norm_20ms = f_20ms/max(f_20ms);
f_norm_n_g_20ms = f_n_g_20ms/max(f_n_g_20ms);
f_norm_n_l_20ms = f_n_l_20ms/max(f_n_l_20ms);

complete_model = plot(theta_20ms,f_norm_20ms,'LineWidth',2,'DisplayName','Complete Model');
no_lateral = plot(theta_n_l_20ms,f_norm_n_l_20ms,'LineWidth',2,'DisplayName','No Lateral Inhibition');
no_global = plot(theta_n_g_20ms,f_norm_n_g_20ms,'LineWidth',2,'DisplayName','No Global Inhibition');

ylabel({'Normalized Population Firing Rate for';'Feedforward Inhibitory Neurons'});
xlabel('Stimulus Angular Size, theta [deg]');        
legend([complete_model,no_lateral,no_global],'Location','northwest');
title('B');

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

savefig('inhibition_plots');
saveas(gcf,'Fig10.pdf');

%% Functions used in Figure 8
function [] = plot_firing(theta_model,f_model,sigma_model,theta_exp,f_exp,title_text,use_figure,position,beta_model)
    %Plots a averaged model outputs and their standard deviations alongside
    %experimental values and standard deviations.  Used for Figure 8 A,B&C
    figure(use_figure);
    subplot(2,2,position);

    hold on;
    
    %Find mean and standard deviation across the 6 locusts in 
    %Wang et al. 2018 Figure S7
    f_exp_mean = mean(f_exp,2);
    sigma_exp = std(f_exp,0,2);
    
    %Create borders for shading the experimental standard deviation region
    X_std_dev_exp_border = [theta_exp;flipud(theta_exp)];
    Y_std_dev_exp_border = [(f_exp_mean-sigma_exp);flipud(f_exp_mean+sigma_exp)];
    
    %Create borders for shading the model standard deviation region
    X_std_dev_border = [theta_model;flipud(theta_model)];
    Y_std_dev_border = [(f_model-sigma_model);flipud(f_model+sigma_model)];

    %plot model standard deviation
    standard_dev = fill(X_std_dev_border,Y_std_dev_border,[0 0.4470 0.7410]);
    set(standard_dev,'facealpha',0.5,'EdgeColor','none');
    
    %plot experimental standard deviation
    standard_dev_exp = fill(X_std_dev_exp_border,Y_std_dev_exp_border,[0.8500 0.3250 0.0980]);
    set(standard_dev_exp,'facealpha',0.5,'EdgeColor','none');

    %plot model output 
    data_model = plot(theta_model,f_model,'Color',[0 0.4470 0.7410],'LineWidth',2,'DisplayName','Model');
    %plot the experimental data
    data_exp = plot(theta_exp,f_exp_mean,'Color',[0.8500 0.3250 0.0980],'LineWidth',2,'DisplayName','Experiment');
    
    %plot a linear fit to model output
    data_fit = plot(theta_model,beta_model(1)*theta_model+beta_model(2),'Color',[0 0 0],'LineWidth',1,'DisplayName','Linear Fit to Model');
    
    ylabel({'Population Firing Rate for';'Feedforward Inhibitory Neurons [spikes/s]'});
    
    xlabel('Stimulus Angular Size, theta [deg]');        
    
    legend([data_model,data_exp,data_fit],'Location','northwest');
    
    title(title_text);
    
    xlim([0 90]);
    ylim([0 250]);
end
function [] = plot_single(theta_model,f_model,theta_exp,f_exp,title_text,use_figure,position,beta_model)
    %Plots a single model run alongside experimental values and standard
    %deviations.  Used for Figure 8 D
    
    figure(use_figure);
    subplot(2,2,position);

    hold on;
    
    %Find mean and standard deviation across the 6 locusts in 
    %Wang et al. 2018 Figure S7
    f_exp_mean = mean(f_exp,2);
    sigma_exp = std(f_exp,0,2);

    %Create borders for shading the experimental standard deviation region
    X_std_dev_exp_border = [theta_exp;flipud(theta_exp)];
    Y_std_dev_exp_border = [(f_exp_mean-sigma_exp);flipud(f_exp_mean+sigma_exp)];

    %plot experimental standard deviation
    standard_dev_exp = fill(X_std_dev_exp_border,Y_std_dev_exp_border,[0.8500 0.3250 0.0980]);
    set(standard_dev_exp,'facealpha',0.5,'EdgeColor','none');

    %plot model output 
    data_model = plot(theta_model,f_model,'Color',[0 0.4470 0.7410],'LineWidth',2,'DisplayName','Model');
    %plot the experimental data
    data_exp = plot(theta_exp,f_exp_mean,'Color',[0.8500 0.3250 0.0980],'LineWidth',2,'DisplayName','Experiment');
    
    %plot a linear fit to model output
    data_fit = plot(theta_model,beta_model(1)*theta_model+beta_model(2),'Color',[0 0 0],'LineWidth',1,'DisplayName','Linear Fit to Model');
    
    ylabel({'Population Firing Rate for';'Feedforward Inhibitory Neurons [spikes/s]'});
    
    xlabel('Stimulus Angular Size, theta [deg]');        
    
    legend([data_model,data_exp,data_fit],'Location','northwest');
    
    title(title_text);
    
    xlim([0 90]);
    ylim([0 250]);
end