% ClvR: Deterministic model of the 2 locus, autosomal element cleaving 
% target gene at separate autosomal locus, plots multiple single runs

% Tobin Ivy (tivy@caltech.edu)
% Last Edit: 08/14/2018

clear all;

% Single run funcion format:
% [Pcb, Ptb, Pt, Pn] = ClvR_2lAA_multirun_support(FC,FCnc,HI,Ca,M,mut,i1,i2,gens)

gens = 2000;

% Plots at showing threshold at 50% haploinsufficiency, with 0.1607 being
% an approximation of the actual threshold release ratio
% [Pcb1, Ptb1, Pt1, Pn1] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.1607,0,gens);
% 
% [Pcb2, Ptb2, Pt2, Pn2] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.21,0,gens);
% [Pcb3, Ptb3, Pt3, Pn3] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.26,0,gens);
% [Pcb4, Ptb4, Pt4, Pn4] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.31,0,gens);
% 
% [Pcb5, Ptb5, Pt5, Pn5] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.11,0,gens);
% [Pcb6, Ptb6, Pt6, Pn6] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.06,0,gens);
% [Pcb7, Ptb7, Pt7, Pn7] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.01,0,gens);


% Plots at showing threshold at 100% haploinsufficiency (haplolethality),
% with 0.315954 being an approximation of the actual threshold release
% ratio
[Pcb1, Ptb1, Pt1, Pn1] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.315954,0,gens);

[Pcb2, Ptb2, Pt2, Pn2] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.36,0,gens);
[Pcb3, Ptb3, Pt3, Pn3] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.41,0,gens);
[Pcb4, Ptb4, Pt4, Pn4] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.46,0,gens);

[Pcb5, Ptb5, Pt5, Pn5] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.26,0,gens);
[Pcb6, Ptb6, Pt6, Pn6] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.21,0,gens);
[Pcb7, Ptb7, Pt7, Pn7] = ClvR_2lAA_multirun_support(0,0,1,0.9,0.9,0,0.16,0,gens);


% g- generations of simulation
g = linspace(0,gens,gens+1);

% Plot multiple runs
figure(4)
plot(g,Ptb1,'k',...
    g,Ptb2,'k',...
    g,Ptb3,'k',...
    g,Ptb4,'k',...
    g,Ptb5,'k',...
    g,Ptb6,'k',...
    g,Ptb7,'k','LineWidth', 3)
xlabel('Generation')
ylabel('Population Frequency (%)')
set(gca,'FontSize',20)
set(gca,'FontWeight','bold')
set(gca,'Fontname','Timesnewroman')
% legend('','')
title('100% Haploinsufficiency')
axis([0 20 0 100])