Why model?

Before building a part or a device in lab, it is probably a good idea to simulate it in a simulation package — such as SimBiology, Cello, or GenoCAD — to check if the part or device one intends to make would work the way one wants it to. Although simulation can't solve all the problems that one might encounter during the actual physical design, it can certainly reduce them, and help in getting a better understanding of how a design behaves under different parameters and conditions.

Our project has two biological ystems: one for detection of corbon monoxide, and another for detection for oxides of nitrogen. Below we show you how we modelled and simulated both these systems.

Modeling of the CO-sensing system

In our CO-sensing system, the CO binds to CooA (an activator). It activates the activator protein, and so the activated activator protein (CooA)2CO binds to promoter and initiates transcription.

The reactions which take place can be given as:

$$\textrm{[CO]}+\textrm{[CooA]} \longleftrightarrow \textrm{[(CooA)2CO]} \label{R1} \tag{R1}$$

The equation above represents the first reaction that takes place in our system in which the CO present in the environment is taken up into the cell where it binds to the CooA protein. The CooA protein is constantly expressed in our system using a constitutive promoter. The coefficient 2 represents the number of CO which bind to CooA; this number was taken from the literature.

To develop a dissociation constant of this equation, we can write it as:

$$KI=\frac{\textrm{[CO].[CooA]}}{\textrm{[(CooA)2CO]}} \label{1} \tag{1}$$

Note that the concentrations of $\textrm{[CO]}$, $\textrm{[CooA]}$, and $\textrm{[(CooA)2CO]}$ are at equilibrium before the start of the reaction R1.

The second reaction that takes place in our system is the binding of the activated activator protein complex of $\textrm{[(CooA)2CO]}$ to the promoter Pr.

It can be described as:

$$\textrm{[(CooA)2CO]} \longleftrightarrow \textrm{[(CooA)2CO.Pr]} \label{R2} \tag{R2}$$

Here, $\textrm{[(CooA)2CO.Pr]}$ represents the activated protein $\textrm{[(CooA)2CO]}$ bound to the promoter.

To develop a dissociation constant of this equation, we can write as:

$$KP=\frac{\textrm{[(CooA)2CO]}}{\textrm{[(CooA)2CO.Pr]}} \label{2} \tag{2}$$

Now, in order to get an equation which describes the change in the concentration of mRNA with time we develop the following differential equation:

$$\frac{d(mRNA)}{dt}=k_{on1}.P_{exp}-\gamma_1.[mRNA] \label{3} \tag{3}$$

where, $k_{on1}$ is represents the formation rate constant of the mRNA, and $\gamma_1$ represents the degradation rate constant of mRNA.

In Eq. (\ref{3}) it can be seen that the change in the $\textrm{[mRNA]}$ with time is equal to the difference of the formation rate of mRNA and its degradation rate.

The initial conditions for the Eq. (\ref{3}) are assumed to be:

$$\textrm{[mRNA(t=0)] = 0 molar}$$

This means that at the time $t=0$ sec (initial time), there are no molecules of mRNA encoding for our desired chromoprotein in the cell.

In Eq. (\ref{3}), $P_{exp}$ represents the probability that the promoter will be activated by the activated activator protein complex $\textrm{[(CooA)2CO]}$. It can be calculated by the following equation:

$$P_{exp} = \frac{\textrm{number of situations in which the promoter will be activated}}{\textrm{total number of situations}} = \frac{w}{z} \label{4} \tag{4}$$

In our modeling of the CO-sensing system, we have assumed that at any moment only one molecule of the activated activator protein complex $\textrm{[(CooA)2CO]}$ can bind to the promoter, and that the \ref{R1} reaches equilibrium before the transcription of the system begins.

To develop our equations further, we say that the number of situations in which a promoter will be activated is equal to:

$$ w=\textrm{[(CooA)2CO.Pr]}=\frac{\textrm{[(CooA)2CO]}}{K_p} \label{5} \tag{5}$$

Eq. (\ref{5}) demonstrates the relation between the binding of the activated activator protein complex $\textrm{[(CooA)2CO]}$, showing that it is equal to the amount of $\textrm{[(CooA)2CO]}$ upon the dissociation constant of the binding of $\textrm{[(CooA)2CO]}$ to the promoter.

To simplify this equation, we assume that the number of situations in which the promoter will not be activated by $\textrm{[(CooA)2CO]}$ is equal to $1$.

Note that this assumption is made because the only situation in which transcription will not take place is when there are no molecules bound to the promoter.

Substituting the assumptions in Eq. (\ref{5}), we get:

$$ z=w+1=\frac{\textrm{[(CooA)2CO]}}{K_p}+1 \label{6} \tag{6}$$

where, $z$ represents the total number of situations, $w$ represents the number of situations in which the promoter will be activated, and $1$ represents the number of situations in which the promoter will not be activated.

Sustituting Eq. (\ref{5}) and (\ref{6}) into Eq. (\ref{4}), we get:

$$P_{exp} = \frac{w}{z} = \frac{\frac{\textrm{[(CooA)2CO]}}{K_p}}{\frac{\textrm{[(CooA)2CO]}}{K_p}+1} \label{7} \tag{7}$$

As we want our $P_{exp}$ to be expressed in the terms of $\textrm{[CO]}$ and $\textrm{[CooA]}$ separately instead of $\textrm{[(CooA)2CO]}$, we repeat these above equations by taking into consideration \ref{R1} and its dissociation constant $\textrm{KI}$.

By doing so, we get:

$$P_{exp} = \frac{\frac{\textrm{[CooA]}.\frac{\frac{\textrm{CO}}{K_I}}{1+\frac{\textrm{CO}}{K_I}}}{K_p}} {1+\frac{\textrm{[CooA]}.\frac{\frac{\textrm{CO}}{K_I}}{1+\frac{\textrm{CO}}{K_I}}}{K_p}} \label{8} \tag{8}$$

As the Eq. (\ref{8}) demonstrates, $P_{exp}$ depends only upon the concentration of $\textrm{[CO]}$ and $\textrm{[CooA]}$.

Recall Eq. (\ref{3})

$$\frac{d(\textrm{mRNA})}{dt}=k_{on1}.P_{exp}-\gamma_1.[\textrm{mRNA}]$$

The initial condition for this equation was $\textrm{[mRNA(t=0)] = 0 molar}$.

The analytical solution of this differential equation is equal to:

$$[mRNA]=\frac{k_{on1}.P_{exp}}{\gamma_1}(1-e^{-\gamma_1t})$$

Now that we have got the equation describing the change in concentration of mRNA, we can use the same method to derive an equation which describes the change in the concentration of the chromoprotein in time as following:

$$\frac{d\textrm{[CP]}}{dt} = k_{on2}.\textrm{[mRNA]}-\gamma_2\textrm{[CP]} \label{9} \tag{9}$$

Here, $k_{on2}$ is the formation rate constant of the chromoprotein and $\gamma_2$ is the degradation rate constant of the chromoprotein.

Again, the initial conditions for this equation are assumed to be:

$$\textrm{[CP(t=0)] = 0 molar}$$

which states that at the initial time, the amount of chromoprotein in the cell is equal to $0$.

Here on forth, we will repeat the steps from Eq. (\ref{4}) to Eq. (\ref{8}) and arrive to the analytical solution of the differential equation of Eq. (\ref{9}), which is:

$$\textrm{[CP]} = k_{on2}.\frac{k_{on1}.P_{exp}}{\gamma_1.\gamma_2}(1-e^{-\gamma_2t}) + k_{on2}.\frac{k_{on1}.P_{exp}}{\gamma_1(\gamma_1-\gamma_2)}(e^{-\gamma_2t}-e^{-\gamma_1t}) \label{10} \tag{10}$$

Eq. (\ref{10}) is the final equation which describes the amount of $\textrm{[CP]}$ along with the amount of $\textrm{[mRNA]}$.

Assumptions

1. 1. Reactions \ref{R1} and \ref{R2} reach equilibrium before the transcription begins.
2. 2. Only one molecule of the complex $\textrm{[(CooA)2CO]}$ can be bound to the promoter at a time.
3. 3. $\textrm{[mRNA(t=0)] = 0 molar}$
4. 4. $\textrm{[CO]}$ and $\textrm{[CooA]}$ concentrations inside the cell are constant.

Simulation of the CO-sensing system

Simulation parameters

Table 1: Parameters, their description, and their values used in MATLAB simulation of CO-sensing system.

Parameter	Description	Value
$\textrm{[CO]}$	Carbon monoxide concentration at equilibrium	$\textrm{10}$
$\textrm{[CooA]}$	CooA protein concentration at equilibrium	$11$
$\textrm{KI}$	Dissociation constant of the binding of CO to CooA	$\textrm{0.011 mM}$
$\textrm{KP}$	Dissociation Constant of the binding of $\textrm{(CO2CooA)}$ to the promoter	$\textrm{0.000009}$ $\textrm{min}^{-1}$
Kon1	Formation rate constant of mRNA	$\textrm{0.001}$ $\textrm{min}^{-1}$
Kon2	Formation rate constant of the chromoprotein	0.006 $\textrm{min}^{-1}$ / 0.0051 $\textrm{min}^{-1}$
Gama1 ($\gamma$)	Degradation rate constant of the mRNA	$\textrm{0.00039}$ $\textrm{min}^{-1}$
Gama2 ($\gamma$)	Degradation rate constant of the chromoprotein	$\textrm{0.001}$ $\textrm{min}^{-1}$

Simulation results

MATLAB Code

function dx = ODEfun_CooM(t,x)
    dx      = zeros(1,2);
    %x is a vector (x(1),x(2))
    %x(1) = [mRNA]
    %x(2) = [Chromoprotein]
    KI      = 0.011;    %Dissociation constant of the binding of CO to CooA
    KP      = 0.1;      %Dissociation constant of the binding of 
                        %(CO2CooA) to the promoter.
    Kon1    = 0.001;    %Formation rate constant of mRNA
    Kon2    = 0.0051;   %Formation rate constant of the CP
    gama1   = 0.00039;  %Degradation rate constant of mRNA
    gama2   = 0.002;    %Degradation rate constant of the CP
    co      = 10;       %Concentration of CO
    cooa    = 11;       %Concentration of CooA
    w       = (cooa*co/(KI*(1+co/KI)))/KP;
    z       = (cooa*co/(KI*(1+co/KI)))/KP+1;
    P_exp   = w/z;
    dx1     = (Kon1*P_exp)-(gama1*x(1));    %Rate equation of mRNA
    dx2     = (Kon2*x(1))-(gama2*x(2));     %Rate equation of CP
    dx      = [dx1;dx2];
end

%Save this code as ODEfun_CooM.m file

function dx = ODEfun_CooF(t,x)
    dx      = zeros(1,2);
    %x is a vector [x(1),x(2)]
    %x(1) = [mRNA]
    %x(2) = [chromoprotein]
    KI      = 0.011;    %Dissociation constant of the binding of CO to CooA
    KP      = 0.1;      %Dissociation constant of the binding of 
                        %(CO2CooA) to the promoter.
    Kon1    = 0.001;    %Formation rate constant of mRNA
    Kon2    = 0.006;    %Formation rate constant of the CP
    gama1   = 0.00039;  %Degradation rate constant of mRNA
    gama2   = 0.002;    %Degradation rate constant of the CP
    co      = 10;       %Concentration of CO
    cooa    = 11;       %Concentration of CooA
    w       = (cooa*co/(KI*(1+co/KI)))/KP;
    z       = (cooa*co/(KI*(1+co/KI)))/KP+1;
    P_exp   = w/z;
    dx1     = (Kon1*P_exp)-(gama1*x(1));    %Rate equation of mRNA
    dx2     = (Kon2*x(1))-(gama2*x(2));     %Rate equation of CP
    dx      = [dx1;dx2];
end

%Save this code as ODEfun_CooF.m file

% Clear worskpace
clc 
clear all
% Code for calling the previous function:
time_range      = 18000; 
[T1,X1]         = ode45(@ODEfun_CooF,[0 time_range],[0 0]); 
[T2,X2]         = ode45(@ODEfun_CooM,[0 time_range],[0 0]); 
figure 
plot(T1,X1(:,2),'r', T2,X2(:,2),'b', T1,X1(:,1),'g')
legend('Conc. of chromoprotein with pCooF promoter ',...
    'Conc. of chromoprotein with pCooM promoter',...
    'mRNA concentration') 
xlabel('Time [sec]'); 
ylabel('Concentration [molar]'); 
title('Chromoprotein concentration with different promoters (pCooF, pCooM)')

%Save this code as co_sensing_system_simulation.m file

Modeling of the NOx-sensing system

To gain an insight into the stability of NO-sensing circuit, we consider the cooperativity of NO binding with NsrR and the law of mass action. The goal of the modelling is to determine the stable expression of chromoprotein by considering a specified concentration of different parameters of the dynamic system.

The following ordinary differential equations are used for determining the expression levels of chromoproteins (CP) and its mRNA (mCP):

$$\frac{d\textrm{[mCP]}}{dt} = k_{transcription}\textrm{[PyeaR]} - \gamma_m \textrm{[mCP]} \label{11} \tag{11}$$ $$\frac{d\textrm{[CP]}}{dt}=k_{translation}\textrm{[mCP]}-\gamma_p\textrm{[CP]} \label{12} \tag{12}$$

Moreover, the cooperativity of $\textrm{NO}$ binding with NsrR can be demonstrated by using the Hill Equation, where free $\textrm{NO}^*$ combines with the $\textrm{n}$ amount of NsrR, leading to the formation of $\textrm{n[NsrR][NO]}$ complex. The total amount of inducer $NO_T$ concentration can be described by the following equation:

$$\textrm{NO}_T = \textrm{n[NsrR][NO]} + [\textrm{NO}^*]$$

Additionally, $\textrm{NsrR}$ dimerization with respect to unbound $\textrm{NsrR}$ can be represented by the following equation by considering the degree of cooperativity of $\textrm{NO}$ binding:

$$\frac{[\textrm{NO}^*]}{[\textrm{NO}_T]} = \frac{1}{1+\left(\frac{NsrR}{K_{NO}} \right)^n} \label{13} \tag{13}$$

The correlation between PyeaR activity and rate of transcription of mCP (\beta) can represented as:

$$\textrm{[PyeaR]}activity = \frac{\beta}{1+ \left( \frac{\textrm{NO}^*}{K_{d(\textrm{NsrR})}} \right)} \label{14} \tag{14}$$

Combining the Eqs. (\ref{13}) and (\ref{14}) into equation (\ref{11}), yields:

$$\frac{d[\textrm{mCP}]}{dt} = \frac{k_{transcription}[\textrm{PyeaR}]}{1 + \left(\frac{[\textrm{NsrR}]}{\left(1+\left[\frac{[\textrm{NO}]}{K_{NO}}\right]^n \right)}\right) k_{d(\textrm{NsrR})} } - \gamma_m [\textrm{mCP}]$$

In order to further simplify the system, let

$$K=\frac{k_{transcription}[\textrm{PyeaR}]}{1 + \left(\frac{[\textrm{NsrR}]}{\left(1+\left[\frac{[\textrm{NO}]}{K_{NO}}\right]^n \right)}\right) k_{d(\textrm{NsrR})} }$$

Additionally, the value of $\textrm{K}$ is inversely proportional to the amount of $\textrm{NsrR}$ and directly proportional to the $\textrm{NO}$ concentration in the system. Substituting $\textrm{K}$ into the differential equation of $\textrm{mCP}$ generates the following equation:

$$\frac{d[\textrm{mCP}]}{dt} = K - \gamma_m [\textrm{mCP}]$$

After solving the equation through homogenous and particular solution, the following equation is obtained for steady—state concentration of mCP generated after PyeaR is triggered:

$$[\textrm{mCP}]=\frac{K}{\gamma_m} - \frac{K}{\gamma_m}e^{-\gamma_mt}$$

The steady—state concentration of mCP is then substituted into Eq. (\ref{12}) and is also solved through homogenous and particular solution. The resulting equation gives the steady—state concentration of CP in the desired system according to the specified formation and degradation rates of chromoprotein in response to promotor stimulation:

$$[\textrm{CP}]=\frac{K k_{translation}}{\gamma_p \gamma_m} + \frac{K k_{translation}}{{\gamma_m}^2-\gamma_p \gamma_m} e^{-\gamma_mt} - \left(\frac{K k_{translation}}{\gamma_p \gamma_m} + \frac{K k_{translation}}{{\gamma_m}^2-\gamma_p \gamma_m} \right)e^{-\gamma_pt} $$

Simulation of the NOx-sensing system

Simulation parameters

Table 2: Parameters, their description, and their values used in MATLAB simulation of the NOx-sensing system.

Parameter	Description	Value
$[\textrm{PyeaR}]$	Concentration of PyeaR	1 copy/cell
$k_{transcription}$	Rate of chromoprotein mRNA synthesis	0.167$\textrm{min}^{-1}$
$k_{translation}$	Rate of chromoprotein synthesis	$0.0011 \textrm{min}^{-1}$
$\gamma_m$	mRNA degradation rate	$0.19 \textrm{min}^{-1}$
$\gamma_p$	Chromoprotein degradation rate	$0.18 \textrm{min}^{-1}$
$K_{\textrm{NO}}$	Dissociation constant of NO	$0.12 \textrm{min}^{-1}$
$K_{\textrm{d(NsrR)}}$	Dissociation constant of NsrR	$\textrm{0.035}$ mM
n	Cooperativity of NO binding to NsrR	2

Simulation results

Matlab Code

    
%Defining the parameters
ktr     = 0.167;     %mRNA synthesis rate of 
                     %Chromoprotein per minute (rate of transcription)
ktrans  = 0.0011;    %Chromoprotein synthesis rate 
                     %per minute (rate of translation)
ycp     = 0.18;      %Chromoproteins degradation rate per minute
ym      = 0.19;      %mRNA degradation rate of chromoproteins per minute 
Kno     = 0.12;      %Association constant for inducer (NO)
Kdn     = 0.035;     %Dissociation constant for repressor (NsrR)
pyear   = 1;         %Number of plasmids per cell
n       = 2;         %Cooperativity of inducer (NO)

%Assume that the repressor is removed 
NsrR    = 0;
NO      = 2;
t       = 0:0.1:240; 

%Calculate mRNA and Chromoprotein values
K       = ktr*pyear/(1+(NsrR/((1+(NO/Kno)^n)*Kdn)));
mCP     = K/ym-K/ym*exp(-ym*t);
CP      = K*ktrans/(ycp*ym)+K*ktrans/(ym^2-ycp*ym)*exp(-ym*t)...
    -(K*ktrans/(ycp*ym)+K*ktrans/(ym^2-ycp*ym))*exp(-ycp*t);

disp('Steady State mCP Concentration:')
disp(mCP(end)) %Steady state mRNA value
disp('Steady State CP Concentration:')
disp(CP(end)) %Steady state chromoprotein value
 
%Plotting mRNA concentration over time
figure(1)
plot(t, mCP,'r','LineWidth',2)
hold on
grid on
xlabel('Time (minutes)')
ylabel('Number of mCP per cell')
title('mCP Expression with \gamma_m=0.19min^-1')

%Plotting the chromoprotein concentrations over time
figure(2)
plot(t,CP,'r', 'LineWidth',2)
grid on
hold on
xlabel('Time (minutes)')
ylabel('Number of chromoproteins molecules per cell')
title('Chromoprotein Expression with \gamma_p=0.18 min^-1')

%Save this as nox_sensing_system_simulation.m file
    

References

1. Li, Juan et al. "A comprehensive kinetic mechanism for CO, CH2O, and CH3OH combustion." International Journal of Chemical Kinetics 39.3 (2007) : 109–136.

2. Chen, H., Shiroguchi, K., Ge, H., & Xie, X. S. (2015). Genome–wide study of mRNA degradation and transcript elongation in Escherichia coli. Molecular Systems Biology, 11(1), 781. doi:10.15252/msb.20145794

3. Bergmann, John E, and Harvey F Lodish. "A kinetic model of protein synthesis. Application to hemoglobin synthesis and translational control." The Journal of biological chemistry 254.23 (1979) : 11927–37.

4. Leduc, Jason et al. "Mapping CooA.RNA polymerase interactions: Identification of activating regions 2 and 3 in CooA, the co–sensing transcriptional activator." Journal of Biological Chemistry 276.43 (2001) : 39968–39973.

5. Nick Csicsery, and Ricky O'Laughlin. "A Mathematical Model of a Synthetically Constructed Genetic Toggle Switch" (2015)