NT1 = 90;   % Majority traA cells

NT2 = 10;   % Minority traA cells

NFB = 10;   % Introduced aggressors

NT = NT1 + NT2;

Nmax = 500;   % Maximal cell number

HW = 5;  % Halfwidth of arena

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

if ~exist('MassComp')

    AggCond = 0;   % Condition of the initial spatial distribution
    % 0: No aggressors
    % 1: Introduced aggressors (added sitAI by HGT)

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

    GT_Cond = 3;   % Genotype condition
    % 1: Same TraA (Control)
    % 2: Similar TraA
    % 3: Different (incompatible) TraA

end

traA_nutri_fluc = 0;

% 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:(NT1+NT2)) = 1;

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

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

Gtype(1:NT1) = 1;
Gtype((NT1+1):NT) = 2;

T1C = [1 0 0]; % Type 1 color
T2C = [0 0 1]; % Type 2 color
T3C = [0 0 0];   % Type 3 color: Type 1's aggressor
T4C = [1 1 1];   % Type 4 color: Type 2's aggressor

TypeCmap = [T1C; T2C; T3C; T4C];
TypeN = 4;

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

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

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 = 1*ones( length(GradY), length(GradX)  );
        traA_nutri_fluc = 1;
        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);

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

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 = (NT2/NT)*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);

for Mcount = 1:NT1

    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(NT2, 1);
minYvals = minY1 + (minY2 - minY1)*rand(NT2, 1);

PosCell((NT1+1):(NT1+NT2), :) = [minXvals, minYvals];
PosCell((NT+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
SAcyto = zeros(Nmax, SitA_FN*SitA_sN);  % Cytoplasmic SitA

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

tempSitA1 = [ SitAIseg, SitAIzero, SitAIseg, SitAIzero, SitAIseg];   % Normal bacterial sitA type
tempSitAFB = [ SitAIseg, SitAIFB,  SitAIseg, SitAIFB, SitAIseg   ];  % Aggressor sitA type

switch GT_Cond

    case 1   % Same traA genotype

        TA_GT(1:NT,1) = 3;
        TA_GT(1:NT,2) = 1;

        p_SitA = 0.2;

        for z0 = 1:NT
            SA_GT(z0,:) = tempSitA1 ;
            SAcyto(z0,:) = tempSitA1 ;
        end

    case 2   % Similar traA genotype  (Same sitAI genotype)

        TA_GT(1:NT1,1) = 3;
        TA_GT(1:NT1,2) = 1;

        TA_GT((NT1+1):NT, 1) = 4;
        TA_GT((NT1+1):NT, 2) = 1;

        for z0 = 1:NT
            SA_GT(z0,:) = tempSitA1 ;
            SAcyto(z0,:) = tempSitA1 ;
        end

    case 3  % Different traA genotype  (Same sitAI genotype)

        TA_GT(1:NT1,1) = 2;
        TA_GT(1:NT1,2) = 1;

        TA_GT((NT1+1):NT, 1) = 5;
        TA_GT((NT1+1):NT, 2) = 4;


        for z0 = 1:NT
            SA_GT(z0,:) = tempSitA1 ;
            SAcyto(z0,:) = tempSitA1 ;
        end

end

if AggCond == 1

    FB_ID = randperm(NT, NFB );   % Aggressor IDs

    FB_ID1 = FB_ID(FB_ID <= NT1); % Type 1 -> FB
    FB_ID2 = FB_ID( (FB_ID > NT1)); % Type 2 -> FB
    Gtype(FB_ID1) = 3;
    Gtype(FB_ID2) = 4;

    for f0 = 1:NFB
        fbID = FB_ID(f0);
        SA_GT(fbID,:) = tempSitAFB ;
        SAcyto(fbID,:) = tempSitAFB ;

    end

end

SA_GT = logical(SA_GT);

SI_GT = SA_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




