function [SABOE,SABOA,SABOL,Larvae,Waiters,OPFM]=TrackDynamics(InitEggs,Larvae,Waiters,tmax,dominance,...
    MaxClutchsize,DensityDependence,ddparameter,hybrid_fitness,larval_weekly_survival,b,T,L,dev_dist,niche_shift,seasonality_option,plotnumber)

% Matrices with 6 rows for 6 genotypes: 3 pure, and 3 hybrids
% Em1 = row1; % FM/FM genotype; red in plots,
% Em2 = row2; % NM/NM genotype; blue in plots,
% Ef = row3; % SL/SL genotype; green in plots,
% row 4,5 and 6 are for Em1/Em2, Em1/Ef and Em2/Ef hybrid genotypes, respectively

nof_seasons=max(1,seasonality_option); % run for a minimum of 1 seasons 
%   When seasonality_option=0, the season always starts in the same moon phase
%   When seasonality_option>0, the season  starts with a random moon phase 

prop_day=[.10 .20 .40 .20 .10]; % Proportion of adult population emerging on each of 5 days

SABOL=nan(tmax,nof_seasons,6); % tracking larvae in each strain over time - abbreviation stands for strain abundance of larvae, right?
SABOE=nan(tmax,nof_seasons,6); % tracking eggs in each strain over time
SABOA=nan(tmax,nof_seasons,6); % tracking adult in each strain over time
OPFM=nan(tmax,nof_seasons); % tracking Overall Probability that Female Mated over time

if dev_dist == "dev_unif" % All results in the paper use this distribution
    prop_mature=[ 0 0 0 0 0 0.166667 .2 .25 .33343 .5 1]; %Probability of being done developing at each age
elseif dev_dist == "dev_norm"
    prop_mature=[ 0 0 0 0 0 0.0500    0.1579    0.3750    0.6000    0.7500    1.0000];
elseif dev_dist == "dev_rate"
    prop_mature=[ 0 0 0 0 0 .8 .8 .8 .8 .8 .8];
end
prop_mature=ones([6 1]).*prop_mature;

% Indexing: row = strain (6 of them)
% Eggs are therefore 6 x 1
% Larvae are 6 x 11, column denotes stage of development; could in principle be something else than 11 of them if prop_mature has fewer or more than 11 numbers
% Waiters are 6 x 4, 1st column contains Larvae that just entered the
% waiting-to-emerge stage, last column contains Larvae that have been waiting for 4 weeks.  
% Adults are 6 x 2, for two sexes 
% Some variables below are 6 x tmax, these collect summary data over time.

%If no initial conditions are specified for these stages, start them at
%zero 
if isempty(Larvae) Larvae=zeros([6 length(prop_mature)]); end;
if isempty(Waiters) Waiters=zeros([6 4]); end; % 6 strains, 4 waiting stages


%dd parameter determines the strength of density dependence (experienced
%equally by all genotypes). The first two entries of ddparameter are parameters alpha_n1 a
% alpha_n2 in the manuscript, they determine the strenght of density-dependence in the two 
% different niches. The last entry of ddparameter is only used in the absence of a niche shift, 
% and determines the strenght of density dependence for that single niche.  
if niche_shift==2
    ddparameter=[1e-8*1 1e-8*2.6 1e-8];
elseif niche_shift==4
    ddparameter=[1e-8*1 1e-7*1.96 1e-8];
elseif niche_shift==6
    ddparameter=[1e-8*1 1e-6 1e-8];
end

% Defining timing phenotypes each genotype in matrix, corresponding to rows (strains)
if dominance == "S"                  % First dominance scenario: synchrony is dominant
    S1_phenotype = [1 4 5];          % Strain 1 plus S1-containing hybrids express S1 phenotype (interpretable as emerging at full moon)
    S2_phenotype = [2 6];            % Strain 2 and the AS2 hybrid express S2 phenotype (interpretable as emerging at new moon)
    A_phenotype =  3;                % Since synchrony is dominant, only AA individuals express the A phenotype
elseif dominance == "A"              % Second dominance scenario: asynchrony is dominant
    S1_phenotype = [1 4];
    S2_phenotype = 2;
    A_phenotype =  [3 5 6];
end

Eggs=InitEggs';

for season=1:nof_seasons % this correctly indicates the number of seasons run, as 1 or more

    if seasonality_option==0
        timeshift=0;
    else % randomized start
        %timeshift=unidrnd(4)-1; % shift from 0 to 3 weeks (4 diffe[rent options)
        %Above commented-out code only works with certain matlab packages, below code works without the packages
        timeshift=randi([1 4],1,1)-1;
    end

    for t=(1:tmax)+timeshift

        % First, compute consequences of density dependence. It modifies the survival rates of larvae
        % (we compute survival of larvae in a way that combines weekly survival with its potential density dependent modifications)
        if DensityDependence=="A"        % Niche shift, resources are shared among other larvae under or over age of niche shift
            clutch_size=MaxClutchsize;
            competition(1)=sum(Eggs(:))+sum(sum(Larvae(:,1:niche_shift-1)));
            competition(2)=sum(sum(Larvae(:,niche_shift:end)))+sum(Waiters, [1 2]);
            larvalsurvivalnow(1)=larval_weekly_survival*exp(-(ddparameter(1))*competition(1).^2);
            larvalsurvivalnow(2)=larval_weekly_survival*exp(-(ddparameter(2))*competition(2).^2);
            larvalsurvivalnow(3)=larval_weekly_survival;
        elseif DensityDependence=="B" % Resources are shared among all larvae.
            clutch_size=MaxClutchsize;
            competition=sum(Eggs(:))+sum(Larvae(:),[1 2])+sum(Waiters,[1 2]);
            larvalsurvivalnow(1:3)=larval_weekly_survival*exp(-ddparameter(3)*competition.^2); %
        end
        %Waiters moving to the next waiting stage
        Waiters(:,2:end)=larvalsurvivalnow(2)*Waiters(:,1:end-1);
        % maturation of larvae into waiters
        Waiters(:,1)=larvalsurvivalnow(2)*sum((Larvae.*prop_mature)')'; %

        % larvae become older
        Larvae(:,niche_shift+1:end)=Larvae(:,niche_shift:end-1).*(larvalsurvivalnow(2).*(1-prop_mature(:,niche_shift:end-1)));
        Larvae(:,2:niche_shift)=Larvae(:,1:niche_shift-1).*(larvalsurvivalnow(1).*(1-prop_mature(:,1:niche_shift-1)));
        %Eggs become larvae
        Larvae(:,1)=Eggs.*larvalsurvivalnow(1);
        %Resetting eggs to zero (will be filled again by reproduction in the code below)
        Eggs=zeros([6 1]);

        % if the time is right, we also have waiting larvae change state to become adults.
        % This is instantaneous (weekly mortalities are taken into account elsewhere) so they're 
        % also immediately removed from the waiter state, which is in other words emptied right now.

        if mod(t,4) == 0  % FM and SL strain waiters emerge
            Adult(S1_phenotype,:,1) = Waiters(S1_phenotype,:)/2;
            Adult(S1_phenotype,:,2) = Waiters(S1_phenotype,:)/2; 
            Waiters(S1_phenotype,:,:) = 0;
            Adult(A_phenotype,:,1) = Waiters(A_phenotype,:)/2;
            Adult(A_phenotype,:,2) = Waiters(A_phenotype,:)/2; 
            Waiters(A_phenotype,:) = 0;
        elseif mod(t, 4) == 2  %  NM and SL strain waiters emerge.
            Adult(S2_phenotype,:,1) = Waiters(S2_phenotype,:)./2;
            Adult(S2_phenotype,:,2) = Waiters(S2_phenotype,:)/2; 
            Waiters(S2_phenotype,:,:) = 0;
            Adult(A_phenotype,:,1) = Waiters(A_phenotype,:)/2;
            Adult(A_phenotype,:,2) = Waiters(A_phenotype,:)/2; 
            Waiters(A_phenotype,:) = 0;
        else
            Adult=zeros([6 4 2]);
        end

        % all the possibilities to emerge have now occurred. If there are adult
        % females & males, then we look at reproduction (and afterwards make
        % sure there's no adults left)

        if sum(Adult(:,:,2),[1 2])>0 % there are females

            % Number of females laying eggs inheriting allele p,q and r
            number_female_p = Adult(1,:,2) + 0.5 * Adult(4,:,2) + 0.5 * Adult(5,:,2);
            number_female_q = Adult(2,:,2) + 0.5 * Adult(4,:,2) + 0.5 * Adult(6,:,2);
            number_female_r = Adult(3,:,2) + 0.5 * Adult(5,:,2) + 0.5 * Adult(6,:,2);

            if sum(Adult(:,:,1),[1 2]) > 0 % there are males

                % Probability of inhereting each allele from father.
                prob_male_p = (sum(Adult(1,:,1)) + 0.5 * sum(Adult(4,:,1)) + 0.5 * sum(Adult(5,:,1)) ) / sum(Adult(:,:,1),[1 2]);
                prob_male_q = (sum(Adult(2,:,1)) + 0.5 * sum(Adult(4,:,1)) + 0.5 * sum(Adult(6,:,1)) ) / sum(Adult(:,:,1),[1 2]);
                prob_male_r = (sum(Adult(3,:,1)) + 0.5 * sum(Adult(5,:,1)) + 0.5 * sum(Adult(6,:,1)) ) / sum(Adult(:,:,1),[1 2]);
        %Here we calculate the mating dynamics if strains only mate with
        %themselves, so as if they're different species). 
                if isnan(hybrid_fitness) % this is used to denote that the strains are kept separate (as if they were different species)

                    % note: here overall_prob_female_mated has to have 6 values, not just
                    % 1, because mating probabilities are strain-specific.
                    % Mostly, this section is intended to be used with just
                    % 'pure' strains present  - note that in this fully assortative
                    % setting, a population that starts with pure strains will
                    % not produce hybrids

                    if isnan(b) % we assume that everyone mates
                        overall_prob_female_mated=ones([6 1]);
                    else % b is defined, therefore mate-finding needs to be modelled properly
                        for d=1:length(prop_day)
                            % per female probability of becoming mated, conditional
                            % on emerging on the dth day of this emergence
                            % week, calculating for each of the 6 genotypes
                            % separately 
                            for i=1:6
                                matedp(i,d)=matingprob(sum(Adult(i,:,2),2)*prop_day(d),sum(Adult(i,:,1),2)*prop_day(d),b,T,L,1);
                            end
                        end
                        matedp(isnan(matedp))=0; % if no-one was there to mate, we define their mating success as 0 rather than nan
                        overall_prob_female_mated=matedp*prop_day'; % this order of matrix multiplication means we'll get a 6x1 matrix out
                    end
                    % if we want OPFM out as a single number (for each time step, week not day), we need to weigh by adult female numbers.
                    OPFM(t)=sum(Adult(:,:,2).*overall_prob_female_mated,[1 2], 'omitnan')/(eps+sum(Adult(:,:,2),[1 2], 'omitnan'));

                    % now actually create the offspring. Here we assume 'strain
                    % begets [same] strain', even in hybrids, but note that
                    % this will really only be run without hybrids
                    Eggs(1) = sum(Adult(1,:,2).* clutch_size * overall_prob_female_mated(1));
                    Eggs(2) = sum(Adult(2,:,2).* clutch_size * overall_prob_female_mated(2));
                    Eggs(3) = sum(Adult(3,:,2).* clutch_size * overall_prob_female_mated(3));
                    Eggs(4) = sum(Adult(4,:,2).* clutch_size * overall_prob_female_mated(4));
                    Eggs(5) = sum(Adult(5,:,2).* clutch_size * overall_prob_female_mated(5));
                    Eggs(6) = sum(Adult(6,:,2).* clutch_size * overall_prob_female_mated(6));

                else % random mating
                    if isnan(b) % we assume that everyone mates
                        overall_prob_female_mated=1;
                    else % b is defined, therefore mate-finding needs to be modelled properly
                        for d=1:length(prop_day)
                            % per female probability of becoming mated, conditional
                            % on emerging on the dth day of this emergence
                            % week, the same for all genotypes \
                            matedp(d)=matingprob(sum(Adult(:,:,2),[1 2])*prop_day(d),sum(Adult(:,:,1),[1 2])*prop_day(d),b,T,L,1);
                        end
                        overall_prob_female_mated=sum(prop_day.*matedp);
                    end
                    OPFM(t-timeshift,season)=overall_prob_female_mated;

                    % now create the offspring
                    Eggs(1) = sum(number_female_p.* clutch_size * overall_prob_female_mated * prob_male_p);
                    Eggs(2) = sum(number_female_q.* clutch_size * overall_prob_female_mated * prob_male_q);
                    Eggs(3) = sum(number_female_r.* clutch_size * overall_prob_female_mated * prob_male_r);
                    Eggs(4) = sum( ( number_female_p * prob_male_q...
                        + number_female_q * prob_male_p ) .* clutch_size * overall_prob_female_mated * hybrid_fitness);

                    Eggs(5) = sum( ( number_female_p * prob_male_r...
                        + number_female_r * prob_male_p ) .* clutch_size * overall_prob_female_mated * hybrid_fitness);

                    Eggs(6) = Eggs(6)+ sum( ( number_female_q * prob_male_r...
                        + number_female_r * prob_male_q ) .* clutch_size * overall_prob_female_mated * hybrid_fitness);
                end
            end
        end
 
        % collect data
        RelevantLarvae=sum(Larvae,2); % Sum over the different developmental stages
        RelevantWaiters=sum(Waiters,2); % Sum over the 4 different waiting lengths
        for k=1:6
            SABOE(t-timeshift,season,k)=Eggs(k);
            SABOL(t-timeshift,season,k)=RelevantLarvae(k);
            %SABOA(t-timeshift,season,k)=Adult(k,2);
            SABOW(t-timeshift,season,k)=RelevantWaiters(k);
        end
    end
    % prepare the next generation
    for k=1:6
       Eggs(k)=sum(SABOE(end-3:end,season,k));
    end
    % same constant number of eggs as was there in the very first season
    Eggs=sum(InitEggs)*Eggs/sum(Eggs);
end

if plotnumber>0 % plotting will happen
    figure(plotnumber);
    EggDynamics1=SABOE(:,:,1); EggDynamics1=EggDynamics1(:);
    EggDynamics2=SABOE(:,:,2); EggDynamics2=EggDynamics2(:);
    EggDynamics3=SABOE(:,:,3); EggDynamics3=EggDynamics3(:);
    EggDynamics4=SABOE(:,:,4); EggDynamics4=EggDynamics4(:);
    EggDynamics5=SABOE(:,:,5); EggDynamics5=EggDynamics5(:);
    EggDynamics6=SABOE(:,:,6); EggDynamics6=EggDynamics6(:);
    tplot=1:length(EggDynamics1);
    plot(tplot,EggDynamics1,'ro',tplot,EggDynamics2,'bx',tplot,EggDynamics3,'gs',tplot,EggDynamics4,'r:',tplot,EggDynamics5,'k:',tplot,EggDynamics6,'g:'); ylabel('Eggs')
    drawnow
end