To quantitatively analyze the workflow of the Glycemic Stabilizer, ensure consistency in the design of steps, validate experimental data, predict product operation, and further improve the product, we have established a dynamic mathematical model. This model is primarily constructed using Delay Differential Equations (DDE), based on existing data. It simulates the gene circuits of engineered cells that sense blood glucose and blue light, followed by insulin secretion and the associated main physiological processes.

We have improved the existing model to make it suitable for our Glycemic Stabilizer design. The DDEs model is mainly composed of three modules:

This module primarily involves the process where glucose enters engineered cells and alters the cellular metabolic state. This change triggers a series of reactions, leading to the initiation of insulin expression and secretion upon blue light exposure. Additionally, the engineered cells have a feedback inhibition mechanism that suppresses insulin expression when extracellular insulin is detected.

This module explains the protective mechanism involving the binding of miRNA to insulin mRNA and the binding of sponge to miRNA.

This module simulates the dynamic changes in blood glucose levels within the metabolic environment of the entire body's blood and tissues, as well as the body's response to insulin.

Due to space limitations, the module selection will include a brief introduction to the model's detailed description rather than an explanation of each equation. Specific parameters can be downloaded and viewed in PDF format. For more information, please refer to the references of our model.

The equations used in delay differential equations primarily include the chemical master equation, Hill equation, and Michaelis equation. The chemical master equation characterizes the transformations between various molecules; the Hill equation is used to describe the activation level of a substance, where the concentration of the substance depends on another substance; the Michaelis equation primarily models the instantaneous rate as compared to the maximum rate under specific conditions.

To begin, we employ ordinary differential equations (ODEs) to describe the process by which this module functions within the system. The glucose concentration-sensitive promoter serves as the input component and initiates the entire system. Therefore, modeling this part is crucial for evaluating the overall design and feasibility of the system. However, due to insufficient parameter data from wet lab experiments, we conducted simulations to supplement and reorganize the model.

We assume that the system is activated when glucose enters the engineered cells. Let Ga represent the glucose concentration in the encapsulation carrier, R1 the glucose uptake rate by a single engineered cell, and N the number of engineered cells in the encapsulation carrier. Thus, the rate of glucose uptake in the encapsulation carrier can be expressed as N·R1. The glucose metabolism rate within a single engineered cell is denoted as R2. To simplify the equations, we reduce the complex glucose metabolism process into two steps: the production of intermediate metabolites and ATP generation. The rate R2 is dependent on the levels of ATP and intracellular glucose.

As a result, the dynamic changes in glucose can be described using the following differential equation:

We established equations (4) to (7) to describe the relationship between the glucose degradation rate and ATP within engineered cells.

The meanings of the important parameters we have defined are as follows:

R1: The rate at which a single cell absorbs glucose;

R2: The rate at which a single cell metabolizes glucose;

R3: The rate at which glucose metabolites are further broken down to produce ATP.

Within the engineered cells, glucose acts as the signal that triggers subsequent processes, specifically the activation of the transcriptional activators Gal4 and VP16. The activation level of these transcriptional activators depends on the intracellular glucose concentration, denoted as Gin.

The meanings of the important parameters we have defined are as follows:

R4: The fold increase in transcriptional expression rate of the Gal4 promoter in the cell;

R5: The fold increase in transcriptional expression rate of the VP16 promoter in the cell.

Once activated, these transcriptional activators initiate the transcription of the Gal4 and VP16 genes, leading to the production of their respective mRNAs and subsequent translation into Gal4 and VP16 proteins.

At the same time, VP16 is initially inactive, but exposure to blue light can rapidly activate it, transforming it into VP16*. Only then can it reversibly bind with Gal4 to become the transcription factor TF for the insulin gene. VP16* itself will undergo hydrolysis and can also revert to an inactive conformation. This regulatory process can thus be described by the following differential equations:

(R6: The fold increase in transcriptional expression rate of the Ins promoter in the cell)

miRNA can specifically recognize and bind to insulin mRNA, thereby inhibiting the translation of insulin mRNA.Therefore, the following differential equations can be used to describe this regulatory process:

The activation level of the transcriptional activator for sponge depends on the intracellular glucose concentration (Gin). Once activated, the transcriptional activator drives the transcription of sponge. sponge can function as a competitive endogenous RNA (ceRNA), sequestering miRNA away from insulin mRNA, thereby promoting the translation of insulin mRNA.

(R7: The fold increase in transcriptional expression rate of the circ promoter in the cell)

In this experiment, we focused on studying the expression levels of the GAL4-VP16 transcription factor under varying glucose concentrations and temporal conditions. Using simulations based on the Hill equation, we proposed the following hypothesis: As glucose concentration increases and time progresses, the expression level of GAL4-VP16 will show a significant upward trend, which conforms to the dynamic behavior predicted by the Hill equation.

We hypothesize that with increasing glucose concentration, the expression levels of GAL4-VP16 will rise significantly. This assumption is based on the biological properties of GAL4-VP16 being sensitive to glucose concentrations. An increase in glucose may activate certain metabolic pathways, thereby enhancing the expression of GAL4-VP16. Therefore, using the Hill equation for simulation, we predict that GAL4-VP16 expression peaks in a high-glucose environment.

We further hypothesize that the expression levels of GAL4-VP16 are regulated not only by glucose concentration but also closely related to time. Specifically, as time progresses, the expression of GAL4-VP16 will gradually increase until it reaches a steady state. This hypothesis is based on the time-dependency of transcription factor expression, where the expression levels increase progressively with continuous signal transduction within the cell, reaching a maximum at a certain time point.

Moreover, we hypothesize that the expression level of GAL4-VP16 is jointly influenced by both time and glucose concentration. Under conditions of high glucose concentration and prolonged duration, GAL4-VP16 expression will exhibit its strongest activation state. In the simulation results of the Hill equation, we anticipate seeing the expression levels of GAL4-VP16 peak in a three-dimensional graph under combinations of higher glucose concentrations and longer time points, aligning closely with the observed trends in experimental fluorescence changes.

Using RStudio to simulate the glucose-induced transcriptional expression process, the simulation results are as follows:

According to literature, normal human blood glucose levels typically range from 3.9 to 6.1 mmol/L, with postprandial levels rising to between 5.0 and 7.8 mmol/L two hours after a meal. As illustrated, simulations conducted with numerical parameters revealed that the expression of GAL4-VP16 increases as glucose levels rise. Moreover, experimental results demonstrate that with increasing glucose levels, the expression of GAL4-VP16 indeed exhibits an upward trend, which closely aligns with predictions made using the Hill equation. The high concordance between the experimental findings and our model validates our hypothesis that the expression of GAL4-VP16 is regulated by glucose concentration and time, and that its expression behavior is biologically plausible. This experiment further confirms that simulations based on the Hill equation can effectively predict transcription factor expression patterns in complex biological systems, providing a reliable model foundation for subsequent biological experimental design and data analysis.

miRNA can specifically bind to the miRNA binding site (miRNA-BS) on insulin mRNA, thereby affecting the translation of insulin mRNA. This process can be represented by the following equation:

Thus,there is

Equation (26) describes the time-dependent changes in insulin concentration, which are influenced by transcriptional levels, miRNA inhibition, and self-hydrolysis. Equation (27) illustrates the changes in insulin transcription levels, while Equation (28) depicts the time-dependent variation in the concentration of the miRNA-insulin RNA complex.

Moreover, sponge can specifically bind to miRNA, thereby influencing the interaction between miRNA and insulin mRNA, regulating the expression of the insulin gene. This process can be represented by the following equations.

Simultaneously, the synthesis of sponge is influenced by glucose concentration. Consequently, the regulatory process can be described using the following differential equations:

Equation (31) describes the time-dependent changes in sponge concentration, influenced by transcriptional levels, miRNA activity, and self-hydrolysis. Equation (32) captures the time-dependent variation in sponge transcription levels as influenced by glucose concentration, while Equation (33) details the changes in the concentration of the miRNA-sponge complex over time.

miRNA plays a central role in the regulatory process described above, with its concentration governed by the following equation:

Considering the lag in transcription regulation due to changes in product concentration, we define

: the time required for the cell to sense changes in glucose concentration;

: the time required for the cell to sense changes in glucose concentration;

Considering that the mechanisms of action may vary across different targets, distinctions are clearly defined.

The expressions in the above equation are then updated as follows:

In this module, we constructed a simulation model designed to investigate the interactions among miRNA, sponge, and insulin mRNA, as well as their effects on the synthesis and regulation of insulin. Building on our understanding of RNA regulatory mechanisms, we formulated the following hypotheses and validated them through simulation results.

We first hypothesize that as blood glucose levels increase, sponge will bind with miRNA, thereby relieving miRNA's suppression of insulin mRNA. miRNA typically inhibits translation or causes mRNA degradation by binding to the 3'UTR of mRNA. However, the presence of sponge can act as a "sponge" for miRNA, competitively binding and thereby freeing insulin mRNA from miRNA suppression. Consequently, over time, the binding between miRNA and sponge is expected to strengthen, while the interaction between miRNA and insulin mRNA diminishes.

We further hypothesize that the expression levels of insulin mRNA will gradually increase over time, especially when the inhibitory effect of miRNA is reduced. With rising blood glucose levels, the synthesis rate of insulin mRNA will markedly increase, ultimately reaching a relatively stable high expression level. We predict that under conditions of persistently high blood glucose levels, insulin mRNA expression will account for more than 60% of total RNA, meeting the body's insulin production needs and effectively regulating blood glucose levels.

We also hypothesize that within the system, the total amounts of miRNA, sponge, and insulin mRNA remain constant (summing to 1), indicating the presence of a dynamic equilibrium mechanism within the system. As the binding between miRNA and sponge strengthens, the inhibitory effect of miRNA on insulin mRNA decreases, thereby promoting an increase in insulin mRNA expression. When the proportion of insulin mRNA rises to about 0.6, the system's regulatory mechanisms will ensure that insulin production is sufficient to lower blood sugar levels but not excessive, to avoid the risk of hypoglycemia.

Using RStudio to simulate the interaction processes of mRNA, miRNA, and sponge, the simulation results are as follows:

As depicted in the model, when blood glucose levels rise, the synthesis of insulin mRNA increases over time, and simultaneously, sponge binds to miRNA, thereby relieving miRNA's suppression of insulin mRNA. This demonstrates the sponge effect of sponge. Concurrently, experimental results indicate that over time, the binding rate between miRNA and sponge progressively increases, while the binding rate between miRNA and insulin mRNA gradually decreases. The expression level of insulin mRNA significantly rises and stabilizes at a higher level in a high glucose environment. This trend is highly consistent with our hypotheses, confirming the regulatory mechanism within the system where sponge alleviates miRNA's suppression of insulin mRNA through the "sponge effect," thus promoting insulin production. The experimental findings not only validate the complex interactions between miRNA, sponge, and mRNA but also provide new insights into insulin production and blood glucose regulation.

This part of the model simulates the changes in blood glucose and insulin levels under both physiological and pathological conditions. When combined with the two modules mentioned earlier, it allows for the analysis of the dynamics of blood glucose and insulin levels throughout the body during the operation of the Glycemic Stabilizer. This analysis, based on existing data, facilitates subsequent improvements. Here, we use a Proportional-Integral-Derivative (PID) controller to model the changes in insulin concentration.

The equations describing the relationships between insulin mRNA, insulin concentration in the blood, and insulin concentration in alginate are established as follows:

After being transcribed by TF and translated within the alginate, insulin diffuses between the alginate and the blood, driven by the concentration gradient.Thus, the process can be described using the following differential equations:

To effectively describe the dynamics of insulin and glucose under experimental conditions and to reflect their response patterns at different time points, we propose the following hypotheses:

We hypothesize that the concentration of insulin will rapidly decrease over time and gradually stabilize. Insulin, a crucial hormone for regulating blood glucose, is secreted rapidly when blood glucose levels rise and then its concentration decreases as blood glucose stabilizes within the body. Our model aims to fit this dynamic process of rapid decline followed by stabilization. By adjusting the HillSlope and the asymmetry parameter S in the model, we expect to capture the asymmetry of the insulin concentration curve, especially during the initial rapid decrease and subsequent slow stabilization.

 For changes in glucose concentration, we hypothesize that it will initially decrease rapidly over time and stabilize after a certain period. Our hypothesis aims to accurately describe this curve characteristic, particularly in the early stages where the rapid decrease in glucose concentration may be influenced by several factors, including the rapid secretion of insulin and cellular uptake of glucose.

We employed the PID model to simulate the dynamic changes in insulin and glucose levels. The simulation results are as follows:

We employed the PID model based on delay differential equations to fit the data, as depicted by the curve. From the fitting results, our hypothesized model appears to have effectively captured the time-dependent changes in insulin and glucose concentrations. The rapid decline and eventual stabilization of insulin levels are highly consistent with the behavior set in our model, and the decreasing trend in glucose concentration also meets our expectations. Additionally, our model demonstrates good consistency with data from the literature, confirming its accuracy. Furthermore, the robustness of our model is also reflected in the R-squared values consistently exceeding 0.97. The analysis data is sourced from reference.

Most of the parameter values used in the model are derived from estimates. All the model application code has been uploaded to gitlab.

[1] Shiffman, Daniel, Shannon Fry, and Zannah Marsh. The nature of code. D. Shiffman, 2012.
[2] Yazawa, M., Sadaghiani, A. M., Hsueh, B. & Dolmetsch, R. E. Induction of protein-protein interactions in live cells using light. Nat Biotechnol 27, 941–945 (2009).
[3] Xie, M. et al. β-cell-mimetic designer cells provide closed-loop glycemic control. Science (New York, N.Y.) 354, 1296–1301 (2016).
[4] Ausländer, D. et al. A synthetic multifunctional mammalian pH sensor and CO2 transgene-control device. Molecular cell 55, 397–408 (2014).
[5] Li J, Kuang Y, Mason CC. Modeling the glucose-insulin regulatory system and ultradian insulin secretory oscillations with two explicit time delays. J Theor Biol. 2006 Oct 7;242(3):722-35. 
[6] Watson, E. M., Chappell, M. J., Ducrozet, F., Poucher, S. M. & Yates, J. W. T. A new general glucose homeostatic model using a proportional-integral-derivative controller. Computer Methods and Programs in Biomedicine 102, 119–129 (2011).
[7] Baggio, L. L. & Drucker, D. J. Biology of Incretins: GLP-1 and GIP. Gastroenterology 132, 2131–2157 (2007).
[8] Cheung, A. T. et al. Glucose-dependent insulin release from genetically engineered K cells. Science 290, 1959–1962 (2000).
[9] Marchetti G, Barolo M, Jovanovic L, Zisser H, Seborg DE. An improved PID switching control strategy for type 1 diabetes. IEEE Trans Biomed Eng. 2008 Mar;55(3):857-65. 
[10] Wang H, Li J, Kuang Y. Mathematical modeling and qualitative analysis of insulin therapies. Math Biosci. 2007 Nov;210(1):17-33. 
[11] Pacini G, Bergman RN. MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsivity from the frequently sampled intravenous glucose tolerance test. Comput Methods Programs Biomed. 1986 Oct;23(2):113-22.