diff --git a/quasiSpecAdjTst.m b/quasiSpecAdjTst.m deleted file mode 100644 index fbc326d..0000000 --- a/quasiSpecAdjTst.m +++ /dev/null @@ -1,296 +0,0 @@ -% 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 -