quasiHam.m

% Quasispecies simulation.  Includes mutation with Hamming distance.

function quasiHam

%close all
h = colormap(lines);

randpop = 0;

B = 8;
N = 2^B;         % size of mutation space (32)
ep = 0.1;      % average mutation rate: 0.1 to 0.01 typical
lam = 1;     % gain of Hamming distance

if randpop == 1
    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);

for yloop = 1:N
    for xloop = 1:N
        H(yloop,xloop) = hamming(yloop,xloop);
    end
end

%keyboard


Qtemp = 1./(1+ep*H);    %Mutation matrix
Qsum = sum(Qtemp,2);  

% Normalize along species
for yloop = 1:N
    for xloop = 1:N
        Q(yloop,xloop) = Qtemp(yloop,xloop)/Qsum(yloop);
    end
end

xp = 1:N;
alpha = 64;
ftemp = exp(-lam*H(alpha,:));   % Fitness landscape
sf = sum(ftemp);
f = ftemp/sf;

% Transition matrix
for yloop = 1:N
    for xloop = 1:N
        W(yloop,xloop) = f(yloop)*Q(yloop,xloop);
    end
end

tspan = [0 500];


[t,x] = ode45(@quasispec,tspan,x0);

Pop0
[sz,dum] = size(t);
Popend = sum(x(sz,:))

phistar = sum(f.*x(sz,:))


figure(1)
plot(f,'o-')
hold on
figure(1)
plot(x(sz,:),'o-r')
hold off
legend('fitness','population')

figure(2)
for loop = 1:N
    semilogx(t,x(:,loop),'Color',h(1+round(loop*63/N),:))
    hold on
end
hold off

figure(3)
for loop = 1:N
    plot(t,x(:,loop),'Color',h(1 + round(loop*63/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)

figure(4)
for loop = 1:N
    loglog(t,x(:,loop),'Color',h(1+round(loop*63/N),:))
    hold on
end
hold off

figure(5)
for loop = 1:N
    semilogy(t,x(:,loop),'Color',h(1+round(loop*63/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)

%keyboard



% Eigenvalues

[V,D] = eig(W);

max(D(:,1))

figure(6)
plot(V(:,2))
title('eigenvalues')


disp(' ')


    function yd = quasispec(t,y)
        
        phi = sum(f'.*y);   % Average fitness of population

        yd = W*y - phi*y;
        

    end     % end quasispec

end     % end quasi