function FitnessModelwithNoiseInnerForSharing(MEAN_vect,VAR_vect, VAR_type, FitFunct,FIT_var1,FIT_var2,iter,EXP_draws,REPORT)

% MODEL FOR EXPLORING THE EFFECT OF EXPRESSION NOISE ON FITNESS. Based on internal version 26. 

%INPUTS

%   MEAN_vect: numeric, a mean value to explore
%           (desinged to be submitted as a batch script)
%   VAR_vect: numeric, a noise value from a vector of noise values to explore 
%           (supplied to cluster)
%   VAR_type: string variable dictating exploration of noise as: "SD","CV", or 
%           "FANO"
%   FitFunct: string variable dictating fitness function form as linear
%           "Linear" or peaked "Normal"
%   FIT_var1: specifies fitness function, 
%                ->for peaked fitness function: 1
%                   specifies the expression level of fastest division rate
%                   (mean)
%                -> for linear function: 40
%                    specifies negative slope
%   FIT_var2: species fitness function,
%                ->for peaked fitness function: 0.3
%                   specifies the standard dev of a normal dist around the 
%                   expression level of fastest division rate
%                -> for linear function: 160
%                    specifies intercept
%   iter: numer of iterations
%   EXP_draws: variable for exploring consequences of more frequent expression
%           fluctuations within a cell cycle, integer describing number of 
%           unique samples from expression distribution used to calculate
%           doubling time, choose 1 unless explicitly testing
%   REPRT: string variable "concise" or "verbose" determines whether
%           instantaneous samples of expression and doubling time for the 
%            population are saved at intervals within run

%OUTPUT is # cells generated by innermost loop--all 
%iterations for one mean and variance combination, and some summaries of
%the doubling time distributions generated by the sampling. 

%

%In this case, dividing cells produce 2 daugheters drawn from the same 
%underlying distribution (no mother-daughter %distinction or %inheritance.)

%Select a seed population, growth interval & pop size to explore
initialize = 1000;
interval = 250;
duration = 1000;
time = interval:interval:duration;
timepoints = duration/interval;
saturation = 5*10^6; %only pre-allocates the starting size of the matrix, does not constrain it
N_T0=10^3;
edges=0:.15:600;

%Version naming distinguishes size of seed population, starting cells, 
%duration, final population matrix size (not limiting), and expression 
%levels sampled per cell cycle,  
vers=strcat('v26_u',num2str(MEAN_vect),'_s',num2str(VAR_vect),'_1Kspop_1Kdur_1KN_1Kfpop_',num2str(EXP_draws),'draws');

%Pre-allocating matrix for counts of doubling times sampled in each
%iteration
DT_out=zeros(length(MEAN_vect)*length(VAR_vect)*iter,length(edges)-1+2);
  

%NOTE: not currently using or reporting this output either, could delete
%from script shared to streamline?

%Pre-allocating summary variables to describe doubling time distributions
%at end points (simplified by calculating over only innermost loop)
DT_mean=zeros(length(MEAN_vect)*length(VAR_vect),iter+2);
DT_std=zeros(length(MEAN_vect)*length(VAR_vect),iter+2);
DT_skew=zeros(length(MEAN_vect)*length(VAR_vect),iter+2);
DT_kurt=zeros(length(MEAN_vect)*length(VAR_vect),iter+2);

%Define a distribution relating expression to doubling time,
if strcmp(FitFunct,'Normal')
    f=@(x) normpdf(x,FIT_var1,FIT_var2);
    c=160/max(f(0.01:0.01:2)); %choosing 80 as arbitrary lowest value
    f=@(x) 240-c*normpdf(x,FIT_var1,FIT_var2); %choosing 160 as arbitrary highest value, can change
elseif strcmp(FitFunct,'Linear')
    %Nice starting points for FIT var 1 & 2:
    %[40,160],[50,180],[30,140],[80,240]
    f=@(x) FIT_var2-FIT_var1*x;
end


%% Perform the experimental sampling
inner_iter=0;

EXP_mu = MEAN_vect;
    
if strcmp(VAR_type,'SD')
    EXP_var = VAR_vect;
elseif strcmp(VAR_type,'CV')
    EXP_var = VAR_vect*EXP_mu;
elseif strcmp(VAR_type,'FANO')
    EXP_var= sqrt(VAR_vect*EXP_mu);
end

if strcmp(REPORT,'verbose') 
    EXP_SLICE=zeros(saturation,2*iter);
    DT_SLICE=zeros(saturation,2*iter);
end

cells = zeros(iter,length(time));

tic
for G = 1:iter
    inner_iter=inner_iter+1;

    DT = zeros(saturation,1);
    EXP = zeros(saturation,EXP_draws);
    TIME = ones(saturation,1)*interval;
    TIME(1:N_T0) = ones(N_T0,1)*initialize;

    %Sample initial population and randomize their ages
    EXP(1:N_T0,1:EXP_draws)=normrnd(EXP_mu,EXP_var,N_T0,EXP_draws);    
    GR_sample = f(EXP(1:N_T0,1:EXP_draws));
    DT(1:N_T0) = mean(GR_sample,2).*rand(N_T0,1);

    loc1 = N_T0; %loc1 is end of currently dividing pop 
    loc2 = loc1+1; %loc2 is index at which to add new cells
    j=1;

    %Initialize with a run in which data is not recorded so that
    %cells are not synchronized at startpoint
        while j < loc1
                while TIME(j) >= DT(j)
                  EXP(j,:)=normrnd(EXP_mu,EXP_var,1,EXP_draws);
                  TIME(j) = TIME(j) - DT(j);
                  TIME(loc2) = TIME(j) - DT(j);
                  GR_sample = f(EXP(j,:));
                  DT(j) =  mean(GR_sample);
                  GR_sample = f(EXP(loc2,:));
                  DT(loc2)= mean(GR_sample);
                  
                end
                j=j+1;
                loc2=loc2+1;
        end
        begin_idx=datasample(1:1:loc2-1,N_T0);

        if strcmp(REPORT,'verbose') 
            EXP_SLICE(1:N_T0,2*G-1)=EXP(begin_idx);
            DT_SLICE(1:N_T0,2*G-1)=DT(begin_idx);
        end


        TIME(1:N_T0)=TIME(begin_idx)+interval;
        TIME(N_T0+1:saturation) = ones(saturation-N_T0,1)*interval;
            
            %Now collect data
            
            loc1 = N_T0; %loc1 is end of currently dividing pop 
            loc2 = loc1+1; %loc2 is index at which to add new cells
            N = histcounts(DT(1:N_T0),edges);
            DT_out(inner_iter,3:end)=DT_out(inner_iter,3:end)+N;

             for i = 1:timepoints;
                 j=1; %start with first cell in matrix
                    while j < loc2  %while you have not yet run through all dividing cells 
                        
                        while TIME(j) >= DT(j) %while the time left in the timepoint is greater than the length of a cell division
                            EXP(j,:)=normrnd(EXP_mu,EXP_var,1,EXP_draws);
                            TIME(j) = TIME(j) - DT(j);
                            TIME(loc2) = TIME(j) - DT(j);
                            GR_sample = f(EXP(j,:));
                            DT(j) =  mean(GR_sample);
                            GR_sample = f(EXP(loc2,:));
                            DT(loc2)= mean(GR_sample);
                            
                            N = histcounts([DT(loc2);DT(j)],edges);
                            DT_out(inner_iter,3:end)=DT_out(inner_iter,3:end)+N;
                            loc2 = loc2 + 1;

                        end
                        TIME(j) = TIME(j)+ interval;
                        j=j+1;
                   end
                
                cells(G,i) = loc2-1;
                %disp(strcat('Cells = ',num2str(cells(G,i))))
                
             end
         if strcmp(REPORT,'verbose') 
            EXP_SLICE(:,2*G)=EXP(:,1);
            DT_SLICE(:,2*G)=DT;
         end
         disp(strcat('G = ', num2str(G),' iterations'));     


        %Extract features of cumulative DT dist
        idx1=find(DT_out(inner_iter,3:end)); %annoying that this version of MATLAB does implement repelem
        n = cumsum(DT_out(inner_iter,idx1+2));
        idx2 = zeros(n(end),1);
        idx2(n-DT_out(inner_iter,idx1+2)+1) = 1;
        bins=edges(idx1);
        DT_dist= bins(cumsum(idx2));

        DT_mean(inner_iter,G+2)=mean(DT_dist);
        DT_std(inner_iter,G+2)= std(DT_dist);
        DT_skew(inner_iter,G+2)=skewness(DT_dist);
        DT_kurt(inner_iter,G+2)=kurtosis(DT_dist);

end
    
stop_time = toc;
   
    if strcmp(REPORT,'verbose')  
        [r,~]=find(EXP_SLICE,1,'last');
        filename = strcat('EXPslice',vers,'_Exp_m',num2str(EXP_mu),'s',num2str(VAR_vect(VAR_iter)),'_',VAR_type,'_',FitFunct,'_m',num2str(FIT_var1),'s',num2str(FIT_var2),'.txt');
        dlmwrite(filename,EXP_SLICE(1:r,:),'delimiter','\t')
        filename = strcat('DTslice',vers,'_Exp_m',num2str(EXP_mu),'s',num2str(VAR_vect(VAR_iter)),'_',VAR_type,'_',FitFunct,'_m',num2str(FIT_var1),'s',num2str(FIT_var2),'.txt');
        dlmwrite(filename,DT_SLICE(1:r,:),'delimiter','\t')
   end
   
DT_std(1:2)= [EXP_mu,EXP_var];
DT_mean(1:2)=[EXP_mu,EXP_var];
DT_skew(1:2)= [EXP_mu,EXP_var];
DT_kurt(1:2)= [EXP_mu,EXP_var];

DT_out(inner_iter-iter+1:inner_iter,1:2)=repmat([EXP_mu,EXP_var],iter,1);
GROWTH_cells=[EXP_mu,EXP_var,reshape(transpose(cells),1,[])];
disp(strcat('Processing time: ', num2str(stop_time/60),' minutes'));

%Critical output: number of cells generated at endpoint
filename = strcat('Cells',vers,'_',FitFunct,'_m',num2str(FIT_var1),'s',num2str(FIT_var2),'_',VAR_type,'_',num2str(iter),'iter.txt');
dlmwrite(filename,[[0,0,repmat(time,1,iter)];GROWTH_cells],'delimiter','\t')


%Variety of summaries of doubling times sampled: hist counts and mean, st
%dev, skew, kurtosis
filename = strcat('DTcounts',vers,'_',FitFunct,'_m',num2str(FIT_var1),'s',num2str(FIT_var2),'_',VAR_type,'_',num2str(iter),'iter.txt');
dlmwrite(filename,[0,0,edges(1:length(edges)-1);DT_out],'delimiter','\t')
filename = strcat('DT_STD',vers,'_',FitFunct,'_m',num2str(FIT_var1),'s',num2str(FIT_var2),'_',VAR_type,'_',num2str(iter),'iter.txt');
dlmwrite(filename,DT_std,'delimiter','\t')
filename = strcat('DT_MEAN',vers,'_',FitFunct,'_m',num2str(FIT_var1),'s',num2str(FIT_var2),'_',VAR_type,'_',num2str(iter),'iter.txt');
dlmwrite(filename,DT_mean,'delimiter','\t')
filename = strcat('DT_SKEW',vers,'_',FitFunct,'_m',num2str(FIT_var1),'s',num2str(FIT_var2),'_',VAR_type,'_',num2str(iter),'iter.txt');
dlmwrite(filename,DT_skew,'delimiter','\t')
filename = strcat('DT_KURT',vers,'_',FitFunct,'_m',num2str(FIT_var1),'s',num2str(FIT_var2),'_',VAR_type,'_',num2str(iter),'iter.txt');
dlmwrite(filename,DT_kurt,'delimiter','\t')
end
