The model of population dynamics

First, we establish the model of population dynamics to study the variation of E. coli population density.

Theory

Let $N$ $N$ be the population density of E. coli. With the Logistic equation, we know that

where $r$ $r$ and $K$ $K$ are the growth rate and the environmental capacity of E. coli respectively. Equation $(1.1)$ $(1.1)$ shows that:

• When $N<{\displaystyle \frac{K}{2}}$ $N<\dfrac{K}{2}$, the population density grows exponentially;
• When $N>{\displaystyle \frac{K}{2}}$ $N>\dfrac{K}{2}$, the environmental resources have a restrictive effect on E. coli, and finally the population density approaches $K$ $K$.

Parameter

The parameters are shown in the table below.

Result

Let the initial value of population density be $0.01\mathrm{\%}$ $0.01\%$ of the environmental capacity. From the figure, we obtain that the population density reach balance after about $33\mathrm{h}$ $33\mathrm{h}$.

Code

main.m

​x
% 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

Next, to verify the effect of toggle switch, we design a genetic circuits as the figure below shows.

The promoter $\lambda P_R$ $\lambda P_{R}$ controls the production of lacI and green fluorescent protein (GFP), and can be promoted by raising temperature. The promoter $P_{lac}$ $P_{lac}$ controls the production of cI857 and red fluorescent protein (RFP), and can be promoted by IPTG. cl857 restricts the promoter $\lambda P_R$ $\lambda P_{R}$, while lacI restricts the promoter $P_{lac}$ $P_{lac}$.

Theory

For $x$ $x$, let $[dx]$ $[dx]$, $[rx]$ $[rx]$ and $[px]$ $[px]$ be the concentration of corresponding DNA, RNA and protein. And let $[IPTG]$ $[IPTG]$ be the concentration of IPTG.

For the transcription of promoter $\lambda P_R$ $\lambda P_R$, we have

where $k_{syn,x}$ $k_{syn,x}$ and $k_{de,x}$ $k_{de,x}$ are the synthesis rate and the degradation rate of $x$ $x$ respectively, ${\alpha }_{\lambda P_R,0}$ $\alpha_{\lambda P_R,0}$ is the background expression rate of $\lambda P_R$ $\lambda P_R$, and ${\alpha }_{\lambda P_R}$ $\alpha_{\lambda P_R}$ is the induced expression rate of $\lambda P_R$ $\lambda P_R$. Equation $(2.1)$ $(2.1)$ shows that the change rates of $[rlacI]$ $[rlacI]$ and $[rGFP]$ $[rGFP]$ decrease as $[pcI857]$ $[pcI857]$ increases.

For the transcription of promoter $P_{lac}$ $P_{lac}$, we have

where ${\alpha }_{P_{lac},0}$ $\alpha_{P_{lac},0}$ is the background expression rate of $P_{lac}$ $P_{lac}$, and ${\alpha }_{P_{lac}}$ $\alpha_{P_{lac}}$ is the induced expression rate of $P_{lac}$ $P_{lac}$. Equation $(2.2)$ $(2.2)$ shows that the change rates of $[rcI857]$ $[rcI857]$ and $[rRFP]$ $[rRFP]$ decrease as $[placI]$ $[placI]$ increases. And $[IPTG]$ $[IPTG]$ can raise the change rates.

For the translation of protein translation, we have

where $k{p}_{syn,x}$ $kp_{syn,x}$ and $k{p}_{de,x}$ $kp_{de,x}$ are the synthesis rate and the degradation rate of the protein of $x$ $x$ respectively.

Parameter

The parameters are shown in the table below.

Meanwhile, the correction is performed according to the CDS length.

Result

When $t=1000\mathrm{s}$ $t=1000\mathrm{s}$, add IPTG. When $t=2000\mathrm{s}$ $t=2000\mathrm{s}$, remove IPTG and raise temperature. The results are shown in the figure below. From the figure, we can see that there are three stable states.

• At first, the concentration of GFP is more than the concentration of RFP;
• After adding IPTG, the concentration of RFP outnumbered GFP;
• After removing IPTG and raising temperature, the rank of RFP and GFP exchanged again.

Code

main.m

​x
% 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];
        

      
        
% add IPTG
        
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;
            
    
            
    % add IPTG
            
    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]