Add files via upload

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
+