Simplification

Although ODEs provide a thorough, precise description of the whole system, they contain too many equations and parameters which would act as a barrier for simulation and further analysis. A simplification of complicated ODEs is necessary. We simplify every single ODE according to certain appropriate assumptions. Finally, we came up with a set of DDE equations.

Assumptions we have made to deduct the original equations involve:

●Relatively faster reactions such as transcription reactions and binding reactions will reach to Quasi-equilibrium.

●Basal expression of protein is so meager that they can be ignored in modeling.

●Two series of Hill Kinetics equations can be estimated through a single Hill Kinetics equation.

●Protein which is translated from mRNA is proportional to corresponding mRNA in a previous time.

Consider gene lasR and luxR which are all constant genes, since AHL has a relatively low concentration compared to other substances, we assume that protein concentrations of PlasR and PluxR are irrelevant to AHLs. Therefore, concentrations of these proteins will reach to constant quantities. Moreover, certain mRNAs will also be constant. That is:

We have:

Then focus on compound LA1 and LA2. The binding process is far quicker than other reactions such as translation reactions, thus we assume that the compounds are Quasi-equilibrium:

We have:

Take the expressions to equation (6) (15), the functions are transformed to:

The transformed functions have nothing to do with feedback factors compared to original ones. Knowing that such feedback is essential for our system, we add additional feedback factors k_a,k_b manipulatively.

Then we want to come up a clear expression of A2_c1,A1_c2.They are generated by protein PluxI and PlasI, whose translation rate is much slower than the transcription rate. So we assume that certain mRNA is Quasi-equilibrium:

We have:

By ignoring the basal expression:

Knowing that protein PtetR is controlled by LA2 through Hill Function, we assume that concentration of mRNA is directly related to LA2:

Parameter KM4 is relevant to both KM1 and KM3, however accurate function depicting the relationship is unknown. We estimate the quantity of KM4 by multiplying KM1 and KM3. Noting LA1 and LA2 have already been expressed, we have:

Then we have trouble in expressing protein PluxI and PlasI. Since concentration of protein is the integrals of its mRNA, we assume that it is proportional to concentration of mRNA in a previous time. Thus we have a DDE function:

Equations concerning environmental AHLs remain unchanged.

The original ODEs can be transformed to a much simple DDEs by above deductions.

Parameters

Some of the parameters are derived from original ones, some of them are created to describe the new equations, and others are set for further testing.

Table 2 Parameters of DDEs

Results

General results

We coded the system in MATLAB. The result shows that every signal AHL is oscillating.

Figure 3 Signal AHL Oscillations

Stability analysis

In a nonlinear system, proper value ranges of the parameters are vital for producing a periodic oscillating solution. It’s not difficult to find the nonlinear part in our model, thus, we sought to analyze the sensitivity of Hill parameters and thoroughly reveal the internal relationship between Hill parameters and the system’s robustness. For simplicity, here we denoteKM1/ρ₁ , KM4/ρ₂ and k_TL2,k_TL4 by k_m1,k_m2and β₁, β₂ respectively. Without loss of generality, we determined to search the affection to the system’s robustness caused by the fluctuation of binary parameter (k_m2,β₂). After observing distinct oscillation mode (we choose two types of concentration of AHLs in environment as our observed object) and its phase trajectory, we depicted the bifurcation of our system.

Set basic parameters as follows:

k_a=0,k_b=0,μ=10min^-1,γ=2.3min^-1,β₁=30nm/min,k_m1=1,n₁=n₂=2,p=0.2,arab=20,t₁=t₂=0,τ₁=6,τ₂=10.

We changed the binary parameter (β₂，k_m2) and got the result below.

Figure 4 C12 and C6 when (β₂,k_m2)=(30,0.1)

Figure 5 C12 and C6 when (β₂,k_m2)=(120,0.1)

Figure 6 C12 and C6 when (β₂,k_m2)=(60,0.3)

After simulating the system at different parameters, we recorded several critical points for oscillation and made a table as follows.

Table 3 Critical points (β₂,k_m2) for oscillation

Depicting those critical points on an axis, we immediately got the bifurcation line of parameters (β₂,k_m2), which indicates the parameters’ value range when our system can oscillate stably, being marked in ‘bistable’.

Figure 7 Bifurcation Analysis on (β₂,k_m2)

Proportion of cell volume

In actual vivo experiment, volumes of two separate cells might not be the same. We simulate the system by changing the proportion of two types of cells. What we know from the simulation is that cell proportions only affect the amplitude of signal molecules’ oscillation, but no influence to the period or stability of oscillation. The following graphics are drawn under cell proportion 1 and 0.7.

Figure 8 oscillations under different cell proportions

We can find that under different cell proportions, our system can oscillate at different amplitude, and varying in cell proportions may also lead to the change of oscillation period.

Period adjustment

As figure 2 indicates, we expected to adjust oscillation period by adding signal molecule aTc into our system. A small molecule as aTc is, it can easily bind to protein TetR and quickly depletes TetR in cell 2, which results in reduction of TetR net production rate, and indirectly, the protein’s inhibition on promoter5 is crippled. Thus, it would take a longer time for our system to reach each threshold, which is equivalent to prolonging the time delay τ₂ in our simplified model. So we can deduce that changing the amount of aTc added into the system in precise model is equivalent to varying the time delay τ₂ in simplified model.

Simulation results under distinct τ₂ are presented as follows.

Figure 9 Oscillation cycle’s regulation

The result is exactly what we expected, which clearly demonstrates our system can truly be controlled by adding in external signal molecules.

Phase adjustment

In cell 2’s gene circuit, we designed a promoter induced by arabinose, marked by promoter 6 (see in figure 2). Promoter 6 is in a suppressed state until being induced by adding arabinose, and after the inhibition is relieved, signal molecule 3O12HSL will be generated extra. We can also analyze the differential equations describing dA1C2/dt, when adding arabinose during period from t₁ to t₂, dA1C2/dt contains an extra item arab?(t>t₁)*(t2). The parameter arab can reflect the rate of adding arabinose nonlinearly. Here we set arab=20，t₁=80，t₂=120 and simulation result is presented as follow.

Figure 10 Oscillation phase’s regulation