quasiSpecAdj.m

% 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