|
| 1 | +function dydt = RingAttractorODESolver(t,y,session) |
| 2 | + |
| 3 | +persistent y_record; % can be used when an algorithm depends on the activity history |
| 4 | + |
| 5 | +persistent vel_backup; |
| 6 | +persistent adjusted_vel; |
| 7 | + |
| 8 | + |
| 9 | +prm_ra = session.parameters.ring_attractor; |
| 10 | +prm_inputs = session.parameters.inputs; |
| 11 | + |
| 12 | +nw = prm_ra.n_wedge_neurons; |
| 13 | +ni = prm_inputs.n_input_nodes; |
| 14 | + |
| 15 | +%% These two variables will be updated |
| 16 | +y1 = y(1:nw); % wedge neuron activity (ring attractor activity) |
| 17 | +y2 = y(nw+1:end); % W_input |
| 18 | +W_input = reshape(y2, nw, ni); |
| 19 | + |
| 20 | + |
| 21 | +%% Take the index of the current time point |
| 22 | +dt = session.sim_conds.dt; |
| 23 | +ind = t/dt+1; |
| 24 | +indf = floor(ind); |
| 25 | +indc = ceil(ind); |
| 26 | + |
| 27 | + |
| 28 | + |
| 29 | + |
| 30 | +%% Record keeping |
| 31 | +if isempty(y_record) |
| 32 | + y_record = zeros(nw, numel(session.sim_conds.t)); |
| 33 | +end |
| 34 | +if indc >0 && indc <= numel(session.sim_conds.t) |
| 35 | + y_record(:,indc) = y1; |
| 36 | +end |
| 37 | + |
| 38 | + |
| 39 | + |
| 40 | + |
| 41 | +%% Obtain interpolated inputs to wedge neurons |
| 42 | + |
| 43 | +% input_neuron_activity: input from ring neurons (negative if ring neurons are inhibitory) |
| 44 | +% vel_signal: velocity input |
| 45 | +% wedge_injection_signal: direct current injection |
| 46 | + |
| 47 | +if indf==indc |
| 48 | + input_neuron_activity = session.sim_conds.visual_input_neurons(:,indf); |
| 49 | + vel_signal = session.sim_conds.vel(indf); |
| 50 | + wedge_injection_signal = session.sim_conds.wedge_current_injection(:,indf); |
| 51 | +else |
| 52 | + input_neuron_activity = (indc-ind)*session.sim_conds.visual_input_neurons(:,indf) + (ind-indf)*session.sim_conds.visual_input_neurons(:,indc); |
| 53 | + vel_signal = (indc-ind)*session.sim_conds.vel(indf) + (ind-indf)*session.sim_conds.vel(indc); |
| 54 | + wedge_injection_signal = (indc-ind)*session.sim_conds.wedge_current_injection(:,indf) + (ind-indf)*session.sim_conds.wedge_current_injection(:,indc); |
| 55 | +end |
| 56 | + |
| 57 | + |
| 58 | +%% Calculate the current from each input type to wedge neurons |
| 59 | + |
| 60 | +%%% |
| 61 | +% 1. Current from visual neurons |
| 62 | +input_current_from_visual_neurons = W_input * (input_neuron_activity * 2*pi/ni) ; |
| 63 | +if strcmp(session.parameters.plasticity.rule, 'Cope') |
| 64 | + [~,innn] = find(input_neuron_activity<max(input_neuron_activity)); |
| 65 | + input_neuron_activity(innn) = 0; |
| 66 | + input_current_from_visual_neurons = session.parameters.plasticity.w_p * W_input * (input_neuron_activity * 2*pi/ni) ; |
| 67 | +end |
| 68 | +if prm_inputs.input_is_excitatory_1_inhibitory_m1 < 0 |
| 69 | + input_current_from_visual_neurons = -input_current_from_visual_neurons; |
| 70 | +end |
| 71 | + |
| 72 | + |
| 73 | +%%% |
| 74 | +% 2. Turning signal |
| 75 | + |
| 76 | +% Scale the vel signal (discretization) |
| 77 | +vel = vel_signal*nw/2/pi; |
| 78 | +% To avoid asymmetricity, I used { [f(t+dt)-f(t)]/dt + [f(t)-f(t-dt)]/dt }/2 |
| 79 | +turning_signal = vel*(y1([end,1:end-1]) - y1([2:end, 1]))/2; % calculate turning_signal at a given moment; |
| 80 | + |
| 81 | + |
| 82 | +%%% |
| 83 | +% 3. Current injection |
| 84 | +wedge_current_injection = wedge_injection_signal * 2*pi/nw; |
| 85 | + |
| 86 | + |
| 87 | + |
| 88 | + |
| 89 | +%% Calculate the delta of the ring attractor |
| 90 | + |
| 91 | +tmp = prm_ra.W_ring_attractor*y1 + 1 + input_current_from_visual_neurons + turning_signal + wedge_current_injection; |
| 92 | +tmp(tmp>prm_ra.membrane_saturation) = prm_ra.membrane_saturation; |
| 93 | +tmp(tmp<prm_ra.membrane_threshold) = prm_ra.membrane_threshold; |
| 94 | + |
| 95 | +y1_tmp = tmp; clear tmp; |
| 96 | +delta_y1 = (y1_tmp - y1) ./ prm_ra.tau_wedge; |
| 97 | + |
| 98 | + |
| 99 | + |
| 100 | + |
| 101 | + |
| 102 | +%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 103 | +%% Update the W_input with a plasticity rule. |
| 104 | + |
| 105 | +ina_rep = repmat(input_neuron_activity', nw,1); |
| 106 | +wedge_rep = repmat(y1,1,ni); |
| 107 | + |
| 108 | +switch session.parameters.plasticity.rule |
| 109 | + |
| 110 | + case 'No learning' |
| 111 | + dW_input_dt = zeros(size(ina_rep)); |
| 112 | + |
| 113 | + |
| 114 | + otherwise |
| 115 | + |
| 116 | + epsilon = session.parameters.plasticity.epsilon_W_input; |
| 117 | + W_max = session.parameters.plasticity.W_max; |
| 118 | + |
| 119 | + % The learning rate is assumed to be velocity dependent. |
| 120 | + % The actual source of velocity dependent modulation is not known. |
| 121 | + |
| 122 | + if numel(vel_backup) ~= numel(session.sim_conds.vel) || ... |
| 123 | + vel_backup(indc) ~= session.sim_conds.vel(indc) || ... |
| 124 | + t == session.sim_conds.t(1) |
| 125 | + vel_backup = session.sim_conds.vel; |
| 126 | + tmp = vel_backup; |
| 127 | + tmp = tmp.^2; |
| 128 | + sss = mean(tmp)+1.5*std(tmp); |
| 129 | + adjusted_vel = tmp/sss; % scaling |
| 130 | + |
| 131 | + end |
| 132 | + |
| 133 | + % interpolate |
| 134 | + if indf==indc |
| 135 | + fv = adjusted_vel(indf); |
| 136 | + else |
| 137 | + fv = (indc-ind)*adjusted_vel(indf) + (ind-indf)*adjusted_vel(indc); |
| 138 | + end |
| 139 | + |
| 140 | + |
| 141 | + % Compute dW |
| 142 | + |
| 143 | + switch session.parameters.plasticity.rule |
| 144 | + case 'SOM inhib, Post-synaptically gated, input profile' |
| 145 | + f_th = 0.04; % about half of the maximum wedge neuron activity. So, this can be dynamically adjustable, but not implemented. |
| 146 | + [PF, NF] = half_wave_rectify(wedge_rep-f_th); % PF is the positive part, NF is the negative part. Both are positive. |
| 147 | + |
| 148 | + dW_input_dt = 3* fv * epsilon * wedge_rep .* ( W_max - ina_rep - W_input ) ; |
| 149 | + |
| 150 | + % In case, the wedge neurons are noisy, it may need to be |
| 151 | + % thresholded. |
| 152 | + % dW_input = fv * epsilon * PF .* ( W_max - ina_rep - W_input ) ; |
| 153 | + |
| 154 | + case 'Hebb inhib, Pre-synaptically gated, wedge profile' |
| 155 | + %%% *** IMPORTANT: The mamp in the "sim_cond.m" should be small. See line 111 of sim_cond.m. |
| 156 | + g_th = 0.1/3; % About a bit less than the median of input activity. |
| 157 | + g_th = g_th*ones(size(ina_rep)); |
| 158 | + [PG, NG] = half_wave_rectify(ina_rep-g_th); % PG is the positive part, NG is the negative part. Both are positive. |
| 159 | + |
| 160 | + %%% adaptive version |
| 161 | + %t_duration = 5; |
| 162 | + %tmp_ind = max(indf-round(t_duration/dt),1):max(indf,1); |
| 163 | + %scale_factor = repmat( max(0.02, max(y_record(:,tmp_ind),[],2) ), 1, ni); |
| 164 | + %scale_factor = repmat( 0.02 + max(y_record(:,tmp_ind),[],2) , 1, ni); |
| 165 | + |
| 166 | + %%% Fixed version |
| 167 | + scale_factor = 0.08 * ones(size(wedge_rep)); |
| 168 | + |
| 169 | + dW_input_dt = 6* fv * epsilon * ( W_max - (wedge_rep./scale_factor)*W_max - W_input ) .* PG; |
| 170 | + |
| 171 | + end |
| 172 | + |
| 173 | + |
| 174 | + |
| 175 | +end |
| 176 | + |
| 177 | +if t>10 |
| 178 | + t = t; |
| 179 | +end |
| 180 | + |
| 181 | +% New W_input state |
| 182 | +W_input = W_input + dW_input_dt; |
| 183 | + |
| 184 | + |
| 185 | +% Cap the value |
| 186 | +W_input(W_input<0) = 0; |
| 187 | +if exist('W_max', 'var') |
| 188 | + W_input(W_input>W_max) = W_max; |
| 189 | +end |
| 190 | + |
| 191 | +% Calculate the delta |
| 192 | +y2_tmp = reshape(W_input, numel(W_input),1); |
| 193 | +delta_y2 = (y2_tmp - y2); |
| 194 | + |
| 195 | + |
| 196 | + |
| 197 | +%% Combine results |
| 198 | +dydt = [delta_y1;delta_y2]; |
| 199 | + |
| 200 | + |
| 201 | + |
| 202 | +%% Occasionally display the simulation time |
| 203 | +if mod(t,20)<0.001 && rand()>0.7 |
| 204 | + disp([' ' num2str(t) 's']); |
| 205 | +end |
| 206 | + |
| 207 | + |
| 208 | +return; |
0 commit comments