===== Complete code for channel encoding model tutorial ===== %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Make simulated data %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % some variables iNeuron = 1; orientations = 0:179; k = 10; % loop over each neuron tuning function for orientPreference = 0:2:179 % compute the neural response as a Von Mises function %Note the 2 here which makes it so that our 0 - 180 orientation % space gets mapped to all 360 degrees neuralResponse(iNeuron,:) = exp(k*cos(2*pi*(orientations-orientPreference)/180)); % normalize to a height of 1 neuralResponse(iNeuron,:) = neuralResponse(iNeuron,:) / max(neuralResponse(iNeuron,:)); % update counter iNeuron = iNeuron + 1; end % plot the response of neuron 45 figure; plot(orientations,neuralResponse(45,:)); xlabel('Orientation'); ylabel('Channel response (normalized units to 1)'); % make a random weighting of neurons on to each voxel nNeurons = size(neuralResponse,1); nVoxels = 50; neuronToVoxelWeights = rand(nNeurons,nVoxels); % make stimulus array nStimuli = 8; % evenly space stimuli stimuli = 0:180/(nStimuli):179; % number of repeats nRepeats = 20; stimuli = repmat(stimuli,1,nRepeats); % round and make a column array stimuli = round(stimuli(:))+1; % compute the voxelResponse nTrials = nStimuli * nRepeats; for iTrial = 1:nTrials % get the neural response to this stimulus, by indexing the correct column of the neuralResponse matrix thisNeuralResponse = neuralResponse(:,stimuli(iTrial)); % multiply this by the neuronToVoxelWeights to get the voxel response on this trial. Note that you need % to get the matrix dimensions right, so transpose is needed on thisNeuralResponse voxelResponse(iTrial,:) = thisNeuralResponse' * neuronToVoxelWeights; end % plot the voxelResponse for the 7th trial figure; plot(voxelResponse(7,:)); xlabel('Voxel (number)'); ylabel('Voxel response (fake measurement units)'); % plot another trial voxel response figure; plot(voxelResponse(7,:),'b-.'); hold on plot(voxelResponse(7+nStimuli,:),'r-o'); xlabel('Voxel (number)'); ylabel('Voxel response (fake measurement units)'); % add noise to the voxel responses noiseStandardDeviation = 0.05; % normalize response voxelResponse = voxelResponse / mean(voxelResponse(:)); % add gaussian noise voxelResponse = voxelResponse + noiseStandardDeviation*randn(size(voxelResponse)); % check the voxelResponses figure; stim1 = 7; stim2 = 3; subplot(1,3,1); plot(voxelResponse(stim1,:),'b-.'); hold on plot(voxelResponse(stim1+nStimuli,:),'r-o'); xlabel('Voxel (number)'); ylabel('Voxel response (fake measurement units)'); subplot(1,3,2); plot(voxelResponse(stim1,:),voxelResponse(stim1+nStimuli,:),'k.'); xlabel('Response to first presentation'); ylabel('Response to second presentation'); axis square subplot(1,3,3); plot(voxelResponse(stim1,:),voxelResponse(stim2,:),'k.'); xlabel(sprintf('Response to stimulus: %i deg',stimuli(stim1))); ylabel(sprintf('Response to stimulus: %i',stimuli(stim2))); axis square %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Make encoding model %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % make channel basis functions nChannels = 8; exponent = 7; orientations = 0:179; prefOrientation = 0:180/nChannels:179; % loop over each channel for iChannel = 1:nChannels % get sinusoid. Note the 2 here which makes it so that our 0 - 180 orientation % space gets mapped to all 360 degrees thisChannelBasis = cos(2*pi*(orientations-prefOrientation(iChannel))/180); % rectify thisChannelBasis(thisChannelBasis<0) = 0; % raise to exponent thisChannelBasis = thisChannelBasis.^exponent; % keep in matrix channelBasis(:,iChannel) = thisChannelBasis; end % plot channel basis functions figure; plot(orientations,channelBasis); xlabel('Preferred orientation (deg)'); ylabel('Ideal channel response (normalized to 1)'); % compute the channelResponse for each trial for iTrial = 1:nTrials channelResponse(iTrial,:) = channelBasis(stimuli(iTrial),:); end % compute estimated weights estimatedWeights = pinv(channelResponse) * voxelResponse; % compute model prediction modelPrediction = channelResponse * estimatedWeights; % compute residual residualResponse = voxelResponse-modelPrediction; % compute r2 r2 = 1-var(residualResponse(:))/var(voxelResponse(:)) %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Inverted encoding model %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % split half into train and test firstHalf = 1:round(nTrials/2); secondHalf = round(nTrials/2)+1:nTrials; trainVoxelResponse = voxelResponse(firstHalf,:); testVoxelResponse = voxelResponse(secondHalf:end,:); % compute weights on train data estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse; % compute channel response from textVoxelResponses estimatedChannelResponse = testVoxelResponse * pinv(estimatedWeights); % plot channel responses figure;colors = hsv(nStimuli); for iStimuli = 1:nStimuli plot(prefOrientation,mean(estimatedChannelResponse(iStimuli:nStimuli:end,:),1),'-','Color',colors(iStimuli,:)); hold on end xlabel('Channel orientation preference (deg)'); ylabel('Estimated channel response (percentile of max)'); title(sprintf('r2 = %0.4f',r2)); % Compute voxel response without noise nTrials = nStimuli * nRepeats; for iTrial = 1:nTrials % get the neural response to this stimulus, by indexing the correct column of the neuralResponse matrix thisNeuralResponse = neuralResponse(:,stimuli(iTrial)); % multiply this by the neuronToVoxelWeights to get the voxel response on this trial. Note that you need % to get the matrix dimensions right, so transpose is needed on thisNeuralResponse voxelResponseNoisy(iTrial,:) = thisNeuralResponse' * neuronToVoxelWeights; end % add noise noiseStandardDeviation = 0.5; % normalize response voxelResponseNoisy = voxelResponseNoisy / mean(voxelResponseNoisy(:)); % add gaussian noise voxelResponseNoisy = voxelResponseNoisy + noiseStandardDeviation*randn(size(voxelResponseNoisy)); % split into train and tes trainVoxelResponseNoisy = voxelResponseNoisy(firstHalf,:); testVoxelResponseNoisy = voxelResponseNoisy(secondHalf:end,:); % compute weights on train data estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponseNoisy; % compute model prediction on test data modelPrediction = channelResponse(secondHalf,:) * estimatedWeights; % compute residual residualResponse = testVoxelResponseNoisy-modelPrediction; % compute r2 r2 = 1-var(residualResponse(:))/var(testVoxelResponseNoisy(:)) % invert model and compute channel response estimatedChannelResponse = testVoxelResponseNoisy * pinv(estimatedWeights); % plot estimated channel profiles figure;colors = hsv(nStimuli); for iStimuli = 1:nStimuli plot(prefOrientation,mean(estimatedChannelResponse(iStimuli:nStimuli:end,:),1),'-','Color',colors(iStimuli,:)); hold on end xlabel('Channel orientation preference (deg)'); ylabel('Estimated channel response (percentile of max)'); title(sprintf('r2 = %0.4f',r2)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % stimulus likelihood function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % split half into train and test firstHalf = 1:round(nTrials/2); secondHalf = round(nTrials/2)+1:nTrials; trainVoxelResponse = voxelResponse(firstHalf,:); testVoxelResponse = voxelResponse(secondHalf:end,:); % compute weights on train data estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse; % compute model prediction on test data modelPrediction = channelResponse(secondHalf,:) * estimatedWeights; % compute residual residualResponse = testVoxelResponseNoisy-modelPrediction; % compute residual variance, note that this is a scalar residualVariance = var(residualResponse(:)); % make this into a covariance matrix in which the diagonal contains the variance for each voxel % and off diagonals (in this case all 0) contain covariance between voxels modelCovar = eye(nVoxels)*residualVariance; % cycle over each trial nTestTrials = size(testVoxelResponse,1); for iTrial = 1:nTestTrials % now cycle over all possible orientation for iOrientation = 1:179 % compute the mean voxel response predicted by the channel encoding model predictedResponse = channelBasis(iOrientation,:)*estimatedWeights; % now use that mean response and the model covariance to estimate the probability % of seeing this orientation given the response on this trial likelihood(iTrial,iOrientation) = mvnpdf(testVoxelResponse(iTrial,:),predictedResponse,modelCovar); end end figure for iStimuli = 1:nStimuli plot(1:179,mean(likelihood(iStimuli:nStimuli:end,:),1),'-','Color',colors(iStimuli,:)); hold on end xlabel('stimulus orientation (deg)'); ylabel('probability given trial response'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Inverted Encoding model with different channel basis functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % reweight the channels channelReweighting = [0 0.8 0.4 0 0 0 0.4 0.8]'; % make into a full matrix xform to transform the original channels for iChannel = 1:nChannels xform(iChannel,:) = circshift(channelReweighting,iChannel-1); end % and get new bimodal channels bimodalChannelBasis = channelBasis * xform; % display a figure with one of the channels figure plot(orientations,bimodalChannelBasis(:,5)); xlabel('orientation (deg)'); ylabel('Channel response (normalized to 1)'); % compute the channelResponse for each trial for iTrial = 1:nTrials channelResponse(iTrial,:) = bimodalChannelBasis(stimuli(iTrial),:); end % compute estimated weights estimatedWeights = pinv(channelResponse) * voxelResponse; % compute model prediction modelPrediction = channelResponse * estimatedWeights; % compute residual residualResponse = voxelResponse-modelPrediction; % compute r2 r2 = 1-var(residualResponse(:))/var(voxelResponse(:)) % compute estimated channel response profiles estimatedChannelResponse = testVoxelResponse * pinv(estimatedWeights); % and plot one of the channels averaged across all trials figure; plot(prefOrientation,mean(estimatedChannelResponse(5:nStimuli:end,:),1)); xlabel('Channel preferred orientation (deg)'); ylabel('Estimated channel response (percentile of full)'); title(sprintf('r2 = %0.4f',r2)); % compute weights on train data estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse; % compute model prediction on test data modelPrediction = channelResponse(secondHalf,:) * estimatedWeights; % compute residual residualResponse = testVoxelResponseNoisy-modelPrediction; % compute residual variance, note that this is a scalar residualVariance = var(residualResponse(:)); % make this into a covariance matrix in which the diagonal contains the variance for each voxel % and off diagonals (in this case all 0) contain covariance between voxels modelCovar = eye(nVoxels)*residualVariance; % cycle over each trial nTestTrials = size(testVoxelResponse,1); for iTrial = 1:nTestTrials % now cycle over all possible orientation for iOrientation = 1:179 % compute the mean voxel response predicted by the channel encoding model predictedResponse = bimodalChannelBasis(iOrientation,:)*estimatedWeights; % now use that mean response and the model covariance to estimate the probability % of seeing this orientation given the response on this trial likelihood(iTrial,iOrientation) = mvnpdf(testVoxelResponse(iTrial,:),predictedResponse,modelCovar); end end % now plot the likelihood function averaged over repeats figure for iStimuli = 1:nStimuli plot(1:179,mean(likelihood(iStimuli:nStimuli:end,:),1)); hold on end xlabel('stimulus orientation (deg)'); ylabel('probability given trial response');