NT1 = 150;   % Starting cell number

NT2 = 150;

NT = NT1 + NT2;

Nmax = 500;   % Maximal cell number

HW = 10;  % Halfwidth of arena

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

Rg = 0.8;   % The radius of the initial distribution

% 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, 

if ~exist('MassComp')

    GeneCond = 2;   % Initial configuration
    % 1: Control: No difference in gene composition (Same traA and sitAI)
    % 2: Compatible traA with different sitAI
    % 3: Partially compaitble traA with different sitAI

end

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 = [0 0 0]; % Type 1 color
T2C = [0 1 1]; % Type 2 color
T3C = [0 0 0];   % Type 3 color
T4C = [1 1 1];

TypeN = 2;

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

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)' ;
NutVal = exp(- GradX.^2  - (GradY).^2);

NutVal = NutVal+0.05;

NutVal = 3*ones(size(NutVal));

NutValOri = NutVal;   % Original nutrient map
Nut_conv_r = 0.03;    % Convergence rate of the nutrient map to the 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
% 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);

Tcount1 = 1;

while (Tcount1 <= NT1)

    X1 = -Rg*rand(1);
    Y1 = 2*Rg*(rand(1)-0.5);

    DistR = sqrt(X1^2 + Y1^2);

    if DistR <= Rg

        PosCell( Tcount1 , :) = [X1, Y1];
        Tcount1 = Tcount1+1;

    end

end

Tcount2 = 1;

while (Tcount2 <= NT2)

    X2 = Rg*rand(1);
    Y2 = 2*Rg*(rand(1)-0.5);

    DistR = sqrt(X2^2 + Y2^2);

    if DistR <= Rg

        PosCell( NT1 +  Tcount2 , :) = [X2, Y2];
        Tcount2 = Tcount2+1;

    end

end


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
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_Neutral = [SitAIseg, SitAIseg, SitAIseg, SitAIseg, SitAIseg];
Config_A = [SitAIone, SitAIone, SitAIseg, SitAIzero, SitAIzero];
Config_B = [SitAIzero, SitAIzero, SitAIseg, SitAIone, SitAIone];

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


Range1 = 1:NT1;
Range2 = (NT1+1):(NT1+NT2);
RangeAll = 1:(NT1+NT2);

switch GeneCond
% 1: Control: No difference in gene composition (Same traA and sitAI)
% 2: Compatible traA with different sitAI
% 3: Partially compaitble traA with different sitAI

    case 1

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

        SA_GT(RangeAll, :) = repmat(Config_Neutral, [(NT1+NT2), 1]);

    case 2

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

        SA_GT(Range1, :) = repmat(Config_A, [NT1, 1]);
        SA_GT(Range2, :) = repmat(Config_B, [NT2, 1]);

    case 3

        TA_GT(Range1, 1) = 5;
        TA_GT(Range1, 2) = 1;
        TA_GT(Range2, 1) = 5;
        TA_GT(Range2, 2) = 6;

        SA_GT(Range1, :) = repmat(Config_A, [NT1, 1]);
        SA_GT(Range2, :) = repmat(Config_B, [NT2, 1]);

end

SAcyto = SA_GT;  % Cytoplasmic SitA 

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:4

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

    InitNum(ct0) = n_init;

end



