tic
VisOn = 2;
% VisOn: 2: Whole view, 1: Update in interval, 0: No view

view_intv = 20;

InsMode = 1;   % 1: Show SitA genotype, 2: Show accumlated damage level (white when no accumulated dam)
% 3: Show motility protein level (white when motility protein is saturated)
% 4: Transferred material (assuming florescent material)
% 5: Nutrient

Voro_Vis = 0;
OME_Vis= 0;
Text_Vis = 0;

Vs = 0.01;  % Self propulsion speed

maxAge = 50;  % the lifespan of a cell, in the units of timesteps
repAge = 10;  % minimal age to replicate
RD = 1;  % Random change in direction (maximal range in rad unit)
CT = 1;  % Chemotaxis property

EngMax = 3;
EngRep = EngMax*0.7;
DamTol = 1;
Flores = 0;

% <Initiation code>

Init_traA_div

% <Result recording>

TimeLength = 300;

RecInds = zeros(TypeN, TimeLength);  % Individual recording, Each row represents type and numbers is the subpopulation size
RecMots = zeros(TypeN, TimeLength);  % Recording motility protein concentration in each type
RecNutri = zeros(TypeN, TimeLength);  % Recording nutrient strorage in each type
RecDam = zeros(TypeN, TimeLength);  % Recording accumulated damage in each type

if Flores == 1
    RecFlores = zeros(TypeN, TimeLength);   % Florescent material transfer
end

DispersionRecording = 0;

if DispersionRecording == 1

    RecDispersion = zeros(2, TimeLength);  % Recording dispersion of each type (Radius of Gyration)

end

% <Nutrient map making>
if exist('traA_nutri_fluc') && traA_nutri_fluc == 1
    FlucDet = 1;   % 0: No fluctuation, 1: Fluctuation
else
    FlucDet = 0;
end

Tdur = 50;  % Durtaion of target nutrient map
MapContain = MapMake(FlucDet, 2, TimeLength, Tdur, NutVal);

if FlucDet == 0

    NutConstant = MapContain{1};

end

NutConsume = zeros(size(NutVal));   % Representing consumed amount of the nutrient

% <Timestep process>

for t = 1:TimeLength

    PosCell = PosCellNext;

    ReproRec = false(Nmax, 1);
    DeathRec = false(Nmax, 1);

    if VisOn~= 0

        hold off

    end

    TrgMapSN =  floor(t/Tdur) + 1;  % Target Nutrient Map Start Number

    TrgMapEN = ceil(t/Tdur) + 1;  % Targent Nutrient Map End number
    NutMapA = MapContain{TrgMapSN};
    NutMapB = MapContain{TrgMapEN};

    MapWeight = rem(t, Tdur)/Tdur;

    BasalNut = NutMapA*(1-MapWeight) + NutMapB*MapWeight;

    NutConsume = NutConsume*(1-Nut_conv_r);

    NutConsume = imgaussfilt(NutConsume, 0.3);

    NutConsume(NutConsume > -10^(-5)) = 0;   % Note that NutConsume is generally negative

    ActualNut = BasalNut + NutConsume;

    ActualNut = imgaussfilt(ActualNut, 0.3);

    ActualNut = max(ActualNut, 0);

    [GradVecX,GradVecY] = gradient(ActualNut);

    PosLivCell = PosCell(LiveMat, :);
    LiveIDs = ID_list(LiveMat);

    mov_ang_Liv = mov_ang(LiveMat);

    PosT = [PosLivCell; BD_mat];  % Total position including boundaries, excluding dead cells

    CONN_DT_Cell_ord = delaunay(PosLivCell);  % Delaunay triangularization connectivity (each row: [ID1ord, ID2ord, ID3ord])

    CONN_DT_Cell = LiveIDs(CONN_DT_Cell_ord);  % % Order is the order among the live individuals. This is the conversion to IDs

    [Vxy, CONN_VR] = voronoin(PosT);   % Vxy: the coordinates of voronoi vertices.
    % CONN_VR: each row show the group of voronoi vertex IDs that constitute an voronoi cell. Convex hull satisfied.

    HeadX = PosLivCell(:,1) + 0.5*BL*cos(mov_ang_Liv);   % Composed of live cells
    HeadY = PosLivCell(:,2) + 0.5*BL*sin(mov_ang_Liv);
    TailX = PosLivCell(:,1) - 0.5*BL*cos(mov_ang_Liv);
    TailY = PosLivCell(:,2) - 0.5*BL*sin(mov_ang_Liv);

    HeadX_All = PosCell(:,1) + 0.5*BL*cos(mov_ang);   % Row size of maximal cell number
    HeadY_All = PosCell(:,2) + 0.5*BL*sin(mov_ang);
    TailX_All = PosCell(:,1) - 0.5*BL*cos(mov_ang);
    TailY_All = PosCell(:,2) - 0.5*BL*sin(mov_ang);

    IDcolor = zeros(Nmax,3);

    for cm = 1:TypeN

        TypeIDs = find(Gtype == cm);
        IDcolor(TypeIDs, :) = repmat(TypeCmap(cm,:), length(TypeIDs),1 );

    end

    LiveIDcolor = IDcolor(LiveMat, :);

    if VisOn == 2

        MapUpdate = 1;

    elseif VisOn == 1

        if mod(t, view_intv) == 1
            MapUpdate=1;
        else
            MapUpdate=0;
        end

    else

        MapUpdate = 0;

    end

    if MapUpdate

        HueData = TA_GT(LiveMat,1)/TA_FN;
        OneMat = ones(size(HueData));

        TraAcolor = hsv2rgb([HueData , 0.7*OneMat, OneMat]);

        imagesc([-HW HW], [-HW HW], ActualNut, clim)
        set(gca,'Ydir','normal')
        colormap pink

        hold on

        H2T_X = TailX' - HeadX'  ; %  Head-to-tail vector's X values
        H2T_Y = TailY' - HeadY'  ; %  Head-to-tail vector's Y values

        TA_NaN = NaN(1, length(HeadX));

        TA_Cmap = ones(3, length(HeadX), 3);   % TraA Colormap

        TA_Cmap(1,:,1) = TraAcolor(:,1)';
        TA_Cmap(1,:,2) = TraAcolor(:,2)';
        TA_Cmap(1,:,3) = TraAcolor(:,3)';

        HeadXE = HeadX' - (0.5/SitA_FN)*H2T_X; % Extended head position (For the marker representing motility)
        TailXE = HeadX' + ((SitA_FN+0.5)/SitA_FN)*H2T_X;

        HeadYE = HeadY' - (0.5/SitA_FN)*H2T_Y;
        TailYE = HeadY'+((SitA_FN+0.5)/SitA_FN)*H2T_Y;

        scatter(HeadXE, HeadYE, 40, LiveIDcolor, 'filled')

        patch([HeadXE; TailXE; TA_NaN], [HeadYE; TailYE; TA_NaN], TA_Cmap, 'edgecolor', 'flat','FaceColor','none' ,   'linewidth',3)

        if InsMode == 1

            SA_GT_Liv = SA_GT(LiveMat, :);

            SitAden = zeros(sum(LiveMat),   SitA_FN );   % The density of SitA alleles in each family

            for s0 = 1:SitA_FN
                start_idx = (s0-1)*SitA_sN + 1;
                end_idx = s0*SitA_sN;
                SitAden(:, s0) = sum(SA_GT_Liv(:, start_idx:end_idx), 2)/SitA_sN;
            end

            for seg = 1:SitA_FN

                SegStartX = HeadX' + (seg - 1)*H2T_X/SitA_FN;
                SegEndX = HeadX' + (seg )*H2T_X/SitA_FN;

                SegStartY = HeadY' + (seg - 1)*H2T_Y/SitA_FN;
                SegEndY = HeadY' + (seg )*H2T_Y/SitA_FN;

                EndCap = NaN(size(SegStartX));

                BarCmap = ones(3, length(HeadX), 3);
                BarCmap(1,:,1) = SitAden(:, seg)';
                BarCmap(1,:,2) = SitAden(:, seg)';
                BarCmap(1,:,3) = SitAden(:, seg)';

                patch([SegStartX; SegEndX; EndCap], [SegStartY; SegEndY; EndCap], BarCmap, 'edgecolor', 'flat','FaceColor','none' ,'linewidth',2)

            end

        elseif  InsMode == 2

            DamLevel = AccDam(LiveMat)/DamTol;
            Dam_Cmap = ones(3, length(HeadX), 3);
            Dam_Cmap(1,:,1) = 1-DamLevel;
            Dam_Cmap(1,:,2) = 1-DamLevel;
            Dam_Cmap(1,:,3) = 1-DamLevel;

            patch([HeadX'; TailX'; TA_NaN], [HeadY'; TailY'; TA_NaN], Dam_Cmap, 'edgecolor', 'flat','FaceColor','none' ,   'linewidth',2)

        elseif InsMode == 3

            MotLevel = Mot(LiveMat);
            Mot_Cmap = ones(3, length(HeadX), 3);
            Mot_Cmap(1,:,1) = MotLevel;
            Mot_Cmap(1,:,2) = MotLevel;
            Mot_Cmap(1,:,3) = MotLevel;

            patch([HeadX'; TailX'; TA_NaN], [HeadY'; TailY'; TA_NaN], Mot_Cmap, 'edgecolor', 'flat','FaceColor','none' ,   'linewidth',2)

        elseif InsMode == 4

            FlorLevel = FLORES(LiveMat);
            Flor_Cmap = ones(3, length(HeadX), 3);
            Flor_Cmap(1,:,1) = FlorLevel;
            Flor_Cmap(1,:,2) = FlorLevel;
            Flor_Cmap(1,:,3) = FlorLevel;

            patch([HeadX'; TailX'; TA_NaN], [HeadY'; TailY'; TA_NaN], Flor_Cmap, 'edgecolor', 'flat','FaceColor','none' ,   'linewidth',2)

        elseif InsMode == 5
            NutriLevel = EngMat(LiveMat)/EngMax;
            Nutri_Cmap = ones(3, length(HeadX), 3);
            Nutri_Cmap(1,:,1) = NutriLevel;
            Nutri_Cmap(1,:,2) = NutriLevel;
            Nutri_Cmap(1,:,3) = NutriLevel;

            patch([HeadX'; TailX'; TA_NaN], [HeadY'; TailY'; TA_NaN], Nutri_Cmap, 'edgecolor', 'flat','FaceColor','none' ,   'linewidth',2)

        end

        patch([TailX'; TailXE; TA_NaN], [TailY'; TailYE; TA_NaN], TA_Cmap, 'edgecolor', 'flat','FaceColor','none' ,   'linewidth',2)   % The linewidth here should be the same with the line representing SitAI barcodes

    end

    PosCellNext =  PosCell ;

    ID_shf_ord = randperm(length(LiveIDs));

    LiveIDs_SHF = LiveIDs(ID_shf_ord);  % Shuffling of the IDs of the alive cells

    PopSize0 = sum(LiveMat);   % Population size before interactions occur

    AreaMat = zeros(1, PopSize0);

    for i0 = 1:PopSize0

        idx0 = LiveIDs(i0);

        FCx = PosCellNext(idx0,1);   % Here, ID should be used

        FCy = PosCellNext(idx0,2);

        VrVtxIDs0 = [CONN_VR{i0} CONN_VR{i0}(1)];   % Here, order should be used
        if ~ismember(1, VrVtxIDs0)

            Vpoly = polyshape(Vxy(VrVtxIDs0,1)-FCx, Vxy(VrVtxIDs0,2)-FCy);   % Voronoi area polyshape

            PolyOverLay = intersect(Vpoly,Cpoly);

            if MapUpdate && Voro_Vis

                PolyOverLay0 = translate(PolyOverLay, [FCx, FCy]);

            end

            AreaMat(i0) = area(PolyOverLay);

        end

    end

    for i = 1:PopSize0

        idx = LiveIDs_SHF(i);  % ID of the focal cell

        % <Acquisition of the nutrient>

        FCx = PosCellNext(idx,1);  % Focal cell x coordinates
        FCy = PosCellNext(idx,2);  % Focal cell y coordinates

        id_ord = ID_shf_ord(i);   % Order among the LiveMat

        % Voronoi vertex IDs, this should use id_ord because Voronoi was performed before shuffling of the IDs

        VrVtxIDs = [CONN_VR{id_ord} CONN_VR{id_ord}(1)];

        [~, grad_idX] = min(abs(FCx - GradX)) ;   % Finding row and column in nutrient map (matrix)
        [~, grad_idY] = min(abs(FCy - GradY)) ;

        if ~ismember(1, VrVtxIDs)

            NutDen = ActualNut(grad_idY, grad_idX);  % Local nutrient density

            NutArea = AreaMat(id_ord);

            NutGet = nutriK*NutDen*NutArea;

            EngMat(idx) = EngMat(idx) + NutGet;

            NutConsume((grad_idY-hwY):(grad_idY+hwY), (grad_idX-hwX):(grad_idX+hwX)) = ...
                NutConsume((grad_idY-hwY):(grad_idY+hwY), (grad_idX-hwX):(grad_idX+hwX)) - 3*NutGet*NutRed;

        end

        % <Cell Death>

        AgeDeathCond = 0;
        if AgeMat(idx) > maxAge

            AgeDeathCond=1;

        elseif AgeMat(idx) >= 0.9*maxAge

            if rand(1) > 0.5
                AgeDeathCond=1;
            end

        end


        if (EngMat(idx) < 0)  || AgeDeathCond  || (AccDam(idx) > DamTol)

            DeathRec(idx) = 1;

        end

        % <Repulsion to too close neighbors or boundary>

        [Row_DT, ~] = find(CONN_DT_Cell == idx);

        con_IDs_Ordered = CONN_DT_Cell(Row_DT, :);  % Each row is pair of IDs

        con_IDs_total = unique(con_IDs_Ordered);  % Total DT neighbor IDs in a single set

        con_IDs_total = con_IDs_total(con_IDs_total~= idx);  % Excluding the self (focal cell)

        DTneighborDX = PosCell(con_IDs_total, 1) - FCx;  % Delta X to the DT neighbors
        DTneighborDY = PosCell(con_IDs_total, 2) - FCy;

        Repul_Logic = ( (DTneighborDX.^2 + DTneighborDY.^2 ) < ContactMin^2);

        if any(Repul_Logic)

            Rep_Ords = find(Repul_Logic);

            RepDX = DTneighborDX(Rep_Ords);
            RepDY = DTneighborDY(Rep_Ords);

            DistInvSqVals = 0.1 * 1./(RepDX.^2 + RepDY.^2);

            DistInvSqVals = min(DistInvSqVals, 2);

            Rep_Vec = [-sum(RepDX.*DistInvSqVals), -sum(RepDY.*DistInvSqVals)];

        else

            Rep_Vec = [0, 0];

        end

        d_bd_x = (HW*BD_r - abs(FCx));  % X Distance (positive) to the boundary
        BD_rep = 0.1;   % repulsion from the boundary

        if  d_bd_x < 0.5*HW*(1-BD_r)   % if FCx is too close to the boundary, HW*(1-BD_r) is the margin width

            bd_rep_x = -sign(FCx)*exp(-1*d_bd_x);
            Rep_Vec(1) = Rep_Vec(1) + BD_rep*bd_rep_x;

        end

        d_bd_y = (HW*BD_r - abs(FCy));  % X Distance (positive) to the boundary
        if  d_bd_y < 0.5*HW*(1-BD_r)   % if FCx is too close to the boundary, HW*(1-BD_r) is the margin width

            bd_rep_y = -sign(FCy)*exp(-1*d_bd_y);
            Rep_Vec(2) = Rep_Vec(2) + BD_rep*bd_rep_y;

        end

        % <Green-beard interactions>

        GBE_Logic = ( (DTneighborDX.^2 + DTneighborDY.^2 ) < (BL)^2);

        G_cands = con_IDs_total(GBE_Logic);   % Candidates for the green-beard interactions (at most a single GBE allowed)

        if i ~= PopSize0

            Rem_IDs = LiveIDs_SHF( (i+1):end);

        else

            Rem_IDs = [];

        end

        remTF = ismember(G_cands, Rem_IDs); % whether candidate IDs are yet to be focal individuals

        Rem_IDs = G_cands(remTF);

        IntExist = 0;   % Check whether intersection exists;

        if ~isempty(Rem_IDs)

            F_head_xy = [HeadX_All(idx), HeadY_All(idx)];
            F_tail_xy = [TailX_All(idx), TailY_All(idx)];

            Rem_IDs_sf = Rem_IDs(randperm(length(Rem_IDs)));   % shuffled

            for g0 = 1:length(Rem_IDs_sf)

                G_cand_ID = Rem_IDs_sf(g0);

                G_head_xy = [HeadX_All(G_cand_ID), HeadY_All(G_cand_ID)];
                G_tail_xy = [TailX_All(G_cand_ID), TailY_All(G_cand_ID)];

                IntersectTF = isIntersecting(F_head_xy, F_tail_xy, G_head_xy, G_tail_xy);

                if IntersectTF
                    IntExist = 1;
                    IntID = G_cand_ID;
                    break

                end

            end

        end

        if IntExist == 1

            p_OME = OMEprob(TA_GT(idx, :), TA_GT(IntID,:), TA_FN, TA_sN);

            if rand(1) < p_OME    % OME occured

                if  MapUpdate && OME_Vis
                    plot( [PosCell(idx, 1),  PosCell(G_cand_ID, 1)], [PosCell(idx, 2),  PosCell(G_cand_ID, 2)], 'b' , 'linewidth',2)

                end

                EngA = EngMat(idx); EngB = EngMat(IntID);
                [EngAr, EngBr] = OME(EngA, EngB, 0.2, 0.1);   % Exchange of nutrients (energy)
                EngMat(idx) = EngAr; EngMat(IntID) = EngBr;

                motA = Mot(idx); motB = Mot(IntID);
                [motAr, motBr] = OME(motA, motB, 0.8, 0.4);   % Exchange of motility proteins
                Mot(idx) = motAr;  Mot(IntID) = motBr;

                toxA = SAcyto(idx,:);   toxB = SAcyto(IntID,:);
                [toxAr, toxBr] = OME(toxA, toxB, 0.9, 0.4);
                SAcyto(idx,:) = toxAr;  SAcyto(IntID,:) = toxBr;


                if Flores == 1
                    florA = FLORES(idx,:);   florB = FLORES(IntID,:);
                    [florAr, florBr] = OME(florA, florB, 0.95, 0.5);
                    FLORES(idx,:) = florAr; FLORES(IntID,:)=florBr;

                end

            end

        end

        if MapUpdate  &&  Text_Vis

            text(FCx, FCy, num2str(idx),'HorizontalAlignment', 'center')

        end

        if MapUpdate && Voro_Vis

            if ~ismember(1, VrVtxIDs)
                plot(Vxy(VrVtxIDs,1), Vxy(VrVtxIDs,2), 'k')  % Plotting Voronoi cell boundaries
            end

        end

        % <Movement of cells>

        Nut_Vec = [GradVecX(grad_idY, grad_idX), GradVecY(grad_idY, grad_idX)];  % The complication arises due to Cartesian and Row*Col system

        rand_ang = RD*(rand(1)-0.5)*exp(-sum((10*Nut_Vec).^2));

        prop_ang = mov_ang(idx) + rand_ang;

        Vec_Prop = INX*[cos(prop_ang), sin(prop_ang)];

        Vec_CTX = CT*Nut_Vec;

        Vec_Mov0 = Vec_Prop + Vec_CTX + Rep_Vec;

        Vec_Mov0 = Vs*Vec_Mov0/(norm(Vec_Mov0)); % Moving vector without any voronoi constraint and motility protein

        mov_ang(idx) = atan2(Vec_Mov0(2), Vec_Mov0(1));  % Moving angle update

        Vec_Mov =  M2V(Mot(idx))*Vec_Mov0; % Movement considering the motility protein

        Mov_pathX = [FCx, FCx + Vec_Mov(1)];
        Mov_pathY = [FCy, FCy + Vec_Mov(2)];

        PosCellNext(idx, :) = [FCx + Vec_Mov(1), FCy + Vec_Mov(2)];


        % <Collision to Voronoi boundary and total boundary>

        for j=1:(length(VrVtxIDs)-1)

            VR1xy = Vxy(VrVtxIDs(j), :);   % The xy coordinates of Voronoi vertex 1

            VR2xy = Vxy(VrVtxIDs(j+1), :);

            if ~any(isinf([VR1xy, VR2xy]))
                P_intsect = IntsecPoint([FCx, FCy], PosCellNext(idx, :), VR1xy,  VR2xy);  % Intersection point

                if ~isempty(P_intsect)

                    PosCellNext(idx, :) = P_intsect;
                    break
                end
            end

        end

        if any(  abs(PosCellNext(idx, :)) >  HW*BD_r)   % Deterred movement toward the boundary

            [rev_x, rev_y] = BoundaryRev(FCx, FCy, Vec_Mov(1), Vec_Mov(2), HW*BD_r);
            PosCellNext(idx, :) = [rev_x, rev_y];

        end

        PosCellNext(idx, 1) = min(PosCellNext(idx, 1), BD_r *HW);
        PosCellNext(idx, 1) = max(PosCellNext(idx, 1), -BD_r *HW);
        PosCellNext(idx, 2) = min(PosCellNext(idx, 2), BD_r *HW);
        PosCellNext(idx, 2) = max(PosCellNext(idx, 2), -BD_r *HW);

        % <Reproduction>

        if (rand(1) < 0.2)  &&  (rand(1) < (1- AccDam(idx)/DamTol ))   % Randomness in reproduction and effect of the accumulated damage

            RC_space = (NutArea > 0.5*Circ_area);   % RC_space: reproduction condition (RC) of space

            RC_nutri = (EngMat(idx)>EngRep);   % Condition on energy

            RC_N = (sum(LiveMat) < Nmax);    % Condition on total cell number

            RC_age = (AgeMat(idx) >= repAge);   % Condition on age

            RC = RC_space*RC_nutri*RC_N*RC_age;

        else

            RC = 0;

        end

        if RC

            ReproRec(idx) = 1;

        end

    end

    % <Material synthesis, decay & cell senescence & energy usage>

    AgeMat(LiveMat) = AgeMat(LiveMat) + 1;   % Cells getting older

    mot_syn = 0.05;   % Motility synthesis in each time step
    mot_dec = 0.002;   % Motility decay in each time step
    Mot(M_GT) = Mot(M_GT) + mot_syn;
    Mot(~M_GT) = Mot(~M_GT) - mot_dec;
    Mot = min(Mot, 1);  Mot = max(Mot, 0);

    tox_syn = 0.1;
    tox_dec = 0.002;
    SAcyto(SA_GT) = SAcyto(SA_GT) + tox_syn;
    SAcyto(~SA_GT) = SAcyto(~SA_GT) - tox_dec;
    SAcyto = min(SAcyto, 1); SAcyto = max(SAcyto,0);

    BaseCost = 0.03;
    CostSitA = 0.03*mean(SA_GT, 2);  % Metabolic cost of SitA alleles
    CostSitI =0.03*mean(SI_GT, 2);  % Metabolic cost of SitI alleles
    CostMot = 0.03*M_GT;  % Metabolic cost of having motility allele

    MetaCost = BaseCost + CostSitA + CostSitI + CostMot;   % factors for higher metabolism: more sitA & sitI alleles, Motility,

    EngMat(LiveMat) = EngMat(LiveMat) - MetaCost(LiveMat);

    EngMat = min(EngMat, EngMax);
    EngMat = max(EngMat, 0);

    if Flores == 1

        flor_syn = 0.2;
        flor_dec = 0.001;
        FLORES(FLR_GT) = FLORES(FLR_GT) + flor_syn;
        FLORES(~FLR_GT) = FLORES(~FLR_GT) - flor_dec;
        FLORES = min(FLORES, 1); FLORES = max(FLORES,0);

    end

    % <Cellular damage due to toxins>

    DamVals = SitAIdam(SAcyto, SI_GT, SitA_FN, SitA_sN, LiveMat);

    AccDam(LiveMat) = AccDam(LiveMat) + DamVals;

    % <Death process>

    LiveMat(DeathRec) = 0;
    EngMat(DeathRec) = 0;
    AgeMat(DeathRec) = 0;

    AccDam(DeathRec) = 0;
    TA_GT(DeathRec,:) = 0;
    SA_GT(DeathRec,:) = 0;
    SAcyto(DeathRec,:) = 0;
    SI_GT(DeathRec,:) =0;
    M_GT(DeathRec) = 0;
    Mot(DeathRec) = 0;
    Gtype(DeathRec) = 0;

    DeathNum = sum(DeathRec);

    % <Birth Process>

    ReproRecActual = ReproRec;

    ReproRecActual(DeathRec == 1) = 0;

    PossibleMaxNumRepro = sum(LiveMat == 0);    % Vacant position for the newborns

    BirthNum = min(sum(ReproRecActual), PossibleMaxNumRepro);

    ReproIDs_all = ID_list(ReproRecActual);

    TrunShf = randperm( sum(ReproRecActual), BirthNum );

    ReproIDs = ReproIDs_all(TrunShf);

    for Bn = 1:BirthNum

        rID = ReproIDs(Bn);

        NewCellID = find(LiveMat == 0, 1);

        LiveMat(NewCellID) = 1;

        EngMat(NewCellID) = EngMat(rID)*0.5;
        EngMat(rID) = 0.5*EngMat(rID);

        % Add position and angle, etc

        CellRndAng = 2*pi*rand(1);  % Random perturbation in position of new cell

        NewPos = 0.5*( PosCellNext(rID, :) + PosCell(rID, :) ) + r_abs*rand(1)*([cos(CellRndAng), sin(CellRndAng)]);

        PosCellNext(NewCellID, :) = NewPos;

        mov_ang(NewCellID) = mov_ang(rID) + 0.5*(rand(1)-0.5);

        AgeMat(rID) = 0;
        AgeMat(NewCellID) = 0;

        TA_GT(NewCellID, :) = TA_GT(rID, :);
        SA_GT(NewCellID, :) = SA_GT(rID, :);
        SI_GT(NewCellID, :) = SI_GT(rID, :);
        M_GT(NewCellID, :) = M_GT(rID, :);

        SAcyto(NewCellID,:)  = 0.5*SAcyto(rID,:);    % Division of cytosolic SitA materials (toxin)
        SAcyto(rID,:)  = 0.5*SAcyto(rID,:);

        Mot(NewCellID,:)  = 0.5*Mot(rID,:);    % Division of motility proteins
        Mot(rID,:)  = 0.5*Mot(rID,:);

        if Flores == 1

            FLR_GT(NewCellID, :) = FLR_GT(rID, :);
            FLORES(NewCellID, :) = 0.5*FLORES(rID,:);   % Florescent material division
            FLORES(rID,:) = 0.5*FLORES(rID,:);

        end

        damReduc = 0.2;   % The reduction of accumulated damage by replication (assumed to be smaller than 0.5)
        AccDam(NewCellID) = damReduc*AccDam(rID);  % Order is important: AccDam(NewCellID) should be first calculated
        AccDam(rID) = damReduc*AccDam(rID);

        Gtype(NewCellID) = Gtype(rID);

    end

    % <Recording>

    for tn = 1:TypeN   % tn represents type

        TypePos = (Gtype == tn);

        TypeLive = (LiveMat & TypePos);

        TypeSize = sum(TypeLive);  % Subpopulation size

        TypeMot = mean(Mot(TypeLive));    % Average motility protein level

        TypeNutri = mean(EngMat(TypeLive))/EngMax;  % Average nutrient storage (proportion to maximal storage)

        TypeDam = mean(AccDam(TypeLive));

        if TypeSize == 0

            TypeMot = 0;
            TypeNutri = 0;
            TypeDam = 0;

        end

        RecInds(tn ,t) = TypeSize;
        RecMots(tn ,t) = TypeMot;
        RecNutri(tn, t) = TypeNutri;
        RecDam(tn, t) = TypeDam;

        if DispersionRecording == 1   % Whether or not recording the dispersion of cells

            TypeCoords = PosCellNext(TypeLive, :);

            centroid = mean(TypeCoords);

            squaredDistances = sum((TypeCoords - centroid).^2, 2);

            R_Gyration = sqrt(mean(squaredDistances));  % Radius of Gyration

            RecDispersion(tn, t) = R_Gyration;

        end

        if Flores == 1
            RecFlores(tn, t) = mean(FLORES(TypeLive));

        end

    end

    if MapUpdate

        xlim([-3 3])
        ylim([-3 3])
        axis square

        pause(0.01)
        t

    else
        t
    end

    if (sum(RecInds(:, t) > 0)  == 1) && (sum(RecInds(:, t) > 0.95*Nmax))   % One type survived and saturated

        break

    end

end

if Flores == 1
    RecFlores(isnan(RecFlores)) = 0;
end

if DispersionRecording == 1
    RecDispersion(isnan(RecDispersion)) = 0;
end

toc

