Spike frequency adaptation in the LGMD (Peron and Gabbiani 2009)

 Download zip file 
Help downloading and running models
Accession:116945
This model is used in the referenced paper to demonstrate that a model of an SK-like calcium-sensitive potassium (KCa) conductance can replicate the spike frequency adaptation (SFA) of the locust lobula giant movement detector (LGMD) neuron. The model simulates current injection experiments with and without KCa block in the LGMD, as well as visual stimulation experiments with and without KCa block.
Reference:
1 . Peron S, Gabbiani F (2009) Spike frequency adaptation mediates looming stimulus selectivity in a collision-detecting neuron. Nat Neurosci 12:318-26 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Locust Lobula Giant Movement Detector (LGMD) neuron;
Channel(s): I Na,t; I K,Ca; I Calcium; I Krp;
Gap Junctions:
Receptor(s): GabaA; Cholinergic Receptors;
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Spike Frequency Adaptation; Vision;
Implementer(s): Gabbiani, F; Peron, Simon [perons at janelia.hhmi.org];
Search NeuronDB for information about:  GabaA; Cholinergic Receptors; I Na,t; I K,Ca; I Calcium; I Krp;
% This will return a conductance-as-function-of-time for treating dendrites
%  as single compartment
%
%   class: 1 = 10x10 box
%          2 = 10x80 bar
%          3 = looming
%          4 = accelarating box
%          5 = loomlike box
%   velocity - in deg/s -- for class = 4, specify start and end velocity in 2-element vector
%   velocity - in deg/s -- for class = 5, specify the angle
%   direction: 1 = AP, -1 = PA, 2=DV, -2=VD
%   dt: time step for *simulation*
%   synaptic_map_path: which map to use? uniform is current dflt
%   n_syns_per_input: how many synapses *per facet*?
%
% The following variables alter timing to account for feedforward inhibition and the observation
%  that larger inputs are more synchronized.  
%
%   dvn: parameters of delay-versus-nfacets function - delay = v(1)./(1+exp((nsyns-v(2))/v(3)))+v(4)
%   sdvn: parameters of sd of delay-versus-nsyns function - delay_jitter = v(1)./(1+exp((nsyns-v(2))/v(3)))+v(4)
%         note that the actual delay is jittered by drawing from a normal distro with the above sd;
%         you read it correctly -- dvn is based on the number of facets, while sdvn is based on the
%         number of synapses
%   vel_sf: this accounts for the observation that there is speed modulation speed_sf = v(1)./(1+exp((nsyns-v(2))/v(3)))+v(4)
%
function g_syn_of_t = get_synaptic_pattern(class, direction, velocity, syn_tau, syn_gmax, ...
                      duration, dt, synaptic_map_path, n_syns_per_input, dvn, sdvn, vel_sf)
  % --- Definitions
  % synaptic map -- to determine how many synapses get triggered
  if (exist('synaptic_map_path') == 0 | synaptic_map_path == -1)
    synaptic_map_path = 'uniform_synmap.mat'; % The synaptic mapping [cmpt frac az el]
  end
	dt_image = 5; % screen inter-step time; monitor specific 
  half_angle = 2.5; % acceptance half-angle in degrees; add

  % for translating stimuli:
  edge_sf = [1 0]; % apply 1 to the leading edge, 2 to the trailing edge [1 0], for instance, means leading only

  % --- Preliminaries
  load(synaptic_map_path);
	n_syns = length(synmap);
	
  % --- derive the timing pattern
	switch (class)
	  case 1 % 10x10 box
			vert_lims = [-5-half_angle 5+half_angle];
			hor_lims = [85-half_angle 95+half_angle];
			bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
		  switch(direction)
			  case 1 % AP
          [synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [50 130], dt_image, bar_thickness, vert_lims, edge_sf);
				case -1 % PA
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [130 50], dt_image, bar_thickness, vert_lims, edge_sf);
				case 2 % DV
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [40 -40], dt_image, bar_thickness, hor_lims, edge_sf);
				case -2 % VD
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [-40 40], dt_image, bar_thickness, hor_lims, edge_sf);
		  end;
	  case 1.1 % 10x10 box, but with bounds that are used for the SFA 'review'/model paper
			vert_lims = [-5-half_angle 5+half_angle];
			hor_lims = [85-half_angle 95+half_angle];
			bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
		  switch(direction)
			  case 1 % AP
          [synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [60 120], dt_image, bar_thickness, vert_lims, edge_sf);
				case -1 % PA
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [120 60], dt_image, bar_thickness, vert_lims, edge_sf);
				case 2 % DV
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [30 -30], dt_image, bar_thickness, hor_lims, edge_sf);
				case -2 % VD
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [-30 30], dt_image, bar_thickness, hor_lims, edge_sf);
		  end;
	  case 2 % 10x80 bar
			vert_lims = [-40-half_angle 40+half_angle];
			hor_lims = [50-half_angle 130+half_angle];
			bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
		  switch(direction)
			  case 1 % AP
          [synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [50 130], dt_image, bar_thickness, vert_lims, edge_sf);
				case -1 % PA
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'v', velocity, [130 50], dt_image, bar_thickness, vert_lims, edge_sf);
				case 2 % DV
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [40 -40], dt_image, bar_thickness, hor_lims, edge_sf);
				case -2 % VD
					[synpos, timing, inst_vel] = pattern_translating_bar(synmap, 'h', velocity, [-40 40], dt_image, bar_thickness, hor_lims, edge_sf);
		  end;
    case 3 % loom l/v = velocity
      % calculate so that you finish with 80 deg, then go for 100 ms more
      end_time = velocity/(tand(80/2));
      start_time = duration - (end_time + 100);
      theta_start = 2*atand(velocity/start_time);
      [synpos, timing, inst_vel] = pattern_circular_looming_stimulus(synmap, velocity, [0 90], theta_start, dt_image);
	  case 4 % 10x10 box, accelarating linearly
			vert_lims = [-5-half_angle 5+half_angle];
			hor_lims = [85-half_angle 95+half_angle];
			bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
		  switch(direction)
			  case 1 % AP
          [synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, 'v', velocity, [60 120], 5, bar_thickness, vert_lims, edge_sf);
				case -1 % PA
					[synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, 'v', velocity, [120 60], 5, bar_thickness, vert_lims, edge_sf);
				case 2 % DV
					[synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, 'h', velocity, [30 -30], 5, bar_thickness, hor_lims, edge_sf);
				case -2 % VD
					[synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, 'h', velocity, [-30 30], 5, bar_thickness, hor_lims, edge_sf);
		  end;
    case 5 % 10x10 box, loomlike linearly
			vert_lims = [-5-half_angle 5+half_angle];
			hor_lims = [85-half_angle 95+half_angle];
			bar_thickness = 10+(2*half_angle); % in dimension parallel to that of travel
		  switch(direction)
			  case 1 % AP
          [synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, 'v', velocity, [60 120], 5, bar_thickness, vert_lims, edge_sf);
				case -1 % PA
					[synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, 'v', velocity, [120 60], 5, bar_thickness, vert_lims, edge_sf);
				case 2 % DV
					[synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, 'h', velocity, [30 -30], 5, bar_thickness, hor_lims, edge_sf);
				case -2 % VD
					[synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, 'h', velocity, [-30 30], 5, bar_thickness, hor_lims, edge_sf);
		  end;
  end;

	% --- the preceding logic returns a timing vector, which tells you the timing of ommatidial
	% --- activation.  however, response at the level of the LGMD is non-linear; specifically,
	% --- synaptic resposne increases superlinearly and occurs faster for larger luminance changes;
	% --- this is true for both inhibitory and excitatory inputs.  We will do this by making the 
	% --- assumption that intensity change is approximated by the fraction of synapses activated per timestep

  % --- for each element of timing, calculate the scaling factor based on velocity ...
	% --- this must be done BEFORE any other timing alterations as it depends on the OMMATIDIAL
	% --- transition time, which is not the same as the timing vetor after considering delay
	vsf = ones(size(timing));
	tv = inst_vel.t_vec;
	for t=0:dt:duration
	  idx = find (timing >= t & timing < t + dt_image);
		ivi = min(find(tv >= t-1));
		if (length(ivi) == 0)  
		  %disp(['warning: ivi not found; t = ' num2str(t)]) ; 
			continue; 
	  end
		vsf(idx) = vel_sf(1)./(1+exp((inst_vel.vel_vec(ivi)-vel_sf(2))/vel_sf(3)))+vel_sf(4);
%if(length(idx) > 0) ;		disp(['vsf: ' num2str(vsf(idx(1))) ' vel: ' num2str(inst_vel.vel_vec(ivi)) ' t: ' num2str(t) ' ivi: ' num2str(ivi)]); end
	end

	% --- first, introduce mean delay via dvn based on number of facets hit in a given timestep
	n_timing = timing; 
	for t=0:dt_image:duration
	  idx = find(timing >= t & timing < t + dt_image); 
	  nfacets = length(idx);
		n_timing(idx) = timing(idx) + dvn(1)./(1+exp((nfacets-dvn(2))/dvn(3)))+dvn(4);
	end
	timing = n_timing;


	% --- Now account for the fact that you do not necessarily have the same number
	% --- of synapses as facests
  n_timing = zeros(1,floor(length(synpos)*n_syns_per_input));
	cf = 1;
	for t=0:dt_image:duration
	  idx = find(timing >= t & timing < t + dt_image); 
	  nfacets = length(idx);
		nsyns = nfacets*n_syns_per_input;

    % For a non-integer number of synapses, scale based on a normal pdf
	  if (nsyns ~= round(nsyns))
		  rem = nsyns - floor(nsyns);
			if (rem > rand) 
			  nsyns = floor(nsyns) + 1;
			else
			  nsyns = floor(nsyns);
			end
    end
	  if (nsyns ~= round(nsyns)) ; disp (['wtf:' num2str(nsyns) ' rem: ' num2str(rem)]);  end

    % Now make nsyns in n_timing have the time
    n_timing(cf:min(length(n_timing),cf+nsyns)) = t;
		if (length(idx) > 0)
			n_vsf(cf:min(length(n_timing),cf+nsyns)) = mean (vsf(idx));
		else
			n_vsf(cf:min(length(n_timing),cf+nsyns)) = 1;
		end
		cf = cf + nsyns;
		if (cf > length(n_timing))
		  disp('Warning: vector too big');
			cf = length(n_timing);
		end
	end
	timing = n_timing(find(n_timing > 0));
	vsf = n_vsf;

  % --- it is still unknown if there is an actual increase in input strength from stronger stimuli;
	% --- what is certain is that jitter declines.  do jitter here. *still based on number of facets*
	n_timing = timing; 
	for t=0:dt_image:duration
	  idx = find(timing >= t & timing < t + dt_image); 
	  nsyns = length(idx);
		jitter = normrnd(0, sdvn(1)./(1+exp(nsyns-sdvn(2))/sdvn(3))+sdvn(4), length(idx), 1);
		if (length(idx) > 0) ; n_timing(idx) = timing(idx) + jitter'; end
	end
	timing = n_timing;

	% --- compute the conductance-as-function-of-time
	g_syn_of_t = zeros(1,round(duration/dt)+1);
	for t=0:dt:duration
	  active_idx = find (timing >= t-syn_tau*20 & timing <= t);
		for a = 1:length(active_idx)
		  T = t-timing(active_idx(a));
			S = vsf(active_idx(a));
		  g_syn_of_t(1+round(t/dt)) = g_syn_of_t(1+round(t/dt)) + S*syn_gmax*(T/syn_tau)*exp(1-T/syn_tau);
		end
	end
	
	% --- sanity check
	if (0 == 1)
	  plot(0:dt:duration, g_syn_of_t);
  end
	
% 
% This should return timing and synapses for loomlike bars of CONSTANT 
%  ANGULAR SIZE AND CHANGING ANGULAR VELOCITY.  Each position is activated
%  TWICE - once  for the ON and once for the OFF.
%
% RETURNS:
%  synpos - [cmpt_idx cmpt_fraction]
%  timing - in ms
%
% ASSIGNED:
%  synmap - [cmpt cmpt_frac az el]
%  direction - 'h' 'v' for horizontal, vertical (DIRECTION OF MOTION, not bar
%              orientation!) - default is vertical
%  velocity - in degrees / second -- l/v ; negative --> receding
%  limits - [start end] where to start; if horizontal motion, give azimuths; 
%           if vertical, elevations. 
%  dt - temporal resolution, in ms
%  thick - thickness of the bar, in degrees (i.e., dimension along axis PARALLEL
%          to motion vector)
%  length_limits - edges of the bar, in degrees for dimension along axis PERPENDICULAR to
%           motion vector
%  
%
function [synpos, timing, inst_vel] = pattern_loomlike_bar(synmap, direction, velocity, ...
                                                 limits, dt, thick, length_limits, ...
                                                 edge_sf)
  % equivalent disc RADIUS at which to start
  theta_start = 5;
 
  % Preliminary
  timing = [];

  % Orientation-dependent component
  if (strcmp('h', direction))
    positional_array = synmap(:,4); 
    length_array = synmap(:,3); 
  else
    positional_array = synmap(:,3); 
    length_array = synmap(:,4); 
  end
  used = zeros(1,length(positional_array));

	% Introduce variation into the eye -- vary each ommatidial view by sd of 5 degrees
  eye_deformation = normrnd(0,10,length(positional_array),1);
  positional_array = positional_array + eye_deformation;
  
  % Determine valid based on length limits
  length_valids = find(length_array > min(length_limits) & length_array < max(length_limits));

  % determine total time based on the limits and the accelaration
  % loom time
  t_o = velocity/tand(theta_start); 
  t_f = velocity/tand(diff(limits) + theta_start);
  if (velocity < 0) ; t_o = -1*t_o; t_f = -1*t_f ; end
  total_time = 5*(1+floor(abs(t_f-t_o)/5));

  % instantaneous velocity vector
  inst_vel.vel_vec = zeros(length(0:dt:total_time),1);
	ivi = 1;

  % loop thru time, and add accordingly
  position = limits(1); % LEADING edge
  for time=0:dt:total_time
    last_position = position;  
    if (velocity < 0) % 'receding'
      theta_t = atand(velocity/(t_f+time));
      theta_tp1 = atand(velocity/(t_f+(time+dt)));
    else % normal
      theta_t = atand(velocity/(t_o-time));
      theta_tp1 = atand(velocity/(t_o-(time+dt)));
    end
    dp = abs(theta_tp1-theta_t);
    
    % change in position/time is velocity!
    int_vel.vel_vec = dp;

    if (time+dt > t_o) ; dp = 0; end

    stimulated_leading = [];
    stimulated_trailing = [];
    
    % Increasing position value (e.g., rightward)
    if (limits(1) < limits(2))
      position = position + dp;
      % Leading edge
      if (position < limits(2))
        stimulated_leading = find(positional_array >= last_position & positional_array < position);
      end
      % Trailing edge
      if (position-thick > limits(1))
        stimulated_trailing = find(positional_array >= last_position-thick & positional_array < position-thick);
      end
    % For decreasing position value
    else 
      position = position - dp;
      % Leading edge
      if (position < limits(1))
        stimulated_leading = find(positional_array < last_position & positional_array >= position);
      end
      % Trailing edge
      if (position+thick > limits(2))
        stimulated_trailing = find(positional_array < last_position+thick & positional_array >= position+thick);
      end
    end
 
    % Constrain by length
    stimulated_trailing = intersect(stimulated_trailing,length_valids);
    stimulated_leading = intersect(stimulated_leading,length_valids);
    if (edge_sf(1) > 0)
      stimulated_leading = stimulated_leading(1:round(length(stimulated_leading)*edge_sf(1)));
    else
      stimulated_leading = [];
    end
    if (edge_sf(2) > 0)
      stimulated_trailing = stimulated_trailing(1:round(length(stimulated_trailing)*edge_sf(2)));
    else
      stimulated_trailing = [];
    end
      
    % And append returned arrays if needbe
    if (length(stimulated_leading) > 0)
      for s=1:length(stimulated_leading)
        synpos(length(timing) + 1, :) = [synmap(stimulated_leading(s),1) synmap(stimulated_leading(s),2)];
        timing(length(timing) + 1) = time;
      end
    end
    if (length(stimulated_trailing) > 0)
      for s=1:length(stimulated_trailing)
        synpos(length(timing) + 1, :) = [synmap(stimulated_trailing(s),1) synmap(stimulated_trailing(s),2)];
        timing(length(timing) + 1) = time;
      end
    end
  end  

  % Scale instantaneous velocity by dt -- in s
  inst_vel.t_vec = 0:dt:total_time;
  inst_vel.vel_vec = abs([inst_vel.vel_vec/(.001*dt)])';
	
% 
% This should return timing and synapses for accelarating bars of CONSTANT 
%  ANGULAR SIZE AND CHANGING ANGULAR VELOCITY.  Each position is activated
%  TWICE - once  for the ON and once for the OFF.
%
% RETURNS:
%  synpos - [cmpt_idx cmpt_fraction]
%  timing - in ms
%
% ASSIGNED:
%  synmap - [cmpt cmpt_frac az el]
%  direction - 'h' 'v' for horizontal, vertical (DIRECTION OF MOTION, not bar
%              orientation!) - default is vertical
%  velocity - in degrees / second (1): start (2): end
%  limits - [start end] where to start; if horizontal motion, give azimuths; 
%           if vertical, elevations. Bar will 'edge in' from one side and edge
%           out the other (i.e., it does not start wholly formed).
%  dt - temporal resolution, in ms
%  thick - thickness of the bar, in degrees (i.e., dimension along axis PARALLEL
%          to motion vector)
%  length_limits - edges of the bar, in degrees for dimension along axis PERPENDICULAR to
%           motion vector
%  
%
function [synpos, timing, inst_vel] = pattern_accelarating_bar(synmap, direction, velocity, ...
                                                    limits, dt, thick, length_limits, ...
                                                    edge_sf)
  % Preliminary
  timing = [];

  % Orientation-dependent component
  if (strcmp('h', direction))
    positional_array = synmap(:,4); 
    length_array = synmap(:,3); 
  else
    positional_array = synmap(:,3); 
    length_array = synmap(:,4); 
  end
  used = zeros(1,length(positional_array));

	% Introduce variation into the eye -- vary each ommatidial view by sd of 5 degrees
  eye_deformation = normrnd(0,10,length(positional_array),1);
  positional_array = positional_array + eye_deformation;
  
  % Determine valid based on length limits
  length_valids = find(length_array > min(length_limits) & length_array < max(length_limits));
  
  % determine total time based on the limits and the accelaration
  % total time = total_distance/mean_velocity [assuming LINEAR accelaration]
  mean_velocity = mean(velocity);
  total_time = abs(abs(diff(limits)))/(mean_velocity/1000);

  % instantaneous velocity vector
  inst_vel.vel_vec = zeros(length(0:dt:total_time),1);
	ivi = 1;
 
  % loop thru time, and add accordingly
  position = limits(1); % LEADING edge
  dv_dt = diff(velocity)/1000/total_time; % again, convert to deg/ms2
  v_o = velocity(1)/1000; % convert to d/ms
  for time=0:dt:total_time
    last_position = position;  
    vel_t = dv_dt*time + v_o;
    dp = vel_t*dt;    
    
    % change in position/time is velocity!
    int_vel.vel_vec = dp;

    stimulated_leading = [];
    stimulated_trailing = [];
    
    % Increasing position value (e.g., rightward)
    if (limits(1) < limits(2))
      position = position + dp;
      % Leading edge
      if (position < limits(2))
        stimulated_leading = find(positional_array >= last_position & positional_array < position);
      end
      % Trailing edge
      if (position-thick > limits(1))
        stimulated_trailing = find(positional_array >= last_position-thick & positional_array < position-thick);
      end
    % For decreasing position value
    else 
      position = position - dp;
      % Leading edge
      if (position < limits(1))
        stimulated_leading = find(positional_array < last_position & positional_array >= position);
      end
      % Trailing edge
      if (position+thick > limits(2))
        stimulated_trailing = find(positional_array < last_position+thick & positional_array >= position+thick);
      end
    end
 
    % Constrain by length
    stimulated_trailing = intersect(stimulated_trailing,length_valids);
    stimulated_leading = intersect(stimulated_leading,length_valids);
    if (edge_sf(1) > 0)
      stimulated_leading = stimulated_leading(1:round(length(stimulated_leading)*edge_sf(1)));
    else
      stimulated_leading = [];
    end
    if (edge_sf(2) > 0)
      stimulated_trailing = stimulated_trailing(1:round(length(stimulated_trailing)*edge_sf(2)));
    else
      stimulated_trailing = [];
    end
      
    % And append returned arrays if needbe
    if (length(stimulated_leading) > 0)
      for s=1:length(stimulated_leading)
        synpos(length(timing) + 1, :) = [synmap(stimulated_leading(s),1) synmap(stimulated_leading(s),2)];
        timing(length(timing) + 1) = time;
      end
    end
    if (length(stimulated_trailing) > 0)
      for s=1:length(stimulated_trailing)
        synpos(length(timing) + 1, :) = [synmap(stimulated_trailing(s),1) synmap(stimulated_trailing(s),2)];
        timing(length(timing) + 1) = time;
      end
    end
  end  

  % Scale instantaneous velocity by dt -- in s
  inst_vel.t_vec = 0:dt:total_time;
  inst_vel.vel_vec = abs([inst_vel.vel_vec/(.001*dt)])';

	
% 
% This should return timing and synapses for translating bars of CONSTANT 
%  ANGULAR SIZE AND ANGULAR VELOCITY.  Each position is activated TWICE - once
%  for the ON and once for the OFF.
%
% RETURNS:
%  synpos - [cmpt_idx cmpt_fraction]
%  timing - in ms
%
% ASSIGNED:
%  synmap - [cmpt cmpt_frac az el]
%  direction - 'h' 'v' for horizontal, vertical (DIRECTION OF MOTION, not bar
%              orientation!) - default is vertical
%  velocity - in degrees / second
%  limits - [start end] where to start; if horizontal motion, give azimuths; 
%           if vertical, elevations. Bar will 'edge in' from one side and edge
%           out the other (i.e., it does not start wholly formed).
%  dt - temporal resolution, in ms
%  thick - thickness of the bar, in degrees (i.e., dimension along axis PARALLEL
%          to motion vector)
%  length_limits - edges of the bar, in degrees for dimension along axis PERPENDICULAR to
%           motion vector
%  edge_sf - scales the input strength of each edge
%  
%
function [synpos, timing, inst_vel] = pattern_translating_bar(synmap, direction, velocity, ...
                                                    limits, dt, thick, length_limits, ...
                                                    edge_sf )
  % Preliminary
  timing = [];

  % Orientation-dependent component
  if (strcmp('h', direction))
    positional_array = synmap(:,4); 
    length_array = synmap(:,3); 
  else
    positional_array = synmap(:,3); 
    length_array = synmap(:,4); 
  end
  used = zeros(1,length(positional_array));

	% Introduce variation into the eye -- vary each ommatidial view by sd of 5 degrees
  eye_deformation = normrnd(0,10,length(positional_array),1);
  positional_array = positional_array + eye_deformation;
  
  % Determine valid based on length limits
  length_valids = find(length_array > min(length_limits) & length_array < max(length_limits));
  
  % determine total time ... limits + thickness scaled by velocity
  total_time = abs(abs(diff(limits)) + thick)/(velocity/1000);

  % loop thru time, and add accordingly
  position = limits(1); % LEADING edge
  dp = (velocity/1000)*dt;
  for time=0:dt:total_time
    last_position = position;

    stimulated_leading = [];
    stimulated_trailing = [];
    
    % Increasing position value (e.g., rightward)
    if (limits(1) < limits(2))
      position = position + dp;
      % Leading edge
      if (position < limits(2))
        stimulated_leading = find(positional_array >= last_position & positional_array < position);
      end
      % Trailing edge
      if (position-thick > limits(1))
        stimulated_trailing = find(positional_array >= last_position-thick & positional_array < position-thick);
      end
    % For decreasing position value
    else 
      position = position - dp;
      % Leading edge
      if (position < limits(1))
        stimulated_leading = find(positional_array < last_position & positional_array >= position);
      end
      % Trailing edge
      if (position+thick > limits(2))
        stimulated_trailing = find(positional_array < last_position+thick & positional_array >= position+thick);
      end
    end
 
    % Constrain by length
    stimulated_trailing = intersect(stimulated_trailing,length_valids);
    stimulated_leading = intersect(stimulated_leading,length_valids);
    if (edge_sf(1) > 0)
      stimulated_leading = stimulated_leading(1:round(length(stimulated_leading)*edge_sf(1)));
    else
      stimulated_leading = [];
    end
    if (edge_sf(2) > 0)
      stimulated_trailing = stimulated_trailing(1:round(length(stimulated_trailing)*edge_sf(2)));
    else
      stimulated_trailing = [];
    end
    
    % And append returned arrays if needbe
    if (length(stimulated_leading) > 0)
      for s=1:length(stimulated_leading)
        synpos(length(timing) + 1, :) = [synmap(stimulated_leading(s),1) synmap(stimulated_leading(s),2)];
        timing(length(timing) + 1) = time;
      end
    end
    if (length(stimulated_trailing) > 0)
      for s=1:length(stimulated_trailing)
        synpos(length(timing) + 1, :) = [synmap(stimulated_trailing(s),1) synmap(stimulated_trailing(s),2)];
        timing(length(timing) + 1) = time;
      end
    end
  end  

  % Scale instantaneous velocity by dt -- in s
  inst_vel.t_vec = 0:dt:total_time;
  inst_vel.vel_vec = abs(velocity*ones(size(inst_vel.t_vec)));

% 
% This should return timing and synapses for looming stimulus. Each position
%  is activated only once.  The stimulus is square.
%  synpos - [cmpt_idx cmpt_fraction]
%  timing - in ms
%  synmap - [cmpt cmpt_frac az el]
%  l_over_v - the l/v parameter
%  center - [el az] where to start
%  theta_start - size at start; zero would take infinite time. in degrees.
%  dt - temporal resolution, in ms
%
function [synpos, timing] = pattern_square_looming_stimulus(synmap, l_over_v, center, theta_start, dt)
  % Preliminary
  timing = [];
  theta_final = 80; % Degrees - final size, as in screen

  % Start, end time - compute it
  start_time = l_over_v/(tand(theta_start/2));
  end_time = l_over_v/(tand(theta_final/2));

  % GEnerate a center-centric coordinate set
  centered_el = synmap(:,4)-center(1);
  centered_az = synmap(:,3)-center(2);

  % determine total time ...
  total_time = start_time - end_time;
 
  % loop thru time, and add accordingly
  theta_0 = theta_start; % Theta at START of step
  for time=0:dt:total_time
    theta_f = 2*atand(l_over_v/(total_time + end_time - time));

    % Find any points between theta_0 and theta_f object size
    stimulated_1 = find(centered_el < theta_f & centered_el >= (-1*theta_f) & centered_az > theta_0 & centered_az <= theta_f);
    stimulated_2 = find(centered_el < theta_f & centered_el >= (-1*theta_f) & centered_az < -1*theta_0 & centered_az >= -1*theta_f);
    stimulated_3 = find(centered_az > -1*theta_0 & centered_az < theta_0 & centered_el >= -1*theta_f & centered_el < -1*theta_0);
    stimulated_4 = find(centered_az > -1*theta_0 & centered_az < theta_0 & centered_el <= theta_f & centered_el > theta_0);

    stimulated = unique([stimulated_1' stimulated_2' stimulated_3' stimulated_4']);
    
    % And append returned arrays if needbe
    if (length(stimulated) > 0)
      for s=1:length(stimulated)
        synpos(length(timing) + 1, :) = [synmap(stimulated(s),1) synmap(stimulated(s),2)];
        timing(length(timing) + 1) = time;
      end
    end
    theta_0 = theta_f;
  end  

% 
% This should return timing and synapses for looming stimulus. Each position
%  is activated only once.  The stimulus is square.
%  synpos - [cmpt_idx cmpt_fraction]
%  timing - in ms
%  inst_vel - the instantaneous velocity of the stimulus
%
%  synmap - [cmpt cmpt_frac az el]
%  l_over_v - the l/v parameter
%  center - [el az] where to start
%  theta_start - size at start; zero would take infinite time. in degrees.
%  dt - temporal resolution, in ms
%
function [synpos, timing, inst_vel] = pattern_circular_looming_stimulus(synmap, l_over_v, center, theta_start, dt)
  % Preliminary
  timing = [];
  theta_final = 80; % Degrees - final size, as in screen

  recede = 0;
  if (l_over_v < 0) % recede
    recede = 1;
    l_over_v = -1*l_over_v;
    theta_start = -1*theta_start;
  end

  % Start, end time - compute it
  start_time = l_over_v/(tand(theta_start/2));
  end_time = l_over_v/(tand(theta_final/2));

  % GEnerate a center-centric coordinate set
  centered_el = synmap(:,4)-center(1);
  centered_az = synmap(:,3)-center(2);

  % determine total time ...
  total_time = start_time - end_time;
 
  % For each point, determine the angular deflection from the center 
  for s=1:length(centered_el)
    alpha(s) = acosd(cosd(centered_el(s))*cosd(centered_az(s)));
  end

  % instantaneous velocity vector
  inst_vel.vel_vec = zeros(length(0:dt:total_time),1);
	ivi = 1;
  
  % loop thru time, and add accordingly
  if (recede == 0)
		theta_0 = theta_start; % Theta at START of step
		for time=0:dt:total_time
			theta_f = 2*atand(l_over_v/(total_time + end_time - time));

			% Find any points between theta_0 and theta_f object size
			stimulated = find(alpha < theta_f & alpha > theta_0);
			
			% And append returned arrays if needbe
			if (length(stimulated) > 0)
				for s=1:length(stimulated)
					synpos(length(timing) + 1, :) = [synmap(stimulated(s),1) synmap(stimulated(s),2)];
					timing(length(timing) + 1) = time;
				end
			end
			inst_vel.vel_vec(ivi) = theta_f - theta_0; ivi = ivi + 1;
			theta_0 = theta_f;
		end  
  else
		theta_0 = theta_final;
		for time=total_time:-1*dt:0
			theta_f = 2*atand(l_over_v/(total_time + end_time - time));
			% Find any points between theta_0 and theta_f object size
			stimulated = find(alpha > theta_f & alpha < theta_0);
			% And append returned arrays if needbe
			if (length(stimulated) > 0)
				for s=1:length(stimulated)
					synpos(length(timing) + 1, :) = [synmap(stimulated(s),1) synmap(stimulated(s),2)];
					timing(length(timing) + 1) = total_time-time;
				end
			end
			inst_vel.vel_vec(ivi) = theta_f - theta_0; ivi = ivi + 1;
			theta_0 = theta_f;
		end 
  end

  % Scale instantaneous velocity by dt -- in s
  inst_vel.t_vec = 0:dt:total_time;
  inst_vel.vel_vec = abs([inst_vel.vel_vec/(.001*dt)])';


Loading data, please wait...