# Team:XJTU-China/Model

Team:XJTU-China/Model

Modelling

# Mathematical Modelling

## Protein Modelling

To predict the experimental result of our project, we constructed a set of mathematical models in five different aspects. Here is the showcase of our math. Additionally, to investigate the mechanism of allosteric inhibition on AroG by Phe and the alleviation in S211F mutant, it is proposed to be quantified and visualized using PyMOL, Gaussian16.0W, GaussView6.0, Swiss, AutoDockTools software. To see the result of our protein model, plese check our wiki of protein modelling.

## Summary

Our modeling includes five steps:

• Establish the model of population dynamics, which displays the population change of E. coli;
• Establish the model of toggle switch, where the production of red fluorescent protein (RFP) and green fluorescent protein (GFP) shows the effect of toggle switch;
• Establish the model of genetic circuits based on the model of toggle switch;
• Establish the model of synthesis of tryptophan based on Michaelis-Menten equation;
• Finally, integrate the above models to establish the model of production.

## Establishment of model

Due to the size limit of iGEM wiki, please click the following buttons to view our models establishment!

## Result and conclusion

### The population density of E. coli

Here, $K=6.08×{10}^{9}\mathrm{CFU}/\mathrm{ml}$ $K=6.08×10^9\mathrm{CFU/ml}$. The figure below shows the population density of E. coli.

• When $N<\frac{K}{2}$ $N<\dfrac{K}{2}$, the population density grows exponentially. And when $N>\frac{K}{2}$ $N>\dfrac{K}{2}$, the environmental resources have a restrictive effect on E. coli. Finally the population density approaches $K$ $K$;
• The population density reach balance at about $33\mathrm{h}$ $33\mathrm{h}$.

### The effect of toggle switch

When $t=1000\mathrm{min}$ $t=1000\mathrm{min}$, add IPTG. When $t=2000\mathrm{min}$ $t=2000\mathrm{min}$, remove IPTG and raise temperature. The results are shown in the figure below.

• The change rates of $\left[rlacI\right]$ $[rlacI]$ and $\left[rGFP\right]$ $[rGFP]$ decrease as $\left[pcI857\right]$ $[pcI857]$ increases. The change rates of $\left[rcI857\right]$ $[rcI857]$ and $\left[rRFP\right]$ $[rRFP]$ decrease as $\left[placI\right]$ $[placI]$ increases. $\left[IPTG\right]$ $[IPTG]$ restrains the production of RFP when its concentration is low, and promotes the production when its concentration is high;
• At first, the concentration of GFP is more than the concentration of RFP, and green fluorescence appears. After adding IPTG, the concentration of RFP outnumbered GFP, and red fluorescence appears. After removing IPTG and raising temperature, the rank of RFP and GFP exchanged again, and green fluorescence appears.

### The product of genetic circuits

Add IPTG at the beginning, and when $t=1170\mathrm{min}$ $t=1170\mathrm{min}$, remove IPTG and raise temperature. The figure below shows the concentration change of gene product.

• There are two stable states during the period of time;
• After raising temperature at $1170\mathrm{min}$ $1170\mathrm{min}$, the concentrations of cl857 and pykA go down, while the concentrations of lacI, aroG, trpA and trpB go up.

### The output of tryptophan

Add IPTG at the beginning, and when $t=1170\mathrm{min}$ $t=1170\mathrm{min}$, remove IPTG and raise temperature. The figure below shows the concentration change during the synthesis of tryptophan.

• When reaction starts, Glc begin to convert to PEP, and PEP immediately turns into Pyr and DAHP;
• The concentration of DAHP reaches maximum at about $1900\mathrm{min}$ $1900\mathrm{min}$, and after that it goes down;
• The product of DAHP is 3IGP, and 3IGP immediately converts to Trp;
• The final products of reactions are Pyr and Trp, whose concentrations are stable at $3000\mathrm{min}$ $3000\mathrm{min}$.

### The best production strategy

Let the initial value of $\left[Glc\right]$ $[Glc]$ be $20$ $20$. Add IPTG at the beginning, and change the time of raising the temperature. Let ${t}_{2}$ $t_2$ be the time of raising temperature. The final outputs of $\left[Pyr\right]$ $[Pyr]$ and $\left[Trp\right]$ $[Trp]$ are shown in the figure below.

• When ${t}_{2}<1170\mathrm{min}$ $t_2<1170\mathrm{min}$, the output of tryptophan increases slowly as ${t}_{2}$ $t_2$ increases;
• When ${t}_{2}=1170\mathrm{min}$ $t_2=1170\mathrm{min}$, the maximun output of tryptophan is $15.43$ $15.43$ ( $77.15\mathrm{%}$ $77.15\%$ of $20$ $20$);
• When $1170\mathrm{min}<{t}_{2}<2000\mathrm{min}$ $1170\mathrm{min}, the output of tryptophan drops sharply as ${t}_{2}$ $t_2$ increases;
• Finally, when ${t}_{2}>2000\mathrm{min}$ $t_2>2000\mathrm{min}$, the output of tryptophan is around $8$ $8$ ( $40\mathrm{%}$ $40\%$ of $20$ $20$).

## Reference

[1] Verhulst, P.-F. "Recherches mathématiques sur la loi d'accroissement de la population." Nouv. mém. de l'Academie Royale des Sci. et Belles-Lettres de Bruxelles 18, 1-41, 1845.

[2] Verhulst, P.-F. "Deuxième mémoire sur la loi d'accroissement de la population." Mém. de l'Academie Royale des Sci., des Lettres et des Beaux-Arts de Belgique 20, 1-32, 1847.

[3] XIAN YIN, HYUN-DONG SHIN, et al. 2017. P gas, a Low-pH-Induced Promoter, as a Tool for Dynamic Control of Gene Expression for Metabolic Engineering of Aspergillus niger. Appl Environ Microbiol. [J/OL], 2;83(6):e03222-16.

# Code

Using MATLAB R2019b.

## The model of population dynamics

main.m

% initialization

clc;
close all;
clear all;

% numerical solution
[t, N] = ode45(@(t, N) odefcn(t, N)', [0 3000], 6.08*10^9*0.0001);

% draw figure
figure();
plot(t, N, 'LineWidth', 2);

xlabel('Time/min');

ylabel('Population Density')
grid on;



odefun.m

xxxxxxxxxx

function dN = odefcn(t, N)

r = 0.0073;
K = 6.08*10^9;
dN = r*N*(1-N/K);



## The model of toggle switch

main.m

xxxxxxxxxx

% initialization

clc;
close all;
clear all;

% initial value
tspan = [0 3000];
c0 = [0.1 0.1 0.1 0.1 0.01 0.01 0.01 0.01];

t1 = 1000;
iptg1 = 10;

% raise temperature
t2 = 2000;
alpha = 1.2;

% numerical solution
[t, c] = ode45(@(t, c) odefcn(t, c, t1, iptg1, t2, alpha)', tspan, c0);

% draw figure
figure();
plot(t, c(:, 6), 'g', t, c(:, 8), 'r', 'LineWidth', 2);
legend('[pGFP]', '[pRFP]');
xlabel('Time/s');
ylabel('Concentration');
grid on;



odefun.m

xxxxxxxxxx

function dc = odefcn(t, c, t1, iptg1, t2, alpha)

% lambdaPR: lacI, GFP
lambdaPR = 0.5;
alphal0 = 0.5;
alphal1 = 1;

% Plac: cl857, RFP
Plac = 1.3;
alphap0 = 0.5;
alphap1 = 1;

% RFP
ksynr = 0.019;
kder = 0.0013*10;
kpsynr = 0.47;
kpder = 0.136*10;

% GFP
kg = 0.877;
ksyng = kg*ksynr;
kdeg = kg*kder;
kpsyng = kg*kpsynr;
kpdeg = kg*kpder;

% cI857
kc = 0.95;
ksync = kc*ksynr;
kdec = kc*kder;
kpsync = kc*kpsynr;
kpdec = kc*kpder;

% lacI
kl = 0.626;
ksynl = kl*ksynr;
kdel = kl*kder;
kpsynl = kl*kpsynr;
kpdel = kl*kpder;

iptg = 0;
if (t>t1)
iptg = iptg1;
end

% raise temperature
if (t>t2)
iptg = 0;
alphal1 = alpha;
end

% differential equations
dc(1) = ksynl*alphal0*lambdaPR + ksynl*alphal1*(1/(c(7)))*lambdaPR - kdel*c(1); % [rlacI]
dc(2) = ksyng*alphal0*lambdaPR + ksyng*alphal1*(1/(c(7)))*lambdaPR - kdeg*c(2); % [rGFP]
dc(3) = ksync*alphap0*Plac + ksync*alphap1*(iptg/(iptg+c(5)))*Plac - kdec*c(3); % [rcI857]
dc(4) = ksynr*alphap0*Plac + ksynr*alphap1*(iptg/(iptg+c(5)))*Plac - kder*c(4); % [rRFP]
dc(5) = kpsynl*c(1) - kpdel*c(5); % [placI]
dc(6) = kpsyng*c(2) - kpdeg*c(6); % [pGFP]
dc(7) = kpsync*c(3) - kpdec*c(7); % [pcI857]
dc(8) = kpsynr*c(4) - kpder*c(8); % [pRFP]



## The model of genetic circuits

main.m

xxxxxxxxxx

% initialization

clc;
close all;
clear all;

% initial value
tspan = [0 3000];
c0 = [0.1 0.1 0.1 0.1 0.1 0.1 0.01 0.01 0.01 0.01 0.01 0.01];

t1 = 1000;
iptg1 = 10;

% raise temperature
t2 = 2000;
alpha = 1.2;

% numerical solution
[t, c] = ode45(@(t, c) odefcn(t, c, t1, iptg1, t2, alpha)', tspan, c0);

% draw figure
figure();
plot(t, c(:, 7), t, c(:, 8), t, c(:, 9), t, c(:, 10), t, c(:, 11), t, c(:, 12), 'LineWidth', 2)
legend('[placI]', '[paroG]', '[ptrpB]', '[ptrpA]', '[pcI857]', '[ppykA]');
grid on;



odefun.m

xxxxxxxxxx

function dc = odefcn(t, c, t1, iptg1, t2, alpha)

% lambdaPR: lacI, aroG, trpB, trpA
lambdaPR = 0.5;
alphal0 = 0.5;
alphal1 = 1;

% Plac: cl857, pykA
Plac = 1.3;
alphap0 = 0.5;
alphap1 = 1;

% RFP
ksynr = 0.019;
kder = 0.0013*10;
kpsynr = 0.47;
kpder = 0.136*10;

% lacI
kl = 0.626;
ksynl = kl*ksynr;
kdel = kl*kder;
kpsynl = kl*kpsynr;
kpdel = kl*kpder;

% aroG
ka = 0.644;
ksyna = ka*ksynr;
kdea = ka*kder;
kpsyna = ka*kpsynr;
kpdea = ka*kpder;

% trpB
kB = 0.568;
ksynB = kB*ksynr;
kdeB = kB*kder;
kpsynB = kB*kpsynr;
kpdeB = kB*kpder;

% trpA
kA = 0.84;
ksynA = kA*ksynr;
kdeA = kA*kder;
kpsynA = kA*kpsynr;
kpdeA = kA*kpder;

% cI857
kc = 0.95;
ksync = kc*ksynr;
kdec = kc*kder;
kpsync = kc*kpsynr;
kpdec = kc*kpder;

% pykA
kp = 0.47;
ksynp = kp*ksynr;
kdep = kp*kder;
kpsynp = kp*kpsynr;
kpdep = kp*kpder;

iptg = 0;
if (t>t1)
iptg = iptg1;
end

% raise temperature
if (t>t2)
iptg = 0;
alphal1 = alpha;
end

% differential equations
dc(1) = ksynl*alphal0*lambdaPR + ksynl*alphal1*(1/(c(11)))*lambdaPR - kdel*c(1); % [rlacI]
dc(2) = ksyna*alphal0*lambdaPR + ksyna*alphal1*(1/(c(11)))*lambdaPR - kdea*c(2); % [raroG]
dc(3) = ksynB*alphal0*lambdaPR + ksynB*alphal1*(1/(c(11)))*lambdaPR - kdeB*c(3); % [rtrpB]
dc(4) = ksynA*alphal0*lambdaPR + ksynA*alphal1*(1/(c(11)))*lambdaPR - kdeB*c(4); % [rtrpA]
dc(5) = ksync*alphap0*Plac + ksync*alphap1*(iptg/(iptg+c(7)))*Plac - kdec*c(5); % [rcI857]
dc(6) = ksynp*alphap0*Plac + ksynp*alphap1*(iptg/(iptg+c(7)))*Plac - kder*c(6); % [rpykA]
dc(7) = kpsynl*c(1) - kpdel*c(7); % [placI]
dc(8) = kpsyna*c(2) - kpdea*c(8); % [paroG]
dc(9) = kpsynB*c(3) - kpdeB*c(9); % [ptrpB]
dc(10) = kpsynA*c(4) - kpdeA*c(10); % [ptrpA]
dc(11) = kpsync*c(5) - kpdec*c(11); % [pcI857]
dc(12) = kpsynp*c(6) - kpdep*c(12); % [ppykA]



## The model of synthesis of tryptophan

main.m

xxxxxxxxxx

% initialization

clc;
close all;

% initial value
tspan = [0 3000];
c0 = [0 0 0 0 0 20];
aroG = 1;
pykA = 0.15;
trpAB = 2.5;

% numerical solution
[t, c] = ode45(@(t, c) odefcn(t, c, aroG, pykA, trpAB)', tspan, c0);

% draw figure
figure();
plot(t, c(:, 6), t, c(:, 1), t, c(:, 2), t, c(:, 3), t, c(:, 4), t, c(:, 5), 'LineWidth', 2)
legend('[Glc]',  '[PEP]', '[DAHP]', '[Pyr]', '[3IGP]', '[Trp]');
xlabel('Time/min');
ylabel('Concetration');
grid on;
axis([0 3000 0 20]);



odefun.m

xxxxxxxxxx

function dc = odefcn(t, c, aroG, pykA, trpAB)

% reaction 1
vmax1 = 0.1;
km1 = 0.005;

% reaction 2
vmax2 = 0.03;
km2 = 0.002;

% reaction aroG
kcataroG = 0.1;
kmaroG = 0.009;

% reaction pykA
kcatpykA = 0.15;
kmpykA = 0.05;

% reaction trpAB
kcattrpAB = 0.009;
kmtrpAB = 0.1;

% differential equations
dc(1) = vmax1*c(6)/(c(6)+km1)-kcataroG*aroG*c(1)/(c(1)+kmaroG)-kcatpykA*pykA*c(1)/(c(1)+kmpykA); % [PEP]
dc(2) = kcataroG*aroG*c(1)/(c(1)+kmaroG)-vmax2*c(2)/(c(2)+km2); % [DAHP]
dc(3) = kcatpykA*pykA*c(1)/(c(1)+kmpykA); % [Pyr]
dc(4) = vmax2*c(2)/(c(2)+km2)-kcattrpAB*trpAB*c(4)/(c(4)+kmtrpAB); % [3IGP]
dc(5) = kcattrpAB*trpAB*c(4)/(c(4)+kmtrpAB); % [Trp]
dc(6) = -vmax1*c(6)/(c(6)+km1); % [Glc]



## The model of production

main.m

xxxxxxxxxx

% initialization

clc;
close all;
clear all;

% initial value
tspan = [0 3000];
c0 = [0.1 0.1 0.1 0.1 0.1 0.1 0.01 0.01 0.01 0.01 0.01 0.01 0 0 0 0 0 20 6.08*10^9*0.001];

n = 10;
maxt = 3000;
mint = 0;
for i = 1:n
clc;
fprintf('%d/%d\n',i,n);
t2 = mint + (maxt-mint)*i/n;
[t, c] = ode45(@(t, c) odefcn(t, c, 0, 10, t2, 1.2)', tspan, c0);
result(i,1) = c(size(t,1),15);
result(i,2) = c(size(t,1),17);
end

% figure 1
tspan2 = mint+(maxt-mint)/n:(maxt-mint)/n:maxt;
plot(tspan2, result(:,1), tspan2, result(:,2), 'LineWidth', 2);
xlabel('Time of raising temperature/min');
ylabel('Output')
legend('[Pyr]', '[Trp]');
grid on;

% find the maximun
clc;
[m, i] = max(result(:,2));
t2 = mint + (maxt-mint)*i/n;
fprintf('When raising temperature at %dmin, the maximun output of Trp is %f. \n', t2, m);
[t, c] = ode45(@(t, c) odefcn(t, c, 0, 10, t2, 1.2)', tspan, c0);

% figure 2
figure();
plot(t, c(:, 19), 'LineWidth', 2)
legend('N');
xlabel('Time/min');
ylabel('Population Density');
grid on;

% figure 3
figure();
plot(t, c(:, 7), t, c(:, 8), t, c(:, 9), t, c(:, 10), t, c(:, 11), t, c(:, 12), 'LineWidth', 2)
legend('[placI]', '[paroG]', '[ptrpB]', '[ptrpA]', '[pcI857]', '[ppykA]');
xlabel('Time/min');
ylabel('Concetration');
grid on;

% figure 4
figure();
plot(t, c(:, 18), t, c(:, 13), t, c(:, 14), t, c(:, 15), t, c(:, 16), t, c(:, 17), 'LineWidth', 2)
legend('[Glc]', '[PEP]', '[DAHP]', '[Pyr]', '[3IGP]', '[Trp]');
xlabel('Time/min');
ylabel('Concetration');
axis([0 3000 0 20]);
grid on;



odefun.m

xxxxxxxxxx

function dc = odefcn(t, c, t1, iptg1, t2, alpha)

% lambdaPR: lacI, aroG, trpB, trpA
lambdaPR = 0.5;
alphal0 = 0.5;
alphal1 = 1;

% Plac: cl857, pykA
Plac = 1.3;
alphap0 = 0.5;
alphap1 = 1;

% RFP
ksynr = 0.019;
kder = 0.0013*10;
kpsynr = 0.47;
kpder = 0.136*10;

% lacI
kl = 0.626;
ksynl = kl*ksynr;
kdel = kl*kder;
kpsynl = kl*kpsynr;
kpdel = kl*kpder;

% aroG
ka = 0.644;
ksyna = ka*ksynr;
kdea = ka*kder;
kpsyna = ka*kpsynr;
kpdea = ka*kpder;

% trpB
kB = 0.568;
ksynB = kB*ksynr;
kdeB = kB*kder;
kpsynB = kB*kpsynr;
kpdeB = kB*kpder;

% trpA
kA = 0.84;
ksynA = kA*ksynr;
kdeA = kA*kder;
kpsynA = kA*kpsynr;
kpdeA = kA*kpder;

% cI857
kc = 0.95;
ksync = kc*ksynr;
kdec = kc*kder;
kpsync = kc*kpsynr;
kpdec = kc*kpder;

% pykA
kp = 0.47;
ksynp = kp*ksynr;
kdep = kp*kder;
kpsynp = kp*kpsynr;
kpdep = kp*kpder;

% reaction 1
vmax1 = 0.1;
km1 = 0.005;

% reaction 2
vmax2 = 0.03;
km2 = 0.002;

% reaction aroG
kcataroG = 0.1;
kmaroG = 0.009;

% reaction pykA
kcatpykA = 0.15;
kmpykA = 0.05;

% reaction trpAB
kcattrpAB = 0.009;
kmtrpAB = 0.1;

% population dynamics
K = 6.08*10^9;
r = (3.1364*10^(-04)*c(15)+0.0052);
k = 5*10^(-11);

iptg = 0;
if (t>t1)
iptg = iptg1;
end

% raise temperature
if (t>t2)
iptg = 0;
alphal1 = alpha;
end

% differential equations
dc(1) = (ksynl*alphal0*lambdaPR + ksynl*alphal1*(1/(c(11)))*lambdaPR - kdel*c(1)); % [rlacI]
dc(2) = (ksyna*alphal0*lambdaPR + ksyna*alphal1*(1/(c(11)))*lambdaPR - kdea*c(2)); % [raroG]
dc(3) = (ksynB*alphal0*lambdaPR + ksynB*alphal1*(1/(c(11)))*lambdaPR - kdeB*c(3)); % [rtrpB]
dc(4) = (ksynA*alphal0*lambdaPR + ksynA*alphal1*(1/(c(11)))*lambdaPR - kdeB*c(4)); % [rtrpA]
dc(5) = (ksync*alphap0*Plac + ksync*alphap1*(iptg/(iptg+c(7)))*Plac - kdec*c(5)); % [rcI857]
dc(6) = (ksynp*alphap0*Plac + ksynp*alphap1*(iptg/(iptg+c(7)))*Plac - kder*c(6)); % [rpykA]
dc(7) = (kpsynl*c(1) - kpdel*c(7)); % [placI]
dc(8) = (kpsyna*c(2) - kpdea*c(8)); % [paroG]
dc(9) = (kpsynB*c(3) - kpdeB*c(9)); % [ptrpB]
dc(10) = (kpsynA*c(4) - kpdeA*c(10)); % [ptrpA]
dc(11) = (kpsync*c(5) - kpdec*c(11)); % [pcI857]
dc(12) = (kpsynp*c(6) - kpdep*c(12)); % [ppykA]
dc(13) = k*c(19)*(vmax1*c(18)/(c(18)+km1)-kcataroG*c(8)*c(13)/(c(13)+kmaroG)-kcatpykA*c(12)*c(13)/(c(13)+kmpykA)); % [PEP]
dc(14) = k*c(19)*(kcataroG*c(8)*c(13)/(c(13)+kmaroG)-vmax2*c(14)/(c(14)+km2)); % [DAHP]
dc(15) = k*c(19)*(kcatpykA*c(12)*c(13)/(c(13)+kmpykA)); % [Pyr]
dc(16) = k*c(19)*(vmax2*c(14)/(c(14)+km2)-kcattrpAB*(c(3)+c(4))*c(16)/(c(16)+kmtrpAB)); % [3IGP]
dc(17) = k*c(19)*(kcattrpAB*(c(3)+c(4))*c(16)/(c(16)+kmtrpAB)); % [Trp]
dc(18) = k*c(19)*(-vmax1*c(18)/(c(18)+km1)); % [Glc]
dc(19) = r*c(19)*(1-c(19)/K); % N



For the labwork of our project, please check:
Engineering Success
Proof Of Concept