13-Oct-2015 18:47:02 function [SpatialCorrelates] = spatialcorrelates(model_name,NeuralDir, BehaviorDir, Trigger, exclusions, options) % SPATIALCORRELATES % % [SpatialCorrelates] = spatialcorrelates(NeuralDir, BehaviorDir, Trigger, exlcusions, options) % % This routine calculates spatial correlates between neural and behavioral % data. These are: 1) Place Fields, 2) Horizontal Head Direction, and 3) % Vertical Head Direction. This processing assumes that the behavioral data % are formatted according to the Motion capture System framework, and that % the Neural Data are formatted according to the Buzsaki lab frame (.res % and .clu files). These specifications are a little fragile, and are % detailed below. % % author: John D. Long II, PhD contact: jlong29@gmail.com % %%%%%%%%%%%%%% %%% INPUTS %%% %%%%%%%%%%%%%% % model_name: model name of the kinematic model associated w/subject % e.g. 'buddy140_km' % NeuralDir: a string designating the absolute path to the neural data % directory. No trailing filesep. % BehaviorDir: a string designating the absolute path to the behavioral % data directory. No trailing filesep. % Trigger: a vector at the lfp sampling rate containing camera % triggers used to align the neural data to the behavioral % data. % exclusions: a structure with the fields: % -triggers: to exclude erroneous trigger blocks e.g. % accidental camera acquisition % -trials: to exclude trials e.g. too short, animal % problem % -shanks: to exclude bad shanks % e.g. [1 3 7 8] i.e. lists of integer indices % options: an optional structure allowing the user to overwrite any % free parameters set in the routine. Just set a field of % "options" to the variables name with an appropriate value. %%%%%%%%%%%%%%% %%% OUTPUTS %%% %%%%%%%%%%%%%%% % SpatialCorrelates: a hierarchical structure with fields: % params: is a structure of input parameters % model_name : the name of the kinematic model used % neuralpath : the directory in which the neural data resides % behaviorpath : the directory in which the behavioral data resides % filebase : the base name of the neural files used % HDres : resolution for binning head direction, in degrees % SpatialScale : sets the interpolation factor for SpatialHistogram % SRspikes : spike sampling rate % SReeg : eeg, trigger channel, sampling rate % thresh : threshold for camera triggers % fps : frame rate of the cameras % excludeshanks : shanks excluded from the analysis % excludetriggers : triggers excluded from the analysis % theta : bins used for rotation about z-axis for HD % phi : bins used for rotation about y-axis for HD % distfactor : for scaling distances in environment to a standard % reward_info : reward locations in (SpatialHistogramBinNumber,X,Y,Z) % % shared: a structure detailing per trial/behavioral variables % Placetime: place occupancy of the subject % HDtime: occupancy of the subject, conditional upon head direction % HDts: time-series of hd direction in 3D (theta,phi) % trajectory: trajectory in the local 3D coordinate frame % location: time-series of location in 2D environment histogram coordinates % speed: speed of the subject in centimeters per second % TrialID: trial id, which might be different from consecutive indices % timeline: per trial time stamps in seconds (this is padded a bit at the end of the trial, but everything is aligned from the start) % CoordFrame: the local coordinate system, relative to the environment % Transl_g: transformation from the local to the reference coord frame % % pershank: a structure detailing per shank variables % ShankID: formatted xx 2-digit [1-99] shankID % UnitIds: formatted xxyyy with xx being shankID and yyy being unit id % UnitInfo: A cell array {unit}{trial} with a bottom level matrix: % [#CameraFramesContainingSpikes x (spike count,timestamp,spatialbin)] % UnitFR: firing rates within each trial and across session % UnitVar: std of unit firing rates within each trial and across session % PlaceFrMap: spike counts as a function of spatial location % This is a structure that contains location and orientation % parameters also used by FovFrMap % HDFrmap: spike counts as a function of head direction, 3D % %%%%%%%%%%%%%%%%%%%% %%% INPUT CHECKS %%% %%%%%%%%%%%%%%%%%%%% if isdir(NeuralDir) % check for and access expected res and clu files resfiles = dir(sprintf('%s%s*.res.*',NeuralDir,filesep)); if isempty(resfiles) fprintf(1,'FileNotFound: there are no res files in the NeuralDir.\n'); fprintf(1,'\tCheck the ''NeuralDir'' directory.\n'); return else clufiles = dir(sprintf('%s%s*.clu.*',NeuralDir,filesep)); %check to make sure the number of res and clue files are the same if length(resfiles) ~= length(clufiles) fprintf(1,'CountError: the number of res and clu files is\n'); fprintf(1,'\tnot equal. Check the directories.\n'); return end end else fprintf(1,'PathError: the path ''NeuralDir'' input by the user does\n'); fprintf(1,'\tnot exist. Confirm path and then try again.\n'); return end [~,filebase,~] = fileparts(NeuralDir); if isdir(BehaviorDir) % Store the directory in which this function was launched CurrentDir = pwd; cd(BehaviorDir) % check for kinematic model frame tmp = dir(sprintf('..%s%s_frame.m',filesep,model_name)); if isempty(tmp) error('Error: No kinematic frame was found for this subject.'); end seps = regexp(tmp.name,'_'); % cutting off _km_frame subject = tmp.name(1:seps(1)-1); % cutting off _frame modelName = tmp.name(1:seps(2)-1); fprintf(1,'\tSubject: %s\n',subject); fprintf(1,'\tModel Name: %s\n',modelName); % check if kinematic data has been loaded and saved before chk = dir([modelName '*_Global_kmData.mat']); if isempty(chk) try kmData = loadKinematicData(modelName,'Data'); catch DamnErr fprintf(1,'\t\tError: %s\n',DamnErr.message); return; end else fprintf(1,'\tLoading %s found in %s.\n',chk(1).name, BehaviorDir); tmp = load(chk(1).name); kmData = tmp.kmData; end % load data_frame for nodeLabels load(sprintf('..%s%s.mat',filesep,model_name)) cd(CurrentDir) % check for and access rend_env files rend_envfiles = dir(sprintf('%s%s*_rend_env.mat',BehaviorDir,filesep)); if isempty(rend_envfiles) fprintf(1,'FileNotFound: there are no rend_env.mat files in the NeuralDir.\n'); fprintf(1,'\tCheck the BehaviorDir path supplied.\n'); return end else fprintf(1,'PathError: the path ''BehaviorDir'' input by the user does\n'); fprintf(1,'\tnot exist. Confirm path and then try again.\n'); return end % check trigger input if nargin < 3 || isempty(Trigger) fprintf(1,'InputError: the user must supply a ''Trigger''\n'); fprintf(1,'\tfor aligning neural and camera data.\n'); return else if isnumeric(Trigger) if size(Trigger,2) > size(Trigger,1) Trigger = Trigger'; if size(Trigger,2) > 1 fprintf(1,'InputError: the intput ''Trigger'' must be a vector.\n'); fprintf(1,'\tThe user input an array.\n'); return end end else fprintf(1,'InputError: ''Triggers'' must be a vector of numbers.\n'); return end end % check for any exclusions if nargin < 4 || isempty(exclusions) excludeshanks = []; excludetriggers = []; excludetrials = []; else if isfield(exclusions,'shanks') excludeshanks = exclusions.shanks; else fprintf(1,'InputWarning: the structure ''exclusions'' does not have\n'); fprintf(1,'\tthe field ''shanks''. Using th Default; [].'\n'); excludeshanks = []; end if isfield(exclusions,'triggers') excludetriggers = exclusions.triggers; else fprintf(1,'InputWarning: the structure ''exclusions'' does not have\n'); fprintf(1,'\tthe field ''triggers''. Using th Default; [].\n'); excludetriggers = []; end if isfield(exclusions,'trials') excludetrials = exclusions.trials; else fprintf(1,'InputWarning: the structure ''exclusions'' does not have\n'); fprintf(1,'\tthe field ''trials''. Using th Default; [].\n'); excludetrials = []; end end %%%%%%%%%%%%%%%%%%%%%%%%%%% %%% SET Free Parameters %%% %%%%%%%%%%%%%%%%%%%%%%%%%%% SpeedThresh = 0; % for place cell firing (cm/second); HDres = 10; % resolution for binning head direction, in degrees SpatialScale = 4; % sets the interpolation factor for SpatialHistogram (Default --> 2x2cm) SRspikes = 20000; % spike sampling rate SReeg = 1250; % eeg, trigger channel, sampling rate thresh = 3000; % threshold for camera triggers fps = 50; % frame rate of the camera (on average) factor = 10; % multiplier to determine trial boundaries OPTIONS MaxSpikeBytes = 5*1024^3;% maximum size used to determine processing branch % Check for options structure for overwriting free parameters if nargin > 5 && ~isempty(options) v2struct(options) end % store parameter settings for this run params.model_name = model_name; params.neuralpath = NeuralDir; params.behaviorpath = BehaviorDir; params.filebase = filebase; params.SpeedThresh = SpeedThresh; params.HDres = HDres; params.SpatialScale = SpatialScale; params.SRspikes = SRspikes; params.SReeg = SReeg; params.thresh = thresh; params.fps = fps; params.excludeshanks = excludeshanks; params.excludetriggers = excludetriggers; %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Camera Triggers Pre-processing %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nsamp = ceil(SReeg/fps); % number of eeg samples per frame out = suprathresh(Trigger,thresh); % find camera triggers interTrial = find(diff(out(:,1))>(Nsamp*factor))+1; Ntrials = length(interTrial)+1; Trials = cell(Ntrials,1); Nframes = zeros(length(Trials),1); %frames per trial % calculate the timestamp of each frame in each trial at the camera frame % rate, plus some padding for the final frame in a trial for ii = 1:Ntrials if ii < 2 Trials{ii} = [out(1:interTrial(ii)-1,1); ... out(interTrial(ii)-1,1)+Nsamp]/SReeg; elseif ii < Ntrials Trials{ii} = [out(interTrial(ii-1)+1:interTrial(ii)-1,1); ... out(interTrial(ii)-1,1)+Nsamp]/SReeg; else Trials{ii} = [out(interTrial(ii-1)+1:end,1); ... out(end)+Nsamp]/SReeg; end % the number of frames this trial Nframes(ii) = length(Trials{ii}); end % Exclusions: first eliminate bad trigger blocks, then trials Trials(excludetriggers) = []; Nframes(excludetriggers) = []; TrialsList = 1:length(Trials); Trials(excludetrials) = []; Nframes(excludetrials) = []; TrialsList(excludetrials) = []; Ntrials = length(TrialsList); clear Trigger out % Use the first environment rendering to construct the 2d grid for planar % mapping cd(BehaviorDir) chk = dir([BehaviorDir filesep 'Data*_rend_env.mat']); tmp = strfind(chk(1).name,'_'); trialName = chk(1).name(1:tmp(1)-1); [table, ~, template, distfactor] = SpatialHistogram(trialName,SpatialScale); params.SpatialRes = 1/SpatialScale; params.EnvTemplate = template; cd(CurrentDir) % use template to set dimensions of all other histograms Dim = size(template,1); %%%% END: Camera Trigger Pre-processing %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Neural Data Pre-processing %%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % exclude bad shanks % load spike data and populate unit IDs Nshanks = length(resfiles); bad_shanks = false(length(Nshanks),1); for ii = 1:Nshanks shankID = strfind(resfiles(ii).name,'.'); shankID = str2double(resfiles(ii).name(shankID(end)+1:end)); if ismember(shankID,excludeshanks) bad_shanks(ii) = true; end end resfiles(bad_shanks) = []; clufiles(bad_shanks) = []; Nshanks = length(resfiles); % file size check for processing SpikeBytes = sum([cat(2,resfiles(:).bytes); cat(2,clufiles(:).bytes)]); %%% END: Neural Data Pre-processing %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Initialize Variables: shared and pershank %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DECLARE: Shared PER TRIAL Variables: % Placetime: place occupancy of the subject % planarGrid: the complete 2D planar grid used for mapping from 3 to 2D % HDtime: occupancy of the subject, conditional upon head direction % HDts: time-series of hd direction in 3D (theta,phi) % trajectory: trajectory in the local 3D coordinate frame % location: time-series of location in 2D environment histogram coordinates % speed: speed of the subject in centimeters per second % TrialID: trial id, which might be different from consecutive indices % timeline: per trial time stamps in seconds at camera frame rate % CoordFrame: the local coordinate system, relative to the environment % Transl_g: transformation from the local to the reference coord frame shared = struct('Placetime', cell(Ntrials,1), ... 'planarGrid', cell(Ntrials,1), ... 'HDtime', cell(Ntrials,1), ... 'HDts', cell(Ntrials,1), ... 'trajectory', cell(Ntrials,1), ... 'location', cell(Ntrials,1), ... 'speed', cell(Ntrials,1), ... 'TrialID', cell(Ntrials,1), ... 'timeline', cell(Ntrials,1), ... 'CoordFrame', cell(Ntrials,1), ... 'Transl_g', cell(Ntrials,1)); %%% INITIALIZE: Shared PER TRIAL Variables: Ntheta = ceil(360/HDres); % quantization of 3D head direction Nphi = ceil(180/HDres); Placetime = zeros(Dim,Dim,Ntrials); % 2D environment histogram planarGrid = cell(Ntrials,1); trajectory = cell(Ntrials,1); HDtime = struct('Theta', cell(Ntrials,1), ... 'Phi', cell(Ntrials,1)); HDts = cell(Ntrials,1); location = cell(Ntrials,1); % 2D environment histogram coordinates speed = cell(Ntrials,1); % cm/s CoordFrame = zeros(3,3,Ntrials); Transl_g = zeros(4,4,Ntrials); params.theta = linspace(-pi,pi,Ntheta); % rotation about z-axis params.phi = linspace(-pi/2,pi/2,Nphi); % out of plane rotation % initialize time-series variables for ii = 1:Ntrials HDts{ii} = nan(Nframes(ii),2); % (theta,phi) location{ii} = nan(Nframes(ii),1); speed{ii} = nan(Nframes(ii),1); % upper bound on model fits trajectory{ii} = nan(Nframes(ii),3); HDtime(ii).Theta = zeros(Ntheta,1); HDtime(ii).Phi = zeros(Nphi,1); end % DECLARE: pershank variables: % ShankID: formatted xx 2-digit [1-99] shankID % UnitIds: formatted xxyyy with xx being shankID and yyy being unit id % UnitInfo: A cell array {unit}{trial} with a bottom level matrix: % FramesContainingUnitSpikes x (spike count, timestamp,pixellocation) % UnitFR: firing rates within each trial and across session % UnitVar: std of unit firing rates within each trial and across session % PlaceFrMap: spike counts as a function of spatial location % This is a structure that contains location and orientation % parameters also used by FovFrMap % HDFrmap: spike counts as a function of head direction, 3D pershank = struct('ShankID', cell(Nshanks,1), ... 'UnitIds', cell(Nshanks,1), ... 'UnitInfo', cell(Nshanks,1), ... 'UnitFR', cell(Nshanks,1), ... 'UnitSTD', cell(Nshanks,1), ... 'PlaceFrMap', cell(Nshanks,1), ... 'HDFrMap', cell(Nshanks,1)); % DECLARE: pershank sub-variables % Note: Field-of-View cells inherent location and orientation parameters % from place field estimator % UnitInfo and associated ancillary variables UnitInfo = cell(Nshanks,1); UnitInfoTemp = cell(Nshanks,1); UnitCounter = cell(Nshanks,1); UnitFR = struct('pertrial', cell(Nshanks,1), ... 'session', cell(Nshanks,1)); UnitSTD = struct('pertrial', cell(Nshanks,1), ... 'session', cell(Nshanks,1)); PlaceFrMap = struct('frmap', cell(Nshanks,1), ... 'location', cell(Nshanks,1),... 'orientation',cell(Nshanks,1)); HDFrMap = struct('frtheta', cell(Nshanks,1), ... 'frphi', cell(Nshanks,1)); %%% END: Initialize Variables: shared and pershank %%% %%% Access kinematic model variables %%%%%%%% nodeLabels = data_frame.nodeLabels; Nnodes = numel(data_frame.nodeLabels); noseID = strcmp(nodeLabels,'nose'); baseskullID = strcmp(nodeLabels,'baseskull'); %neckID = strcmp(nodeLabels,'neck'); %%% END: access kinematic model variables %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Main Processing Loop %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%% % check processing branch based upon file size if SpikeBytes < MaxSpikeBytes % access Neural Data cd(NeuralDir) % SpikesData is a cell array with an Nx2 entry [res clu] for each shank % By logging all the spikes in advance we can count spikes within a % frame across all shanks at once. SpikeData = struct('SpikeMat',cell(length(resfiles),1)); % ids for matrix indexing Ids = cell(length(resfiles),1); Nneurons = zeros(length(resfiles),1); % load spike data and populate unit IDs for ii = 1:Nshanks % log shank ID shankID = strfind(resfiles(ii).name,'.'); shankID = str2double(resfiles(ii).name(shankID(end)+1:end)); % log shank id pershank(ii).ShankID = shankID; % load res and clu file clu = load(clufiles(ii).name); res = load(resfiles(ii).name); % convert spike times to seconds res = res/SRspikes; % The number of units and a vector of ids Nneurons(ii) = clu(1); ids = 2:Nneurons(ii)-1; % Update output UnitIds: to find index: e.g. round(strfind(UnitIds,'02010')/6)+1 pershank(ii).UnitIds = ... sprintf('%02d%03d\n',[shankID*ones(1,length(ids));ids]); % get rid of the artifact and multiunit clu and generate SpikeMat Ids{ii} = ids - 1; Nneurons(ii) = length(Ids{ii}); SpikeMat = [res clu(2:end)-1]; SpikeMat((SpikeMat(:,2) == -1) | (SpikeMat(:,2) == 0),:) = []; SpikeData(ii).SpikeMat = SpikeMat; % get rid of spiking data outside the period in which we are using % the camera system tooEarly = SpikeData(ii).SpikeMat(:,1) < Trials{1}(1); tooLate = SpikeData(ii).SpikeMat(:,1) > Trials{end}(end); excludeTs = tooEarly | tooLate; SpikeData(ii).SpikeMat(excludeTs,:) = []; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Initialize parameters requiring NEURON COUNT S%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Unit information cell arrays: {shank}{unit}{trial} UnitInfo{ii} = cell(Nneurons(ii),1); for jj = 1:Nneurons(ii) UnitInfo{ii}{jj} = cell(Ntrials,1); end % Unit firing rates and stds will be estimated recursively % over the trials and session UnitFR(ii).pertrial = zeros(Nneurons(ii),Ntrials); UnitFR(ii).session = zeros(Nneurons(ii),1); UnitSTD(ii).pertrial = zeros(Nneurons(ii),Ntrials); UnitSTD(ii).session = zeros(Nneurons(ii),1); % Place field information PlaceFrMap(ii).frmap = zeros(Dim,Dim,Ntrials,Nneurons(ii)); PlaceFrMap(ii).location = zeros(Dim,Dim,Ntrials,Nneurons(ii)); PlaceFrMap(ii).orientation.x = zeros(Dim,Dim,Ntrials,Nneurons(ii)); PlaceFrMap(ii).orientation.y = zeros(Dim,Dim,Ntrials,Nneurons(ii)); % Head direction information HDFrMap(ii).frtheta = zeros(Ntheta,Ntrials,Nneurons(ii)); HDFrMap(ii).frphi = zeros(Nphi,Ntrials,Nneurons(ii)); end clear clu res SpikeMat %%%%%%%%%%%%%%%%%%%%%%%% %%% LOOP OVER TRIALS %%% %%%%%%%%%%%%%%%%%%%%%%%% % access behavioral data cd(BehaviorDir) % Load reward configuration, if available reward_ids = dir('Reward_configuration*'); if isempty(reward_ids) params.reward_info = 'No reward configuration available.'; else reward_ids = load(reward_ids.name); reward_ids = reward_ids.reward_ids; % Generate array for logging reward locations in 2 and 3D. % (linear histogram index, X,Y,Z)' for each reward point RewardLocs = zeros(numel(reward_ids),1+3,Ntrials); end TotalFrameCounter = 0; % for session firing rate and std for ii = 1:Ntrials fprintf(1,'Trial index %d\n',ii); % INITIALIZE within trial UnitInfo variables for jj = 1:Nshanks UnitInfoTemp{jj} = zeros(Nframes(ii),3,length(Ids{jj})); UnitCounter{jj} = zeros(Nneurons(jj),1); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% LOAD and pre-process environment data %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % find the rend_envfiles corresponding to the current trial for jj = 1:length(rend_envfiles) % access environment rendering for this trial env_data = rend_envfiles(jj).name; % access trial name and number chk = strfind(env_data,'_'); trialname = env_data(1:chk(1)-1); trialnumber = str2double(trialname(regexp(trialname,'[0-9]'))); if trialnumber == TrialsList(ii) break end end fprintf(1,'\tTrial number %d\n',trialnumber); % log trial id shared(ii).TrialID = trialnumber; % Load the environment rendering for this trial load(env_data) env_pts = rend_env.env_pts; neighbors = rend_env.neighbors; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% create 3D environment to 2D histogram mapping %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This is always relative to the current position of the % cheeseboard maze with respect to the cameras planarMap = cheeseboardMaze2planarMap(rend_env); planarPts = planarMap.coords; coordFrames = planarMap.coordFrames; % Log reward locations in 2 and 3D if ~isempty(reward_ids) RewardLocs(:,2:4,ii) = env_pts(:,reward_ids)'; % project point into plane planarPos = planarPts(reward_ids,:); % find linear indices corresponding to these planar locations RewardLocs(:,1,ii) = table(dsearchn(table(:,1:2),planarPos),3); end % The right-handed coordinate frame is set to follow the % conventions in label_circ_grid_rewards.m % positive x-axis: 12 o'clock % positive y-axis: 9 o'clock CTR1 = env_pts(:,89); Xaxis = env_pts(:,neighbors{89,1}(4))-CTR1; Xaxis = Xaxis/norm(Xaxis); Yaxis = env_pts(:,neighbors{89,1}(3))-CTR1; Yaxis = Yaxis/norm(Yaxis); % planar correction to make as close to orthogonal as possible Yaxis = Yaxis - (Yaxis'*Xaxis)*Xaxis; Yaxis = Yaxis/norm(Yaxis); Zaxis = cross(Xaxis,Yaxis); Zaxis = Zaxis/norm(Zaxis); % log coordinate frame for this trial CoordFrame(:,:,ii) = [Xaxis Yaxis Zaxis]; if ii == 1 % if this is the first, reference trial, the transformation to % the reference frame is trivial CTR0 = CTR1; Transl_g(:,:,ii) = eye(4); else % Otherwise, we need to construct the transformation that will % take us from the current coordinate frame back to the % reference frame [~,Rl_g,Tl_g] = estsimt(... [CoordFrame(:,:,ii) -CoordFrame(:,1:2,ii)] + repmat(CTR1,[1 5]),... [CoordFrame(:,:,1) -CoordFrame(:,1:2,1)] + repmat(CTR0,[1 5])); Transl_g(:,:,ii) = [Rl_g Tl_g;zeros(1,3) 1]; end % Update the number of frames per trial according to the km folder % size. There can be a few more triggers than model fits. Nkms = size(dir(sprintf('%s%s%s_km%skm*',trialname,filesep,trialname,filesep)),1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% LOOP OVER CAMERA FRAMES %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The subject's speed for frame j-1 is estimated as the position % difference from frame j-2 to frame j-1 over the interval of the % camera exposures. Thus, its speed in frame j-1 is our estimate of % its speed at the beginning of the frame. Though, the spikes % associated with frame j-1 are counted from j-1 -> j. This is how % I captured the notion of "within frame" spike counts. first = 2; % for handling possible missing km_modelfits TrialFrameCounter = 0; % for within trial firing rate and std fprintf(1,'\t\tProcessing frames for this trial...\n'); for jj = 2:Nkms %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Generate kinematic model pose %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nodes = reshape(kmData(ii).Data(jj-1,:),[3 Nnodes]); offset = kmData(ii).offsets(jj-1,:); nodes = nodes + repmat(offset',[1 Nnodes]); % assuming kinematic model conventions for node placement nose = nodes(:,noseID); baseofskull = nodes(:,baseskullID); % continue if missing data for this frame if isnan(nodes(1)) % increment first frame is necessary if first first = first + 1; end continue else if jj == first ctrm0 = offset; trajectory{ii}(jj-1,:) = ctrm0; else ctrm1 = offset; trajectory{ii}(jj-1,:) = ctrm1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% LOG SPEED FOR THIS FRAME %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% speed{ii}(jj-1) = norm(ctrm1-ctrm0)*distfactor*fps; ctrm0 = ctrm1; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% GENERATE PLACE FIELD OCCUPANCY FOR THIS POSE %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % the place field will be determined relative to the nose of % the model and mapped into the 2D histogram via nearest % neighbor search %project into the plane ind = dsearchn(env_pts',nose'); chk0 = nose-env_pts(:,ind); planarPos = planarPts(ind,:) + ... [coordFrames(:,1,ind)'*chk0, ... coordFrames(:,2,ind)'*chk0]; placeInd = table(dsearchn(table(:,1:2),planarPos),3); Placemask = zeros(Dim,Dim); Placemask(placeInd) = 1; % log location time-series data location{ii}(jj-1) = placeInd; %%% PLACE estimate: end %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% GENERATE 3D HEAD DIRECTION FOR THIS POSE %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % planar head dir is set by vector from base of skull to nose hd_vectorH = nose-baseofskull; hd_vectorH = hd_vectorH/norm(hd_vectorH); % vertical head dir is set by vector from base of skull to nose hd_vectorV = hd_vectorH; % Calculate theta and phi xH = hd_vectorH'*Xaxis; yH = hd_vectorH'*Yaxis; xV = hd_vectorV'*Xaxis; yV = hd_vectorV'*Yaxis; zd = hd_vectorV'*Zaxis; h = sqrt(xV^2+yV^2); dtheta = atan2(yH,xH); dphi = atan(zd/h); % scale xH and yH such that the resulting vector is unit % length in the local xy-plane (for rate estimation) hd_scale = (1/norm(Xaxis*xH+Yaxis*yH)); xH = hd_scale*xH; yH = hd_scale*yH; % log head direction time-series data HDts{ii}(jj-1,1) = dtheta; HDts{ii}(jj-1,2) = dphi; % now bin the 3D head directions dtheta = histc(dtheta,params.theta,1); dphi = histc(dphi,params.phi,1); % log head direction occupancy HDtime(ii).Theta = HDtime(ii).Theta + dtheta; HDtime(ii).Phi = HDtime(ii).Phi + dphi; %%% Head Direction estimate: end %%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% LOOP OVER ALL NEURAL DATA SETS %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for kk = 1:Nshanks ind = find(SpikeData(kk).SpikeMat(:,1)>Trials{ii}(jj-1) & SpikeData(kk).SpikeMat(:,1)0); Nids = length(ids); % Increment UnitCounter UnitCounter{kk}(ids) = UnitCounter{kk}(ids) + 1; % Update UnitInfo (intense indexing to avoid forloop) UnitInfoTemp{kk}(... repmat(Nframes(ii)*3*(ids-1),[1 3]) + ... repmat(UnitCounter{kk}(ids),[1 3]) + ... repmat([0 Nframes(ii) 2*Nframes(ii)],[Nids 1])) = ... [spikes(ids) repmat(Trials{ii}(jj-1),[Nids 1]) repmat(placeInd,[Nids 1])]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% UPDATE FIRING RATE AND VARIANCE ESTIMATES %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Update means recursively UnitFR(kk).pertrial(:,ii) = ... 1/(TrialFrameCounter+1)*... (TrialFrameCounter*UnitFR(kk).pertrial(:,ii) + spikes); UnitFR(kk).session(:) = ... 1/(TotalFrameCounter+1)* ... (TotalFrameCounter*UnitFR(kk).session(:) + spikes); % Update variances recursively if TrialFrameCounter > 1 UnitSTD(kk).pertrial(:,ii) = ... (TrialFrameCounter-1)/TrialFrameCounter* ... UnitSTD(kk).pertrial(:,ii) + ... (1/(TrialFrameCounter-1)) * ... (spikes-UnitFR(kk).pertrial(:,ii)).^2; end if TotalFrameCounter > 1 UnitSTD(kk).session(:) = ... (TotalFrameCounter-1)/TotalFrameCounter* ... UnitSTD(kk).session(:) + ... (1/(TotalFrameCounter-1)) * ... (spikes-UnitFR(kk).session(:)).^2; end % advance to next frame if no spikes to log if isempty(ind) continue end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% UPDATE PLACE FIELD INFORMATION %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % UPDATE: place cell information % increment place occupancy Placetime(:,:,ii) = Placetime(:,:,ii) + Placemask; tilePlacemask = repmat(Placemask,[1 1 1 Nids]); PlaceFrMap(kk).frmap(:,:,ii,ids) = ... PlaceFrMap(kk).frmap(:,:,ii,ids) + ... reshape(repmat(spikes(ids)',[Dim^2 1]),[Dim Dim 1 Nids]).*tilePlacemask; % Update place location array PlaceFrMap(kk).location(:,:,ii,ids) = ... PlaceFrMap(kk).location(:,:,ii,ids) + tilePlacemask; % Update place orientation array (weighted by spike count) PlaceFrMap(kk).orientation.x(:,:,ii,ids) = ... PlaceFrMap(kk).orientation.x(:,:,ii,ids) + ... reshape(repmat((xH*spikes(ids))',[Dim^2 1]),[Dim Dim 1 Nids]).*tilePlacemask; PlaceFrMap(kk).orientation.y(:,:,ii,ids) = ... PlaceFrMap(kk).orientation.y(:,:,ii,ids) + ... reshape(repmat((yH*spikes(ids))',[Dim^2 1]),[Dim Dim 1 Nids]).*tilePlacemask; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% UPDATE HEAD DIRECTION INFORMATION %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HDFrMap(kk).frtheta(:,ii,ids) = ... HDFrMap(kk).frtheta(:,ii,ids) + ... reshape(repmat((spikes(ids))',[Ntheta 1]),[Ntheta 1 Nids]).*... repmat(dtheta,[1 1 Nids]); HDFrMap(kk).frphi(:,ii,ids) = ... HDFrMap(kk).frphi(:,ii,ids) + ... reshape(repmat((spikes(ids))',[Nphi 1]),[Nphi 1 Nids]).*... repmat(dphi,[1 1 Nids]); end % For firing rates and stds, I'm only counting those % frames in which I have behavioral data. This results in a % minimal distortion and keeps the statistics consistent. % increment trial frame counter TrialFrameCounter = TrialFrameCounter + 1; % increment session frame counter TotalFrameCounter = TotalFrameCounter + 1; end % Deal out UnitInfo on a per neuron basis for jj = 1:Nshanks for kk = 1:Nneurons(jj) if UnitCounter{jj}(kk) > 0 UnitInfo{jj}{kk}{ii} = ... squeeze(UnitInfoTemp{jj}(1:UnitCounter{jj}(kk),:,kk)); end end end end else fprintf(1,'There is a maximum file size of %6.3 Gb accepted by this code.\n',MaxSpikeByes/1024^3); fprintf(1,'\tChange the internal variable "MaxSpikeByes" if this is too low'); return end %%%%%%%%%%%%%%%%%%%%%%% %%% GENERATE OUTPUT %%% %%%%%%%%%%%%%%%%%%%%%%% % populate shared information for ii = 1:Ntrials shared(ii).Placetime = Placetime(:,:,ii); shared(ii).planarGrid = planarGrid{ii}; shared(ii).trajectory = trajectory{ii}; shared(ii).HDtime = HDtime(ii); shared(ii).HDts = HDts{ii}; shared(ii).location = location{ii}; shared(ii).speed = speed{ii}; shared(ii).timeline = Trials{ii}; shared(ii).CoordFrame = CoordFrame(:,:,ii); shared(ii).Transl_g = Transl_g(:,:,ii); end % populated per shank information for ii = 1:Nshanks % Convert firing rates and standard deviations from per frame % to per second UnitFR(ii).pertrial = UnitFR(ii).pertrial*fps; UnitFR(ii).session = UnitFR(ii).session*fps; UnitSTD(ii).pertrial = sqrt(UnitSTD(ii).pertrial)*fps; UnitSTD(ii).session = sqrt(UnitSTD(ii).session)*fps; pershank(ii).UnitInfo = UnitInfo{ii}; pershank(ii).UnitFR = UnitFR(ii); pershank(ii).UnitSTD = UnitSTD(ii); pershank(ii).PlaceFrMap = PlaceFrMap(ii); pershank(ii).HDFrMap = HDFrMap(ii); end % update params params.distfactor = distfactor; if ~isempty(reward_ids) params.reward_info = RewardLocs; end params.Ntrials = Ntrials; % put it all together SpatialCorrelates.shared = shared; SpatialCorrelates.pershank = pershank; SpatialCorrelates.params = params; % save output to NeuralDir save(sprintf('%s%s%s_spcorrs.mat',NeuralDir,filesep,filebase),'SpatialCorrelates') %return to caller directory cd(CurrentDir) % write out log file logAnalysisFile(mfilename('fullpath')); end