12-Mar-2014 12:52:52 function [SpatialCorrelates] = spatialcorrelates2_2(model_name,NeuralDir, BehaviorDir, Trigger, exclusions, options) % [SpatialCorrelates] = spatialcorrelates(NeuralDir, BehaviorDir, Trigger, exlcusions, options) % % This routine calculates spatial correlates between neural and behavioral % data. These correlates are: 1) Place Cells, 2) Field of View Cells, and 3) % Head Direction Cells. 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). This specifications are a little fragile, and are % detailed within. This routine assumes the user have already included the % eye locations and norms in the kinematic model. % % Note: This code can we run anywhere, because the user must supply % complete paths to both the Neural and Behavioral datas. Currently, this % routine excludes the artifact and multi-unit/noise clusters. % % 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' % NeuralData: a string designating the absolute path to the neural data % directory. No trailing filesep. % BehaviorData: a string designating the absolute path to the behavioral % data directory. No trailing filesep. % Trigger: a Nx1 vector of 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 ... %%% 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 if isdir(BehaviorDir) % copy over the kinematic model, but first check if it exists in this % folder if isempty(dir(sprintf('%s%s%s.m',BehaviorDir,filesep,model_name))) % Store the directory in which this function was launched CurrentDir = pwd; cd(BehaviorDir) % copy in the model fprintf(1,'Copying in model.\n'); copyfile(sprintf('..%s%s.m',filesep,model_name),pwd); % load data_frame for nodeLabels load(sprintf('..%s%s.mat',filesep,model_name)) cd(CurrentDir) else % Store the directory in which this function was launched CurrentDir = pwd; cd(BehaviorDir) fprintf(1,'Model found in diretory...'); % Check whether the model in BehaviorDir is the same as the one in % the subject root directory M1 = dir(sprintf('%s.m',model_name)); M2 = dir(sprintf('..%s%s.m',filesep,model_name)); if ~strcmp(M1.date,M2.date) % delete existing file fprintf(1,'deleting.\n'); delete(sprintf('%s.m',model_name)); % copy in the new file fprintf(1,'Copying in model.\n'); copyfile(sprintf('..%s%s.m',filesep,model_name),pwd); else fprintf(1,'up to date.\n'); end % load data_frame for nodeLabels load(sprintf('..%s%s.mat',filesep,model_name)) cd(CurrentDir) end kmodel = eval(sprintf('@(w,x,y,z)%s(w,x,y,z)',model_name)); kmodel_frame = eval(sprintf('@(x,y,z)%s_frame(x,y,z)',model_name)); % 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 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.SpeedThresh = SpeedThresh; params.HDres = HDres; params.SpatialScale = SpatialScale; params.SRspikes = SRspikes; params.SReeg = SReeg; params.thresh = thresh; params.fps = fps; params.factor = factor; params.MaxSpikeBytes = MaxSpikeBytes; 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 tiemstamp of each frame in each trial, 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 % Exclude bad triggers i.e. accidental camera acquisition Trials(excludetriggers) = []; Nframes(excludetriggers) = []; TrialsList = 1:length(Trials); Trials(excludetrials) = []; Nframes(excludetrials) = []; TrialsList(excludetrials) = []; Ntrials = length(TrialsList); clear Trigger out % for interpolating FOV and position estimates range = linspace(0,1,SpatialScale+1); range = range(1:SpatialScale); %%%% 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: occupancy of the subject, conditional upon Place Cell firing % trajectory: trajectory in the local 3D coordinate frame % FOVtime: occupancy of the subject, conditional upon View Cell firing % FOVts: time-series of FOV, stored as pixel coordinates on boundary % HDtime: occupancy of the subject, conditional upon head direction % HDts: time-series of hd direction in 3D (theta,phi) % 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 % 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), ... 'trajectory', cell(Ntrials,1), ... 'FOVtime', cell(Ntrials,1), ... 'FOVts', cell(Ntrials,1), ... 'HDtime', cell(Ntrials,1), ... 'HDts', 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: Dim = 14*SpatialScale; % 14 is a parameter specific to the cheeseboard maze Ntheta = ceil(360/HDres); % quantization of 3D head direction Nphi = ceil(180/HDres); Placetime = zeros(Dim,Dim,Ntrials); % 2D environment histogram trajectory = cell(Ntrials,1); FOVtime = zeros(Dim,Dim,Ntrials); % 2D environment histogram FOVts = 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); % rotate about y-axis % initialize time-series variables for ii = 1:Ntrials FOVts{ii} = cell(Nframes(ii),1); % FOV verts 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: % 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 % FovFrMap: spike counts as function of postulated binocular view % 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), ... 'FovFrMap', cell(Nshanks,1), ... 'HDFrMap', cell(Nshanks,1)); % DECLARE: pershank sub-variables % For both Place Cell and View Cell estimation, I will also be keeping % track of the subject's spike count weighted, average head orientation % in the LOCAL XY-plane. Hence, the field 'orientation' is a structure % with fields x and y. % 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)); FovFrMap = 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; noseID = strcmp(nodeLabels,'nose'); baseskullID = strcmp(nodeLabels,'baseskull'); %%% 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 at a given % time 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 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 cells and FOV cells have their own location and orientation % arrays to along for disjoint exclusion criteria when allocating % spikes to different functional types 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)); FovFrMap(ii).frmap = zeros(Dim,Dim,Ntrials,Nneurons(ii)); FovFrMap(ii).location = zeros(Dim,Dim,Ntrials,Nneurons(ii)); FovFrMap(ii).orientation.x = zeros(Dim,Dim,Ntrials,Nneurons(ii)); FovFrMap(ii).orientation.y = zeros(Dim,Dim,Ntrials,Nneurons(ii)); 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 index, X,Y,Z)' for each reward point RewardLocs = zeros(1 + numel(reward_ids),3,Ntrials); end TotalFrameCounter = 0; % for session firing rate and std distfactor = 0; % for converting distances in 3D to cm for ii = 1:Ntrials disp('ii') disp(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 disp('trialnumber') disp(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 table = SpatialHistogram(trialname,SpatialScale); % Log reward locations in 2 and 3D if ~isempty(reward_ids) RewardLocs(2:end,:,ii) = env_pts(:,reward_ids); temp_subs = table(dsearchn(table(:,3:end),env_pts(:,reward_ids)'),1:2); RewardLocs(1,:,ii) = sub2ind([Dim Dim],temp_subs(:,1),temp_subs(:,2)); end % generate scale conversion factor to map from 3D frame to % centrimeters distances = zeros(648,1); conv_ind = 1; for jj = 1:size(env_pts,2) accum = sqrt(sqdist(env_pts(:,neighbors{jj,1}),env_pts(:,jj))); distances(conv_ind:conv_ind+length(accum)-1) = accum; conv_ind = conv_ind + length(accum); end % factor to convert to cm (8cm grid spacing) % Update recursively across trials distfactor = (1/ii)*((ii-1)*distfactor + 8/median(distances)); % The right-handed coordinate frame is set to follow the % conventions in label_circ_grid_rewards.m % positive x-axis: 3 o'clock % positive y-axis: 12 o'clock CTR1 = env_pts(:,89); Xaxis = env_pts(:,neighbors{89,1}(2))-CTR1; Xaxis = Xaxis/norm(Xaxis); Yaxis = env_pts(:,neighbors{89,1}(4))-CTR1; 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%s',trialname,filesep,trialname,filesep)),1)-2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 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 for jj = 2:Nkms % access kinematic model data try load(sprintf('%s%s%s_km%skm_%05d.mat',trialname,filesep,trialname,filesep,jj-1)) catch % end of trial in terms of model fits break end % continue if missing data for this frame if isempty(model_fit) % increment first frame is necessary if first first = first + 1; end continue else if jj == first ctrm0 = model_fit.offset; trajectory{ii}(jj-1,:) = ctrm0'; else ctrm1 = model_fit.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 kinematic model pose %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [~,~,~,eyes,eyeN] = kmodel([0 0 0],model_fit.sp_coords,model_fit.theta,[]); nodes = kmodel_frame([0 0 0],model_fit.sp_coords,model_fit.theta); % assuming kinematic model conventions for node placement nose = nodes(:,noseID); baseofskull = nodes(:,baseskullID); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 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 temp_subs = table(dsearchn(table(:,3:end),nose'),1:2); placeInd = sub2ind([Dim Dim],temp_subs(:,1),temp_subs(:,2)); Placemask = zeros(Dim,Dim); Placemask(placeInd) = 1; % log location time-series data location{ii}(jj-1) = placeInd; %%% PLACE estimate: end %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% GENERATE FIELD OF VIEW FOR THIS POSE %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The intersection of the binocular field of view with the % environment is defined by the union of the positive % half-spaces determined by the normals defined the monocular % fields of view FOVind = (sign(eyeN(:,1)'*(env_pts-repmat(eyes(:,1),[1 size(env_pts,2)]))) > 0) & ... (sign(eyeN(:,2)'*(env_pts-repmat(eyes(:,2),[1 size(env_pts,2)]))) > 0); % Interpolate along the boundaries of our FOV FOV = find(FOVind); % intialize array of boundary points bounds = zeros(1000,3); counter = 0; for kk = 1:length(FOV) % check if FOV(kk) is on the boundary if size(neighbors{FOV(kk),1}) < 4 counter = counter + 1; bounds(counter,:) = env_pts(:,FOV(kk))'; end % if not all neighbors are within the FOV, interpolate between and find % the last point within the FOV if sum(FOVind(neighbors{FOV(kk),1})) < length(neighbors{FOV(kk),1}) edge = FOV(kk); nearby = neighbors{FOV(kk),1}(~FOVind(neighbors{FOV(kk),1})); for ll = 1:length(nearby) % generate candidate points candidates = repmat(env_pts(:,edge),[1 SpatialScale]) + ... repmat(range,[3 1]).*... (repmat(env_pts(:,nearby(ll)),[1 SpatialScale]) - ... repmat(env_pts(:,edge),[1 SpatialScale])); % find all those that fall within the FOV indc = (sign(eyeN(:,1)'*(candidates-repmat(eyes(:,1),[1 SpatialScale]))) > 0) & ... (sign(eyeN(:,2)'*(candidates-repmat(eyes(:,2),[1 SpatialScale])))>0); % find that point that is within the FOV but furtherest from % edge inside = find(indc,1,'last'); counter = counter +1; bounds(counter,:) = candidates(:,inside)'; end end end % truncate bounds bounds = bounds(1:counter,:); bounds = unique(bounds,'rows','stable'); % Create image mask representing FOV pixels = table(dsearchn(table(:,3:end),bounds),1:2); try k = convhull(pixels(:,1),pixels(:,2)); pixels = pixels(k,:); % add 0.5 pixels to border with mean as center to follow % matlab standard and guarantee the correct FOV pixels = pixels + ... 0.5*sign(pixels-... repmat(mean(pixels),[size(pixels,1) 1])); % log FOV time-series information FOVts{ii}{jj-1} = pixels; % create mask within bounds and increment FOVtime FOVmask = poly2mask(pixels(:,2),pixels(:,1),Dim,Dim); catch % handle the case of just 1 point FOVinds = sub2ind([Dim Dim],pixels(:,1),pixels(:,2)); % log FOV time-series information FOVts{ii}{jj-1} = pixels; FOVmask = zeros(Dim,Dim); FOVmask(FOVinds) = 1; end %%% FOV estimate: END %%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% GENERATE 3D HEAD DIRECTION FOR THIS POSE %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% hd_vector = nose-baseofskull; hd_vector = hd_vector/norm(hd_vector); % Calculate theta and phi xd = hd_vector'*Xaxis; yd = hd_vector'*Yaxis; zd = hd_vector'*Zaxis; h = sqrt(xd^2+yd^2); dtheta = atan2(yd,xd); dphi = atan(zd/h); % scale xd and yd such that the resulting vector is of unit % length in the local xy-plane (for rate estimation) hd_scale = (1/norm(Xaxis*xd+Yaxis*yd)); xd = hd_scale*xd; yd = hd_scale*yd; % 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 AND FIELD OF VIEW INFORMATION %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% CONDITIONAL INCLUSION OF SPIKES TO PLACE/VIEW CELLS %%% % check speed condition for assigning spikes to place cells % or FOV cells if speed{ii}(jj-1) > SpeedThresh % UPDATE: field of view information % increment FOV occupancy FOVtime(:,:,ii) = FOVtime(:,:,ii) + FOVmask; tileFOVmask = repmat(FOVmask,[1 1 1 Nids]); tilePlacemask = repmat(Placemask,[1 1 1 Nids]); % increment FOV spike counts FovFrMap(kk).frmap(:,:,ii,ids) = ... FovFrMap(kk).frmap(:,:,ii,ids) + ... reshape(repmat(spikes(ids)',[Dim^2 1]),[Dim Dim 1 Nids]).*tileFOVmask; % Update FOV location array FovFrMap(kk).location(:,:,ii,ids) = ... FovFrMap(kk).location(:,:,ii,ids) + tilePlacemask; % Update FOV orientation array (weighted by count) FovFrMap(kk).orientation.x(:,:,ii,ids) = ... FovFrMap(kk).orientation.x(:,:,ii,ids) + ... reshape(repmat((xd*spikes(ids))',[Dim^2 1]),[Dim Dim 1 Nids]).*tilePlacemask; FovFrMap(kk).orientation.y(:,:,ii,ids) = ... FovFrMap(kk).orientation.y(:,:,ii,ids) + ... reshape(repmat((yd*spikes(ids))',[Dim^2 1]),[Dim Dim 1 Nids]).*tilePlacemask; % 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((xd*spikes(ids))',[Dim^2 1]),[Dim Dim 1 Nids]).*tilePlacemask; PlaceFrMap(kk).orientation.y(:,:,ii,ids) = ... PlaceFrMap(kk).orientation.y(:,:,ii,ids) + ... reshape(repmat((yd*spikes(ids))',[Dim^2 1]),[Dim Dim 1 Nids]).*tilePlacemask; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 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).trajectory = trajectory{ii}; shared(ii).FOVtime = FOVtime(:,:,ii); shared(ii).FOVts = FOVts{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).FovFrMap = FovFrMap(ii); pershank(ii).HDFrMap = HDFrMap(ii); end % update params params.distfactor = distfactor; if ~isempty(reward_ids) params.reward_info = RewardLocs; end % put it all together SpatialCorrelates.shared = shared; SpatialCorrelates.pershank = pershank; SpatialCorrelates.params = params; cd(CurrentDir) % write out log file logAnalysisFile(mfilename('fullpath')); end