function y=evol_rate(y,n) % this function calculates mutations for the phages and hosts %parameters mu=5.4*1e-8;%10^-6; %mutation rate for hosts nu=1.25*1e-4;%10^-5; %mutation rate for phages phi=5.4*10^-7;%10^-7; %back mutation for hosts sigma=0;%10^-8; %back mutation for phages % setting the final values of the integration as the intial values of monte carlo for i=1:(2^n) hj(i)=y(i); %only import y(end,:) hk(i)=y(2^n+i); hl(i)=y(2*(2^n)+i); hm(i)=y(3*(2^n)+i); hn(i)=y(4*(2^n)+i); end pj=y(n*(2^n)+1); pk=y(n*(2^n)+2); pl=y(n*(2^n)+3); pm=y(n*(2^n)+4); pn=y(n*(2^n)+5); %hosts %without back-mutation % evol_matH=[0 0 0 0 0; % mu 0 0 0 0; % mu^2 mu 0 0 0; % mu^3 mu^2 mu 0 0; % mu^4 mu^3 mu^2 mu 0]; %with back-mutation evol_matH=[0 phi phi^2 phi^3 phi^4; mu 0 phi phi^2 phi^3; mu^2 mu 0 phi phi^2; mu^3 mu^2 mu 0 phi; mu^4 mu^3 mu^2 mu 0]; %phages evol_matP=[0 sigma sigma^2 sigma^3 sigma^4; nu 0 sigma sigma^2 sigma^3; nu^2 nu 0 sigma sigma^2; nu^3 nu^2 nu 0 sigma; nu^4 nu^3 nu^2 nu 0]; %% Hosts % Calculating mutations for i=1:2^n H=[hj(i) hk(i) hl(i) hm(i) hn(i)]; M=[H;H;H;H;H]; M=evol_matH.*M; for j=1:n for k=1:n if M(j,k)>=5 M(j,k)=round(M(j,k)); else M(j,k)=poissrnd2(M(j,k)); end end end %%to apply the mutations for l=1:n for m=1:n H(l)=H(l)+M(l,m)-M(m,l); end if H(l)<0 H(l)=0; end end %for setting values back for hosts for m=1:n y((m-1)*(2^n)+i)=H(m); end end %% Phages P=[pj pk pl pm pn]; M2=[P;P;P;P;P]; M2=M2.*evol_matP; %calculating mutations for j=1:n for k=1:n if M2(j,k)>=5 M2(j,k)=round(M2(j,k)); else M2(j,k)=poissrnd2(M2(j,k)); end end end %applying mutations for l=1:n for m=1:n P(l)=P(l)+M2(l,m)-M2(m,l); end if P(l)<0 P(l)=0; end end %returning the new values for phages y(n*(2^n)+1)=P(1); y(n*(2^n)+2)=P(2); y(n*(2^n)+3)=P(3); y(n*(2^n)+4)=P(4); y(n*(2^n)+5)=P(5);