function [ScriptVar_new, FRAPDataVG, FRAPDataNonVG, FRAPDataVG_bkg, ...
    FRAPDataNonVG_bkg] = InSilicoFRAP_sepC(ScriptVar, AllCoords, ...
    xyFRAPRegion, xyFRAPBkgRegion, FRAP_bleaching_frame)
    % Outputs -
    %    ScriptVar_new - modified ScriptVar with two added keys: xyFRAPRegion,
    %                    FRAP_bleaching_frame 
    %    FRAPDataVG - 1xk array of number of unbleached particles inside the
    %               viscous granule portions of the ROI; k = Total number of frames
    %    FRAPDataNonVG - 1xk array of number of unbleached particles outside
    %               viscous granules but inside the ROI; k = Total number of frames
    %    FRAPDataVG_bkg - 1xk array of number of unbleached particles inside the
    %               viscous granule portions of the background region; k = Total number of frames
    %    FRAPDataNonVG_bkg - 1xk array of number of unbleached particles outside 
    %               viscous granule but inside the background region; k = Total number of frames
    % Inputs - 
    %    ScriptVar - Parameters container map
    %    AllCoords - Particle coordinates
    %    xyFRAPRegion - 1x1 cell whose 1x3 array
    %                   has the coordinates of the center of the ROI,
    %                   radius
    %    FRAP_bleaching_frame - the frame in which in silico photobleaching occurs
    % 
    % Example - 
    % [ScriptVar_new, FRAPDataVG, FRAPDataNonVG, FRAPDataVG_bkg, ...
    % FRAPDataNonVG_bkg] = InSilicoFRAP_sepC(ScriptVar, AllCoords, ...
    % {[0.3 0.45/2 0.2]}, {[0.7 0.45/2 0.2]}, 4000)
    
    tic
    fprintf('Checkpoint -5\n')
    toc
    fprintf('\n')
    AllCoords = single(AllCoords); % for reducing memory usage
    N_Particles = ScriptVar('N_Particles');
    Total_Particles = sum(N_Particles);
    
    %% Assign photobleaching status for all particles, ...
    %% in the photobleaching frame and assign whether inside or outside granules
    fprintf('Checkpoint -4\n')
    toc
    fprintf('\n')
    FRAP_indices = IsParticleInsideFigure(AllCoords(:,1,FRAP_bleaching_frame),...
        AllCoords(:,2,FRAP_bleaching_frame), xyFRAPRegion)';
    % FRAP_indices is a kx1 array where k = total number of particles
    % It is an array of logical values where 1 => photobleached during and
    % beyond FRAP_bleaching_frame, 0 => otherwise
    fprintf('Checkpoint -3\n')
    toc
    fprintf('\n')
    
    %% Assign particle indices
    Particle_index = zeros(size(FRAP_indices));
    for i = 1:Total_Particles
       Particle_index(i,1) = ParticleTypeCalculator(ScriptVar,i);
    end
    fprintf('Checkpoint -2\n')
    toc
    fprintf('\n')
    
    %% Carry over the photobleaching status values to subsequent frames
    %% Carry over VG index (inside or outside glyco granules) to all frames
    % For frames before the photobleaching frame, the status value = 0 for
    % all particles
    N_Frames = size(AllCoords,3);
    AllCoords = cat(2, AllCoords, zeros(Total_Particles, 2, N_Frames));
    for i = 1:N_Frames
        AllCoords(:,5,i) = AllCoords(:,3,i); % Make fifth column the 
        % viscous granule indices (1 if particle inside granule, else 0)
        AllCoords(:,4,i) = Particle_index; % Fourth column = which particle type
        % For default case, fourth column = 1 throughout
        if i >= FRAP_bleaching_frame
            AllCoords(:,3,i) = FRAP_indices; % FRAP indices assigned for all frames
            % >= in silico photobleaching frame
        else
            AllCoords(:,3,i) = 0; % unbleached as in silico photobleaching not yet performed
        end
    end
    clear FRAP_indices Particle_index
    fprintf('Checkpoint -1\n')
    toc
    fprintf('\n')
    
    %% Count the number of unbleached particles within the ROI
    n_inside_VG = zeros(length(N_Particles),N_Frames);
    n_inside_NonVG = zeros(length(N_Particles),N_Frames);
    n_inside_bkg_VG = zeros(size(n_inside_VG));
    n_inside_bkg_NonVG = zeros(size(n_inside_VG));
    xyFRAPRegion = xyFRAPRegion{1};
    xyFRAPBkgRegion = xyFRAPBkgRegion{1};
    fprintf('Checkpoint 0\n')
    toc
    fprintf('\n')
    
    %% ROI
    ind = AllCoords(:,3,:)==1;
    ind = cat(2,ind,ind,ind,ind,ind);
    AllCoords(ind) = NaN;
    fprintf('Checkpoint 1\n')
    toc
    fprintf('\n')
    
    if length(xyFRAPRegion) == 3 % circular ROI
        ind = (AllCoords(:,1,:) - xyFRAPRegion(1)).^2 + ...
            (AllCoords(:,2,:) - xyFRAPRegion(2)).^2;
        ind = ind <= xyFRAPRegion(3)^2;
    elseif length(xyFRAPRegion) == 4 % rect ROI
        ind = AllCoords(:,1,:) <= xyFRAPRegion(1)+xyFRAPRegion(3) ...
            & AllCoords(:,1,:) >= xyFRAPRegion(1) ...
            & AllCoords(:,2,:) <= xyFRAPRegion(2)+xyFRAPRegion(4) ...
            & AllCoords(:,2,:) >= xyFRAPRegion(2);
    end
    
    ind = cat(2,ind,ind,ind,ind,ind);
    x_ROI = AllCoords;
    x_ROI(~ind) = NaN;
    fprintf('Checkpoint 2\n')
    toc
    fprintf('\n')
    
    %% Bkg
    if length(xyFRAPBkgRegion) == 3 % circular Background ROI
        ind = (AllCoords(:,1,:) - xyFRAPBkgRegion(1)).^2 + ...
            (AllCoords(:,2,:) - xyFRAPBkgRegion(2)).^2;
        ind = ind <= xyFRAPBkgRegion(3)^2;
    elseif length(xyFRAPBkgRegion) == 4 % rect Background ROI
        ind = AllCoords(:,1,:) <= xyFRAPBkgRegion(1)+xyFRAPBkgRegion(3) ...
        & AllCoords(:,1,:) >= xyFRAPBkgRegion(1) ...
        & AllCoords(:,2,:) <= xyFRAPBkgRegion(2)+xyFRAPBkgRegion(4) ...
        & AllCoords(:,2,:) >= xyFRAPBkgRegion(2);
    end
    
    ind = cat(2,ind,ind,ind,ind,ind);
    x_bkg = AllCoords;
    x_bkg(~ind) = NaN;
    fprintf('Checkpoint 3\n')
    toc
    fprintf('\n')
    clear AllCoords
    
    %% ROI VG and NonVG   
    x_ROI_VG = x_ROI;
    ind = x_ROI(:,5,:)==0;
    ind = cat(2,ind,ind,ind,ind,ind);
    x_ROI_VG(ind) = NaN;
    fprintf('Checkpoint 4\n')
    toc
    fprintf('\n')
    
    x_ROI_NonVG = x_ROI;
    ind = x_ROI(:,5,:)==1;
    ind = cat(2,ind,ind,ind,ind,ind);
    x_ROI_NonVG(ind) = NaN;
    fprintf('Checkpoint 5\n')
    toc
    fprintf('\n')
    clear x_ROI
    
    %% Bkg VG and NonVG
    x_bkg_VG = x_bkg;
    ind = x_bkg(:,5,:)==0;
    ind = cat(2,ind,ind,ind,ind,ind);
    x_bkg_VG(ind) = NaN;
    fprintf('Checkpoint 6\n')
    toc
    fprintf('\n')
    
    x_bkg_NonVG = x_bkg;
    ind = x_bkg(:,5,:)==1;
    ind = cat(2,ind,ind,ind,ind,ind);
    x_bkg_NonVG(ind) = NaN;
    fprintf('Checkpoint 7\n')
    toc
    fprintf('\n')
    clear x_bkg
    
    %% Make final data
    for i = 1:length(N_Particles)
        x_ROI_VG = sum(~isnan(x_ROI_VG),1);
        x_ROI_VG = squeeze(x_ROI_VG(1,1,:));
        n_inside_VG(i,:) = x_ROI_VG.';
        
        x_ROI_NonVG = sum(~isnan(x_ROI_NonVG),1);
        x_ROI_NonVG = squeeze(x_ROI_NonVG(1,1,:));
        n_inside_NonVG(i,:) = x_ROI_NonVG.';
        
        x_bkg_VG = sum(~isnan(x_bkg_VG),1);
        x_bkg_VG = squeeze(x_bkg_VG(1,1,:));
        n_inside_bkg_VG(i,:) = x_bkg_VG.';

        x_bkg_NonVG = sum(~isnan(x_bkg_NonVG),1);
        x_bkg_NonVG = squeeze(x_bkg_NonVG(1,1,:));
        n_inside_bkg_NonVG(i,:) = x_bkg_NonVG.';
    end
    fprintf('Checkpoint 8\n')
    toc
    fprintf('\n')
    %% Make output vars
    
    FRAPDataVG = n_inside_VG;
    FRAPDataNonVG = n_inside_NonVG;
    
    FRAPDataVG_bkg = n_inside_bkg_VG;
    FRAPDataNonVG_bkg = n_inside_bkg_NonVG;
    
    ScriptVar_new = containers.Map(ScriptVar.keys, ScriptVar.values);
    ScriptVar_new('xyFRAPRegion') = {xyFRAPRegion};
    ScriptVar_new('xyFRAPBkgRegion') = {xyFRAPBkgRegion};
    ScriptVar_new('FRAP_bleaching_frame') = FRAP_bleaching_frame;
    toc
    % Note that the division of the number of unbleached particles in the
    % ROI and the background regions is not done here, as for
    % normalization, the background region is chosen from a uniform
    % viscosity simulation for standardization purposes.
end