Simulation one

Module one

For each model in the system, the concentration of the extracellular carrier is set to be [Pi'], and the intracellular carrier concentration is [Pi]. The transfection efficiency by fluorescence identification is ni.We have formula (1) for them.

Module1 has the following conditions:

Wherein,[P1] represents the carrier concentration and NP0 is the corresponding concentration threshold.

Gene expression process

At first the gene expression process is modeled. Gene expression satisfies the central rule[1], and its expression model is shown in Figure 2.

In the picture, all parameter units are seconds and V is set to unit volume. For the model in the figure, there are six reaction channels. The single trend function and the state change vector are omitted here, and the chemical master equation of the corresponding model is obtained.

For the equilibrium state of gene expression:

The result is used to express the average of the number of corresponding molecules in a cell at a given time. From Equation (2), the equation that the average value is satisfied can be obtained:

Solve the equations:

The final concentration of GAD protein can be obtained:

[GAD] is the concentration of GAD protein.

To simplify the model, we introduced the comprehensive expression efficiency of DNA. The comprehensive expression efficiency of module one is:

Finally, GAD protein concentration is:

Module three

We hypothesized that in the cell, the input 5-FC drug concentration can be kept consistent with the outside world[2]. Then we can get the intracellular drug concentration [5-FC].

Binding protein expression

After GAD is combined with G8p, if the conditions are met:

There will be expression of HIF-1a ODDD and yCD binding protein.

We think that the number of GAD proteins is less than the number of DNA vectors of module three and it can achieve complete binding to G8p. Then its expression is satisfied:

[H-y] is the concentration of the HIF-1a ODDD and yCD binding protein. Thus, the final expression of the protein expression amount can be obtained:

After protein expression, the following reactions will be catalyzed:

When the amount of 5-FU reaches a certain level, the cells will be killed. The final concentration can be measured experimentally. We use F as the vitality:

Result

Cells can be divided into normal cell (called N) and cancer cell (called C). According to the sensitivity to the P1 system, they can be classified as Nor1, Nor2, Can2 and Can2. They satisfy the relationship in terms of the proportion:

The sensitivity of different cells is essentially positively correlated with the ayn1 of cell. For Nor1, Nor2, Can1 and Can2 cells, their relationship is:

The final system inactivates a certain number of cells, which is mainly determined by the amount of enzyme, which is HIF-1a ODDD and yCD binding protein. Then we set kcell as the coefficient corresponding to the enzyme-activated cells. The probability PNor1 that the cells are inactivated can be expressed as:

Therefore, in the experiment, in the case of the same number of normal cells and cancer cells, the number of normal cells and cancer cell inactivated cells obtained meets the following conditions:

The comparison between the simulation results and the experimental results is shown in Figure 3:(The left picture is the simulation data of dead cell number. The right picture is the experiment data of dead cell number.)

Simulation two

In order to reduce the killing of normal cells with high sensitivity to P1, we introduced miRNA regulation[3,4]. miRNAs are more commonly found in normal cells and a small fraction are present in cancer cells[5,6]. So in terms of quantity, the relationship that different cells satisfy is:

Module one

miRNA regulation process

If the influence of other factors on expression is not considered, because P1, GAD and BS are on the same component, for gene expression processes that have miRNA regulation, the following reactions exist:

For mRNA regulation process, we constructed the ODE model[7]:

In the model k is the chemical equilibrium constant which represents the rate of self-degradation of RNA or protein.

When the reaction reaches equilibrium, the effects of mRNA and miRNA self-degradation are negligible, and the final concentration of GAD protein can be obtained:：

[GAD] is the concentration that GAD is expressed, and [mRG] indicates the mRNA concentration of GAD. mRG*miR means mRNAGAD*miRNA.

Finally, GAD protein concentration is updated to be:

Module three

We also have the ODE model for the module three:

The concentration of binding protein was changed, too:

In the formula (27),[mRH-y] is the mRNA concentration of the binding protein,mRH-y*miR means mRNAH-y*miRNA.

Finally, the expression formula of the protein expression amount can be obtained:

Compared with simulation one, the number of inactivated cells in both normal cells and cancer cells was decreased. According to formula (19), the number of inactivated cells decreased was different, and the second simulation result was obtained. The comparison between the simulation results and the experimental results is shown in Figure 4:(The left picture is the simulation data of dead cell number. The right picture is the experiment data of dead cell number.)

Simulation three

After the addition of miRNA, the system has reduced the killing of cancer cells, so we introduced module two to reduce the effect of miRNA on the system in cancer cells.

Module two

First we set the threshold of miRNA regulation for the system to be N0. Then if the miRNA does not affect the system, the following condition must be met:

For a cancer cell in which miRNA regulation is present, it is only necessary to reduce its concentration to N0.

For the module two where the P2 promoter is located, [P2] is set to the P2 concentration which can express Sponge. The expressed Sponge concentration satisfies:

Among them, β is the expression coefficient which indicates the comprehensive expression efficiency of expressing Sponge in P2. [Sponge] is the Sponge concentration. Sponge is able to adsorb miRNAs, effectively reducing miRNA levels to level N0. If the adsorption rate is ξ, the condition of [miR] which is not added to module2 (remarked as [miR]0) is:

By increasing the concentration of P2, the miRNA will never reach the condition that is effective, and the number of cancer cells inactivated will reach the level of simulation one. From this we get the third simulation result. The comparison between the simulation results and the experimental results is shown in Figure 5:(The left picture is the simulation data of dead cell number. The right picture is the experiment data of dead cell number.)

Finally, according to formula (13), hypoxia can be used to control the killing effect on cells. The simulation results are shown in the following figure:(The left picture is the simulation data of hypoxia situation. The right picture is the simulation data of normoxia situation.)

Conclusion

Our model mathematically describes the process of gene expression and regulation of the system. Three simulation predictions were carried out for the experiment, which provided valuable theoretical guidance for the experiment. With the support of mathematical models, the construction of the system was more objective. The experiment was also more directional, and all material collection and experimental steps were carried out in the modeling approach. The three predictions with experimental results have proved that our mathematical model is accurate. At the same time, we also have the fourth forecast, which will be verified later.

References

[1]Lei J. Z, systems biology[M].Shanghai:Shanghai Science and Technology Press,2010:24-27.

[2]Lee, I. et al. New class of microRNA targets containing simultaneous 5′-UTR and 3′-UTR interaction sites. Genome Res. 19, 1175–1183 (2009).

[3]Cloonan, N. Re-thinking miRNA-mRNA interactions: intertwining issues confound target discovery. Bioessays 37, 379–388 (2015).

[4]Jonas, S. & Izaurralde, E. Towards a molecular understanding of microRNAmediated gene silencing. Nat. Rev. Genet. 16, 421–433 (2015).

[5]Hon, L. S. & Zhang, Z. The roles of binding site arrangement and combinatorial targeting in microRNA repression of gene expression. Genome Biol. 8, R166 (2007).

[6]Chou, T. C. & Talalay, P. Generalized equations for the analysis of inhibitions of Michaelis-Menten and higher-order kinetic systems with two or more mutually exclusive and nonexclusive inhibitors. Eur. J. Biochem. 252, 6438–6442 (1981).

[7]Jeremy J. G, JONATHAN B & RON W．A mixed antagonistic/synergistic miRNA repression model enables accurate predictions of multi-input miRNA sensor activity[J]．Nature Communications，2018，9:2430 ．