diff --git a/quasiSpecAdj.m b/quasiSpecAdj.m new file mode 100644 index 0000000..fbc326d --- /dev/null +++ b/quasiSpecAdj.m @@ -0,0 +1,296 @@ +% Quasispecies and Replicator-Mutator simulation with adjacency matrix. + +function node = quasiSpecAdj + +%close all +h = colormap(lines); + +randpop = 1; % 0) = spike population; 1) = random population +mutype = 2; % 0) = Hamming; 1) = rand; 2) = network distance +fitype = 5; % 0) = Hamming; 1) = 2-peak; 2) = rand+gauss; 3) = freq-dep; 4) = Network distance 5) = degree + +B = 7; +N = 2^B; % size of mutation space (128) +lam = 1; % Hamming fitness only +gamma = 1; % freq-dep fitness (payoff matrix only) +relran = 0.025; % relative random contrib to fitness +time_expand = 50; + +ep = 0.05; % average mutation rate: 0.1 to 0.01 typical (0.4835) (0.0290) + +% Set up adjacency matrix + +node = makeER(N,0.05); +%node = makeSW(N,3,0.2); +%node = makeSF(N,3); +%node = makeglobal(N); + +[N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2,meandistance,diam] = clusterstats(node); +deg = degreedist(node); + +disp(' ') +displine('avg degree = ',avgdegree') +displine('number of clusters =', numclus) +displine('diameter = ',diam) +disp(' ') + +figure(1) +drawnet(node) +%DrawNetC(node) + +[A,degree,Lap] = adjacency(node); +dist = node2distance(node,100); +figure(2) +imagesc(dist) +colormap(jet) +colorbar +caxis([0 10]) +title('Distance Matrix') + +%keyboard + +%%%%% Set original population +if randpop == 1 + %rng(0); + x0temp = rand(1,N); % Initial population + sx = sum(x0temp); + x0 = x0temp/sx; +else + x0 = zeros(1,N); + x0(1) = 0.667; x0(2) = 0.333; +end + +Pop0 = sum(x0); + + +%%%%%% Set Hamming distance +for yloop = 1:N + for xloop = 1:N + H(yloop,xloop) = hamming(yloop-1,xloop-1); + end +end + + + +%%%%%%% Set Mutation matrix +if mutype == 0 % Hamming + Qtemp = 1./(1+H/ep); %Mutation matrix on Hamming + %Qtemp = exp(-H/(ep*50)); + Qsum = sum(Qtemp,2); + + % Normalize mutation among species + for yloop = 1:N + for xloop = 1:N + Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(xloop); + end + end + +end +if mutype == 1 % Random mutation + %rng(0); + S = stochasticmatrix(N); + Stemp = S - diag(diag(S)); + Qtemp = ep*Stemp; + sm = sum(Qtemp,2)'; + Q = Qtemp + diag(ones(1,N) - sm); +end +if mutype == 2 % Network distance + Qtemp = 1./(1+dist/ep); %Mutation matrix on Hamming + Qsum = sum(Qtemp,2); + + % Normalize mutation among species + for yloop = 1:N + for xloop = 1:N + Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(xloop); + end + end + +end + +%keyboard + + +%%%%%%% Set fitness landscape + +if fitype == 0 % Hamming + x = 1:N; + alpha = 84; + ftemp = exp(-lam*H(alpha,:)); % Fitness landscape + sf = sum(ftemp); + f = ftemp/sf; +end +if fitype == 1 % double peak and rand + %rng(1); + f = rand(1,N); + x = 1:N; + delg = 20; + sig1 = 1; + sig2 = 4; + g1 = gaussprob(x,(N/2 - delg),sig1); + g2 = 3*gaussprob(x,(N/2 + delg),sig2); + ftemp = relran*f + g1 + g2; + f = ftemp/sum(ftemp); +end +if fitype == 2 % rand + Gauss + %rng(0); + f = rand(1,N); + x = 1:N; + ftemp = relran*f + gauss((x-N/2)/2); % Fitness landscape + f = ftemp/sum(ftemp); +end +if fitype == 3 % frequency-dependent Hamming + avgdis = mean(mean(H)); + %payoff = exp(-gamma*(H - avgdis)); % payoff matrix + %payoff = H.^2; + %payoff = ones(size(H)); + payoff = exp(-gamma*H); +end +if fitype == 4 % Network Distance from node 64 + x = 1:N; + alpha = 64; + ftemp = exp(-lam*dist(alpha,:)); % Fitness landscape + sf = sum(ftemp); + f = ftemp/sf; +end +if fitype == 5 + ftemp = exp(-lam*deg); + sf = sum(ftemp); + f = ftemp/sf; +end + +%keyboard + +% Run time evolution +tspan = [0 1000]; +[t,x] = ode45(@quasispec,tspan,x0); + +Pop0 +[sz,dum] = size(t); +Popend = sum(x(sz,:)) + +phistar = sum(f.*x(sz,:)) % final average fitness + + +figure(3) +plot(f,'-') +hold on +figure(3) +plot(x(sz,:),'r') +hold off +legend('fitness','population') + +figure(4) +loglog(f,x(sz,:),'or') +xlabel('Fitness') +ylabel('Population') + +%keyboard + + + +% figure(4) +% for loop = 1:N +% semilogx(t,x(:,loop),'Color',h(round(loop*64/N),:),'LineWidth',1.25) +% hold on +% end +% hold off +% set(gcf,'Color','white') +% xlabel('Time','FontSize',14) +% ylabel('Population','FontSize',14) +% hh = gca; +% set(hh,'FontSize',14) +% title('Semilogx') +% +% figure(5) +% for loop = 1:N +% plot(t,x(:,loop),'Color',h(round(loop*64/N),:)) +% hold on +% end +% hold off +% title('Linear') +% +% figure(6) +% for loop = 1:N +% loglog(t,x(:,loop),'Color',h(round(loop*64/N),:)) +% hold on +% end +% hold off +% title('Log-Log') + +figure(7) +for loop = 1:N + semilogy(t,x(:,loop),'Color',h(round(loop*64/N),:)) + hold on +end +hold off +title('Semilogy') + + +% Eigenvalues + +[V,D] = eig(W); + +Lyap = max(D); + +minD = min(Lyap) +mxD = max(D(:,1)) +mnD = mean(Lyap) + + +figure(8) +%semilogy(abs(V(:,1))) +plot((V(:,1))) + +%histfixplot(Lyap,20,0,1.1*mxD); +figure(9) +plot(real(Lyap)) +title('Eigenvalue Spectrum') + +%keyboard + +disp(' ') + + +if fitype == 1 + xlo = N/2 - delg - 2*sig1; + xhi = N/2 - delg + 2*sig1; + fit44 = sum(f(xlo:xhi)) + pop44 = sum(x(sz,xlo:xhi))/Popend + xlo = N/2 + delg - 2*sig2; + xhi = N/2 + delg + 2*sig2; + fit84 = sum(f(xlo:xhi)) + pop84 = sum(x(sz,xlo:xhi))/Popend + +end + + + + function yd = quasispec(~,y) + + if fitype == 3 % frequency-dependent Hamming + for loop = 1:N + ftemp(loop) = sum(payoff(:,loop).*y); + end + f = time_expand*ftemp/sum(ftemp); + end + + + % Transition matrix + for yloop = 1:N + for xloop = 1:N + W(yloop,xloop) = f(yloop)*Q(yloop,xloop); + end + end + + phi = sum(f'.*y); % Average fitness of population + + yd = W*y - phi*y; + + + end % end quasispec + + + + +end % end quasiSpec +