NAT = 100;   % Group A size (including minority)

NAm = 10;   % Group A minority size (traA minority)

NBT = 100;   % Group B size (including minority)

NBm = 10;   % Group B minority size (traA minority)

NT = NAT + NBT;

TypeN = 6;

Nmax = 500;   % Maximal cell number

HW = 5;  % Halfwidth of arena

ID_list = 1:Nmax;   % For designating IDs to cells

if ~exist('MassComp')

    NutriCond = 2;   % Condition of the initial nutrient distribution
    % 1: High
    % 2: Low
    % 3: Nutrition source

    minExist = 1;  % Whether minority group exists
    % 0: no minority group (Some traA_Cond are overlappin in this case)
    % 1: minority group exists

    traA_Cond = 4;   % Matching condition
    % 1: Majority: same traA, minority: partially compatible traA to majority
    % 2: Majority: same traA, minority: incompatible to majority
    % 3: Majority: partlly compatible traA, minority: partially compatible traA to majority
    % 4: Majority: partlly compatible traA, minority: incompatible traA to majority
    % 5: (Control) Majority: same traA, minority has same traA with majority
    % 6: (Control) Majority: partlly compatible traA, minority has same traA with majority

    SpaCond = 1;  % Spatial distribution of the two groups
    % 1: Well mixed
    % 2: Separated

    RPS_cond = 3;
    % 1: RP (rock left, paper right)
    % 2: PS (paper left, scissors right)
    % 3: RS (rock left, scissors right)

end

% Color schematics
% Rock: red, sitI / incompatible rock (traA minority): yellow (Type 1, 2)
% Paper: blue, no sitAI / incompatible paper (traA minority): light blue (Type 3, 4)
% Scissors: white, sitA and sitI / incompatible scissors (traA minority): gray  (Type 5, 6)

% Possessions of bacteria: Energy (nutrient), Motility protein, SitA accumulation
% Properties of bacteria: Alive, Age, Genotype of TraA and SitAI, Genotype for motility, Accumulated damage, Position, Moving angle,

LiveMat = logical(zeros(Nmax,1));  % Alive or dead: 0 or 1
LiveMat(1:(NAT+NBT)) = 1;

InitNum = zeros(TypeN,1);   % Subpopulation at the initial condition

EngMat =  zeros(Nmax,1);  % Energy storage
EngMat(1:NT) = 1;

T1C = [1, 0, 0]; % Rock major
T2C = [1, 1, 0]; % Rock minor
T3C = [0, 0, 1];   % Paper major
T4C = [0.5, 0.5, 1];  % Paper minor
T5C = [1, 1, 1]; % Scissors major
T6C = [0.5, 0.5, 0.5]; % Scissors minor

TypeCmap = [T1C; T2C; T3C; T4C; T5C; T6C];

AgeMat =  zeros(Nmax,1);  % Age of the cells

AgeMat(1:NT) = maxAge*rand(NT,1);

BD_r = 0.85;

AccDam = zeros(Nmax,1);   % Accmulated damage due to SitA

Gtype = zeros(Nmax,1);  % Types for represent on the visual map (different marker color)

switch RPS_cond
    case 1   % RP
        Gtype(1:NAT) = 1;
        Gtype((NAT+1):(NAT+NBT)) = 3;

    case 2   % PS

        Gtype(1:NAT) = 3;
        Gtype((NAT+1):(NAT+NBT)) = 5;

    case 3   % RS
        Gtype(1:NAT) = 1;
        Gtype((NAT+1):(NAT+NBT)) = 5;
end


if minExist == 1
    Gtype((NAT - NAm + 1):NAT) = Gtype((NAT - NAm + 1):NAT) + 1;
    Gtype((NAT+NBT - NBm + 1):(NAT+NBT)) = Gtype((NAT+NBT - NBm + 1):(NAT+NBT)) + 1;
end

INX = 0.05; % Inertia to sustain the previous moving direction
nutriK = 10;   % Nutrient absorption coefficient
r_abs = 0.1; % Radius of nutrien absorbtion
C_theta = 0:0.4:2*pi;

CircleX =  r_abs.*cos(C_theta);
CircleY =  r_abs.*sin(C_theta);

Cpoly = polyshape(CircleX, CircleY);   % Circle polyshape
Circ_area = area(Cpoly);

BL = 0.15;   % body length of myxobacteria

ContactMin = 0.05;  % Minimal contact range

n_boundary = 30;  % should consider the inner boundary size

bd_series0 = linspace(-HW, HW, n_boundary+2);

bd_series = bd_series0(2:(n_boundary+1))';
bd_constant = HW*ones(size(bd_series));

R_bd = [bd_constant, bd_series];  % Right boundary
L_bd = [-bd_constant, bd_series];  % Left boundary
U_bd = [bd_series, bd_constant]; % Upper boundary
B_bd = [bd_series, -bd_constant]; % Bottom boundary
TestBedVtx = [HW, -HW; -HW, HW; HW, HW; -HW, -HW];

GradX = -HW:0.02:HW;  % Row of NutVal
GradY = (-HW:0.02:HW)' ;

switch NutriCond

    case 1
        NutVal = 2*ones( length(GradY), length(GradX)  );
        Nut_conv_r = 0.03;    % Convergence rate of the nutrient map to the original nutrient map, modified in this code


    case 2
        NutVal = 1*ones( length(GradY), length(GradX)  );
        Nut_conv_r = 0.01;    % Convergence rate of the nutrient map to the original nutrient map, modified in this code

    case 3

        NutVal = 3*exp(- GradX.^2  - (GradY).^2);
        NutVal = max(NutVal, 0.05);
        Nut_conv_r = 0.05;    % Convergence rate of the nutrient map to the original nutrient map, modified in this code

end

NutValOri = NutVal;   % Original nutrient map

clim = [0, 4];
NutValZero = zeros(size(NutVal));

R_reduc = 0.1;  % Radius of the reduction in nutrient map
[NutRed, hwX, hwY] = BasalNutriReduc(R_reduc, GradX, GradY);
% NutRed is the nutrient reduction map by nutrient absorption of a cell
% hwX is the half width of the nutrient reduction map

BD_mat = [R_bd; L_bd; U_bd; B_bd; TestBedVtx];   % Boundary xy values

mov_ang = zeros(Nmax, 1);
mov_ang(1:NT) = 2*pi*rand(NT,1);

PosCell = zeros(Nmax, 2);

if minExist == 0

    for TP = 1:2   % for type A and type B

        if SpaCond == 1 % Well mixed

            MajRangeX = [-1 1]; % Range of X position for majority group
            MajRangeY = [-1 1];

        elseif SpaCond == 2  % separated

            if TP == 1
                MajRangeX = [-1,  0];
                MajRangeY = [-1, 1];

            elseif TP == 2
                MajRangeX = [0, 1];
                MajRangeY = [-1, 1];

            end
        end

        if TP == 1
            PosRange0 = 1:NAT;
        elseif TP == 2
            PosRange0 = (NAT+1):(NAT+NBT);
        end

        Lx = max(MajRangeX) - min(MajRangeX);
        Ly = max(MajRangeY) - min(MajRangeY);

        MajX = min(MajRangeX) + rand(length(PosRange0), 1)*Lx;
        MajY = min(MajRangeY) + rand(length(PosRange0), 1)*Ly;
        PosCell(PosRange0, :) = [MajX, MajY];

    end

end

if minExist == 1
    for TP = 1:2   % for type A and type B

        if SpaCond == 1 % Well mixed

            MajRangeX = [-1 1]; % Range of X position for majority group
            MajRangeY = [-1 1];

        elseif SpaCond == 2  % separated

            if TP == 1
                MajRangeX = [-1,  0];
                MajRangeY = [-1, 1];

            elseif TP == 2
                MajRangeX = [0, 1];
                MajRangeY = [-1, 1];

            end

        end

        Lx = abs(MajRangeX(2) - MajRangeX(1));
        Ly = abs(MajRangeY(2) - MajRangeY(1));

        MajArea = Lx*Ly;  % Area of the Majority group including the minority group area

        minArea = (NAm/NAT)*MajArea;  % To maintain same cell density for major and minor groups

        minWidth = sqrt(minArea);  % Width of a square for minority group

        minMidX = min(MajRangeX) + (minWidth/2) + rand(1)*(Lx - minWidth);   % Mid X point of minority region
        minMidY = min(MajRangeY) + (minWidth/2) + rand(1)*(Ly - minWidth);

        minRangeX = [minMidX - minWidth/2 ,  minMidX + minWidth/2];
        minRangeY = [minMidY - minWidth/2 ,  minMidY + minWidth/2];

        minX1 = min(minRangeX);
        minX2 = max(minRangeX);

        minY1 = min(minRangeY);
        minY2 = max(minRangeY);

        if TP == 1
            McountRange = 1:(NAT - NAm);
            min_range = (NAT-NAm+1):(NAT);

        elseif TP == 2
            McountRange = (NAT+1):(NAT+NBT - NBm);
            min_range = (NAT+NBT-NBm+1):(NAT + NBT);
        end


        for Mcount = McountRange

            MajX = min(MajRangeX) + rand(1)*Lx;
            MajY = min(MajRangeY) + rand(1)*Ly;

            while (MajX > minX1) && (MajX < minX2) && (MajY > minY1) && (MajY < minY2)    % included in the minority group region

                MajX = min(MajRangeX) + rand(1)*Lx;
                MajY = min(MajRangeY) + rand(1)*Ly;

            end

            PosCell(Mcount, :) = [MajX, MajY];

        end

        minXvals = minX1 + (minX2 - minX1)*rand(NAm, 1);
        minYvals = minY1 + (minY2 - minY1)*rand(NBm, 1);

        PosCell(min_range, :) = [minXvals, minYvals];

    end

end

PosCell((NAT+NBT+1):end, :) = NaN;

PosCellNext = PosCell;   % PosCellNext is the position in the next time step

% <Genotype Determination>

TA_FN = 6;   % TraA family number
TA_sN = 10;   % Subtype number in each family of TraA

TA_GT = zeros(Nmax, 2);

SitA_FN = 5;    % SitA (and cognate SitI) family number
SitA_sN = 7;    % Subtype number in each family of SitA (and cognate SitI)

SA_GT = zeros(Nmax, SitA_FN*SitA_sN); % Genotype for SitA
SI_GT = zeros(Nmax, SitA_FN*SitA_sN);  % Cytoplasmic SitA

SitAIseg = [0 1 0 1 0 1 0];
SitAIzero = [0 0 0 0 0 0 0];
SitAIone = [1 1 1 1 1 1 1];

Config_saturated = [SitAIone, SitAIone, SitAIone, SitAIseg, SitAIseg];  % First 3 segment filled, rest 4 are half
Config_deficient = [SitAIzero, SitAIzero, SitAIzero, SitAIseg, SitAIseg]; % First 3 segment zero, rest 4 are half

% Rock
R_sitA = Config_deficient;
R_sitI = Config_saturated;

% Paper
P_sitA = Config_deficient;
P_sitI = Config_deficient;

% Scissors
S_sitA = Config_saturated;
S_sitI = Config_saturated;

switch RPS_cond
    % 1: RP (rock left, paper right)
    % 2: PS (paper left, scissors right)
    % 3: RS (rock left, scissors right)

    case 1

        A_sitA = R_sitA;
        A_sitI = R_sitI;

        B_sitA = P_sitA;
        B_sitI = P_sitI;

    case 2

        A_sitA = P_sitA;
        A_sitI = P_sitI;

        B_sitA = S_sitA;
        B_sitI = S_sitI;

    case 3

        A_sitA = R_sitA;
        A_sitI = R_sitI;

        B_sitA = S_sitA;
        B_sitI = S_sitI;

end

A_MAJ_range = 1:(NAT - NAm);
A_min_range = (NAT-NAm+1):(NAT);
A_all_range = 1:NAT;

B_MAJ_range = (NAT+1):(NAT+NBT - NBm);
B_min_range = (NAT+NBT-NBm+1):(NAT + NBT);
B_all_range = (NAT+1):(NAT+NBT);

if minExist == 0
    A_MAJ_range = 1:NAT;
    A_min_range = [];

    B_MAJ_range = (NAT+1):(NAT+NBT);
    B_min_range = [];

end

switch traA_Cond
    % 1: Majority: same traA, minority: partially compatible traA to majority
    % 2: Majority: same traA, minority: incompatible to majority
    % 3: Majority: partlly compatible traA, minority: partially compatible traA to majority
    % 4: Majority: partlly compatible traA, minority: incompatible traA to majority

    % Majority partial compatibility is determined by allele difference in the same traA family

    case 1

        TA_GT(A_MAJ_range,1) = 3;
        TA_GT(A_MAJ_range,2) = 1;

        TA_GT(A_min_range,1) = 2;
        TA_GT(A_min_range,2) = 1;

        TA_GT(B_MAJ_range,1) = 3;
        TA_GT(B_MAJ_range,2) = 1;

        TA_GT(B_min_range,1) = 4;
        TA_GT(B_min_range,2) = 1;

        SA_GT(A_all_range, :) = repmat(A_sitA, [(NAT), 1]);
        SI_GT(A_all_range, :) =  repmat(A_sitI, [(NAT), 1]);

        SA_GT(B_all_range, :) = repmat(B_sitA, [(NBT), 1])  ;
        SI_GT(B_all_range, :) = repmat(B_sitI, [(NBT), 1])  ;

    case 2

        TA_GT(A_MAJ_range,1) = 3;
        TA_GT(A_MAJ_range,2) = 1;

        TA_GT(A_min_range,1) = 1;
        TA_GT(A_min_range,2) = 1;

        TA_GT(B_MAJ_range,1) = 3;
        TA_GT(B_MAJ_range,2) = 1;

        TA_GT(B_min_range,1) = 5;
        TA_GT(B_min_range,2) = 1;

        SA_GT(A_all_range, :) =  repmat(A_sitA, [(NAT), 1]);
        SI_GT(A_all_range, :) =  repmat(A_sitI, [(NAT), 1]);

        SA_GT(B_all_range, :) = repmat(B_sitA, [(NBT), 1])  ;
        SI_GT(B_all_range, :) = repmat(B_sitI, [(NBT), 1])  ;


    case 3

        TA_GT(A_MAJ_range,1) = 3;
        TA_GT(A_MAJ_range,2) = 1;

        TA_GT(A_min_range,1) = 2;
        TA_GT(A_min_range,2) = 1;

        TA_GT(B_MAJ_range,1) = 3;
        TA_GT(B_MAJ_range,2) = 6;

        TA_GT(B_min_range,1) = 4;
        TA_GT(B_min_range,2) = 6;

        SA_GT(A_all_range, :) =  repmat(A_sitA, [(NAT), 1]);
        SI_GT(A_all_range, :) =  repmat(A_sitI, [(NAT), 1]);

        SA_GT(B_all_range, :) = repmat(B_sitA, [(NBT), 1])  ;
        SI_GT(B_all_range, :) = repmat(B_sitI, [(NBT), 1])  ;

    case 4

        TA_GT(A_MAJ_range,1) = 3;
        TA_GT(A_MAJ_range,2) = 1;

        TA_GT(A_min_range,1) = 1;
        TA_GT(A_min_range,2) = 1;

        TA_GT(B_MAJ_range,1) = 3;
        TA_GT(B_MAJ_range,2) = 6;

        TA_GT(B_min_range,1) = 5;
        TA_GT(B_min_range,2) = 6;

        SA_GT(A_all_range, :) =  repmat(A_sitA, [(NAT), 1]);
        SI_GT(A_all_range, :) =  repmat(A_sitI, [(NAT), 1]);

        SA_GT(B_all_range, :) = repmat(B_sitA, [(NBT), 1])  ;
        SI_GT(B_all_range, :) = repmat(B_sitI, [(NBT), 1])  ;

    case 5

        TA_GT(A_MAJ_range,1) = 3;
        TA_GT(A_MAJ_range,2) = 1;

        TA_GT(A_min_range,1) = 3;
        TA_GT(A_min_range,2) = 1;

        TA_GT(B_MAJ_range,1) = 3;
        TA_GT(B_MAJ_range,2) = 1;

        TA_GT(B_min_range,1) = 3;
        TA_GT(B_min_range,2) = 1;

        SA_GT(A_all_range, :) = repmat(A_sitA, [(NAT), 1]);
        SI_GT(A_all_range, :) =  repmat(A_sitI, [(NAT), 1]);

        SA_GT(B_all_range, :) = repmat(B_sitA, [(NBT), 1])  ;
        SI_GT(B_all_range, :) = repmat(B_sitI, [(NBT), 1])  ;

    case 6

        TA_GT(A_MAJ_range,1) = 3;
        TA_GT(A_MAJ_range,2) = 1;

        TA_GT(A_min_range,1) = 3;
        TA_GT(A_min_range,2) = 1;

        TA_GT(B_MAJ_range,1) = 3;
        TA_GT(B_MAJ_range,2) = 6;

        TA_GT(B_min_range,1) = 3;
        TA_GT(B_min_range,2) = 6;

        SA_GT(A_all_range, :) =  repmat(A_sitA, [(NAT), 1]);
        SI_GT(A_all_range, :) =  repmat(A_sitI, [(NAT), 1]);

        SA_GT(B_all_range, :) = repmat(B_sitA, [(NBT), 1])  ;
        SI_GT(B_all_range, :) = repmat(B_sitI, [(NBT), 1])  ;

end

SAcyto = SA_GT;

SA_GT = logical(SA_GT);
SI_GT = logical(SI_GT);

M_GT = logical(zeros(Nmax, 1));  % Motility genotype
M_GT(1:NT) = 1;

Mot =  zeros(Nmax, 1);
Mot(1:NT) = 1;

for ct0 = 1:TypeN

    n_init = sum((Gtype == ct0) & (LiveMat));

    InitNum(ct0) = n_init;

end
