function dydt = RingAttractorODESolver(t,y,session)
persistent y_record; % can be used when an algorithm depends on the activity history
persistent vel_backup;
persistent adjusted_vel;
prm_ra = session.parameters.ring_attractor;
prm_inputs = session.parameters.inputs;
nw = prm_ra.n_wedge_neurons;
ni = prm_inputs.n_input_nodes;
%% These two variables will be updated
y1 = y(1:nw); % wedge neuron activity (ring attractor activity)
y2 = y(nw+1:end); % W_input
W_input = reshape(y2, nw, ni);
%% Take the index of the current time point
dt = session.sim_conds.dt;
ind = t/dt+1;
indf = floor(ind);
indc = ceil(ind);
%% Record keeping
if isempty(y_record)
y_record = zeros(nw, numel(session.sim_conds.t));
end
if indc >0 && indc <= numel(session.sim_conds.t)
y_record(:,indc) = y1;
end
%% Obtain interpolated inputs to wedge neurons
% input_neuron_activity: input from ring neurons (negative if ring neurons are inhibitory)
% vel_signal: velocity input
% wedge_injection_signal: direct current injection
if indf==indc
input_neuron_activity = session.sim_conds.visual_input_neurons(:,indf);
vel_signal = session.sim_conds.vel(indf);
wedge_injection_signal = session.sim_conds.wedge_current_injection(:,indf);
else
input_neuron_activity = (indc-ind)*session.sim_conds.visual_input_neurons(:,indf) + (ind-indf)*session.sim_conds.visual_input_neurons(:,indc);
vel_signal = (indc-ind)*session.sim_conds.vel(indf) + (ind-indf)*session.sim_conds.vel(indc);
wedge_injection_signal = (indc-ind)*session.sim_conds.wedge_current_injection(:,indf) + (ind-indf)*session.sim_conds.wedge_current_injection(:,indc);
end
%% Calculate the current from each input type to wedge neurons
%%%
% 1. Current from visual neurons
input_current_from_visual_neurons = W_input * (input_neuron_activity * 2*pi/ni) ;
if strcmp(session.parameters.plasticity.rule, 'Cope')
[~,innn] = find(input_neuron_activity<max(input_neuron_activity));
input_neuron_activity(innn) = 0;
input_current_from_visual_neurons = session.parameters.plasticity.w_p * W_input * (input_neuron_activity * 2*pi/ni) ;
end
if prm_inputs.input_is_excitatory_1_inhibitory_m1 < 0
input_current_from_visual_neurons = -input_current_from_visual_neurons;
end
%%%
% 2. Turning signal
% Scale the vel signal (discretization)
vel = vel_signal*nw/2/pi;
% To avoid asymmetricity, I used { [f(t+dt)-f(t)]/dt + [f(t)-f(t-dt)]/dt }/2
turning_signal = vel*(y1([end,1:end-1]) - y1([2:end, 1]))/2; % calculate turning_signal at a given moment;
%%%
% 3. Current injection
wedge_current_injection = wedge_injection_signal * 2*pi/nw;
%% Calculate the delta of the ring attractor
tmp = prm_ra.W_ring_attractor*y1 + 1 + input_current_from_visual_neurons + turning_signal + wedge_current_injection;
tmp(tmp>prm_ra.membrane_saturation) = prm_ra.membrane_saturation;
tmp(tmp<prm_ra.membrane_threshold) = prm_ra.membrane_threshold;
y1_tmp = tmp; clear tmp;
delta_y1 = (y1_tmp - y1) ./ prm_ra.tau_wedge;
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Update the W_input with a plasticity rule.
ina_rep = repmat(input_neuron_activity', nw,1);
wedge_rep = repmat(y1,1,ni);
switch session.parameters.plasticity.rule
case 'No learning'
dW_input_dt = zeros(size(ina_rep));
otherwise
epsilon = session.parameters.plasticity.epsilon_W_input;
W_max = session.parameters.plasticity.W_max;
% The learning rate is assumed to be velocity dependent.
% The actual source of velocity dependent modulation is not known.
if numel(vel_backup) ~= numel(session.sim_conds.vel) || ...
vel_backup(indc) ~= session.sim_conds.vel(indc) || ...
t == session.sim_conds.t(1)
vel_backup = session.sim_conds.vel;
tmp = vel_backup;
tmp = tmp.^2;
sss = mean(tmp)+1.5*std(tmp);
adjusted_vel = tmp/sss; % scaling
end
% interpolate
if indf==indc
fv = adjusted_vel(indf);
else
fv = (indc-ind)*adjusted_vel(indf) + (ind-indf)*adjusted_vel(indc);
end
% Compute dW
switch session.parameters.plasticity.rule
case 'SOM inhib, Post-synaptically gated, input profile'
f_th = 0.04; % about half of the maximum wedge neuron activity. So, this can be dynamically adjustable, but not implemented.
[PF, NF] = half_wave_rectify(wedge_rep-f_th); % PF is the positive part, NF is the negative part. Both are positive.
dW_input_dt = 3* fv * epsilon * wedge_rep .* ( W_max - ina_rep - W_input ) ;
% In case, the wedge neurons are noisy, it may need to be
% thresholded.
% dW_input = fv * epsilon * PF .* ( W_max - ina_rep - W_input ) ;
case 'Hebb inhib, Pre-synaptically gated, wedge profile'
%%% *** IMPORTANT: The mamp in the "sim_cond.m" should be small. See line 111 of sim_cond.m.
g_th = 0.1/3; % About a bit less than the median of input activity.
g_th = g_th*ones(size(ina_rep));
[PG, NG] = half_wave_rectify(ina_rep-g_th); % PG is the positive part, NG is the negative part. Both are positive.
%%% adaptive version
%t_duration = 5;
%tmp_ind = max(indf-round(t_duration/dt),1):max(indf,1);
%scale_factor = repmat( max(0.02, max(y_record(:,tmp_ind),[],2) ), 1, ni);
%scale_factor = repmat( 0.02 + max(y_record(:,tmp_ind),[],2) , 1, ni);
%%% Fixed version
scale_factor = 0.08 * ones(size(wedge_rep));
dW_input_dt = 6* fv * epsilon * ( W_max - (wedge_rep./scale_factor)*W_max - W_input ) .* PG;
end
end
if t>10
t = t;
end
% New W_input state
W_input = W_input + dW_input_dt;
% Cap the value
W_input(W_input<0) = 0;
if exist('W_max', 'var')
W_input(W_input>W_max) = W_max;
end
% Calculate the delta
y2_tmp = reshape(W_input, numel(W_input),1);
delta_y2 = (y2_tmp - y2);
%% Combine results
dydt = [delta_y1;delta_y2];
%% Occasionally display the simulation time
if mod(t,20)<0.001 && rand()>0.7
disp([' ' num2str(t) 's']);
end
return;
