function [MT,MST]=MODE_motion(Stim,time,noise_sig_MT,noise_sig_MST,K3,K4,InhWeights_MTMST)
%
% MODE_motion :: Function that performs the computations in the motion
% pathway from Retina to medial superior temporal area in each trial
%
%% Input variables
% Stim :: Random dot motion stimulus (3D matrix: 2 spatial dimensions and 1
% temporal) (Scale of pixel [luminance] intensities is 0-255)
%
% time :: Duration (in sec) for which middle temporal and medial superior
% temporal recordings are computed in response to the random dot motion stimulus
%
% noise_sig_MT :: Parameter that scales the standard deviation of the
% cellular noise that can be injected into middle temporal dynamics
%
% noise_sig_MST :: Parameter that scales the standard deviation of the
% cellular noise that can be injected into medial superior temporal dynamics
%
% K3 :: Parameter that scales inhibition in the directional inhibitory interneuron layer
%
% K4 :: Parameter that scales inhibition in the directional transient cell layer
%
% InhWeights_MTMST :: Parameter vector that scales the inter-directional competition
% in the motion capture circuit of middle temporal and medial superior temporal areas
% [w0 w1 w2 w3 w4], where wj modulates the amount of competition between
% global motion cells whose preferred directions differ by j*pi/4
% NOTE: MODE model is sensitive only to eight equally-spaced directions
%
%% Output variables
% MT :: Dynamic population activity of various middle temporal pools of
% neuorns tuned to different directions
%
% MST :: Dynamic population activity of various medial superior temporal
% pools of neurons tuned to different directions
%
%% Reference
% Grossberg, S. and Pilly, P. K. (2008). Temporal dyanamics of decision-making during motion perception in the visual cortex. Vision Research, 48(12), 1345-1373.
%
%% Author
% Praveen K. Pilly (advaitp@gmail.com)
%
%% NOTE
% The motion processing stages of the MODE model are adapted from
% Grossberg et al. (2001) and Berzhanskaya et al. (2007)
%
%% License policy
% Written by Praveen K. Pilly, Department of Cognitive and Neural Systems, Boston University
% Copyright 2009, Trustees of Boston University
%
% Permission to use, copy, modify, distribute, and sell this software and its documentation for any purpose is hereby granted
% without fee, provided that the above copyright notice and this permission notice appear in all copies, derivative works and
% associated documentation, and that neither the name of Boston University nor that of the author(s) be used in advertising or
% publicity pertaining to the distribution or sale of the software without specific, prior written permission. Neither Boston
% University nor its agents make any representations about the suitability of this software for any purpose. It is provided "as
% is" without warranty of any kind, either express or implied. Neither Boston University nor the author indemnify any
% infringement of copyright, patent, trademark, or trade secret resulting from the use, modification, distribution or sale of
% this software.
%
%% Last modified
% June 25, 2009
%%
% To set the pseudo-random number generator to a random state to begin with
rand('state',sum(100*clock))
randn('state',sum(100*clock))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The forward Euler's method is used to numerically integrate the system of
% model equations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fixed time step of numerical integration
dt=0.001; % (in sec)
%%% In the experiments, the stimulus was presented at 60 Hz
% frame duration
Fdur=1/60; % (in sec)
% Number of equally-spaced directions to which the model is sensitive
numdir=8;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Model variable initializations
[rows,cols]=size(Stim(:,:,1)); % size of the stimulus array
x=zeros(rows,cols);
z=ones(rows,cols);
b=zeros(rows,cols);
c=zeros(rows,cols,numdir);
e=zeros(rows,cols,numdir);
f=zeros(rows,cols,numdir);
h=zeros(rows,cols,numdir);
m=zeros(rows,cols,numdir);
y=zeros(rows,cols,numdir);
t=0; % time variable
tn=1; % counter for the time steps
% Outputs
MT=zeros(1,numdir); % at t=0
MST=zeros(1,numdir); % at t=0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Model parameters
% Stage 1
A1=1;
B1=10;
Tb=0.1;
% Stage 2
A2=1;
K2=50;
% Stage 3
A3=5;
C3=5;
% K3=20;
% Stage 4
A4=50;
C4=5;
% K4=20;
Te=0.2;
% Stage 5
A5=0.5;
Tf=0.5;
% Stage 6
A6=5;
D6=1;
% Stage 7
A7=10;
% Stage 8
A8=10;
C8=10;
%%% Various kernel generations
% Kernel G in Stage 5
sigG1=2;
sigG2=0.5;
Ga=100;
% Kernel J in Stage 6
sigJ1=3;
sigJ2=1;
Ja=75;
% Kernel K in Stage 6
sigK=4;
Ka=75;
% Kernel L in Stage 7
sigL1=10;
sigL2=3;
La=10;
% Kernel P in Stages 7 and 8
sigP=8;
Pa=1;
% Anisotropic kernels [G, J, and L] and offset isotropic kernel [K]
for d=1:numdir
theta=(d-1)*2*pi/numdir;
G(:,:,d)=anisot(Ga,theta,sigG1,sigG2);
J(:,:,d)=anisot(Ja,theta,sigJ1,sigJ2);
L(:,:,d)=anisot(La,theta,sigL1,sigL2);
K(:,:,d)=isot(Ka,sigK,2,theta);
end
% Isotropic kernel
P=isot(Pa,sigP,1,0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tic
for t=dt:dt:time
t=t+dt;
tn=tn+1;
Fnum=ceil(t/Fdur); % frame number in the motion stimulus
In=Stim(:,:,Fnum);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change-sensitive receptors: Non-directional transient cells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=x+dt*A1*(-B1*x+(1-x).*In);
z=z+dt*A2*(1-z-K2*x.*z);
b=x.*z;
B=rect(b-Tb); % thresholding and half-wave rectification
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Directional tansient cell network
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Directional inhibitory interneurons
c_old=c;
for d=1:numdir
D=modD(d+numdir/2,numdir); % opponent direction
theta=(d-1)*2*pi/numdir;
c(:,:,d)=c(:,:,d)+dt*A3*(-c(:,:,d)+C3*B-K3*rect(shift(c_old(:,:,D),theta)));
end
% Directional transient cells
for d=1:numdir
D=modD(d+numdir/2,numdir);
theta=(d-1)*2*pi/numdir;
e(:,:,d)=e(:,:,d)+dt*A4*(-e(:,:,d)+C4*B-K4*rect(shift(c_old(:,:,D),theta)));
end
E=rect(e-Te);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Directional short-range filters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fE=zeros(rows,cols,numdir); % excitation
for d=1:numdir
% Accumulation of directional evidence across a small spatial range
fE(:,:,d)=imfilter(E(:,:,d),G(:,:,d),'conv');
end
f=f+dt*A5*(-f+fE);
F=rect(f-Tf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spatial and opponent directional competition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hE=zeros(rows,cols,numdir); % excitation
hI=zeros(rows,cols,numdir); % inhibition
for d=1:numdir
D=modD(d+numdir/2,numdir);
% Opponent directional competition
h_D(:,:,d)=F(:,:,D);
% Spatial competition
hE(:,:,d)=imfilter(F(:,:,d),J(:,:,d),'conv'); % on-center
hI(:,:,d)=imfilter(F(:,:,d),K(:,:,d),'conv'); % off-surround
end
h=h+dt*A6*(-h+(1-h).*hE-(h+0.5).*(hI+D6*h_D));
H=rect(h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Directional long-range filters and Directional grouping
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mE=zeros(rows,cols,numdir); % excitation
for d=1:numdir
mE(:,:,d)=imfilter(H(:,:,d),L(:,:,d),'conv');
end
for d=1:numdir
yP=imfilter(rect(y),P,'conv');
FB_d=zeros(rows,cols);
for dd=1:numdir
FB_d=FB_d+InhWeights_MTMST(min(abs(d-dd),8-abs(d-dd))+1)*yP(:,:,dd);
end
FB(:,:,d)=FB_d; % Inhibitory feedback signals
end
% Model MT
m=m+dt*A7*(-m+(1-m).*mE-(m+0.5).*FB)+(sqrt(dt))*noise_sig_MT*randn(rows,cols,numdir);
M=rect(m);
% Model MST
y=y+dt*A8*(-y+C8*(1-y).*M-(y+0.5).*FB)+(sqrt(dt))*noise_sig_MST*randn(rows,cols,numdir);
Y=rect(y);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Saving the global directional estimates at each time instant from
% model areas MT and MST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
est_mt=[];
est_mst=[];
for d=1:numdir
est_mt=[est_mt sum(sum(M(:,:,d)))];
est_mst=[est_mst sum(sum(Y(:,:,d)))];
end
MT=[MT; est_mt];
MST=[MST; est_mst];
end
% toc
return
