%%%%%load model%%%%%%%%
sbioloadproject IGEM-selection-system

%%%%set SolverOtions%%%%%%%%%%
configsetObj = m1.getconfigset('active');
configsetObj.SolverType = 'ssa';
configsetObj.SolverOptions.LogDecimation = 100;

%%%%%load Species%%%
C=m1.Species(1);
V=m1.Species(2);
M=m1.Species(3);
CV=m1.Species(4);
MV=m1.Species(5);

%%%%%%%%load parameters%%%%%%%%%%%%
k_death=m1.Reactions(1).KineticLaw.Parameters;
k_rep=m1.Reactions(2).KineticLaw.Parameters;
k_mut=m1.Reactions(3).KineticLaw.Parameters;
k_intC=m1.Reactions(4).KineticLaw.Parameters;
k_intM=m1.Reactions(5).KineticLaw.Parameters;
k_rep1=m1.Reactions(6).KineticLaw.Parameters;
k_death1=m1.Reactions(7).KineticLaw.Parameters;
k_repstim=m1.Reactions(8).KineticLaw.Parameters;
k_deathstim=m1.Reactions(9).KineticLaw.Parameters;
k_repstim1=m1.Reactions(10).KineticLaw.Parameters;
k_deathstim1=m1.Reactions(11).KineticLaw.Parameters;
k_loss=m1.Reactions(12).KineticLaw.Parameters;
k_loss1=m1.Reactions(13).KineticLaw.Parameters;

%%%%set patameter values%%%%%
set(C,'InitialAmount',1000);
set(V,'InitialAmount',10000);
set(k_rep,'Value',0.02);
set(k_rep1,'Value',0.02);
set(k_death,'Value',0.03);
set(k_death1,'Value',0.03);
set(k_mut,'Value',1e-5);
set(k_intC,'Value',1e-6);
set(k_intM,'Value',1e-4);
set(k_repstim,'Value',0.1);
set(k_repstim1,'Value',0.1);
set(k_loss,'Value',0.2);
set(k_loss1,'Value',0.2);
set(k_deathstim,'Value',0.03);
set(k_deathstim1,'Value',0.03);

%%%%%%%
numRuns = 100; %number of runs

y1=[];
num=1;
z1=[];
e1=[];
for l=2:1:5
    y2=[];
    x1=[];
    set(k_death,'Value',0.01*l);
    set(k_death1,'Value',0.01*l);
for k=2:0.5:4
    set(C,'InitialAmount',1*10^k);
    %set(k_rep1,'Value',0.01*k);
    allpos=[];
    for j=1:3
        simdata = sbioensemblerun(m1, numRuns);
        [t, M1] = selectbyname(simdata, 'M'); %select M
        [t, C1] = selectbyname(simdata, 'C'); %select C
        pos=0;
        for i=1:numRuns
            %figure(j);
            %hold on;
            %plot(t{i}, M1{i});
            %plot(t{i}, C1{i},'r');
            %figure(j+3);
            %plot(t{i}, M1{i}./C1{i});
            if (M1{i}(length(M1{i}))./C1{i}(length(C1{i}))>1) pos=pos+1;end %count positiv results (successful selections)
        end
        allpos=[allpos pos];
        %ylim([0 3e10]);
        %legend('[M]','[C]');
        %title(['initial amount of V:' num2str(get(V,'InitialAmount'))],'HorizontalAlignment','right','Fontsize',14);
        %saveas(figure(j),['initialV00' num2str(num) '.jpg'])
        %close all;
        %num=num+1;
    end
    x1=[x1 get(C,'InitialAmount')]
    y2=[y2 mean(allpos)]
    e1=[e1 std(allpos)];
end
z1=[z1 get(k_death,'Value')]
y1=[y1; y2]
end


%%%3d plot%%%%%
figure;
surf(x1,z1,y1);
xlabel('initial amount WT cell');
ylabel('k_death');
zlabel('% of succcess');

%%%%plot%%%%%
figure;
hold on
h=errorbar(x1,y1,e1,'r');
xlabel('k_rep');
ylabel('% of success');
