function runRTtask(Motion,coherences,period,A9,B9,g_f,g_h,lambda,sigLIP, ...
pre_y,muu,n_f,n_h,Tdec,Tsacc,g_BG,TD,varargin)
% runRTtask :: Function, which is called by runTask.m, to run the reaction time task (Roitman & Shadlen, 2002)
%
%% Input variables
% Motion :: 5D matrix that stores the dynamics of population activity of
% middle temporal and medial superior temporal neurons
% [1: time steps, 2: preferred direction, 3: coherence level, 4: trial, 5:
% middle temporal or medial superior temporal]
%
% coherences :: Coherence levels between 0 and 100%
%
% period :: Duration for which recordings were made in model middle
% temporal and medial superior temporal areas
%
% A9 :: Parameter that scales passive decay in lateral intraparietal
% neuronal responses
%
% B9 :: Maximal rate of the model lateral intraparietal neurons
%
% g_f :: Parameter that scales the gain of self-excitation in lateral
% intraparietal recurrent competitive field
%
% g_h :: Parameter that scales the gain of recurrent inhition in lateral
% intraparietal recurrent competitive field
%
% lambda :: Parameter that scales the bottom-up excitation to a lateral
% intraparietal cellfrom the media superior temporal pool tuned to the
% preferred direction
%
% sigLIP :: Standard deviation of the Brownian motion process that injects
% noise into the dynamics of lateral intraparietal area neurons
%
% pre_y :: Pre-motion onset period rate of each model lateral intraparietal
% neuron whose receptive field contains a choice target
%
% muu :: Threshold parameter of the recurrent signal functions in the
% lateral intraparietal recurrent competitive field
%
% n_f :: Steepness parameter of the self-excitatory signal function in the
% lateral intraparietal recurrent competitive field
%
% n_h :: Steepness parameter of the recurrent inhibitory signal function in
% the lateral intraparietal recurrent competitive field
%
% Tdec :: Threshold lateral intraparietal activity for making decision in
% the reaction time task
%
% Tsacc :: Rate of the selected lateral intraparietal cell when the saccade
% is initiated
%
% g_BG :: Parameter that scales the gain of self-excitation for the selected
% lateral intraparietal call after decision threshold-crossing
%
% TD :: 1 if top-down, volitional 'forced choice' mechanism is in place; 0
% otherwise
%
%% 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)
%
%% 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
%%
% fixed time step of numerical integration
dt=0.001; % (in sec)
cohL=size(Motion,3);
numtrials=size(Motion,4);
num=length(0:dt:period); % number of time steps
RTc=zeros(numtrials,cohL); % correct trial RTs
Y1c=zeros(num,numtrials,cohL); % correct trial LIP dynamics of recorded population
Y2c=zeros(num,numtrials,cohL); % correct trial LIP dynamics of the other population
RTe=zeros(numtrials,cohL); % error trial RTs
Y1e=zeros(num,numtrials,cohL); % error trial LIP dynamics of recorded population
Y2e=zeros(num,numtrials,cohL); % error trial LIP dynamics of the other population
Acc=zeros(1,cohL); % tracks the number of correct responses
FoCh=zeros(2,cohL); % tracks the number of forced choices (if any)
for coh=1:cohL
coh % to track the progress
for trial=1:numtrials
% trial % to track the progress
% To set the pseudo-random number generator to a random state to
% begin with
rand('state',sum(100*clock))
randn('state',sum(100*clock))
% Note stimulus grid is 60 x 60
x1=Motion(:,1,coh,trial,2)/(60*60); % average rightward MST cell activity
x2=Motion(:,5,coh,trial,2)/(60*60); % average leftward MST cell activity
% run model LIP module for RT task
[Y1,Y2,RT,Ch,FC]=MODE_lip_RT(x1,x2,A9,B9,g_f,g_h,lambda,sigLIP,pre_y,muu,n_f,n_h,Tdec,Tsacc,g_BG,TD);
if Ch==1
RTc(trial,coh)=RT;
Y1c(:,trial,coh)=Y1;
Y2c(:,trial,coh)=Y2;
% Remember that the correct stimulus direction in the simulations
% is always rightward (d=1), for simplicity
Acc(coh)=Acc(coh)+1;
elseif Ch==2
RTe(trial,coh)=RT;
Y1e(:,trial,coh)=Y1;
Y2e(:,trial,coh)=Y2;
end
if FC==1
FoCh(1,coh)=FoCh(1,coh)+1;
elseif FC==2
FoCh(2,coh)=FoCh(2,coh)+1;
end
if nargin > 0
hwait = varargin{1};
waitbar((1+(trial+(coh-1)*numtrials))/(cohL*numtrials + 1), ...
hwait, 'Run in progress . . .');
end
end
end
% First the behavioral data is analyzed
[Beh]=analyzeRT(RTc,RTe);
% Convert coherence levels into log10 scale for usage in figures
log10cohs=max(log10(coherences+eps),-0.1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Chronometric function
figure
hold on
errorbar(log10cohs,Beh(:,3),Beh(:,4),'b.')
errorbar(log10cohs,Beh(:,5),Beh(:,6),'ro')
axis([-0.15 log10(100) 400 1100])
hold off
xlabel('Motion strength (% coherence)','Fontsize',15)
ylabel('Reaction time (ms)','Fontsize',15)
title('Chronometric function (RT as a function of coherence and trial type: correct or error) for the RT task','Fontsize',15)
set(gca,'Fontsize',12)
box off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Psychometric function
% clear global x_ml t_ml n_xl
%
% [w_ml]=ml_e(coherences',Acc'/numtrials,numtrials);
% % The first two inputs need to be column vectors
%
% Xcoh=0.7:0.1:10; Xcoh=[Xcoh 11:100]; % in units of percentage
% Yperf=1-0.5*exp(-(Xcoh/w_ml(1,1)).^w_ml(2,1));
%
% % Fitted parameters of the psychometric function
% alpha=w_ml(1,1) % coherence threshold
% beta=w_ml(2,1) % parameter that controls the steepness of the psychometric function
figure
% plot(log10(Xcoh),100*Yperf,'k')
hold on
plot(log10cohs,Acc*100/numtrials,'ko-')
hold off
xlabel('Motion strength (% coherence)','Fontsize',15)
ylabel('Accuracy (% correct)','Fontsize',15)
title('Psychometric function (Accuracy as a function of coherence)','Fontsize',15)
set(gca,'Fontsize',12)
box off
axis([-0.15 log10(52) 40 100])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Average LIP responses during correct trials at all coherences
% 'tn' signifies time-step; 'c': correct; 'e': error; 'F': forward; 'B':
% backward
[Y1c_lhs,Y2c_lhs,tnc_Fextent]=computeLIP_lhs(RTc,Y1c,Y2c,cohL,numtrials,dt);
[Y1e_lhs,Y2e_lhs,tne_Fextent]=computeLIP_lhs(RTe,Y1e,Y2e,cohL,numtrials,dt);
RTc_max=max(max(RTc));
RTe_max=max(max(RTe));
tn_extent=ceil(max(RTc_max,RTe_max)+100);
[Y1c_rhs,Y2c_rhs,tnc_Bextent]=computeLIP_rhs(RTc,Y1c,Y2c,cohL,numtrials,dt,tn_extent);
[Y1e_rhs,Y2e_rhs,tne_Bextent]=computeLIP_rhs(RTe,Y1e,Y2e,cohL,numtrials,dt,tn_extent);
span=1+60/(dt*1000); % As in the data, the average responses are smoothed with a 60 ms running mean
colors=['bgrcmy'];
figure
hold on
for coh=1:cohL
if max(RTc(:,coh))~=0 % to check if at least one correct trial at a coherence level has occurred
% Left portion
plot(smooth(Y1c_lhs(1:tnc_Fextent(coh),coh),span),colors(coh)),plot(smooth(Y2c_lhs(1:tnc_Fextent(coh),coh),span),[colors(coh) '--'])
% Right portion
plot(800+[tn_extent-tnc_Bextent(coh):tn_extent],smooth(Y1c_rhs(tn_extent-tnc_Bextent(coh):tn_extent,coh),span),colors(coh))
plot(800+[tn_extent-tnc_Bextent(coh):tn_extent],smooth(Y2c_rhs(tn_extent-tnc_Bextent(coh):tn_extent,coh),span),[colors(coh) '--'])
end
end
box off
hold off
axis([0 ceil((tn_extent+800)/(dt*1000)) 20 70])
set(gca,'Fontsize',12)
xlabel('Time (ms)','Fontsize',15)
ylabel('Activity','Fontsize',15)
title('Average LIP responses during correct trials at all coherence levels','Fontsize',15)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Correct and error trial average LIP responses at a given coherence
coh=3; % select coherence level here
% error trial average responses are shown in black
figure
hold on
% Correct trials
if max(RTc(:,coh))~=0
% Left portion
plot(smooth(Y1c_lhs(1:tnc_Fextent(coh),coh),span),colors(coh)),plot(smooth(Y2c_lhs(1:tnc_Fextent(coh),coh),span),[colors(coh) '--'])
% Right portion
plot(800+[tn_extent-tnc_Bextent(coh):tn_extent],smooth(Y1c_rhs(tn_extent-tnc_Bextent(coh):tn_extent,coh),span),colors(coh))
plot(800+[tn_extent-tnc_Bextent(coh):tn_extent],smooth(Y2c_rhs(tn_extent-tnc_Bextent(coh):tn_extent,coh),span),[colors(coh) '--'])
end
% Error trials
if max(RTe(:,coh))~=0 % to check if at least one error trial at the coherence level has occurred
% Left portion
plot(smooth(Y1e_lhs(1:tne_Fextent(coh),coh),span),'k-'),plot(smooth(Y2e_lhs(1:tne_Fextent(coh),coh),span),'k--')
% Right portion
plot(800+[tn_extent-tne_Bextent(coh):tn_extent],smooth(Y1e_rhs(tn_extent-tne_Bextent(coh):tn_extent,coh),span),'k-')
plot(800+[tn_extent-tne_Bextent(coh):tn_extent],smooth(Y2e_rhs(tn_extent-tne_Bextent(coh):tn_extent,coh),span),'k--')
end
box off
hold off
axis([0 ceil((tn_extent+800)/(dt*1000)) 20 70])
set(gca,'Fontsize',12)
xlabel('Time (ms)','Fontsize',15)
ylabel('Activity','Fontsize',15)
title('Correct and error trial average LIP responses at 6.4% coherence for which the model monkey makes errors','Fontsize',15)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% RT distributions
bincenters=400:50:2500;
% Correct trials
figure
% title('Correct trial RT distributions','Fontsize',12)
for coh=1:cohL
if max(RTc(:,coh))~=0
subplot(cohL,2,2*(coh-1)+1),hist(nonzeros(RTc(:,coh)),bincenters,'k'), box off,axis([400 2300 0 numtrials/2])
if coh==1
title('Correct trial RT distributions at various coherences for the RT task','Fontsize',15)
end
box off
set(gca,'Fontsize',12)
xlabel('RT (ms)','Fontsize',10)
ylabel('No. of trials','Fontsize',10)
end
end
% Error trials
figure
% title('Error trial RT distributions','Fontsize',12)
for coh=1:cohL
if max(RTe(:,coh))~=0
subplot(cohL,2,2*(coh-1)+1),hist(nonzeros(RTc(:,coh)),bincenters,'k'), box off,axis([400 2300 0 numtrials/2])
if coh==1
title('Error trial RT distributions at various coherences for the RT task','Fontsize',15)
end
box off
set(gca,'Fontsize',12)
xlabel('RT (ms)','Fontsize',10)
ylabel('No. of trials','Fontsize',10)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Code to compute the range of reaction times in both correct and error
%%% trials at each coherence level
RT_range=NaN*ones(cohL,4); % each row represents a coherence: 1) minRT (correct trial) 2) maxRT (correct trial) 3) minRT (error trial) 4) maxRT (error trial)
for coh=1:cohL
if max(RTc(:,coh))~=0
RT_range(coh,1:2)=[min(nonzeros(RTc(:,coh))) max(nonzeros(RTc(:,coh)))];
end
if max(RTe(:,coh))~=0
RT_range(coh,3:4)=[min(nonzeros(RTe(:,coh))) max(nonzeros(RTe(:,coh)))];
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Figure that plots the % of correct trials in which RT is > 1 s
RTc_more_1s=100*sum(RTc>1000)./sum(RTc~=0);
figure
plot(log10cohs,RTc_more_1s,'o')
axis([-0.15 log10(100) 0 35])
xlabel('Motion strength (% coherence)','Fontsize',15)
ylabel('% of correct trials with RT > 1s','Fontsize',15)
box off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Figure that plots the % of trials in which 'forced choices' were made,
%%% and if they turned out to be correct or not
% FoCh
figure
plot(log10cohs,FoCh(1,:)*100/numtrials,'s')
hold on
plot(log10cohs,FoCh(2,:)*100/numtrials,'d')
plot(log10cohs,sum(FoCh,1)*100/numtrials,'o')
axis([-0.15 log10(100) 0 15])
xlabel('Motion strength (% coherence)','Fontsize',15)
ylabel('% of forced trials','Fontsize',15)
box off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% To plot individual trial LIP responses
trial=1; % choose trial here
coh=2; % choose coherence here
figure
hold on
if RTc(trial,coh)~=0
plot(Y1c(:,trial,coh),'r')
plot(Y2c(:,trial,coh),'r--')
% plots decision threshold line
plot([0 period*1000],[55 55],'k--') % period is in secs
else
plot(Y1e(:,trial,coh),'r')
plot(Y2e(:,trial,coh),'r--')
% plots decision threshold line
plot([0 period*1000],[55 55],'k--') % period is in secs
end
hold off
xlabel('Time (ms)','Fontsize',15)
ylabel('Activity','Fontsize',15)
title('An individual trial LIP responses for 3.2% coherence','Fontsize',15)
set(gca,'Fontsize',12)
box off
axis([0 period/dt 0 70])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Figure 10 in the article
%%% Relationship between LIP response and reaction time
% Warning: Figure 10 makes sense only if there are atleast 2 correct trials
% at a coherence level
LongShort
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return
