% ClvR: Deterministic model of the 2 locus, autosomal element cleaving 
% target gene at separate autosomal locus across a range of fitness costs
% and release ratios

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

clear all;
tic

% T = transgene bearing autosmal locus 1
% N = null Cas9 copy of transgene
% W = wildtype bearing autosmal locus 1
% C = autosmal locus 2 with essential gene cleaved
% A = wildtype autosmal locus 2

% G1 = TTCC
% G2 = TTCA
% G3 = TTAA
% G4 = TNCC
% G5 = TNCA
% G6 = TNAA
% G7 = TWCC
% G8 = TWCA
% G9 = TWAA
% G10 = NNCC
% G11 = NNCA
% G12 = NNAA
% G13 = NWCC
% G14 = NWCA
% G15 = NWAA
% G17 = WWCA
% G18 = WWAA

%% Set parameters:

% Range of tested conditions:
% Nfc- number of fitness costs
% Nif- number of introduction frequencies
Nfc = 1001;
Nif = 801;

% Fitness Costs:
% FC- maximum fitness cost of T homozygote
% FCnc- proportion of fitness cost carried by the null cas9 mutant
% HI- maximum haploinsufficiency cost of T
% s(genotype)- fitness cost for a given genotype
FC = 0;
FCnc = 0;
HI = 1;

sG01M = linspace(0,FC,Nfc);
sG02M = linspace(0,FC,Nfc);
sG03M = linspace(0,FC,Nfc);
sG04M = linspace(0,FC,Nfc)/2+linspace(0,FCnc,Nfc)/2+linspace(0,HI,Nfc);
sG05M = linspace(0,FC,Nfc)/2+linspace(0,FCnc,Nfc)/2;
sG06M = linspace(0,FC,Nfc)/2+linspace(0,FCnc,Nfc)/2;
sG07M = linspace(0,FC,Nfc)/2+linspace(0,HI,Nfc);
sG08M = linspace(0,FC,Nfc)/2;
sG09M = linspace(0,FC,Nfc)/2;
sG10M = 0+linspace(0,FCnc,Nfc);
sG11M = 0+linspace(0,FCnc,Nfc);
sG12M = 0+linspace(0,FCnc,Nfc);
sG13M = 0+linspace(0,FCnc,Nfc)/2+linspace(0,HI,Nfc);
sG14M = 0+linspace(0,FCnc,Nfc)/2;
sG15M = 0+linspace(0,FCnc,Nfc)/2;
sG17M = linspace(0,0,Nfc)+linspace(0,HI,Nfc);
sG01F = linspace(0,FC,Nfc);
sG02F = linspace(0,FC,Nfc);
sG03F = linspace(0,FC,Nfc);
sG04F = linspace(0,FC,Nfc)/2+linspace(0,FCnc,Nfc)/2+linspace(0,HI,Nfc);
sG05F = linspace(0,FC,Nfc)/2+linspace(0,FCnc,Nfc)/2;
sG06F = linspace(0,FC,Nfc)/2+linspace(0,FCnc,Nfc)/2;
sG07F = linspace(0,FC,Nfc)/2+linspace(0,HI,Nfc);
sG08F = linspace(0,FC,Nfc)/2;
sG09F = linspace(0,FC,Nfc)/2;
sG10F = 0+linspace(0,FCnc,Nfc);
sG11F = 0+linspace(0,FCnc,Nfc);
sG12F = 0+linspace(0,FCnc,Nfc);
sG13F = 0+linspace(0,FCnc,Nfc)/2+linspace(0,HI,Nfc);
sG14F = 0+linspace(0,FCnc,Nfc)/2;
sG15F = 0+linspace(0,FCnc,Nfc)/2;
sG17F = linspace(0,0,Nfc)+linspace(0,HI,Nfc);

% Cas9 mediated events (a for activity in adults):
% CLEa- % chance of a single copy of T cleaving X into C
% MC- % chance of maternal carryover of Cas9 and gRNAs from a single copy
% TMC- % chance of maternal carryover of Cas9 and gRNAs from two copies

CLEa = 0.9;
MC = CLEa;
TMC = MC;


% Mutation:
% mutCas9- % chance of Cas9 mutating into inactivity
mutCas9 = 1*10^(-5);
mutCas9 = 0;

% Releases, run time, and matrix:
% IF- introduction frequency of desired genotype, as fraction of new total
% population
% IF2- introduction frequency of TNCY heterozygotes as a fraction of IF
% (logic is fraction of factory organisms that have a mutated drive)
% numGensRelease- number of sequential generational releases
% numGens- number of simulated generations
% FC_IF- matrix for storing the number of generations it took for the T to
% be present in >99% of individuals across all tested fitness costs and
% introduction frequencies
% sinput- vector of fitness costs used as input for support function for a
% single run

IF = linspace(0,0.8,Nif);
IF2 = 0.0*IF;
numGensRelease = 1;
numGens = 100;
FC_IF = 10000*ones(Nif,Nfc);

sinput = zeros(1, 32);

% Loop through fitness costs:
for f = 1:Nfc
    
    % Set run specific genotype fitness costs:
    sinput(1) = sG01M(f);
    sinput(2) = sG02M(f);
    sinput(3) = sG03M(f);
    sinput(4) = sG04M(f);
    sinput(5) = sG05M(f);
    sinput(6) = sG06M(f);
    sinput(7) = sG07M(f);
    sinput(8) = sG08M(f);
    sinput(9) = sG09M(f);
    sinput(10) = sG10M(f);
    sinput(11) = sG11M(f);
    sinput(12) = sG12M(f);
    sinput(13) = sG13M(f);
    sinput(14) = sG14M(f);
    sinput(15) = sG15M(f);
    sinput(16) = sG17M(f);
    sinput(17) = sG01F(f);
    sinput(18) = sG02F(f);
    sinput(19) = sG03F(f);
    sinput(20) = sG04F(f);
    sinput(21) = sG05F(f);
    sinput(22) = sG06F(f);
    sinput(23) = sG07F(f);
    sinput(24) = sG08F(f);
    sinput(25) = sG09F(f);
    sinput(26) = sG10F(f);
    sinput(27) = sG11F(f);
    sinput(28) = sG12F(f);
    sinput(29) = sG13F(f);
    sinput(30) = sG14F(f);
    sinput(31) = sG15F(f);
    sinput(32) = sG17F(f);
        
    % Loop through introduction frequencies:
    for i = 1:Nif
        FC_IF(i,f) = ClvR_2lAA_FCvIF_support(sinput,CLEa,MC,mutCas9,IF(i),IF2(i));
    end
end
toc

FC_IF2 = FC_IF;

for i = 1:Nif
    for f = 4:Nfc
        if FC_IF2(i,f) == 10000
        elseif FC_IF2(i,f) == FC_IF2(i,f-3)
            FC_IF2(i,f-3) = 10000;
        end
    end
end


% Heatmap thresholds
% Plot haploinsufficiency cost of T vs. release ratio of transgenics
figure(3)
h = image(sG04M*100,IF*100,FC_IF2);
xlabel('Haploinsufficiency (%)')
ylabel('Release Ratio (%)')
set(gca,'FontSize',20)
set(gca,'FontWeight','bold')
set(gca,'Fontname','Arial')
set(gca,'YDir','normal')
axis([0 100 0 50])

% Adjust colormap (so that no takeover is white) 
z = 1:numGens;
color = [jet(length(z));1 1 1];
set(gcf, 'ColorMap', color);
colorbar


% Heatmap full
% Plot haploinsufficiency cost of T vs. release ratio of transgenics
figure(1)
h = image(sG04M*100,IF*100,FC_IF);
xlabel('Haploinsufficiency (%)')
ylabel('Release Ratio (%)')
set(gca,'FontSize',20)
set(gca,'FontWeight','bold')
set(gca,'Fontname','Arial')
set(gca,'YDir','normal')
axis([0 100 0 50])

% Adjust colormap (so that no takeover is white) 
z = 1:numGens;
color = [jet(length(z));1 1 1];
set(gcf, 'ColorMap', color);
colorbar



load gong.mat;
soundsc(y);