clearvars

frames = 600; %number of frames in the movie
StoN = 1; %signal to noise
taus = 3; % number of tau subunits in the clamp loader
replisomes = 300; % number of replisomes
mol_per_replisome = 3; % number of replilsomes per cell
fracbndt0 = 0.5; %initial fraction bound

randomise = 0; % off = 0; on = 1

Kon_repl = 0.01; % Kon for Pol III*
Koff_repl = 0.01; % Koff for Pol III*

Tbleach_core = 250; % Photobleaching time for the core
Kon_core = 1; % Kon for the core
Koff_core = 0.0000001; % Koff for the core

Tbleach_cl = 250; %Photobleading time for the CLC
Kon_cl = 1; %Kon for the CLC
Koff_cl = 0.0000001; %Koff for the CLC

Nmolec = replisomes*mol_per_replisome; %total number of molecules

core_signal = zeros(frames,Nmolec); %matrix for storing core signal
cl_signal = zeros(frames,Nmolec); %matrix for storing CLC signal

% Determine raw signal for all (also unbound) molecules in the cell
for m = 1:Nmolec
    core_signal(:,m) = Core_binding(taus,frames,1,Kon_core,Koff_core,Tbleach_core,0);
    cl_signal(:,m) = Clamploader_bleaching(taus,frames,1,Kon_cl,Koff_cl,Tbleach_cl,0);
end


% initialize input parameters
Nsites = replisomes;                % number of trajectories to obtain
Nframes = frames;              % length of the trajectories to obtain
Signal = 1;                 % 'signal' strength 
KonVal =     [1 0.1]*  Kon_repl;   % units are 'frames'
KoffVal =    [1 0.1]*  Koff_repl;

% determine the kon and koff rates for every frame at every site
Ton = norminv(rand(Nframes,Nsites),KonVal(1), KonVal(2));
Toff = norminv(rand(Nframes,Nsites),KoffVal(1),KoffVal(2));

% set random selection array
SelArray = rand(Nframes,Nsites);

% define an ID for each molecule having a bleaching traj
moleculeID = [zeros(Nmolec,1), rand(Nmolec,1)]; % initialize with all molecules unbleached
			% in-use state  , associated random number
% moleculeID = sortrows(moleculeID,2);


% \\ Set binding trajectories
bound = zeros(Nframes,Nsites);

% initialize the bound state
f = 1;
bound(f,:) = rand(1,Nsites);
bound(f,bound(f,:) > fracbndt0) = 1; 
bound(f,bound(f,:) <= fracbndt0) = 0;
idxbd = find(bound(f,:)==1);

if ~isempty(idxbd)
    [~,idx] =  min(abs(rand()-moleculeID(:,2)));
    bound(f,idxbd(1)) = idx;
    moleculeID(idx(1),1) = 1;
    
    for i = 2:size(idxbd,2)
        
        s = idxbd(i);
        molsel = 0;
        while molsel == 0
            
            [~,idx] = min(abs(rand()-moleculeID(:,2)));
            if moleculeID(idx,1) == 0
                
                molsel = 1;
                bound(f,s) = idx;
                moleculeID(idx,1) = 1;
                
            end
            
        end
        
    end
end

% Set binding state for all subsequent frames

for f = 2:Nframes
    
    idxbd = find(bound(f-1,:) ~= 0); % bound sites
    idxun = find(bound(f-1,:) == 0); % unbound sites
    
    % with bound sites, determine if they stay bound, or unbind in the
    % current frame
    for i = 1:size(idxbd,2)
        
        s = idxbd(i);
        
        Punbind = Toff(f,s)*exp(-Toff(f,s));
        %         Pstay = 1-Punbind;
        if SelArray(f,s) <= Punbind % molecule unbinds
            bound(f,s) = 0;
            moleculeID(bound(f-1,s),1) = 0;
        else % molecule stays bound
            bound(f,s) = bound(f-1,s);
        end
        
    end
    
    % with unbound sites, determine if they stay unbound, or bind a new
    % molecule
    for i = 1:size(idxun,2)
        
        s = idxun(i);
        Pbind = Ton(f,s)*exp(-Ton(f,s));
        if SelArray(f,s) <= Pbind
            
            molsel = 0;
            while molsel == 0
                
                [~,idx] = min(abs(rand()-moleculeID(:,2)));
                if moleculeID(idx,1) == 0
                    
                    molsel = 1;
                    bound(f,s) = idx;
                    moleculeID(idx,1) = 1;
                    
                end
                
            end
            
        end
        
    end
    
    
end
	
	
traj_bound = bound;
traj_bound(traj_bound>0)=Signal;

% Combine bound state and molecule bleaching to define signal trajectory

% initialize the signal trajectory as only noise around zero
traj_signal_core = norminv(rand(Nframes,Nsites),0,Signal/StoN);
traj_signal_cl = norminv(rand(Nframes,Nsites),0,Signal/StoN);

% if a site has a non-bleached molecule bound, increase the 
for f = 1:Nframes
    
    for s = 1:Nsites
        
        molsel = bound(f,s);
        if molsel~=0 
            
            traj_signal_core(f,s) = core_signal(f,s) + traj_signal_core(f,s);        
            traj_signal_cl(f,s) = cl_signal(f,s) + traj_signal_cl(f,s);
        end
        
    end
    
end

if randomise ==1
    traj_signal_core = traj_signal_core(:,randperm(size(traj_signal_core,2)));
else
end

x = 1:frames;
x2= -frames+1:frames-1;

% traj_signal_core=traj_signal_core(501:600,:);
% traj_signal_cl=traj_signal_cl(501:600,:);

x = 1:frames;
x2= -frames+1:frames-1;

% Calculate correlation functions
for r = 1:replisomes
    acf_core(:,r) = xcorr(traj_signal_core(:,r));
    acf_cl(:,r) = xcorr(traj_signal_cl(:,r));
    ccf(:,r) = xcorr(traj_signal_core(:,r),traj_signal_cl(:,r));
end

acf_core = acf_core(frames:end,:);
acf_cl = acf_cl(frames:end,:);

avg_acf_core = mean(acf_core,2);
avg_acf_cl = mean(acf_cl,2);
avg_ccf = mean(ccf,2);
flip_avg_ccf=flipud(avg_ccf);

% Find photobleaching trajectories and correlation functions
photobleach_cl = mean(traj_signal_core,2);
photobleach_core = mean(traj_signal_cl,2);

acf_photobleach_cl = xcorr(photobleach_cl);
acf_photobleach_core = xcorr(photobleach_core);

acf_photobleach_cl = acf_photobleach_cl(frames:end,:);
acf_photobleach_core = acf_photobleach_core(frames:end,:);

ccf_photobleach = xcorr(photobleach_cl,photobleach_core);
flip_ccf_photobleach = flipud(ccf_photobleach);

%Normalize the average correlation functions to the correlation functions
%of the photobleaching.
norm_mean_acf_cl = avg_acf_cl./acf_photobleach_cl;
norm_mean_acf_core = avg_acf_core./acf_photobleach_core;
norm_mean_ccf = avg_ccf./ccf_photobleach;
flip_norm_ccf = flipud(norm_mean_ccf);
norm_ccf = (norm_mean_ccf + flip_norm_ccf)./2;