endtime = 10;

% Data
dataayeast = zeros(5, 5, 1);
dataalphayeast = zeros(5, 5, 1);
datadiploidyeast = zeros(5, 5, 1);
dataashmoo = zeros(5, 5, 1);
dataalphashmoo = zeros(5, 5, 1);
dataafactor = zeros(5, 5, 1);
dataalphafactor = zeros(5, 5, 1);
datatime = [];

% Initial location of yeast cell
ayeast = zeros(5, 5);
alphayeast = zeros(5, 5);
diploidyeast = zeros(5, 5);
ashmoo = zeros(5, 5, 4);
alphashmoo = zeros(5, 5, 4);

ayeast(1, 5) = 1;
alphayeast(5, 1) = 1;

% Initial pheromone concentration
afactor = zeros(5, 5);
alphafactor = zeros(5, 5);

% Initial time
time = 0;
run = 1;

% Initial Data
datatime = [datatime time];
dataayeast(:, :, 1) = ayeast;
dataalphayeast(:, :, 1) = alphayeast;
datadiploidyeast(:, :, 1) = diploidyeast;
dataashmoo(:, :, 1) = sum(ashmoo, 3);
dataalphashmoo(:, :, 1) = sum(alphashmoo, 3);
dataafactor(:, :, 1) = afactor;
dataalphafactor(:, :, 1) = alphafactor;

% Random drift of yeast in solution

% Drift rate
drift = 0.1;

% Drift distance
displacement = 1;

% Mitosis rate
mitosis = 0.1;

% Death rate
death = 0.001;

% Pheromone production rate
production = 5;

% Sex switch rate
change = 0.5;

% Pheromone

% Diffusion Rate
diffusion = 10;

% Decay rate
decay = 0.1;

% Shmoo
shmoo = 10;
unshmoo = 0.01;

% Diploid

% Meiosis Rate
meiosis = 0.01;



% Pheromone

% Diffusion distance
diffusiondisplacement = 1;

while time < endtime

    % Total rate
    populationa = sum(sum(ayeast));
    populationalpha = sum(sum(alphayeast));
    populationshmooa = sum(sum(sum(ashmoo)));
    populationshmooalpha = sum(sum(sum(alphashmoo)));
    populationshmoo = populationshmooa + populationshmooalpha;
    particlesa = sum(sum(afactor));
    particlesalpha = sum(sum(alphafactor));
    populationdiploid = sum(sum(diploidyeast));

    haploidrate = (populationa + populationalpha)*(drift + mitosis + death + change + production);
    shmoorate = (populationshmooa + populationshmooalpha)*(drift + death + production + unshmoo + shmoo);
    particlerate = (particlesa + particlesalpha)*(diffusion + decay);
    diploidrate = populationdiploid*(drift + mitosis + death + meiosis);

    totalrate = haploidrate + shmoorate + particlerate + diploidrate;




    % Action
    random = 1 - rand(1);
    if random <= populationa*drift/totalrate
        pick = Pick(ayeast);
        % Brownian motion
        newcell = Direction(pick, displacement);

        if ayeast(pick(1, 1), pick(1, 2)) >= 0
            ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) - 1;
        else
            ayeast(pick(1, 1), pick(1, 2)) = 0;
        end
        if alphafactor(newcell(1, 1), newcell(1, 2)) == 0
            ayeast(newcell(1, 1), newcell(1, 2)) = ayeast(newcell(1, 1), newcell(1, 2)) + 1;
        else
            ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
            alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
        end

    elseif random <= (populationa + populationalpha)*drift/totalrate
        pick = Pick(alphayeast);
        newcell = Direction(pick, displacement);
        if alphayeast(pick(1, 1), pick(1, 2)) >= 0
            alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) - 1;
        else
            alphayeast(pick(1, 1), pick(1, 2)) = 0;
        end
        if afactor(newcell(1, 1), newcell(1, 2)) == 0
            alphayeast(newcell(1, 1), newcell(1, 2)) = alphayeast(newcell(1, 1), newcell(1, 2)) + 1;
        else
            alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
            afactor(newcell(1, 1), newcell(1, 2)) = 0;
        end


    elseif random <= (populationa*(drift + mitosis) + populationalpha*drift)/totalrate
        pick = Pick(ayeast);
        % Mitosis

        ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) + 1;
    elseif random <= (populationa + populationalpha)*(drift + mitosis)/totalrate
        pick = Pick(alphayeast);
        alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) + 1;

    elseif random <= (populationa*(drift + mitosis + death) + populationalpha*(drift + mitosis))/totalrate
        pick = Pick(ayeast);
        % Cell death

        ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) - 1;
    elseif random <= (populationa + populationalpha)*(drift + mitosis + death)/totalrate
        pick = Pick(alphayeast);
        alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) - 1;


    elseif random <= (populationa*(drift + mitosis + death + change) + populationalpha*(drift + mitosis + death))/totalrate
        pick = Pick(ayeast);
        % Sex switching

        ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) - 1;
        alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) + 1;
    elseif random <= (populationa + populationalpha)*(drift + mitosis + death + change)/totalrate
        pick = Pick(alphayeast);
        ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) + 1;
        alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) - 1;


    elseif random <= (haploidrate - populationalpha*production)/totalrate
        pick = Pick(ayeast);
        % Pheromone production

        afactor(pick(1, 1), pick(1, 2)) = afactor(pick(1, 1), pick(1, 2)) + 1;
    elseif random <= haploidrate/totalrate
        pick = Pick(alphayeast);
        alphafactor(pick(1, 1), pick(1, 2)) = alphafactor(pick(1, 1), pick(1, 2)) + 1;

    elseif  random <= (haploidrate + particlesa*diffusion)/totalrate;
        % Pheromone action
        % Diffusion
        particlepick = Pick(afactor);
        newparticle = Direction(particlepick, diffusiondisplacement);
        newshmoo = newparticle(1, 3) - (-1)^newparticle(1, 3);

        %             r1 = afactor(particlepick(1, 1), particlepick(1, 2))
        if afactor(particlepick(1, 1), particlepick(1, 2)) >= 0
            afactor(particlepick(1, 1), particlepick(1, 2)) = afactor(particlepick(1, 1), particlepick(1, 2)) - 1;
        else
            afactor(particlepick(1, 1), particlepick(1, 2)) = 0;
        end
        %             r2 = afactor(particlepick(1, 1), particlepick(1, 2))
        if sum(alphashmoo(newparticle(1, 1), newparticle(1, 2), :)) + alphayeast(newparticle(1, 1), newparticle(1, 2)) == 0
            afactor(newparticle(1, 1), newparticle(1, 2)) = afactor(newparticle(1, 1), newparticle(1, 2)) + 1;
        else
            drandom  = 1 - rand(1);
            alphashmool = alphashmoo(newparticle(1, 1), newparticle(1, 2), 1);
            alphashmoor = alphashmoo(newparticle(1, 1), newparticle(1, 2), 2);
            alphashmood = alphashmoo(newparticle(1, 1), newparticle(1, 2), 3);
            alphashmoou = alphashmoo(newparticle(1, 1), newparticle(1, 2), 4);
            nparticles = alphashmool + alphashmoor + alphashmood + alphashmoou + alphayeast(newparticle(1, 1), newparticle(1, 2));
            if drandom <= alphashmool/nparticles;
                alphashmoo(newparticle(1, 1), newparticle(1, 2), 1) = alphashmoo(newparticle(1, 1), newparticle(1, 2), 1) - 1;
                alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) = alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) + 1;
            elseif drandom <= (alphashmool + alphashmoor)/nparticles
                alphashmoo(newparticle(1, 1), newparticle(1, 2), 2) = alphashmoo(newparticle(1, 1), newparticle(1, 2), 2) - 1;
                alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) = alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) + 1;
            elseif drandom <= (alphashmool + alphashmoor + alphashmood)/nparticles
                alphashmoo(newparticle(1, 1), newparticle(1, 2), 3) = alphashmoo(newparticle(1, 1), newparticle(1, 2), 3) - 1;
                alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) = alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) + 1;
            elseif drandom <= (alphashmool + alphashmoor + alphashmood + alphashmoou)/nparticles
                alphashmoo(newparticle(1, 1), newparticle(1, 2), 4) = alphashmoo(newparticle(1, 1), newparticle(1, 2), 4) - 1;
                alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) = alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) + 1;
            else
                alphayeast(newparticle(1, 1), newparticle(1, 2)) = alphayeast(newparticle(1, 1), newparticle(1, 2)) - 1;
                alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) = alphashmoo(newparticle(1, 1), newparticle(1, 2), newshmoo) + 1;
            end
        end

    elseif random <= (haploidrate + (particlesa +  particlesalpha)*diffusion)/totalrate;
        particlepick = Pick(alphafactor);
                newparticle = Direction(particlepick, diffusiondisplacement);
        newshmoo = newparticle(1, 3) - (-1)^newparticle(1, 3);
        if alphafactor(particlepick(1, 1), particlepick(1, 2)) >= 0
            alphafactor(particlepick(1, 1), particlepick(1, 2)) = alphafactor(particlepick(1, 1), particlepick(1, 2)) - 1;
        else
            alphafactor(particlepick(1, 1), particlepick(1, 2)) = 0;
        end
        if sum(ashmoo(newparticle(1, 1), newparticle(1, 2), :)) + ayeast(newparticle(1, 1), newparticle(1, 2)) == 0
            alphafactor(newparticle(1, 1), newparticle(1, 2)) = alphafactor(newparticle(1, 1), newparticle(1, 2)) + 1;
        else
            drandom  = 1 - rand(1);
            ashmool = ashmoo(newparticle(1, 1), newparticle(1, 2), 1);
            ashmoor = ashmoo(newparticle(1, 1), newparticle(1, 2), 2);
            ashmood = ashmoo(newparticle(1, 1), newparticle(1, 2), 3);
            ashmoou = ashmoo(newparticle(1, 1), newparticle(1, 2), 4);
            nparticles = ashmool + ashmoor + ashmood + ashmoou + ayeast(newparticle(1, 1), newparticle(1, 2));
            if drandom <= ashmool/nparticles;
                ashmoo(newparticle(1, 1), newparticle(1, 2), 1) = ashmoo(newparticle(1, 1), newparticle(1, 2), 1) - 1;
                ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) = ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) + 1;
            elseif drandom <= (ashmool + ashmoor)/nparticles
                ashmoo(newparticle(1, 1), newparticle(1, 2), 2) = ashmoo(newparticle(1, 1), newparticle(1, 2), 2) - 1;
                ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) = ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) + 1;
            elseif drandom <= (ashmool + ashmoor + ashmood)/nparticles
                ashmoo(newparticle(1, 1), newparticle(1, 2), 3) = ashmoo(newparticle(1, 1), newparticle(1, 2), 3) - 1;
                ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) = ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) + 1;
            elseif drandom <= (ashmool + ashmoor + ashmood + ashmoou)/nparticles
                ashmoo(newparticle(1, 1), newparticle(1, 2), 4) = ashmoo(newparticle(1, 1), newparticle(1, 2), 4) - 1;
                ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) = ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) + 1;
            else
                ayeast(newparticle(1, 1), newparticle(1, 2)) = ayeast(newparticle(1, 1), newparticle(1, 2)) - 1;
                ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) = ashmoo(newparticle(1, 1), newparticle(1, 2), newparticle(1, 3)) + 1;
            end
        end



    elseif random <= (haploidrate + particlerate - particlesalpha*decay)/totalrate;
        % Decay
        particlepick = Pick(afactor);

        if afactor(particlepick(1, 1), particlepick(1, 2)) >= 0
            afactor(particlepick(1, 1), particlepick(1, 2)) = afactor(particlepick(1, 1), particlepick(1, 2)) - 1;
        else
            afactor(particlepick(1, 1), particlepick(1, 2)) = 0;
        end
        %                 r3 = afactor(particlepick(1, 1), particlepick(1, 2))
        %                 afactor(particlepick(1, 1), particlepick(1, 2)) = afactor(particlepick(1, 1), particlepick(1, 2)) - 1
        %                 r4 = afactor(particlepick(1, 1), particlepick(1, 2))
    elseif random <= (haploidrate + particlerate)/totalrate
        particlepick = Pick(alphafactor);
        if alphafactor(particlepick(1, 1), particlepick(1, 2)) >= 0
            alphafactor(particlepick(1, 1), particlepick(1, 2)) = alphafactor(particlepick(1, 1), particlepick(1, 2)) - 1;
        else
            alphafactor(particlepick(1, 1), particlepick(1, 2)) = 0;
        end

        % Shmoo
    elseif random <= (haploidrate + particlerate + populationshmoo*drift)/totalrate
        % Drift
        pick = PickShmoo(ashmoo, alphashmoo);
        % Brownian motion
        newcell = Direction(pick, displacement);
        if pick(1, 3) == 1
            ashmoo(pick(1, 1), pick(1, 2), 1) = ashmoo(pick(1, 1), pick(1, 2), 1) - 1;
            if alphafactor(newcell(1, 1), newcell(1, 2)) == 0
                ashmoo(newcell(1, 1), newcell(1, 2), 1) = ashmoo(newcell(1, 1), newcell(1, 2), 1) + 1;
            else
                ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
                alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
            end

        elseif pick(1, 3) == 2
            ashmoo(pick(1, 1), pick(1, 2), 2) = ashmoo(pick(1, 1), pick(1, 2), 2) - 1;
            if alphafactor(newcell(1, 1), newcell(1, 2)) == 0
                ashmoo(newcell(1, 1), newcell(1, 2), 2) = ashmoo(newcell(1, 1), newcell(1, 2), 2) + 1;
            else
                ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
                alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
            end
        elseif pick(1, 3) == 3
            ashmoo(pick(1, 1), pick(1, 2), 3) = ashmoo(pick(1, 3), pick(1, 2), 1) - 1;
            if alphafactor(newcell(1, 1), newcell(1, 2)) == 0
                ashmoo(newcell(1, 1), newcell(1, 2), 3) = ashmoo(newcell(1, 1), newcell(1, 2), 3) + 1;
            else
                ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
                alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
            end

        elseif pick(1, 3) == 4
            ashmoo(pick(1, 1), pick(1, 2), 4) = ashmoo(pick(1, 1), pick(1, 2), 4) - 1;
            if alphafactor(newcell(1, 1), newcell(1, 2)) == 0
                ashmoo(newcell(1, 1), newcell(1, 2), 4) = ashmoo(newcell(1, 1), newcell(1, 2), 4) + 1;
            else
                ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = ashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
                alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
            end

        elseif pick(1, 3) == 5
            alphashmoo(pick(1, 1), pick(1, 2), 1) = alphashmoo(pick(1, 1), pick(1, 2), 1) - 1;
            if afactor(newcell(1, 1), newcell(1, 2)) == 0
                alphashmoo(newcell(1, 1), newcell(1, 2), 1) = alphashmoo(newcell(1, 1), newcell(1, 2), 1) + 1;
            else
                alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
                afactor(newcell(1, 1), newcell(1, 2)) = 0;
            end
        elseif pick(1, 3) == 6
            alphashmoo(pick(1, 1), pick(1, 2), 2) = alphashmoo(pick(1, 1), pick(1, 2), 2) - 1;
            if afactor(newcell(1, 1), newcell(1, 2)) == 0
                alphashmoo(newcell(1, 1), newcell(1, 2), 2) = alphashmoo(newcell(1, 1), newcell(1, 2), 2) + 1;
            else
                alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
                afactor(newcell(1, 1), newcell(1, 2)) = 0;
            end
        elseif pick(1, 3) == 7
            alphashmoo(pick(1, 1), pick(1, 2), 3) = alphashmoo(pick(1, 1), pick(1, 2), 3) - 1;
            if afactor(newcell(1, 1), newcell(1, 2)) == 0
                alphashmoo(newcell(1, 1), newcell(1, 2), 3) = alphashmoo(newcell(1, 1), newcell(1, 2), 3) + 1;
            else
                alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
                afactor(newcell(1, 1), newcell(1, 2)) = 0;
            end
        elseif pick(1, 3) == 8
            alphashmoo(pick(1, 1), pick(1, 2), 4) = alphashmoo(pick(1, 1), pick(1, 2), 4) - 1;
            if afactor(newcell(1, 1), newcell(1, 2)) == 0
                alphashmoo(newcell(1, 1), newcell(1, 2), 4) = alphashmoo(newcell(1, 1), newcell(1, 2), 4) + 1;
            else
                alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) = alphashmoo(newcell(1, 1), newcell(1, 2), newcell(1, 3)) + 1;
                afactor(newcell(1, 1), newcell(1, 2)) = 0;
            end
        end

    elseif random > (haploidrate + particlerate + populationshmoo*drift)/totalrate && random <= (haploidrate + particlerate + populationshmoo*(drift + death))/totalrate
        % Death
        pick = PickShmoo(ashmoo, alphashmoo);

        if pick(1, 3) == 1
            ashmoo(pick(1, 1), pick(1, 2), 1) = ashmoo(pick(1, 1), pick(1, 2), 1) - 1;
        elseif pick(1, 3) == 2
            ashmoo(pick(1, 1), pick(1, 2), 2) = ashmoo(pick(1, 1), pick(1, 2), 2) - 1;
        elseif pick(1, 3) == 3
            ashmoo(pick(1, 1), pick(1, 2), 3) = ashmoo(pick(1, 3), pick(1, 2), 1) - 1;
        elseif pick(1, 3) == 4
            ashmoo(pick(1, 1), pick(1, 2), 4) = ashmoo(pick(1, 1), pick(1, 2), 4) - 1;
        elseif pick(1, 3) == 5
            alphashmoo(pick(1, 1), pick(1, 2), 1) = alphashmoo(pick(1, 1), pick(1, 2), 1) - 1;
        elseif pick(1, 3) == 6
            alphashmoo(pick(1, 1), pick(1, 2), 2) = alphashmoo(pick(1, 1), pick(1, 2), 2) - 1;
        elseif pick(1, 3) == 7
            alphashmoo(pick(1, 1), pick(1, 2), 3) = alphashmoo(pick(1, 1), pick(1, 2), 3) - 1;
        elseif pick(1, 3) == 8
            alphashmoo(pick(1, 1), pick(1, 2), 4) = alphashmoo(pick(1, 1), pick(1, 2), 4) - 1;

        elseif random  > (haploidrate + particlerate + populationshmoo*(drift + death))/totalrate && random <= (haploidrate + particlerate + populationshmoo*(drift + death + production))/totalrate
            % Production
            pick = PickShmoo(ashmoo, alphashmoo);

            if pick(1, 3) <= 4
                afactor(pick(1, 1), pick(1, 2)) = afactor(pick(1, 1), pick(1, 2)) + 1;
            else
                alphafactor(pick(1, 1), pick(1, 2)) = alphafactor(pick(1, 1), pick(1, 2)) + 1;
            end
        elseif random > (haploidrate + particlerate + populationshmoo*(drift + death + production))/totalrate && random <= (haploidrate + particlerate + populationshmoo*(drift + death + production + unshmoo))/totalrate
            % Unshmoo
            pick = PickShmoo(ashmoo, alphashmoo);

            if pick(1, 3) == 1
                ashmoo(pick(1, 1), pick(1, 2), 1) = ashmoo(pick(1, 1), pick(1, 2), 1) - 1;
                ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) + 1;
            elseif pick(1, 3) == 2
                ashmoo(pick(1, 1), pick(1, 2), 2) = ashmoo(pick(1, 1), pick(1, 2), 2) - 1;
                ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) + 1;
            elseif pick(1, 3) == 3
                ashmoo(pick(1, 1), pick(1, 2), 3) = ashmoo(pick(1, 3), pick(1, 2), 1) - 1;
                ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) + 1;
            elseif pick(1, 3) == 4
                ashmoo(pick(1, 1), pick(1, 2), 4) = ashmoo(pick(1, 1), pick(1, 2), 4) - 1;
                ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) + 1;
            elseif pick(1, 3) == 5
                alphashmoo(pick(1, 1), pick(1, 2), 1) = alphashmoo(pick(1, 1), pick(1, 2), 1) - 1;
                alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) + 1;
            elseif pick(1, 3) == 6
                alphashmoo(pick(1, 1), pick(1, 2), 2) = alphashmoo(pick(1, 1), pick(1, 2), 2) - 1;
                alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) + 1;
            elseif pick(1, 3) == 7
                alphashmoo(pick(1, 1), pick(1, 2), 3) = alphashmoo(pick(1, 1), pick(1, 2), 3) - 1;
                alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) + 1;
            elseif pick(1, 3) == 8
                alphashmoo(pick(1, 1), pick(1, 2), 4) = alphashmoo(pick(1, 1), pick(1, 2), 4) - 1;
                alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) + 1;

            elseif random > (haploidrate + particlerate + populationshmoo*(drift + death + production + unshmoo))/totalrate && random <= (haploidrate + particlerate + shmoorate)/totalrate
                % Shmoo
                pick = PickShmoo(ashmoo, alphashmoo);
                newcell = ShmooMove(pick);
                if newcell(1, 3) == 1
                    newcell(1, 1) = newcell(1, 1) - displacement;
                    ashmoo(pick(1, 1), pick(1, 2), 1) = ashmoo(pick(1, 1), pick(1, 2), 1) - 1;
                    alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
                    if sum(alphashmoo(newcell(1, 1), newcell(1, 2), :)) == 0
                        ashmoo(newcell(1, 1), newcell(1, 2), 1) = ashmoo(newcell(1, 1), newcell(1, 2), 1) + 1;
                    else
                        shmate = PickMoreShmoo(alphashmoo(newcell(1, 1), newcell(1, 2), :));
                        alphashmoo(newcell(1, 1), newcell(1, 2), shmate) = alphashmoo(newcell(1, 1), newcell(1, 2), shmate) - 1;
                        diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
                    end


                elseif newcell(1, 3) == 2
                    newcell(1, 1) = newcell(1, 1) + displacement;
                    ashmoo(pick(1, 1), pick(1, 2), 2) = ashmoo(pick(1, 1), pick(1, 2), 2) - 1;
                    alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
                    if sum(alphashmoo(newcell(1, 1), newcell(1, 2), :)) == 0
                        ashmoo(newcell(1, 1), newcell(1, 2), 2) = ashmoo(newcell(1, 1), newcell(1, 2), 2) + 1;
                    else
                        shmate = PickMoreShmoo(alphashmoo(newcell(1, 1), newcell(1, 2), :));
                        alphashmoo(newcell(1, 1), newcell(1, 2), shmate) = alphashmoo(newcell(1, 1), newcell(1, 2), shmate) - 1;
                        diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
                    end

                elseif newcell(1, 3) == 3
                    newcell(1, 2) = newcell(1, 2) - displacement;
                    ashmoo(pick(1, 1), pick(1, 2), 3) = ashmoo(pick(1, 1), pick(1, 2), 3) - 1;
                    alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
                    if sum(alphashmoo(newcell(1, 1), newcell(1, 2), :)) == 0
                        ashmoo(newcell(1, 1), newcell(1, 2), 3) = ashmoo(newcell(1, 1), newcell(1, 2), 3) + 1;
                    else
                        shmate = PickMoreShmoo(alphashmoo(newcell(1, 1), newcell(1, 2), :));
                        alphashmoo(newcell(1, 1), newcell(1, 2), shmate) = alphashmoo(newcell(1, 1), newcell(1, 2), shmate) - 1;
                        diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
                    end

                elseif newcell(1, 3) == 4
                    newcell(1, 2) = newcell(1, 2) + displacement;
                    ashmoo(pick(1, 1), pick(1, 2), 4) = ashmoo(pick(1, 1), pick(1, 2), 4) - 1;
                    alphafactor(newcell(1, 1), newcell(1, 2)) = 0;
                    if sum(alphashmoo(newcell(1, 1), newcell(1, 2), :)) == 0
                        ashmoo(newcell(1, 1), newcell(1, 2), 4) = ashmoo(newcell(1, 1), newcell(1, 2), 4) + 1;
                    else
                        shmate = PickMoreShmoo(alphashmoo(newcell(1, 1), newcell(1, 2), :));
                        alphashmoo(newcell(1, 1), newcell(1, 2), shmate) = alphashmoo(newcell(1, 1), newcell(1, 2), shmate) - 1;
                        diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
                    end

                elseif newcell(1, 3) == 5
                    newcell(1, 1) = newcell(1, 1) - displacement;
                    alphashmoo(pick(1, 1), pick(1, 2), 1) = alphashmoo(pick(1, 1), pick(1, 2), 1) - 1;
                    afactor(newcell(1, 1), newcell(1, 2)) = 0;
                    if sum(ashmoo(newcell(1, 1), newcell(1, 2), :)) == 0
                        alphashmoo(newcell(1, 1), newcell(1, 2), 1) = alphashmoo(newcell(1, 1), newcell(1, 2), 1) + 1;
                    else
                        shmate = PickMoreShmoo(ashmoo(newcell(1, 1), newcell(1, 2), :));
                        ashmoo(newcell(1, 1), newcell(1, 2), shmate) = ashmoo(newcell(1, 1), newcell(1, 2), shmate) - 1;
                        diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
                    end

                elseif newcell(1, 3) == 6
                    newcell(1, 1) = newcell(1, 1) + displacement;
                    alphashmoo(pick(1, 1), pick(1, 2), 2) = alphashmoo(pick(1, 1), pick(1, 2), 2) - 1;
                    afactor(newcell(1, 1), newcell(1, 2)) = 0;
                    if sum(ashmoo(newcell(1, 1), newcell(1, 2), :)) == 0
                        alphashmoo(newcell(1, 1), newcell(1, 2), 2) = alphashmoo(newcell(1, 1), newcell(1, 2), 2) + 1;
                    else
                        shmate = PickMoreShmoo(ashmoo(newcell(1, 1), newcell(1, 2), :));
                        ashmoo(newcell(1, 1), newcell(1, 2), shmate) = ashmoo(newcell(1, 1), newcell(1, 2), shmate) - 1;
                        diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
                    end

                elseif newcell(1, 3) == 7
                    newcell(1, 2) = newcell(1, 2) - displacement;
                    alphashmoo(pick(1, 1), pick(1, 2), 3) = alphashmoo(pick(1, 1), pick(1, 2), 3) - 1;
                    afactor(newcell(1, 1), newcell(1, 2)) = 0;
                    if sum(ashmoo(newcell(1, 1), newcell(1, 2), :)) == 0
                        alphashmoo(newcell(1, 1), newcell(1, 2), 3) = alphashmoo(newcell(1, 1), newcell(1, 2), 3) + 1;
                    else
                        shmate = PickMoreShmoo(ashmoo(newcell(1, 1), newcell(1, 2), :));
                        ashmoo(newcell(1, 1), newcell(1, 2), shmate) = ashmoo(newcell(1, 1), newcell(1, 2), shmate) - 1;
                        diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
                    end

                elseif newcell(1, 3) == 8
                    newcell(1, 2) = newcell(1, 2) + displacement;
                    alphashmoo(pick(1, 1), pick(1, 2), 4) = alphashmoo(pick(1, 1), pick(1, 2), 4) - 1;
                    afactor(newcell(1, 1), newcell(1, 2)) = 0;
                    if sum(ashmoo(newcell(1, 1), newcell(1, 2), :)) == 0
                        alphashmoo(newcell(1, 1), newcell(1, 2), 4) = alphashmoo(newcell(1, 1), newcell(1, 2), 4) + 1;
                    else
                        shmate = PickMoreShmoo(ashmoo(newcell(1, 1), newcell(1, 2), :));
                        ashmoo(newcell(1, 1), newcell(1, 2), shmate) = ashmoo(newcell(1, 1), newcell(1, 2), shmate) - 1;
                        diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
                    end
                elseif newcell(1, 3) == 0
                    ashmoo(pick(1, 1), pick(1, 2), pick(1, 3)) = ashmoo(pick(1, 1), pick(1, 2), pick(1, 3)) - 1;
                    ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) + 1;
                elseif newcell(1, 3) == 9
                    alphashmoo(pick(1, 1), pick(1, 2), (pick(1, 3) - 4)) = alphashmoo(pick(1, 1), pick(1, 2), (pick(1, 3) - 4)) - 1;
                    alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) + 1;
                end
            end


        elseif random > (haploidrate + particlerate + shmoorate)/totalrate && random <= (haploidrate + particlerate + shmoorate + populationdiploid*drift)/totalrate
            pick = PickDiploid(diploidyeast);
            % Brownian motion
            newcell = Direction(pick, displacement);
            if pick(1, 3) == 0
                diploidyeast(pick(1, 1), pick(1, 2)) = diploidyeast(pick(1, 1), pick(1, 2)) - 1;
                diploidyeast(newcell(1, 1), newcell(1, 2)) = diploidyeast(newcell(1, 1), newcell(1, 2)) + 1;
            end
        elseif random > (haploidrate + particlerate + shmoorate + populationdiploid*drift)/totalrate && random <= (haploidrate + particlerate + shmoorate + populationdiploid*(drift + mitosis))/totalrate
            pick = PickDiploid(diploidyeast);
            % Mitosis
            diploidyeast(pick(1, 1), pick(1, 2)) = diploidyeast(pick(1, 1), pick(1, 2)) + 1;
        elseif random > (haploidrate + particlerate + shmoorate + populationdiploid*(drift + mitosis))/totalrate && random <= (haploidrate + particlerate + shmoorate + populationdiploid*(drift + mitosis + death))/totalrate
            pick = PickDiploid(diploidyeast);
            % Mitosis
            diploidyeast(pick(1, 1), pick(1, 2)) = diploidyeast(pick(1, 1), pick(1, 2)) - 1;
        elseif random > (haploidrate + particlerate + shmoorate + populationdiploid*(drift + mitosis + death))/totalrate && random <= random < (haploidrate + particlerate + shmoorate + diploidrate)/totalrate
            pick = PickDiploid(diploidyeast);
            % Meiosis
            diploidyeast(pick(1, 1), pick(1, 2)) = diploidyeast(pick(1, 1), pick(1, 2)) - 1;
            ayeast(pick(1, 1), pick(1, 2)) = ayeast(pick(1, 1), pick(1, 2)) + 2;
            alphayeast(pick(1, 1), pick(1, 2)) = alphayeast(pick(1, 1), pick(1, 2)) + 2;
        end
    end





    % Time
    time = time + Time(totalrate);
    run = run + 1;
if mod(run,100) == 0
    % New Data
    datatime = [datatime time];
    dataayeast(:, :, run) = ayeast;
    dataalphayeast(:, :, run) = alphayeast;
    datadiploidyeast(:, :, run) = diploidyeast;
    dataashmoo(:, :, run) = sum(ashmoo, 3);
    dataalphashmoo(:, :, run) = sum(alphashmoo, 3);
    dataafactor(:, :, run) = afactor;
    dataalphafactor(:, :, run) = alphafactor;
end

end