Overview

To comprehensively evaluate the feasibility and controllability of our oral medicine, we established a two-part mathematical framework. Model 1 focuses on the pharmacokinetics of OMVs, simulating their in vivo transport, absorption, and distribution to predict the concentration of therapeutic proteins at hemorrhoidal sites. Building upon this foundation, Model 2 introduces an expiry date circuit within the engineered bacteria to ensure it can survive in the human body for a sufficient length of time while maintaining human safety.

Together, these models provide a dynamic of drug delivery and bacterial self-regulation: the first quantifies how effectively the active protein reaches the target site; the second ensures that the engineered bacteria only die off after surviving for a sufficient period of time. This integrated modeling approach offers both mechanistic validation and safety assurance for our oral medicine.

Model 1: Pharmacokinetic Model of OMVs

In order to simulate the in vivo transport efficiency of OMVs (Outer Membrane Vesicles) and to preliminarily estimate the final drug concentration targeting hemorrhoids, we adopted a pharmacokinetic model to validate the effectiveness of the approach at the modeling level.

1. Target of Final Simulation

This study ultimately focuses on the active proteins that exert therapeutic effects at the wound site following oral administration:

a. CAR: Anchored on the surface of OMVs, which mediates wound-targeting activity;

b. bioPROTAC: Encapsulated within OMVs, which exerts anti-angiogenic effects.

2. Target of Process Simulation

In the pharmacokinetic process, OMVs are considered the primary simulation subject. Since OMVs serve as the delivery vehicles for active proteins, the entire process—from protein production to their therapeutic action on hemorrhoids—relies on OMV-mediated transport. Therefore, parameters such as the degradation rate and transport rate of the active proteins in the bloodstream are assumed to be equivalent to the corresponding pharmacokinetic characteristics of OMVs.

3. Action Target of Drug

Considering the differences in venous return between types of hemorrhoids: For IH (internal hemorrhoids), venous blood must pass through the liver before entering the systemic circulation, thereby undergoing significant first-pass metabolism. In contrast, for EH (external hemorrhoids), venous blood drains directly into the systemic circulation without hepatic passage. To capture this pharmacokinetic distinction, the IH module incorporates hepatic metabolism pathways and related parameters, whereas in the EH module, drug absorption flux is directly coupled to the systemic circulation[1].

To further improve the overall model and to verify the theoretical feasibility of "the endocytosis process of OMVs crossing the small intestinal villus epithelial cell membrane can proceed spontaneously", we conducted the following thermodynamic analysis:

From a thermodynamic perspective, the spontaneity of substance transmembrane processes is determined by the ΔG (Gibbs free energy change): if ΔG = 0, the process can proceed spontaneously. Therefore, it is necessary to derive the spontaneous condition for OMV transmembrane endocytosis by analyzing the equilibrium relationship between interfacial tension and membrane bending energy, namely the competition between the energy driving membrane invagination and the energy resisting membrane bending.

1. Total Energy of the Initial State: The initial state involves the surface energy of three types of interfaces

$$\displaystyle E_1 = \sigma_1 \cdot 2\pi R_1^2 + \sigma_2 \cdot \pi R_1^2 + \sigma_3 \cdot \left(S_{\text{tot}} - \pi R_1^2\right)$$

2. Total Energy of the Final State: Based on the Helfrich membrane elasticity theory, the final state involves vesicle interfacial energy and membrane bending elastic energy

$$\displaystyle E_2 = \sigma_2 \cdot 4\pi R_2^2 + \sigma_3 \cdot \left(S_{\text{tot}} - 4\pi R_2^2\right) + k \cdot \left(\frac{1}{R_2}\right)^2 \cdot 4\pi R_2^2$$

Since the wrapped volume remains unchanged before and after endocytosis:

$$\displaystyle \frac{2}{3}\pi R_1^3 = \frac{4}{3}\pi R_2^3$$

Simplification gives:

$$\displaystyle R_2^2 \approx 0.630 R_1^2$$

The thermodynamic criterion for a spontaneous process is

$$\displaystyle E_2 < E_1$$

the total energy decreases as

$$\displaystyle \Delta G = E_2 - E_1 < 0$$

By substituting the parameters into the above thermodynamic criterion and simplifying, we obtain:

$$\displaystyle \frac{R_1^2 \cdot \left(0.5\sigma_1 - 0.38\sigma_2 + 0.38\sigma_3\right)}{k} = 15.3 > 1$$

To further verify the model's reliability, a robustness analysis was conducted by varying the model parameters: If the OMV size is small (with a diameter of only 10 nm) , the inequality becomes 2.46>1, indicating that the spontaneous condition is still satisfied under physiological conditions relying on a high interfacial tension difference; If the interfacial tension difference is negative, the driving force is insufficient, and the process is non-spontaneous. However, this scenario contradicts the high affinity between OMV surface components and the membrane, making it rare under physiological conditions; If the cell membrane rigidity increases ( k = 2 × 10 ^-19 J ), the spontaneous tendency is weakened, but the inequality 4.83>1 still holds, meaning the spontaneous condition is satisfied under typical parameter settings.

In conclusion, this derivation thermodynamically proves that OMVs can spontaneously cross the small intestinal villus epithelial cell membrane. This provides a physical theoretical basis for OMVs to be used as nanocarriers for drug delivery and explains the energy rationality of the interaction between intestinal flora and host cells via OMVs.

Table 1. Parameters and values in thermodynamic analysis of OMVs

Parameters

Parameter	Physical meaning	Value	Unit	Source
R2	Initial characteristic radius of OMVs	5 × 10 -8	m	[2]
σ1	Interfacial tension between OMVs and extracellular fluid	10 -3	N·m -1	[4]
σ2	Interfacial tension between OMVs and the cell membrane	8 × 10 -4	N·m -1	[3]
σ3	Interfacial tension between the cell membrane and intracellular fluid	5 × 10 -4	N·m -1	[4]
k	Bending modulus of the cell membrane	10 -19	J	[5]

Our main governing equations include Growth Rate Equation, Small Intestine Transport Equation, and Blood Circulation Equation. Through the model, these three components are overall coupled to represent the complete process of oral administration, OMV transport and absorption in the small intestine, and delivery to the systemic circulation and the hemorrhoid site. The following sections present the detailed design, equations, and corresponding results of these regulatory equations.

The experimental group employed a logistic growth curve for fitting(click here). Based on this, we incorporating the active protein content corresponding to a single bacterial cell, the initial parameter values for the model were further obtained.

$$\displaystyle \frac{df}{dt} = \frac{A k e^{-k(t - T)}}{\left(1 + e^{-k(t - T)}\right)^2}$$

In the above equation:

t: time, in h

A: maximum OD600 value, 0.3911

k: growth rate constant, 0.008227 h^-1

T: inflection point of the system, 2.79 h

The ADAM model (Advanced Dissolution and Absorption and Metabolism Model) is a physiology-based pharmacokinetic model widely used to simulate the dissolution, absorption, and metabolism of orally administered drugs in vivo[6]. Considering the mechanistic differences between OMVs-based delivery systems and conventional small-molecule oral drugs, we appropriately simplified and modified the original ADAM model. In particular, the core "transport-absorption" module was retained and optimized to better reflect the characteristics of OMVs-mediated delivery.

a. For the transport module, the small intestine was partitioned into 7 functional segments[7] (remove the colon part). The OMVs transport rate within each segment was computed iteratively, allowing the dynamic profile of active protein concentration along the intestinal tract to be determined.

The following equations describe the transport processes within the small intestine module:

$$\displaystyle \frac{dA_D[n]}{dt} = G + I - \left(k_{\mathrm{deg}}[n] + k_{t}[n]\right) A_D[n]$$

b. For the absorption module, the trans-epithelial transport rate of OMVs into intestinal cells was calculated, allowing the temporal profile of active protein concentrations within the epithelial cells to be determined.

The following equations describe the absorption processes across the intestinal epithelium:

$$\displaystyle \frac{dC_{\mathrm{ent}}[n]}{dt} = \frac{1}{V_{\mathrm{ent}}[n]} \left(A_D[n] k_a[n] - C_{\mathrm{ent}}[n] Q_{\mathrm{ent}}[n]\right)$$

In the above equation:

A_D[n]: protein amount in the n-th intestinal segment, in ng

C_ent[n]: protein concentration in the epithelial cells of the n-th segment, in ng/ml

G(protein_generation): protein generation rate, namely the input derived from the growth rate

I(input_prev): input from the previous segment, in ng/min

k_t[n]: protein transport rate constant in the lumen of the n-th segment, in min^-1

Table 2. Parameters and values in transport and absorption equations

Parameters

Parameter	Description	Value	Unit	Source
kdeg	protein degradation rate constant in the lumen of the n-th segment	2.08 × 10-3	min -1	[9]
Vent	volume of the n-th small intestine section	79	ml	[8]
ka	absorption rate constant from lumen to epithelium	0.1	min -1	[-]
Qent	clearance rate of epithelial tissue in the n-th segment	129.08	ml/min	[8]

This section is based on the mass balance equation of exosomes from iGEM NJU 2013. Considering that in our model, OMVs enter the bloodstream from the intestinal epithelium as a continuous input rather than a single pulse, the original input term was modified to a continuous form proportional to the epithelial cell concentration \(C_\mathrm{cell}(t)\). At the same time, since there is no clear distinction between EH and IH in medical concepts, we consider them as rapidly perfused tissues. Accordingly, the mass balance equations for IH were derived:

a. Equation for EH:

$$\displaystyle \frac{dC_\mathrm{EH}}{dt} = K_\mathrm{in} \cdot C_\mathrm{cell}(t,n) - K_\mathrm{out} \cdot C_\mathrm{EH}(t)$$

Where:

$$\displaystyle K_\mathrm{in} = \frac{\frac{Q_\mathrm{EH} P_\mathrm{EH}}{Q_c}}{\sum_{\text{tissue}} \frac{Q_i P_i}{Q_c} + k_\mathrm{mblood} - k_{m,\mathrm{EH}}}, \quad K_\mathrm{out} = k_{m,\mathrm{EH}}$$

b. Equation for IH:

$$\displaystyle \frac{dC_\mathrm{IH}}{dt} = K_\mathrm{in} \cdot C_\mathrm{cell}(t,n) - K_\mathrm{out} \cdot C_\mathrm{IH}(t)$$

Where:

$$\displaystyle K_\mathrm{in} = \frac{\frac{Q_\mathrm{IH} P_\mathrm{IH} - Q_\mathrm{liver} P_\mathrm{liver}}{Q_c}}{\sum_{\text{tissue}} \frac{Q_i P_i}{Q_c} + k_\mathrm{mblood} + k_\mathrm{m,liver} - k_{m,\mathrm{IH}}}, \quad K_\mathrm{out} = k_{m,\mathrm{IH}}$$

In the above equation:

C_cell[t,n]: concentration of active protein generated in the intestinal epithelium, serving as the input source, in μg/ml

C_EH[t],C_IH[t]: concentrations of active protein in the blood for external hemorrhoids (EH) and internal hemorrhoids (IH), respectively, in μg/ml

K_in: effective input rate from the intestinal epithelium to the bloodstream, in min^-1

K_out: clearance rate of active protein in the bloodstream, in min^-1

Table 3. Parameters and values in blood circulation equation

Parameters

Parameter	Description	Value	Unit	Source
Qc	cardiac output, namely total systemic blood flow	372	L/h	[10]
Qliver	blood flow of the liver	83.75	L/h	[11]
Pliver	liver tissue partition coefficient	0.12	[-]	[12]
km,liver	metabolism rate constant of the liver	2.625×10-3	min-1	[11]
QEH,QIH	blood flow rates for EH and IH tissues	105.86	L/h	[10]
km,EH,km,IH	clearance rate constants for EH and IH tissues, respectively	0.1	min-1	[-]
PEH,PIH	tissue partition coefficients for EH and IH	0.66	[-]	[12]
kmblood	nonspecific degradation rate constant of OMVs in the bloodstream	0.05	min-1	[-]

Based on the experimental data, we obtained the initial concentrations of CAR and bioPROTAC. Since OMVs function as a continuous input system, the concentration of active protein does not decline over time but instead approaches a steady state after a certain period.

First, we fitted the concentrations of active proteins in the intestinal epithelial tissue:

Counting from the time when the engineered bacteria have colonized, the results showed that both CAR and bioPROTAC reached stability at approximately 10 minutes. The steady-state concentration of CAR was about 5.5 μg/ml, while that of bioPROTAC was about 2.5 μg/ml.

Next, by applying the blood circulation equation, we further simulated the concentrations of CAR and bioPROTAC in EH and IH. Considering that hemorrhoidal tissue is not purely a rapidly perfused tissue, the actual concentration of active protein may be lower than the simulated values.

In EH, at around 60 minutes, the steady-state concentration of CAR was about 47.5 μg/ml, while that of bioPROTAC was about 21.5 μg/ml.

Figure 6. Concentration of active protein in IH.

In IH, at around 60 minutes, the steady-state concentration of CAR was about 40 μg/ml, while that of bioPROTAC was about 18 μg/ml.

1. The steady-state concentration of CAR was approximately 47.5 μg/ml in EH and about 40 μg/ml in IH. Although the concentration in internal hemorrhoids was lower, this level can still achieve a certain degree of targeted action;

2. The steady-state concentration of bioPROTAC was approximately 21.5 μg/ml in EH and about 18 μg/ml in IH. Both concentrations exhibit therapeutic potential, though the effect in internal hemorrhoids is expected to be weaker than in EH;

3. In addition to IH and EH, mixed hemorrhoids are also common in clinical practice. The corresponding concentration of active proteins is expected to fall between EH and IH, suggesting that oral administration may also exert therapeutic effects in this condition.

Model 2: Expiry-Date Circuit

In this study, we designed a expiry-date circuit (Figure 7). Conceptually, upon food intake, patients ingest a defined amount of small molecules, which serve as ligands binding to the riboswitch and thereby continuously suppressing the synthesis of the toxin protein CcdB. Once food intake ceases, the small molecules gradually degrade and are consumed, allowing the production of CcdB to resume and finally leading to the elimination of the engineered bacteria.

During the modeling process, our objective was to predict how variations in the type and dosage of small molecules affect the duration of CcdB suppression. Based on these predictions, we estimated the optimal dosing interval for small-molecule ingestion, thereby providing further validation of the effectiveness and feasibility of the proposed oral drug delivery strategy.

1. We constructed an ODE (ordinary differential equation) model to describe the gene expression process. The inhibitory effect of the riboswitch on downstream genes, as well as its persistence, was modeled using Hill function;

2. In addition to predicting the duration of CcdB suppression under different ligand dosages, we selected an initial small-molecule concentration of 20nM for simulation. This concentration represents a relatively low dosage level, consistent with the plausible range of small-molecule concentrations in food reaching the small intestine, and is considered representative in our model;

3. The threshold was defined as 10% of the maximum concentration of CcdB. When expression levels reach this threshold, the engineered bacteria are assumed to die. At the corresponding time point, the patient is required to supplement with the appropriate ligand dosage to maintain effective control over the engineered bacteria.

$$\displaystyle {Gene}_{MerR} \stackrel{\alpha_{MerR} \cdot H_{rep,X}}{\longrightarrow} mRNA_{MerR}\quad\quad(r1.1)$$

$$\displaystyle mRNA_{MerR} \stackrel{k_{\text{tln}, MerR}}{\longrightarrow} MerR \quad\quad(r1.2)$$

$$\displaystyle {Gene}_{cI} \stackrel{\alpha_{cI} \cdot H_{rep,MerR}}{\longrightarrow} mRNA_{cI}\quad\quad(r1.3)$$

$$\displaystyle mRNA_{cI} \stackrel{k_{\text{tln}, cI}}{\longrightarrow} cI \quad\quad(r1.4)$$

$$\displaystyle {Gene}_{ccdB} \stackrel{\alpha_{ccdB} \cdot H_{rep,cI}}{\longrightarrow} mRNA_{ccdB}\quad\quad(r1.5)$$

$$\displaystyle mRNA_{ccdB} \stackrel{k_{\text{tln}, ccdB}}{\longrightarrow} CcdB \quad\quad(r1.6)$$

where,

$$\displaystyle \quad H_{rep,X} = \frac{(K_X)^{n_X}}{(K_X)^{n_X} + (X)^{n_X}},$$

$$\displaystyle \quad H_{rep,MerR} = \frac{(K_{MerR})^{n_{MerR}}}{(K_{MerR})^{n_{MerR}} + ([MerR])^{n_{MerR}}},$$

$$\displaystyle \quad H_{rep,cI} = \frac{(K_{cI})^{n_{cI}}}{(K_{cI})^{n_{cI}} + ([cI])^{n_{cI}}}$$

Degradation processes as follows:

$$\displaystyle mRNA_{MerR} \stackrel{\delta_{m_{MerR}}}{\longrightarrow} \varnothing\quad\quad(r1.7)$$

$$\displaystyle MerR \stackrel{\delta_{MerR}}{\longrightarrow} \varnothing\quad\quad(r1.8)$$

$$\displaystyle mRNA_{cI} \stackrel{\delta_{m_{cI}}}{\longrightarrow} \varnothing\quad\quad(r1.9)$$

$$\displaystyle cI \stackrel{\delta_{cI}}{\longrightarrow} \varnothing\quad\quad(r1.10)$$

$$\displaystyle mRNA_{ccdB} \stackrel{\delta_{m_{ccdB}}}{\longrightarrow} \varnothing\quad\quad(r1.11)$$

$$\displaystyle CcdB \stackrel{\delta_{CcdB}}{\longrightarrow} \varnothing\quad\quad(r1.12)$$

$$\displaystyle \frac{d[mRNA_{MerR}]}{dt} = \alpha_{MerR} \cdot \frac{(K_X)^{n_X}}{(K_X)^{n_X} + (X)^{n_X}} - \delta_{m_{MerR}} \cdot [mRNA_{MerR}]\quad\quad(f1.1)$$

$$\displaystyle \frac{d[MerR]}{dt} = k_{\text{tln}, MerR} \cdot [mRNA_{MerR}] - \delta_{MerR} \cdot [MerR]\quad\quad(f1.2)$$

$$\displaystyle \frac{d[mRNA_{cI}]}{dt} = \alpha_{cI} \cdot \frac{(K_{MerR})^{n_{MerR}}}{(K_{MerR})^{n_{MerR}} + ([MerR])^{n_{MerR}}} - \delta_{m_{cI}} \cdot [mRNA_{cI}]\quad\quad(f1.3)$$

$$\displaystyle \frac{d[cI]}{dt} = k_{\text{tln, cI}} \cdot [mRNA_{cI}] - \delta_{cI} \cdot [cI]\quad\quad(f1.4)$$

$$\displaystyle \frac{d[mRNA_{ccdB}]}{dt} = \alpha_{ccdB} \cdot \frac{(K_{cI})^{n_{cI}}}{(K_{cI})^{n_{cI}} + ([cI])^{n_{cI}}} - \delta_{m_{ccdB}} \cdot [mRNA_{ccdB}]\quad\quad(f1.5)$$

$$\displaystyle \frac{d[CcdB]}{dt} = k_{\text{tln, ccdB}} \cdot [mRNA_{ccdB}] - \delta_{CcdB} \cdot [CcdB]\quad\quad(f1.6)$$

The model is initialized at the onset of small-molecule release, assuming that its concentration decreases exponentially over time：

$$\displaystyle X(t) =\begin{cases}X_{\text{init}}, & t < 0 \\X_{\text{init}} \cdot e^{-k_{\text{decayX}} \cdot t}, & t \ge 0\end{cases}\quad\quad(f1.7)$$

In the riboswitch design module, we calculated the binding free energy between the aptamer and its ligand to predict their interaction characteristics. The results indicated that, under our designed riboswitch configuration, IC₅₀ (the half-inhibitory concentration) of hippuric acid is 7507.83 nM(click here). Based on this parameter, we simulated the system dynamics and obtained the following results.

The results showed that although IC₅₀ exhibits a relatively high half-inhibitory concentration, the time required for CcdB to reach the predefined threshold remained stable at 41.42 h. This suggests that patients would only need to supplement the small molecule approximately every 41.42 h. Such an interval falls within a reasonable range and represents a patient-friendly dosing schedule.

Furthermore, considering the diversity of riboswitch-ligand pairs in oral drug delivery system, we extended the simulations to examine CcdB expression dynamics under different IC₅₀.

Given that the IC₅₀ may vary across a wide range, we set broad hypothetical values in the model. The results demonstrated that when IC₅₀ = 35 nM and IC₅₀ = 1×10⁸ nM, the resulting CcdB expression curves exhibited minimal differences, indicating that beyond a certain IC₅₀ threshold, the time for CcdB to reach the predefined level becomes only marginally affected. This phenomenon is consistent with the properties of the Hill equation.

We also evaluated the influence of different X_init (initial ligand concentrations) on the model outcome. Under the assumption of IC₅₀ = 50 nM, we simulated CcdB expression curves at various initial concentrations.

The results indicated that even under extreme conditions, when the initial ligand concentration was as low as 1 nM, the time required for CcdB to reach the threshold remained around 41.45 h. This finding further supports the effectiveness and robustness of the expiry date circuit in the context of oral drug delivery system.

Table 4. Parameters and values in expiry date circuit

Parameters

Parameter	Description	Value	Unit	Source
kdecayX	Degradation rate of X	8.23 × 10-2	min-1	[13]
αMerR	Transcription rate of MerR mRNA	38.88	nM/min	[14]
nX	Hill coefficient of X to PJ23100	2.23	[-]	[15]
ktln,MerR	Translation rate of MerR mRNA	6.93	min-1	[17]
δmMerR	Degradation rate of MerR mRNA	5.78 × 10-3	min -1	[17]
δMerR	Degradation rate of MerR	0.046	min-1	[15]
KMerR	IC50 of MerR	10	nM	[15]
nMerR	Hill coefficient of MerR to PMer	2.5	[-]	[15]
ktln,cI	Translation rate of cI mRNA	6.93	min-1	[17]
δmcI	Degradation rate of cI mRNA	5.78 × 10-3	min-1	[17]
δcI	Degradation rate of cI	2.5 × 10-3	min-1	[17]
αccdB	Transcription rate of ccdB mRNA	38.98	nM/min	[14]
KcI	IC50 of cI	30	nM	[17]
ncI	Hill coefficient of cI to PL	2	[-]	[-]
ktln,ccdB	Translation rate of ccdB mRNA	1.98	min-1	[16]
δmccdB	Transcription rate of cI mRNA	38.275	nM/min	[16]
δCcdB	Degradation rate of CcdB	0.0693	min-1	[16]

References

[1] Rathi, R., Sanshita, Kumar, A., Vishvakarma, V., Huanbutta, K., Singh, I., & Sangnim, T. (2022). Advancements in rectal drug delivery systems: Clinical trials, and patents perspective. Pharmaceutics, 14(10), 2210.

[2] Deatherage, B. L., & Cookson, B. T. (2009). Bacterial outer membrane vesicles: A rising tide of interest. Nature Reviews Microbiology, 7(11), 803-811.

[3] Kulp, A. L., & Kuehn, M. J. (2010). Biogenesis and functions of bacterial outer membrane vesicles. Microbiology and Molecular Biology Reviews, 74(1), 14-34.

[4] Evans, E., & Skalak, R. (1980). Biomechanics of soft tissues. Biophysical Journal, 30(1), 25-54.

[5] Boal, D. (2012). Mechanics of the Cell. Cambridge University Press.

[6] Jamei, M., Turner, D., Yang, J., Neuhoff, S., Polak, S., Rostami-Hodjegan, A., & Tucker, G. (2009). Population-based mechanistic prediction of oral drug absorption. AAPS Journal, 11(2), 225–237.

[7] Maurer, A. H. (2015). Gastrointestinal motility, part 2: Small-bowel and colon transit. Journal of Nuclear Medicine, 56(9), 1395–1400.

[8] Olivares-Morales, A., Ghosh, A., Aarons, L., & Rostami-Hodjegan, A. (2016). Development of a novel simplified PBPK absorption model to explain the higher relative bioavailability of the OROS® formulation of oxybutynin. AAPS Journal, 18(6), 1532–1549.

[9] Takita, H., Darwich, A. S., Ahmad, A., & Rostami-Hodjegan, A. (2020). Application of the nested enzyme-within-enterocyte (NEWE) turnover model for predicting the time course of pharmacodynamic effects. CPT: Pharmacometrics & Systems Pharmacology, 9(11), 617–627.

[10]https://github.com/Open-Systems-Pharmacology/PK-Sim/releases.

[11]https://2013.igem.org/Team:NJU_China/Modeling.

[12] Wang, J., Hu, L., Huang, H., Yu, Y., Wang, J., Yu, Y., Li, K., Li, Y., Tian, T., & Chen, F. (2020). CAR (CARSKNKDC) peptide modified ReNcell-derived extracellular vesicles as a novel therapeutic agent for targeted pulmonary hypertension therapy. Hypertension, 76(4), 1147–1160.

[13]Nitzan Rosenfeld et al. (2005).Gene Regulation at the Single-Cell Level.Science,307,1962-1965.

[14]https://2009.igem.org/File:PKU_Sensi.PNG_

[15]https://2011.igem.org/Team:Grenoble

[16]https://2020.igem.org/Team:BITSPilani-Goa_India/Model/Kill_Switch#DifferentialEquations

[17]Tran KM, Nong NT, Ren J, Lee K, Lee D, Gsponer J, Lee HM, Na D.(2025).Genetic "expiry-date" circuits control lifespan of synthetic scavenger bacteria for safe bioremediation. Nucleic Acids Res.Jul 19;53(14):gkaf703.