function []= IGEM()

%parameter analysis
MtoC=[];
Xvar=[];


for i=3:1:6
  Yvar=[];
  MtoCone=[];
    for j=0:1:3

par = loadparameters();

%initial state
V0=10^i;
C0=10^j;
CV0=0;
M0=0;
MV0=0;

y0=[V0; C0; CV0; M0; MV0];

%intervall of integration
t_span = [0 450];

%compute solution
[t,y]=ode23s(@(t,y) SelectionModel(t,y,par),t_span,y0);

      MtoCone=[MtoCone;y(450,4)./y(450,2)]; %calculate ratio at the end

      Yvar=[Yvar;C0]; %detection of various parameter values (here e.g. C0)
    end
  MtoC=[MtoC,MtoCone];
  Xvar=[Xvar;V0]; %detection of various parameter values (here e.g. V0)
end

    function dydt = SelectionModel(t,y,par)
        
        %load parameters for end of selection pressure after time point t
        if (t>250) par=loadparametersB();
        end
        
        %get actual values
        V=y(1);
        C=y(2);
        CV=y(3);
        M=y(4);
        MV=y(5);
        
        
        %compute derivations
        dVdt=-par.k_infM*V*M-par.k_infV*C*V-par.k_dest*V;
        dCdt=par.k_div*C-par.k_death*C-par.k_infV*C*V-par.k_mut*C+par.k_loss*CV;
        dCVdt=par.k_infV*C*V+par.k_divAinf*CV-par.k_deathAinf*CV-par.k_loss*CV;
        dMdt=par.k_mut*C+par.k_div*M-par.k_death*M-par.k_infM*M*V+par.k_loss*MV;
        dMVdt=par.k_infM*M*V+par.k_divAinf*MV-par.k_deathAinf*MV-par.k_loss*MV;
        
        dydt=[dVdt; dCdt; dCVdt; dMdt; dMVdt];
    end

    function par = loadparameters()
       par.k_dest=0; %virus destruction
       par.k_div=0.02;      %cell division
       par.k_divAinf=0.1;  %cell division after virus infection
       par.k_death=0.03;    %cell deatch rate
       par.k_deathAinf=par.k_death;  %cell death rate after virus infection
       par.k_mut=0.00001;      %mutation rate to get a better antibody
       par.k_infV=0.000001;     %virus infection rate for wt cell
       par.k_infM=0.0001;     %virus infection rate for mutated cell   
       par.k_loss=0.2;  %loss of viral genes (loss of stimulus)
    end

    function par = loadparametersB()  %parameters for end of selection pressure
       par=loadparameters();
       par.k_rep=0.02;      %cell division
       par.k_death=0.01;    %cell death 
    end

% get solution
        a=t
        V=y(:,1);
        C=y(:,2);
        CV=y(:,3);
        M=y(:,4);
        MV=y(:,5);

        
% plot results
% altogether:
figure(2); clf;
plot(t,y)
xlabel('time');
hleg1=legend('# V','# C','# CV','# M','# MV');

%detached:
figure(1); clf;

subplot(1,5,1);
plot(t,V,'-b'); xlim(t_span); yl = ylim; ylim([0 yl(2)]);
xlabel('time'); ylabel('# V');

subplot(1,5,2);
plot(t,C,'-g'); xlim(t_span); yl = ylim; ylim([0 yl(2)]);
xlabel('time'); ylabel('# C');

subplot(1,5,3);
plot(t,M,'-r'); xlim(t_span); yl = ylim; ylim([0 yl(2)]);
xlabel('time'); ylabel('# M');

subplot(1,5,4);
plot(t,CV,'-b'); xlim(t_span); yl = ylim; ylim([0 yl(2)]);
xlabel('time'); ylabel('# CV');

subplot(1,5,5);
plot(t,MV,'-r'); xlim(t_span); yl = ylim; ylim([0 yl(2)]);
xlabel('time'); ylabel('# MV');

%plot one dimensional parameter analysis:
figure(3); clf;
plot(Xvar,MtoCone);
xlabel('initial virus concentration');
ylabel('ratio mutated to WT cells');

%plot two dimensional parameter analysis
figure; surf(Xvar,Yvar,MtoC);
xlabel('initial virus concentration');
ylabel('initial wt cell concentration');
zlabel('ratio mutated to wt cells');


end