%%Singh, S Supplementary file 5: Monte Carlo simulations updated to include %third intermediate "S" phase like expression state %One parameter set module extracted from "StochProfGUI" from Wang, Lixin, and Kevin A. Janes. Nature Protocols 8.2 (2013): 282. %% Parameters %Currently set to simulate 3-state cell cycle with the parameters listed in %Supplementary Figure S10 clear i = 1; %counter for first variable sweep j = 1; %counter for second variable sweep k = 1; %counter for number of Monte Carlo simulations CI = 0.90; %confidence interval for p values maxiters=50; %number of simulations numexpts=[5:2:25]; %number of 10-cell samples collected (variable sweep) CVtest=[0:0.05:0.5]; %range of heterogeneities to assess (variable sweep) numcells=10; %number of cells per sample %Parameter F decides two state vs 3-state simulation (used by multinomial %function) %average frequency of cells in different cell states:3 state model,cell-cycle distribution F=[0.65,0.1,0.25]; %First value is G1 proportion, then S proportion, then G2/M %F=[0.7,0.0,0.3]; %Sphase frequency set to 0 to re-capitulate 2 state model D=3; % fold-change of expression between high/low or G1/G2-M states CVref=0.2; %biological variance in samples %% MC sims for k = 1:maxiters i=1; %numexpts counter for numexptstemp = numexpts j=1; %CVtest counter for CVtesttemp = CVtest clear genemeas refmeas controlmeas for m=1:numexptstemp genemeastemp=0; refmeastemp=0; controlmeastemp=0; hititer = mnrnd(numcells, F); %changed this from binomial to multinomial sampling to pick 3 frequencies for n=1:hititer(3) %this is the high expression state or G2/M genemeastemp=genemeastemp+lognrnd(log(D),sqrt(log(CVtesttemp^2+1))); refmeastemp=refmeastemp+lognrnd(0,sqrt(log(CVref^2+1))); controlmeastemp=controlmeastemp+lognrnd(0,sqrt(log(CVtesttemp^2+1))); end for n=1:hititer(2) %this is the S phase distribution - bounds set according to the other two states genemeastemp=genemeastemp+unifrnd((0+(sqrt(log(CVtesttemp^2+1)))),(log(D)-(sqrt(log(CVtesttemp^2+1))))); refmeastemp=refmeastemp+unifrnd((0-(sqrt(log(CVref^2+1)))),0+(sqrt(log(CVref^2+1)))); controlmeastemp=controlmeastemp+unifrnd((0-(sqrt(log(CVtesttemp^2+1)))),0+(sqrt(log(CVtesttemp^2+1)))); end for n=1:hititer(1) %this is the high expression state or G1 genemeastemp=genemeastemp+lognrnd(0,sqrt(log(CVtesttemp^2+1))); refmeastemp=refmeastemp+lognrnd(0,sqrt(log(CVref^2+1))); controlmeastemp=controlmeastemp+lognrnd(0,sqrt(log(CVtesttemp^2+1))); end genemeas(m) = genemeastemp; refmeas(m) = refmeastemp; controlmeas(m) = controlmeastemp; end [h,p(k,i,j),kstat(k,i,j)]=kstest(log(genemeas)',[log(genemeas)' normcdf(log(genemeas)', ... mean(log(genemeas)),std(log(refmeas)))],0.05,'unequal'); [h,pcontrol(k,i,j),kstatcontrol(k,i,j)]=kstest(log(controlmeas)',[log(controlmeas)' ... normcdf(log(controlmeas)',mean(log(controlmeas)),std(log(refmeas)))], ... 0.05,'unequal'); j=j+1; end i=i+1; end sprintf('%d percent complete.',round(k/maxiters*100)) end %% Summarize for i=1:length(numexpts) for j=1:length(CVtest) if median(p(:,i,j))<0.05 && median(pcontrol(:,i,j))>=0.05 summary(i,j)=0; %effective stochastic sampling elseif median(p(:,i,j))>=0.05 summary(i,j)=-1; %false negative elseif median(pcontrol(:,i,j))<0.05 && CVtest(j)