%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MEDUSA POPULATION DYNAMICS:            %%
%% John Marshall, 5/26/2014               %%
%% john.marshall@imperial.ac.uk           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. DETERMINISTIC MODEL:                %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% Set default parameters:

sAB = 0;            % fitness cost of XAYB individuals
sAa = 0;            % fitness cost of XAXa individuals
releasePr = 0.5;    % release proportion (XAYB individuals)
numGensRelease = 2; % number of generations over which release occurs
numGens = 20;       % number of generations simulated

% Initial conditions (first geneneration):

pAB(1) = releasePr;         % initial proportion of released XAYB males
pab(1) = 0.5*(1-releasePr); % initial proportion of XaYb (wild-type) males
pAa(1) = 0;                 % initial proportion of XAXa females
paa(1) = 0.5*(1-releasePr); % initial proportion of XaXa (wild-type) females

% Difference equations (second generation onwards):
% For more details, see: Marshall JM, Hay BA (2014) Medusa: A novel gene 
% drive system for confined suppression of insect populations. PLoS ONE

for k = 1:numGens
    % Un-normalized genotype frequencies:
    pABTemp(k+1) = 0.25*pAB(k)*pAa(k);
    pabTemp(k+1) = 0.5 *pab(k)*paa(k);
    pAaTemp(k+1) = 0.5 *pAB(k)*paa(k);
    paaTemp(k+1) = 0.5 *pab(k)*paa(k);
    
    % Normalizing constant (takes into account fitness costs):
    sigma(k+1) = pABTemp(k+1)*(1-sAB) + pabTemp(k+1) + pAaTemp(k+1)*(1-sAa) + paaTemp(k+1);
    
    % Normalized genotype frequencies:
    if k < numGensRelease
        % Taking into account new releases:
        pAB(k+1) = releasePr + (1-releasePr)*pABTemp(k+1)*(1-sAB)/sigma(k+1);
        pab(k+1) = (1-releasePr)*pabTemp(k+1)/sigma(k+1);
        pAa(k+1) = (1-releasePr)*pAaTemp(k+1)*(1-sAa)/sigma(k+1);
        paa(k+1) = (1-releasePr)*paaTemp(k+1)/sigma(k+1);
    else
        % In the absence of new releases:
        pAB(k+1) = pABTemp(k+1)*(1-sAB)/sigma(k+1);        
        pab(k+1) = pabTemp(k+1)/sigma(k+1);
        pAa(k+1) = pAaTemp(k+1)*(1-sAa)/sigma(k+1);
        paa(k+1) = paaTemp(k+1)/sigma(k+1);
    end
end

% Derived quantities (for plotting):

gen = linspace(1,numGens+1,numGens+1); % Generations of simulation
PrTransgenic = 100*(pAa+pAB);          % Percentage of population that are transgenic over time
PrFemale = 100*(paa+pAa);              % Percentage of population that are female over time

% Plot of Medusa spread over time:

figure(1)
plot(gen,PrTransgenic,'g',gen,PrFemale,'b')
xlabel('Generation')
ylabel('Population frequency (%)')
title('Medusa spread, deterministic model, s=0, two 50% releases')
legend('Transgenic','Female')
axis([1 20 0 100])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. SINGLE POPULATION STOCHASTIC MODEL: %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% Mosquito life cycle parameters:

beta    = 32;        % eggs laid per day by female mosquito
beta_Aa = beta;      % eggs laid per day by XAXa mosquito
beta_aa = beta;      % eggs laid per day by XaXa mosquito
mu_m    = 0.123;     % daily mortality of adult mosquito
mu_m_AB = mu_m;      % daily mortality of adult XAYB mosquito
mu_m_Aa = mu_m;      % daily mortality of adult XAXa mosquito
mu_m_aa = mu_m;      % daily mortality of adult XaXa mosquito
mu_m_ab = mu_m;      % daily mortality of adult XaYb mosquito
rm      = 1.096;     % daily mosquito population growth rate
T0      = 1;         % duration of egg life stage (days)
TL      = 14;        % duration of larval life stage (days)
TP      = 1;         % duration of pupal life stage (days)
M_eq    = 10000;     % equilibrium adult population size

% Derived mosquito life cycle parameters:

g       = T0 + TL + TP + 1/mu_m;  % average generation time (days)
Rm      = rm^g;                   % mosquito population growth rate (per generation)
mu_l    = 1 - (Rm*mu_m/((beta/2)*(1-mu_m)))^(1/(T0+TL+TP));        
                                  % daily mortality of larvae
theta0  = (1 - mu_l)^T0;          % probability of surviving egg stage
thetaL  = (1 - mu_l)^TL;          % probability of surviving larval stage
thetaP  = (1 - mu_l)^TP;          % probability of surviving pupal stage
alpha   = (beta*theta0*(M_eq/2)/(Rm-1))*((1-thetaL/Rm)/(1 - (thetaL/Rm)^(1/TL)));
                                  % parameter influencing the strenth of larval density-dependence
L_eq    = round(alpha*(Rm-1));    % equilibrium larval population size

% Initial conditions:

LAB(1:50)    = 0;       % initial number of XAYB male larvae 
LAa(1:50)    = 0;       % initial number of XAXa female larvae
Lab(1:50)    = L_eq/2;  % initial number of XaYb (wild-type) male larvae
Laa(1:50)    = L_eq/2;  % initial number of XaXa (wild-type) female larvae
L(1:50)      = L_eq;    % initial number of larvae

MAB(1:50)    = 0;       % initial number of XAYB male adults
Mab(1:50)    = M_eq/2;  % initial number of XaYb (wild-type) male adults

MAa_AB(1:50) = 0;       % initial number of XAXa female adults mated to XAYB males
MAa_ab(1:50) = 0;       % initial number of XAXa female adults mated to XaYb males
Maa_AB(1:50) = 0;       % initial number of XaXa female adults mated to XAYB males
Maa_ab(1:50) = M_eq/2;  % initial number of XaXa female adults mated to XaYb males

max_t = 600;            % Days of simulation

% Difference equations (second day onwards):
% For more details, see: Marshall JM, Hay BA (2014) Medusa: A novel gene 
% drive system for confined suppression of insect populations. PLoS ONE

for t = 51:max_t

    % Number of larvae having each genotype the next day:
    thetaLA = thetaL;
    for i=1:TL
        thetaLA = thetaLA*((alpha/(alpha+L(t-i)))^(1/TL));
    end
    
    Y1 = poissrnd(beta_aa*Maa_ab(t-T0));
    X1 = binornd(Y1,(1/2));

    LAB(t) = max(0,binornd(LAB(t-1),(1-mu_l)*((alpha/(alpha+L(t-1)))^(1/TL))) + binornd(binornd(poissrnd(beta_Aa*MAa_AB(t-T0)),(1/4)),theta0) - binornd(binornd(poissrnd(beta_Aa*MAa_AB(t-T0-TL)),(1/4)),theta0*thetaLA));
    LAa(t) = max(0,binornd(LAa(t-1),(1-mu_l)*((alpha/(alpha+L(t-1)))^(1/TL))) + binornd(binornd(poissrnd(beta_aa*Maa_AB(t-T0)),(1/2)),theta0) - binornd(binornd(poissrnd(beta_aa*Maa_AB(t-T0-TL)),(1/2)),theta0*thetaLA));
    Lab(t) = max(0,binornd(Lab(t-1),(1-mu_l)*((alpha/(alpha+L(t-1)))^(1/TL))) + binornd(X1,theta0) - binornd(binornd(poissrnd(beta_aa*Maa_ab(t-T0-TL)),(1/2)),theta0*thetaLA));
    Laa(t) = max(0,binornd(Laa(t-1),(1-mu_l)*((alpha/(alpha+L(t-1)))^(1/TL))) + binornd(Y1-X1,theta0) - binornd(binornd(poissrnd(beta_aa*Maa_ab(t-T0-TL)),(1/2)),theta0*thetaLA));
    L(t) = LAB(t-1) + LAa(t-1) + Lab(t-1) + Laa(t-1);
    
    % Number of adults having each genotype the next day:
    thetaLA = thetaL;
    for i=1:TL
        thetaLA = thetaLA*((alpha/(alpha+L(t-i-TP)))^(1/TL));
    end

    Y2 = poissrnd(beta_aa*Maa_ab(t-T0-TL-TP));
    X2 = binornd(Y2,(1/2));
    
    MAB(t) = max(0,binornd(MAB(t-1),(1-mu_m_AB)) + binornd(binornd(poissrnd(beta_Aa*MAa_AB(t-T0-TL-TP)),(1/4)),theta0*thetaLA*thetaP*(1-mu_m_AB)));
    Mab(t) = max(0,binornd(Mab(t-1),(1-mu_m_ab)) + binornd(X2,theta0*thetaLA*thetaP*(1-mu_m_ab)));

    Z1 = binornd(binornd(poissrnd(beta_aa*Maa_AB(t-T0-TL-TP)),(1/2)),theta0*thetaLA*thetaP*(1-mu_m_Aa));
    Z2 = binornd(Y2-X2,theta0*thetaLA*thetaP*(1-mu_m_aa));
    
    MAa_AB(t) = max(0,binornd(MAa_AB(t-1),(1-mu_m_Aa)) + binornd(Z1,max(0,MAB(t-1)/(MAB(t-1)+Mab(t-1)))));
    MAa_ab(t) = max(0,binornd(MAa_ab(t-1),(1-mu_m_Aa)) + binornd(Z1,max(0,Mab(t-1)/(MAB(t-1)+Mab(t-1)))));
    Maa_AB(t) = max(0,binornd(Maa_AB(t-1),(1-mu_m_aa)) + binornd(Z2,max(0,MAB(t-1)/(MAB(t-1)+Mab(t-1)))));
    Maa_ab(t) = max(0,binornd(Maa_ab(t-1),(1-mu_m_aa)) + binornd(Z2,max(0,Mab(t-1)/(MAB(t-1)+Mab(t-1)))));
    
    % 6 consecutive releases of 10,000 adult XAYB mosquitoes every half generation:
    if ((t==51)|(t==51+14*1)|(t==51+14*2)|(t==51+14*3)|(t==51+14*4)|(t==51+14*5))
        MAB(t) = MAB(t) + M_eq;
    end    
end

% Derived quantities (for plotting):

t             = linspace(1,max_t,max_t);           % Days of simulation
Mm            = Mab + MAB;                         % Number of male adult mosquitoes over time
Mf            = MAa_AB + MAa_ab + Maa_AB + Maa_ab; % Number of female adult mosquitoes over time
M             = Mm + Mf;                           % Total number of adult mosquitoes over time
Mtransgenic   = MAa_AB + MAa_ab + MAB;             % Number of transgenic adult mosquitoes over time
Mf_transgenic = MAa_AB + MAa_ab;                   % Number of female transgenic adult mosquitoes over time

% Plot of Medusa spread over time:

figure(2)
plot(t(51:551)-51,M(51:551),'r', t(51:551)-51,Mtransgenic(51:551),'g', t(51:551)-51,Mf(51:551),'b');
xlim([0 500])
ylim([0 22500])
xlabel('Day')
ylabel('Population size')
title('Medusa spread, stochastic model, 6 releases of Medusa males')
legend('Total pop','Transgenic','Female')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3. TWO-POPULATION STOCHASTIC MODEL:    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% Mosquito life cycle parameters:

beta    = 32;        % eggs laid per day by female mosquito
beta_Aa = beta;      % eggs laid per day by XAXa mosquito
beta_aa = beta;      % eggs laid per day by XaXa mosquito
mu_m    = 0.123;     % daily mortality of adult mosquito
mu_m_AB = mu_m;      % daily mortality of adult XAYB mosquito
mu_m_Aa = mu_m;      % daily mortality of adult XAXa mosquito
mu_m_aa = mu_m;      % daily mortality of adult XaXa mosquito
mu_m_ab = mu_m;      % daily mortality of adult XaYb mosquito
rm      = 1.096;     % daily mosquito population growth rate
T0      = 1;         % duration of egg life stage (days)
TL      = 14;        % duration of larval life stage (days)
TP      = 1;         % duration of pupal life stage (days)
M_eq    = 10000;     % equilibrium adult population size

% Derived mosquito life cycle parameters:

g       = T0 + TL + TP + 1/mu_m;  % average generation time (days)
Rm      = rm^g;                   % mosquito population growth rate (per generation)
mu_l    = 1 - (Rm*mu_m/((beta/2)*(1-mu_m)))^(1/(T0+TL+TP));        
                                  % daily mortality of larvae
theta0  = (1 - mu_l)^T0;          % probability of surviving egg stage
thetaL  = (1 - mu_l)^TL;          % probability of surviving larval stage
thetaP  = (1 - mu_l)^TP;          % probability of surviving pupal stage
alpha   = (beta*theta0*(M_eq/2)/(Rm-1))*((1-thetaL/Rm)/(1 - (thetaL/Rm)^(1/TL)));
                                  % parameter influencing the strenth of larval density-dependence
L_eq    = round(alpha*(Rm-1));    % equilibrium larval population size
migrate = 0.01/g;                 % migration rate between populations

% Initial conditions (population 1):

LAB1(1:50)    = 0;       % initial number of XAYB male larvae 
LAa1(1:50)    = 0;       % initial number of XAXa female larvae
Lab1(1:50)    = L_eq/2;  % initial number of XaYb (wild-type) male larvae
Laa1(1:50)    = L_eq/2;  % initial number of XaXa (wild-type) female larvae
L1(1:50)      = L_eq;    % initial number of larvae

MAB1(1:50)    = 0;       % initial number of XAYB male adults
Mab1(1:50)    = M_eq/2;  % initial number of XaYb (wild-type) male adults

MAa_AB1(1:50) = 0;       % initial number of XAXa female adults mated to XAYB males
MAa_ab1(1:50) = 0;       % initial number of XAXa female adults mated to XaYb males
Maa_AB1(1:50) = 0;       % initial number of XaXa female adults mated to XAYB males
Maa_ab1(1:50) = M_eq/2;  % initial number of XaXa female adults mated to XaYb males

% Initial conditions (population 2):

LAB2(1:50)    = 0;       % initial number of XAYB male larvae 
LAa2(1:50)    = 0;       % initial number of XAXa female larvae
Lab2(1:50)    = L_eq/2;  % initial number of XaYb (wild-type) male larvae
Laa2(1:50)    = L_eq/2;  % initial number of XaXa (wild-type) female larvae
L2(1:50)      = L_eq;    % initial number of larvae

MAB2(1:50)    = 0;       % initial number of XAYB male adults
Mab2(1:50)    = M_eq/2;  % initial number of XaYb (wild-type) male adults

MAa_AB2(1:50) = 0;       % initial number of XAXa female adults mated to XAYB males
MAa_ab2(1:50) = 0;       % initial number of XAXa female adults mated to XaYb males
Maa_AB2(1:50) = 0;       % initial number of XaXa female adults mated to XAYB males
Maa_ab2(1:50) = M_eq/2;  % initial number of XaXa female adults mated to XaYb males

max_t = 2100;           % Days of simulation

% Difference equations (second day onwards):
% For more details, see: Marshall JM, Hay BA (2014) Medusa: A novel gene 
% drive system for confined suppression of insect populations. PLoS ONE

for t = 51:max_t

    % Population 1:

    % Number of larvae having each genotype the next day:
    thetaLA = thetaL;
    for i=1:TL
        thetaLA = thetaLA*((alpha/(alpha+L1(t-i)))^(1/TL));
    end
    
    Y1 = poissrnd(beta_aa*Maa_ab1(t-T0));
    X1 = binornd(Y1,(1/2));

    LAB1(t) = max(0,binornd(LAB1(t-1),(1-mu_l)*((alpha/(alpha+L1(t-1)))^(1/TL))) + binornd(binornd(poissrnd(beta_Aa*MAa_AB1(t-T0)),(1/4)),theta0) - binornd(binornd(poissrnd(beta_Aa*MAa_AB1(t-T0-TL)),(1/4)),theta0*thetaLA));
    LAa1(t) = max(0,binornd(LAa1(t-1),(1-mu_l)*((alpha/(alpha+L1(t-1)))^(1/TL))) + binornd(binornd(poissrnd(beta_aa*Maa_AB1(t-T0)),(1/2)),theta0) - binornd(binornd(poissrnd(beta_aa*Maa_AB1(t-T0-TL)),(1/2)),theta0*thetaLA));
    Lab1(t) = max(0,binornd(Lab1(t-1),(1-mu_l)*((alpha/(alpha+L1(t-1)))^(1/TL))) + binornd(X1,theta0) - binornd(binornd(poissrnd(beta_aa*Maa_ab1(t-T0-TL)),(1/2)),theta0*thetaLA));
    Laa1(t) = max(0,binornd(Laa1(t-1),(1-mu_l)*((alpha/(alpha+L1(t-1)))^(1/TL))) + binornd(Y1-X1,theta0) - binornd(binornd(poissrnd(beta_aa*Maa_ab1(t-T0-TL)),(1/2)),theta0*thetaLA));
    L1(t) = LAB1(t-1) + LAa1(t-1) + Lab1(t-1) + Laa1(t-1);
    
    % Number of adults having each genotype the next day:
    thetaLA = thetaL;
    for i=1:TL
        thetaLA = thetaLA*((alpha/(alpha+L1(t-i-TP)))^(1/TL));
    end

    Y2 = poissrnd(beta_aa*Maa_ab1(t-T0-TL-TP));
    X2 = binornd(Y2,(1/2));
    
    MAB1(t) = max(0,binornd(MAB1(t-1),(1-mu_m_AB)) + binornd(binornd(poissrnd(beta_Aa*MAa_AB1(t-T0-TL-TP)),(1/4)),theta0*thetaLA*thetaP*(1-mu_m_AB)));
    Mab1(t) = max(0,binornd(Mab1(t-1),(1-mu_m_ab)) + binornd(X2,theta0*thetaLA*thetaP*(1-mu_m_ab)));

    Z1 = binornd(binornd(poissrnd(beta_aa*Maa_AB1(t-T0-TL-TP)),(1/2)),theta0*thetaLA*thetaP*(1-mu_m_Aa));
    Z2 = binornd(Y2-X2,theta0*thetaLA*thetaP*(1-mu_m_aa));
    
    MAa_AB1(t) = max(0,binornd(MAa_AB1(t-1),(1-mu_m_Aa)) + binornd(Z1,max(0,MAB1(t-1)/(MAB1(t-1)+Mab1(t-1)))));
    MAa_ab1(t) = max(0,binornd(MAa_ab1(t-1),(1-mu_m_Aa)) + binornd(Z1,max(0,Mab1(t-1)/(MAB1(t-1)+Mab1(t-1)))));
    Maa_AB1(t) = max(0,binornd(Maa_AB1(t-1),(1-mu_m_aa)) + binornd(Z2,max(0,MAB1(t-1)/(MAB1(t-1)+Mab1(t-1)))));
    Maa_ab1(t) = max(0,binornd(Maa_ab1(t-1),(1-mu_m_aa)) + binornd(Z2,max(0,Mab1(t-1)/(MAB1(t-1)+Mab1(t-1)))));
    
    % 6 consecutive releases of 10,000 adult XAYB mosquitoes every half generation:
    if ((t==51)|(t==51+14*1)|(t==51+14*2)|(t==51+14*3)|(t==51+14*4)|(t==51+14*5))
        MAB1(t) = MAB1(t) + M_eq;
    end
    
    % Population 2:

    % Number of larvae having each genotype the next day:    
    thetaLA = thetaL;
    for i=1:TL
        thetaLA = thetaLA*((alpha/(alpha+L2(t-i)))^(1/TL));
    end
    
    Y1 = poissrnd(beta_aa*Maa_ab2(t-T0));
    X1 = binornd(Y1,(1/2));

    LAB2(t) = max(0,binornd(LAB2(t-1),(1-mu_l)*((alpha/(alpha+L2(t-1)))^(1/TL))) + binornd(binornd(poissrnd(beta_Aa*MAa_AB2(t-T0)),(1/4)),theta0) - binornd(binornd(poissrnd(beta_Aa*MAa_AB2(t-T0-TL)),(1/4)),theta0*thetaLA));
    LAa2(t) = max(0,binornd(LAa2(t-1),(1-mu_l)*((alpha/(alpha+L2(t-1)))^(1/TL))) + binornd(binornd(poissrnd(beta_aa*Maa_AB2(t-T0)),(1/2)),theta0) - binornd(binornd(poissrnd(beta_aa*Maa_AB2(t-T0-TL)),(1/2)),theta0*thetaLA));
    Lab2(t) = max(0,binornd(Lab2(t-1),(1-mu_l)*((alpha/(alpha+L2(t-1)))^(1/TL))) + binornd(X1,theta0) - binornd(binornd(poissrnd(beta_aa*Maa_ab2(t-T0-TL)),(1/2)),theta0*thetaLA));
    Laa2(t) = max(0,binornd(Laa2(t-1),(1-mu_l)*((alpha/(alpha+L2(t-1)))^(1/TL))) + binornd(Y1-X1,theta0) - binornd(binornd(poissrnd(beta_aa*Maa_ab2(t-T0-TL)),(1/2)),theta0*thetaLA));
    L2(t) = LAB2(t-1) + LAa2(t-1) + Lab2(t-1) + Laa2(t-1);
    
    % Number of adults having each genotype the next day:
    thetaLA = thetaL;
    for i=1:TL
        thetaLA = thetaLA*((alpha/(alpha+L2(t-i-TP)))^(1/TL));
    end

    Y2 = poissrnd(beta_aa*Maa_ab2(t-T0-TL-TP));
    X2 = binornd(Y2,(1/2));
    
    MAB2(t) = max(0,binornd(MAB2(t-1),(1-mu_m_AB)) + binornd(binornd(poissrnd(beta_Aa*MAa_AB2(t-T0-TL-TP)),(1/4)),theta0*thetaLA*thetaP*(1-mu_m_AB)));
    Mab2(t) = max(0,binornd(Mab2(t-1),(1-mu_m_ab)) + binornd(X2,theta0*thetaLA*thetaP*(1-mu_m_ab)));

    Z1 = binornd(binornd(poissrnd(beta_aa*Maa_AB2(t-T0-TL-TP)),(1/2)),theta0*thetaLA*thetaP*(1-mu_m_Aa));
    Z2 = binornd(Y2-X2,theta0*thetaLA*thetaP*(1-mu_m_aa));
    
    MAa_AB2(t) = max(0,binornd(MAa_AB2(t-1),(1-mu_m_Aa)) + binornd(Z1,max(0,MAB2(t-1)/(MAB2(t-1)+Mab2(t-1)))));
    MAa_ab2(t) = max(0,binornd(MAa_ab2(t-1),(1-mu_m_Aa)) + binornd(Z1,max(0,Mab2(t-1)/(MAB2(t-1)+Mab2(t-1)))));
    Maa_AB2(t) = max(0,binornd(Maa_AB2(t-1),(1-mu_m_aa)) + binornd(Z2,max(0,MAB2(t-1)/(MAB2(t-1)+Mab2(t-1)))));
    Maa_ab2(t) = max(0,binornd(Maa_ab2(t-1),(1-mu_m_aa)) + binornd(Z2,max(0,Mab2(t-1)/(MAB2(t-1)+Mab2(t-1)))));
    
    % Migrants pop 1 -> pop 2:
    MAB1to2(t)    = poissrnd(MAB1(t)*migrate);
    Mab1to2(t)    = poissrnd(Mab1(t)*migrate);
    MAa_AB1to2(t) = poissrnd(MAa_AB1(t)*migrate);
    MAa_ab1to2(t) = poissrnd(MAa_ab1(t)*migrate);
    Maa_AB1to2(t) = poissrnd(Maa_AB1(t)*migrate);
    Maa_ab1to2(t) = poissrnd(Maa_ab1(t)*migrate);

    % Migrants pop 2 -> pop 1:
    MAB2to1(t)    = poissrnd(MAB2(t)*migrate);
    Mab2to1(t)    = poissrnd(Mab2(t)*migrate);
    MAa_AB2to1(t) = poissrnd(MAa_AB2(t)*migrate);
    MAa_ab2to1(t) = poissrnd(MAa_ab2(t)*migrate);
    Maa_AB2to1(t) = poissrnd(Maa_AB2(t)*migrate);
    Maa_ab2to1(t) = poissrnd(Maa_ab2(t)*migrate);

    % Population sizes following migration events:
    MAB1(t)    = MAB1(t) - MAB1to2(t) + MAB2to1(t);
    Mab1(t)    = Mab1(t) - Mab1to2(t) + Mab2to1(t);
    MAa_AB1(t) = MAa_AB1(t) - MAa_AB1to2(t) + MAa_AB2to1(t);
    MAa_ab1(t) = MAa_ab1(t) - MAa_ab1to2(t) + MAa_ab2to1(t);
    Maa_AB1(t) = Maa_AB1(t) - Maa_AB1to2(t) + Maa_AB2to1(t);
    Maa_ab1(t) = Maa_ab1(t) - Maa_ab1to2(t) + Maa_ab2to1(t);

    MAB2(t)    = MAB2(t) + MAB1to2(t) - MAB2to1(t);
    Mab2(t)    = Mab2(t) + Mab1to2(t) - Mab2to1(t);
    MAa_AB2(t) = MAa_AB2(t) + MAa_AB1to2(t) - MAa_AB2to1(t);
    MAa_ab2(t) = MAa_ab2(t) + MAa_ab1to2(t) - MAa_ab2to1(t);
    Maa_AB2(t) = Maa_AB2(t) + Maa_AB1to2(t) - Maa_AB2to1(t);
    Maa_ab2(t) = Maa_ab2(t) + Maa_ab1to2(t) - Maa_ab2to1(t);
end

% Derived quantities (for plotting):

t              = linspace(1,max_t,max_t);               % Days of simulation
Mm1            = Mab1 + MAB1;                           % Number of male adult mosquitoes over time
Mf1            = MAa_AB1 + MAa_ab1 + Maa_AB1 + Maa_ab1; % Number of female adult mosquitoes over time
M1             = Mm1 + Mf1;                             % Total number of adult mosquitoes over time
Mtransgenic1   = MAa_AB1 + MAa_ab1 + MAB1;              % Number of transgenic adult mosquitoes over time
Mf_transgenic1 = MAa_AB1 + MAa_ab1;                     % Number of female transgenic adult mosquitoes over time

Mm2            = Mab2 + MAB2;                           % Number of male adult mosquitoes over time
Mf2            = MAa_AB2 + MAa_ab2 + Maa_AB2 + Maa_ab2; % Number of female adult mosquitoes over time
M2             = Mm2 + Mf2;                             % Total number of adult mosquitoes over time
Mtransgenic2   = MAa_AB2 + MAa_ab2 + MAB2;              % Number of transgenic adult mosquitoes over time
Mf_transgenic2 = MAa_AB2 + MAa_ab2;                     % Number of female transgenic adult mosquitoes over time

% Plot of Medusa spread over time:

figure(3)
plot(t(51:2051)-51,M1(51:2051),'r', t(51:2051)-51,Mtransgenic1(51:2051),'g', t(51:2051)-51,M2(51:2051),'k', t(51:2051)-51,Mtransgenic2(51:2051),'b');
xlim([0 2000])
ylim([0 22500])
xlabel('Day')
ylabel('Population size')
title('Medusa spread, stochastic model, 2 pops, migration rate = 1%/gen')
legend('Total pop (pop A)','Transgenic (pop A)','Total pop (pop B)','Transgenic (pop B)')