%% 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;