fzosim_withstartdist_5categories.m

function [n_nucleoids_all,std_nucleoids_all] = fzosim_withstartdist_5categories(nsim,cycles,ncells,lognormal_mu,lognormal_sigma,synbase,syn4,syn7,syn10,syn13,initnorm)

n_nucleoids_all=[];
std_nucleoids_all=[];

for s=1:nsim %repeat simulation

if initnorm == 0
%initialize cells
numberofnucleoids=16;   %initial number of nucleoids at division; following Roussou et al., 2024 (32 for diploids)
nucleoids = numberofnucleoids*ones(ncells,1); %list of number of nucleoids at division for each cell

elseif initnorm == 1
%initialize with normal distribution
init_mu=16;
init_sigma=9.1;
nucleoids_tmp=round(normrnd(init_mu,init_sigma,ncells,1));  %draw starting values from normal distribution
nucleoids=nucleoids_tmp.*(nucleoids_tmp>0); %if nucleoid number should be negative, set it to 0

end

nucleoids_arxiv=nucleoids; %nucleoids_arxiv keeps track of nucleoid numbers for each cell cycle throughout the simulation

for i=1:cycles %loop over all simulated cell cycles

%simulate cell division
ratio_nucleoids = lognrnd(lognormal_mu,lognormal_sigma,ncells,1);   %sampling from the lognormal distribution to get the inheritance ratio for each cell
mother_nucleoids=min(nucleoids,round(nucleoids./(ratio_nucleoids+1)));  %determine the number of nucleoids the mother cell gets at division; draw from distribution of ratios, but mother number cannot be more than total
daughter_nucleoids=nucleoids-mother_nucleoids;  %determine daughter nucleoids given total and mother nucleoids

nucleoids_G1start = randsample([mother_nucleoids;daughter_nucleoids],ncells);   %reduce number of cells after division back to ncells


%simulate mitochondria synthesis during the cell cycle

nucleoids_new=round(nucleoids_G1start.*(synbase+syn4*(nucleoids_G1start<4)+syn7*(nucleoids_G1start<7)+syn10*(nucleoids_G1start<10)+syn13*(nucleoids_G1start<13))); %relative synthesis rate depens on number of nucleoids, described with three categories
nucleoids_arxiv=[nucleoids_arxiv,nucleoids_new];    %append nucleoids of this cell cycle to arxiv
nucleoids=nucleoids_new;   
end

n_nucleoids_all=[n_nucleoids_all;sum(nucleoids_arxiv)./ncells];
std_nucleoids_all=[std_nucleoids_all;std(nucleoids_arxiv)];
    
end

end