NT1 = 100;   % Cells with florescent materials

NT2 = 100;   % Intermediary cell or compatible TraA to Type 1

NT3 = 100;   % Incompatible TraA to Type 1

NT = NT1 + NT2 + NT3;

Nmax = 500;   % Maximal cell number

HW = 3;  % 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 (Non motile cells can survive, nutrient is highly regenerative (high Nut_conv_r))
    % 2: Low (Non motile cells cannot survive)
    % 3: Nutrition source

    SR_Cond = 3;   % Serial transfer condition
    % 1: Type 1 and 2: Fully compatible
    % 2: Type 1 and 2: Partially compatible
    % 3: Type 1, 2 and 3: Type 2 is an intermediate strain

end

Flores = 1;   % Existence of florescent material (for MyxoSimul code)

% 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

FLR_GT = zeros(Nmax,1); % Florescent genotype
FLR_GT(1:NT1) = 1;
FLR_GT = logical(FLR_GT);

FLORES = zeros(Nmax,1);  % Florescent amount
FLORES(1:NT1) = 1;

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

EngMat =  zeros(Nmax,1);  % Energy storage
EngMat(1:(NT1+NT2)) = 1;

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

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

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

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

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 = 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);

PosCell(1:NT, :) = [ 2*(rand(NT,1)-0.5), 2*(rand(NT,1)-0.5) ];

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

p_SitA = 0.2;

tempSitA1 = rand(1, SitA_FN*SitA_sN);
tempSitA1 = (tempSitA1<p_SitA);

switch SR_Cond

    case 1   % Two types, compatible

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

        for z0 = 1:(NT1+NT2)
            SA_GT(z0,:) = tempSitA1 ;
            SAcyto(z0,:) = tempSitA1 ;
        end

        LiveMat(1:(NT1+NT2)) = 1;

    case 2   % Two types, partially compatible

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

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

        for z0 = 1:(NT1+NT2)
            SA_GT(z0,:) = tempSitA1 ;
            SAcyto(z0,:) = tempSitA1 ;
        end

        LiveMat(1:(NT1+NT2)) = 1;

    case 3  % Intermediary

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

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

        TA_GT((NT1+NT2+1):(NT1+NT2+NT3), 1) = 5;
        TA_GT((NT1+NT2+1):(NT1+NT2+NT3), 2) = 1;

        Gtype((NT1+NT2+1):(NT1+NT2+NT3)) = 3;

        for z0 = 1:(NT1+NT2+NT3)
            SA_GT(z0,:) = tempSitA1 ;
            SAcyto(z0,:) = tempSitA1 ;
        end

        EngMat(1:(NT1+NT2+NT3)) = 1;
        LiveMat(1:(NT1+NT2+NT3)) = 1;

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




