%% Main code P_initial= 10^(4); % no. of phage particles P_increment= 10*P_initial; P_final= 10^(6)*P_initial; H=10^(8); % no. of host particles r=0.8; % reproductive potential of phage B=70; % burst size t=5; % threshold for cooperation to start penalty=0.985; % penalty before cooperation starts p = 1-penalty; TotalCount = 300; % total no. of different density data-points D_final = zeros(1,TotalCount); % Final Density of phage D_initial= zeros(1,TotalCount); % Initial Density of phage D_initial=D_initial-1; growthRate=zeros(1, TotalCount); for n=0:1:TotalCount temp_P=P_initial*10^(n*(log10(P_final)-log10(P_initial))/TotalCount); count=n+1; MOI=temp_P/H; % Phage To Host ratio, MOI_min ~ 0, MOI_max = 100 i=0; D=zeros(1,100000); L=zeros(1,100000); while i=t && iMOI+1000 % The poisson distribution falls sharply after the mean MOI, so breaking the loop after the probability of infection class i bbecomes greater than lambda and its conntribution goes to zero breakLimit=L(i+1)*P_final*H; if breakLimit<0.1 break; end end i=i+1; end D_final(count)=sum(D); D_initial(count)=temp_P; growthRate(count)=D_final/D_initial; end logGrowthRate=log(growthRate); %% Final Plots figure(1) semilogx(D_initial,logGrowthRate,'-r'); % For ancestor: p=0.985, B=70, r=0.8, t=5 %semilogx(D_initial,logGrowthRate,'-b'); % For DNAJ1: p=0.73, B=50, r=2, t=6 %semilogx(D_initial,logGrowthRate,'-.g'); % For DNAJ2: p=0.73, B=35, r=3, t=4 title(['p=',num2str(p), ', B=', num2str(B), ', b=', num2str(r), ', T=', num2str(t)]); xlabel('Initial \lambda density per ml'); ylabel('Growth rate'); ax1=gca; ax1.XTick=[10^5 10^6 10^7 10^8 10^9 10^10]; ax1.XAxisLocation='bottom'; ax1.YLim=[-5 2]; yPos=0; hold on plot(ax1.XLim,[yPos yPos],'-k','LineWidth',0.25);