FFIModelParams_16_02_2021.m

%% 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
%Clear the workspace before saving and setting parameters
clear;
%% Compound eye properties

% Results in an eye with 7057 ommatidia, with a rough arc of 1.85 deg per
% ommatidium.  This broadly matches anatomical values from Horridge (1978)
% and Krapp & Gabbiani (2005)
radius = 49; %radius of the eye in ommatidia, including the center ommatidium
theta_max = ((radius-0.5)/48.5)*pi; %Maximum visual subtense angle [rad]
retina_size = 3*radius^2-3*radius+1; %number of retinotopic units (or ommatidia)     

%% Photoreceptors and LMCs

%Photoreceptor properties - obtained by fitting response to data from Jones
%& Gabbiani (2010)
K_photoD = 1.1536; %gain for derivative part of transfer function [V/s]
K_photoP = 6.8945; %gain for proportional part of transfer function [V/s^2]
a_1 = 66.5332; %first pole for 2nd-order transfer function [1/s]
a_2 = 11.3822; %second pole for 2nd-order transfer function [1/s]

%Lamina monopolar cell properties - similarly obtained by fitting response
%to data from Jones & Gabbiani (2010)
K_LMC = 1.0341; %gain for derivative part of transfer function
b = 3.9573; %first pole for 1st-order transfer function [1/s]

%% Lateral and Global Inhibition

%Lateral inhibition properties
tau_lat = 0.3; %time constant for low-pass filter used to delay lateral
                   %inhibition [s].  Adjusted manually based on model response
                   %to both grating and looming stimuli
                   
%The following matrix is used to distribute low-pass filtered LMC outputs
%to adjacent retinotopic units.  Weighting is inversely proportional to
%distance from the center unit.  This, along with the maximum range of six
%ommatidia (11.1 deg), is based on O'Shea & Rowell (1975)
inhibitory_kernel = [[1/36, 1/36, 1/36, 1/36, 1/36, 1/36, 1/36, 0, 0, 0, 0, 0, 0];...
                          [1/36, 1/30, 1/30, 1/30, 1/30, 1/30, 1/30, 1/36, 0, 0, 0, 0, 0];...
                          [1/36, 1/30, 1/24, 1/24, 1/24, 1/24, 1/24, 1/30, 1/36, 0, 0, 0, 0];...
                          [1/36, 1/30, 1/24, 1/18, 1/18, 1/18, 1/18, 1/24, 1/30, 1/36, 0, 0, 0];...
                          [1/36, 1/30, 1/24, 1/18, 1/12, 1/12, 1/12, 1/18, 1/24, 1/30, 1/36, 0, 0];...
                          [1/36, 1/30, 1/24, 1/18, 1/12, 1/6, 1/6, 1/12, 1/18, 1/24, 1/30, 1/36, 0];...
                          [1/36, 1/30, 1/24, 1/18, 1/12, 1/6, 0, 1/6, 1/12, 1/18, 1/24, 1/30, 1/36];...
                          [0, 1/36, 1/30, 1/24, 1/18, 1/12, 1/6, 1/6, 1/12, 1/18, 1/24, 1/30, 1/36];...
                          [0, 0, 1/36, 1/30, 1/24, 1/18, 1/12, 1/12, 1/12, 1/18, 1/24, 1/30, 1/36];...
                          [0, 0, 0, 1/36, 1/30, 1/24, 1/18, 1/18, 1/18, 1/18, 1/24, 1/30, 1/36];...
                          [0, 0, 0, 0, 1/36, 1/30, 1/24, 1/24, 1/24, 1/24, 1/24, 1/30, 1/36];...
                          [0, 0, 0, 0, 0, 1/36, 1/30, 1/30, 1/30, 1/30, 1/30, 1/30, 1/36];...
                          [0, 0, 0, 0, 0, 0, 1/36, 1/36, 1/36, 1/36, 1/36, 1/36, 1/36]];

%Global inhibition properties
K_GI = 0.2625; %value determining the degree of saturation in the sub-linear
               %global inhibition function.  Adjusted manually based on
               %comparison of model response to looming stimuli and grating
               %stimuli with recordings in Wang et al. (2018).  A starting 
               %point was obtained via particle swarm optimization of the
               %model's looming response to match data from Wang et al.
               %(2018)

%% General Spiking Neuron Properties

%Reversal potentials and threshold voltages
E_E = 0; %Reversal potential for excitatory currents [V].  Based on values 
         %used by Jones & Gabbiani (2012a) and Bazhenov et al. (2001)
E_I = -0.08; %Reversal potential for inhibitory currents [V].  Based on
             %K+ channel reversal potentials in other locust vision
             %modelling work (e.g. Peron & Gabbiani 2009)
E_A = -0.08; %Reversal potential for adaptation currents [V].  Based on
             %K+ channel reversal potentials in other locust vision
             %modelling work (e.g. Peron & Gabbiani 2009)
E_L = -0.065; %Reversal potential for leakage currents [V]
              %Used to give resting potential similar to that
              %given for the LGMD in works such as Jones & Gabbiani (2010)
              %and Dewell & Gabbiani (2018)
V_t = -0.055; %Voltage threshold for exponential function governing
              %spike initiation [V].  Difference from resting voltage E_L 
              %based approximately on EPSPs induced by single-facet stimuli
              %in Jones & Gabbiani (2010) which reliably
              %induced spikes in the LGMD.
V_fire = -0.045; %Voltage at which neuron fires and voltage resets [V]
                 %Chosen for computational reasons (to limit excessively
                 %high derivatives from the exponential spike initiation
                 %function)

%Other properties
c = 0.015; %Specific capacitance of the membrane [F/m^2], based on value
             %for the LGMD in Peron et al. (2007).  Standardized across
             %neurons to facilitate easier comparison.
T_noise = 0.00167; %Sample time for random excitatory noise [s], chosen to
                   %limit the maximum frequency component in random noise
                   %to 300 Hz

%% TmA Properties 

%The following parameters were determined via particle swarm optimization,
%using the comparison of a simulated, voltage-clamped LGMD's response to a
%single-facet OFF stimulus with recording data from Jones & Gabbiani (2010)
K_E_m = 1.1485e+02; %Gain for excitatory synapses [S/(m^2*V)]
g_L_m = 2.4075; %Specific conductance of the membrane [S/m^2]
Delta_t_m = 0.00289; %Governs sharpness of exponential spike initiation [V] 
Delta_g_A_m = 0.010742; %Increment in adaptation conductance upon firing [S/m^2]
tau_A_m = 0.14882; %Time constant for adaptation [s];
g_N_m = 0.25030; %Excitatory noise conductance strength [S/m^2]

%Neuron parameters determined via particle swarm optimization of response
%to a looming stimulus with l/|v| = 20 ms, compared with data from Wang et
%al. (2018).  Note that these wallowed to vary alongside FFI neuron
%parameters.
K_lat_m = 168.45; %Gain for lateral inhibitory synapses [S/(m^2*V)]
K_gl_m = 2.0985; %Gain for global inhibitory synapses [S/(m^2)]

%Other neuron parameters:
tau_alpha_m = 0.0003; %time constant for alpha function governing synaptic output [s], 

%% FFI Neuron Properties 

%Receptive field parameters, based on data from Wang et al. (2018)
%For an arc per ommatidium of 1.85 deg
r_field = 13; %Radius of the FFI receptive fields in retinotopic units (or ommatidia).
        %The center unit IS counted.  Corresponds to 46.4 deg width.
r_spacing = 12; %Spacing of the receptive fields in retinotopic units, which determines 
         %the total number of FFI neurons.  Corresponds to spacing of 22.3
         %deg.
r_pop = floor((radius-1)/r_spacing)+1; %Radius of the distribution of FFI neurons.
                                 %The center unit is NOT counted.
pop_f = 3*r_pop^2-3*r_pop+1; %Total number of inhibitory neurons

%Neuron parameters determined via particle swarm optimization and manual
%adjustment, based on comparison of the response to looming and grating 
%stimuli with data from Wang et al. (2018)
K_E_f = 0.0233; %Gain value applied to summed input to determine excitatory 
                  %conductance [S/m^2]
g_L_f = 0.50129; %Membrane leakage conductance [S/m^2]
Delta_t_f = 6.1186e-04; %Governs sharpness of exponential spike initiation [V]
tau_A_f = 0.031297; %Time constant for adaptation [s];
Delta_g_A_f = 3.5069; %Increment in adaptation conductance upon firing [S/m^2]
g_N_f = 0.084727; %Excitatory noise conductance strength [S/m^2]

%Other neuron parameters
tau_alpha_f = 0.003; %time constant for alpha function governing synaptic output
                 %of FFI neurons [s]

%% LGMD Properties

%LGMD geometry values, based on Peron et al. (2007)
L_b = 2e-4; %Dendritic branch length [m]
L_j = 1e-4; %Dendritic junction length [m]
L_c = 1e-4; %Calcium compartment length [m]
r_b = 1.8e-6; %Dendritic branch radius [m]
r_j = 2.5e-6; %Dendritic junction radius [m]
r_c = 2.5e-6; %Calcium compartment radius [m]

%LGMD general properties
R_ax = 0.6; %Axial resistivity [Ohm*m] from Peron et al. (2007)
g_L_LGMD = 1.1; %Leakage conductance for the membrane [S/m^2], based on 
                %Jones & Gabbiani (2012a) and Peron & Gabbiani (2009) 

%Calculated resistances, using Equation 6.30 in Dayan and Abbott (2001)
R_b = R_ax*L_b/(2*pi*r_b^2);
R_j = R_ax*L_j/(2*pi*r_j^2);
R_c = R_ax*L_c/(2*pi*r_c^2);

%LGMD dendritic branch properties
K_E_b = 3.7913e-08; %Gain value applied to inputs to determine excitatory 
                %conductances [S], determined during simulated voltage clamp
A_m_b = 2*pi*L_b*r_b; %Membrane surface area for each dendritic branch [m^2]
G_b_j = 1/(R_j+R_b); %Conductance from each dendritic branch to the dendritic
              %junction [S]

%LGMD dendritic junction properties
A_m_j = 2*pi*L_j*r_j; %Membrane surface area for the dendritic junction [m^2]
G_j_c = 1/(R_c+R_j); %Conductance from the dendritic junction to the calcium
                    %compartment [S]                                              

%% Looming stimulus properties
%These correspond to a looming stimulus with l/|v| of 20 ms
%For other looming stimuli, l should be specified as needed
l = 0.06; %stimulus half-height [m]
v = 3; %stimulus velocity [m/s]
d = 4.5; %stimulus starting distance [m]
T_max = (d-l/tan(theta_max/2))/v;%the end time for a looming simulation [s]
theta_0 = 2*atan(l/d); %initial visual subtense angle [rad]

%% Grating Stimulus Properties

%Chosen to match the grating stimulus presented in Wang et al. (2018)
theta_grating = 40*pi/180; %Width of the grating stimulus [rad]
phi_grating = 28.4*pi/180; %Angular speed of the grating stimulus [rad/s]
t_app = 2; %Time at which the stimulus appears [s]
t_move = 12; %Time at which the stimulus begins to move [s]
t_end = 22; %Time at which the stimulus disappears [s]
T_final = 22; %End time for the simulation [s]

save('FFIModelParams_16_02_2021.mat');
clear;