%% 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
function [I_out,t_out,max_f,f_base,delay_final] = VoltageClampTrialFunction(V_clamp,parFile,smoothing)
%Performs a simulated voltage clamp experiment on the LGMD, showing the
%response to a single-facet stimulus (i.e. input from one TmA).  The
%recording site is proximal to the SIZ.  Results are generated for four
%different luminance change durations: 0 ms, 83 ms, 183 ms and 350 ms.
%Current traces are offset temporally based on cross-correlation with
%experimental data from Jones & Gabbiani (2010).
%Inputs:
%parfile = name of the parameter file containing model parameters
%V_clamp = the clamp voltage [V]
%smoothing = the standard deviation of the Gaussian filter used to
%determine maximum firing rate [ms], which must be given as an integer
%Outputs:
%I_out = the current traces for each stimulus duration [A], averaged across
%ten trials. Each column corresponds to a single stimulus duration, going
%from shortest to longest, left to right
%t_out = the corresponding time values [s]
%max_f = the maximum firing rate of the TmA [spikes/s]
%f_base = the baseline firing rate of an unstimulated TmA, due to noise
%[Hz]
%delay_final = the delay added to current outputs [s]

load('Jones_and_Gabbiani_2010_Voltage_Clamp_Data');
load(parFile);

options = simset('srcworkspace','current','ReturnWorkspaceOutputs','on');
results = sim("VoltageClampTrial.slx",0.449,options);
t = results.yout{1}.Values.Time; %Access time data
y1 = results.yout{1}.Values.Data; %Access dendritic branch current data
y2 = results.yout{2}.Values.Data;
y3 = results.yout{3}.Values.Data;
y4 = results.yout{4}.Values.Data;
spiketrain = results.yout{5}.Values.Data;
Simulink.sdi.clear;

for i = 1:9
    options = simset('srcworkspace','current','ReturnWorkspaceOutputs','on');
    results = sim("VoltageClampTrial.slx",0.449,options);
    t = results.yout{1}.Values.Time; %Access time data
    y1 = y1+results.yout{1}.Values.Data; %Access dendritic branch current data
    y2 = y2+results.yout{2}.Values.Data;
    y3 = y3+results.yout{3}.Values.Data;
    y4 = y4+results.yout{4}.Values.Data;
    spiketrain = spiketrain+results.yout{5}.Values.Data;
    Simulink.sdi.clear;        
end

I_model1 = 0.1*y1;
I_model2 = 0.1*y2;
I_model3 = 0.1*y3;
I_model4 = 0.1*y4;

I_exp1 = interp1(t_exp,I_exp1,t);
I_exp2 = interp1(t_exp,I_exp2,t);
I_exp3 = interp1(t_exp,I_exp3,t);
I_exp4 = interp1(t_exp,I_exp4,t);

cross_corr1 = xcorr(I_model1,I_exp1);
cross_corr2 = xcorr(I_model2,I_exp2);
cross_corr3 = xcorr(I_model3,I_exp3);
cross_corr4 = xcorr(I_model4,I_exp4);
best_corr1 = max(cross_corr1);
best_corr2 = max(cross_corr2);
best_corr3 = max(cross_corr3);
best_corr4 = max(cross_corr4);

indices1 = find(cross_corr1 == best_corr1);
indices2 = find(cross_corr2 == best_corr2);
indices3 = find(cross_corr3 == best_corr3);
indices4 = find(cross_corr4 == best_corr4);
delays1 = indices1-max(size(I_model1));
delays2 = indices2-max(size(I_model2));
delays3 = indices3-max(size(I_model3));
delays4 = indices4-max(size(I_model4));

residuals_sq1 = [];
residuals_sq2 = [];
residuals_sq3 = [];
residuals_sq4 = [];

for i = 1:max(size(indices1))
    add_on = zeros(abs(delays1(i)),1);
    if delays1(i) < 0
        I_model_temp = [add_on;I_model1];
        I_model_temp = I_model_temp(1:max(size(I_exp1)));

        residual = I_model_temp-I_exp1;
        residuals_sq1 = [residuals_sq1,residual'*residual];
    else
        I_model_temp = [I_model1;add_on];
        I_model_temp = I_model_temp(1:max(size(I_exp1)));

        residual = I_model_temp-I_exp1;
        residuals_sq1 = [residuals_sq1,residual'*residual];
    end
end

for i = 1:max(size(indices2))
    add_on = zeros(abs(delays2(i)),1);
    if delays2(i) < 0
        I_model_temp = [add_on;I_model2];
        I_model_temp = I_model_temp(1:max(size(I_exp2)));

        residual = I_model_temp-I_exp2;
        residuals_sq2 = [residuals_sq2,residual'*residual];
    else
        I_model_temp = [I_model2;add_on];
        I_model_temp = I_model_temp(1:max(size(I_exp2)));

        residual = I_model_temp-I_exp2;
        residuals_sq2 = [residuals_sq2,residual'*residual];
    end
end

for i = 1:max(size(indices3))
    add_on = zeros(abs(delays3(i)),1);
    if delays3(i) < 0
        I_model_temp = [add_on;I_model3];
        I_model_temp = I_model_temp(1:max(size(I_exp3)));

        residual = I_model_temp-I_exp3;
        residuals_sq3 = [residuals_sq3,residual'*residual];
    else
        I_model_temp = [I_model3;add_on];
        I_model_temp = I_model_temp(1:max(size(I_exp3)));

        residual = I_model_temp-I_exp3;
        residuals_sq3 = [residuals_sq3,residual'*residual];
    end
end

for i = 1:max(size(indices4))
    add_on = zeros(abs(delays4(i)),1);
    if delays4(i) < 0
        I_model_temp = [add_on;I_model4];
        I_model_temp = I_model_temp(1:max(size(I_exp4)));

        residual = I_model_temp-I_exp4;
        residuals_sq4 = [residuals_sq4,residual'*residual];
    else
        I_model_temp = [I_model4;add_on];
        I_model_temp = I_model_temp(1:max(size(I_exp4)));

        residual = I_model_temp-I_exp4;
        residuals_sq4 = [residuals_sq4,residual'*residual];
    end
end

[~,index1] = min(residuals_sq1);
[~,index2] = min(residuals_sq1);
[~,index3] = min(residuals_sq1);
[~,index4] = min(residuals_sq1);
delay1 = delays1(index1);
delay2 = delays2(index2);
delay3 = delays3(index3);
delay4 = delays4(index4);

delay_int = round(mean([delay1,delay2,delay3,delay4]));

delay_final = -1*delay_int*(t(2)-t(1));

add_on = zeros(abs(delay_int),1);
if delay_int < 0
    I_model_temp = [add_on;I_model1];
    I_model1 = I_model_temp(1:max(size(I_exp1)));
    I_model_temp = [add_on;I_model2];
    I_model2 = I_model_temp(1:max(size(I_exp1)));
    I_model_temp = [add_on;I_model3];
    I_model3 = I_model_temp(1:max(size(I_exp1)));
    I_model_temp = [add_on;I_model4];
    I_model4 = I_model_temp(1:max(size(I_exp1)));
else
    I_model_temp = [I_model1;add_on];
    I_model1 = I_model_temp(1:max(size(I_exp1)));
    I_model_temp = [I_model2;add_on];
    I_model2 = I_model_temp(1:max(size(I_exp1)));
    I_model_temp = [I_model3;add_on];
    I_model3 = I_model_temp(1:max(size(I_exp1)));
    I_model_temp = [I_model4;add_on];
    I_model4 = I_model_temp(1:max(size(I_exp1)));
end

t_out = t;
I_out = [I_model1,I_model2,I_model3,I_model4];

max_f = 0.1*max(histoFiringRate(spiketrain,t,0.001,smoothing));

spikecount = 0;

for i = 1:10
    options = simset('srcworkspace','current','ReturnWorkspaceOutputs','on');
    T_final = 500;
    results = sim("TmABaselineModel.slx",T_final,options);
    y1 = results.yout{1}.Values.Data;
    spikecount = spikecount+sum(y1);
end

f_base = spikecount/(T_final*10);
end