Modeling of GLP-1 protein

---

I.Overview

Glucagon-like peptide-1 (GLP-1) is a peptide hormone primarily secreted by intestinal L cells. Its core physiological functions include promoting insulin secretion from pancreatic β-cells and inhibiting glucagon secretion from α-cells in a glucose-dependent manner, as well as regulating metabolic homeostasis by delaying gastric emptying and suppressing appetite. However, GLP-1 exhibits an extremely short half-life of approximately 2 minutes in blood circulation. Currently, the mechanism of action and structural characteristics of GLP-1’s primary degrading enzyme dipeptidyl peptidase-4 (DPP-4) are relatively well elucidated[1]. In contrast, research on its secondary degrading enzyme, neprilysin (NEP-24.11, hereafter abbreviated as NEP), remains relatively scarce[2]. Herein, we aim to optimize GLP-1 sequence via protein engineering to enhance its resistance to NEP-mediated degradation.

To enhance the plasma stability of GLP-1 against DPP-4 and NEP-24.11, we rationally adopted a novel variant, GLP-1Z[3], based on a systematic literature analysis. This variant incorporates key substitutions at positions A8H, E27K, and K34R within the native human-derived GLP-1(7-36) sequence.

To acquire more comprehensive site optimization data while avoiding the prolonged timelines and high costs associated with traditional high-throughput wet-lab screening, we employed a computational modeling and screening strategy. This approach focused on two goals: rationally selecting mutation sites to identify potential binding sites between GLP-1(Z) and NEP by analyzing their intermolecular hydrogen-bond interaction network, and performing systematic computational scanning of all GLP-1(Z) residues to calculate free energy changes upon non-canonical amino acid substitutions, thereby identifying positions that enhance protein performance and stability.

For structural analyses, we utilized the structure of wild-type GLP-1 (7–36) in complex with its receptor GLP-1R (PDB ID: 9A3F)[4], as well as the crystal structure of NEP (PDB ID: 2QPJ)[5]. The three-dimensional structure of GLP-1Z (7–35) was predicted using AlphaFold3.

II. Molecular Docking

Given the absence of a cocrystal structure for NEP-24.11 in complex with GLP-1 in the RCSB PDB database, we employed multiple molecular docking approaches to obtain optimal binding conformations. The three-dimensional structure corresponding to the GLP-1(Z) sequence was modeled and subjected to molecular docking analyses with both NEP-24.11 and GLP-1R, respectively.

1. Docking of GLP-1 with NEP

Based on the NEP crystal structure obtained from the RCSB PDB database, we initially performed molecular docking. However, no reliable binding conformation was obtained. This may be attributed to the "closed" or inactive state of the active site in this crystal structure, whose conformation is unsuitable for direct docking. We subsequently employed molecular dynamics simulations to sample more open conformations, which may be more suitable for docking analyses.

Figure 1. Structure of NEP (PDB: 2QPJ). The orange region indicates the predicted binding site for docking.

Based on the open conformational state of the protein screened through 300 ns molecular dynamics simulations, we reperformed molecular docking and systematically evaluated multiple platforms, including AutoDock, ClusPro 2.0, ZDOCK, HADDOCK, and HDOCK.

Figure 2. Maximum open conformation of NEP obtained from 300 ns molecular dynamics simulation.

Following multi-platform docking screening, two ideal docking conformations (Model I and Model II) were obtained from ClusPro 2.0 and HADDOCK. Subsequent 300 ns molecular dynamics simulations revealed that only Model I reached a stable state, while Model II exhibited continuous conformational fluctuations. Therefore, the stable Model I was selected for further analysis.

Figure 3. Predicted complex structures of NEP-24.11 and GLP-1 from different docking platforms.

A. Result from HADDOCK. NEP-24.11 is shown in blue, GLP-1 in pink, and the active site is highlighted in red. B. Result from ClusPro 2.0. NEP-24.11 is shown in yellow-green, GLP-1 in blue, and the active site is highlighted in green.

2. Docking of GLP-1Z with NEP

For the engineered GLP-1Z sequence, we performed molecular docking to analyze its complex structure with NEP-24.11, aiming to identify key interaction interfaces and provide theoretical guidance for selecting modification sites using non-canonical amino acids.

Figure 4. Molecular docking conformation of GLP-1Z with NEP-24.11. The optimal model was obtained using the ClusPro 2.0 platform. The active site of NEP-24.11 (pink) is highlighted in purple, and the GLP-1Z peptide chain is shown in green.

3. Docking of GLP-1Z with GLP-1R

To validate the biological function of the engineered GLP-1Z sequence, we predicted its three-dimensional structure using AlphaFold3 and simulated its interaction with the receptor GLP-1R through molecular docking. The results show that GLP-1Z forms a stable binding conformation with GLP-1R, providing key computational evidence that it retains the physiological activity of wild GLP-1.

Figure 5. Comparison of binding conformations between GLP-1/GLP-1Z and GLP-1R.

A. Crystal structure of the GLP-1–GLP-1R complex (PDB ID: 3C5T[6]). B. Predicted model of the GLP-1Z–GLP-1R complex (predicted by AlphaFold3 and optimized with ClusPro 2.0). GLP-1R, GLP-1, and GLP-1Z are colored light purple, deep purple, and gold, respectively.

III. Molecular Dynamics Simulations

Using the optimized NEP enzyme structure as the receptor, we performed molecular docking attempts with wild-type GLP-1 and GLP-1Z across multiple platforms. Successful docking simulations for both complexes were ultimately achieved using the HADDOCK platform[7].

Molecular dynamics simulations were performed using GROMACS 2023.3 on the docking complexes of both ligands with the NEP. Through analysis of structural stability and intermolecular interactions from the simulation results, optimal mutation sites capable of reducing the binding activity between GLP-1 and NEP were identified, with the ultimate goal of effectively prolonging the in vivo half-life of GLP-1.

The docking simulation system of NEP with GLP-1(Z) was simulated at 300 K. To begin with, the system is minimized by 50,000 steps using the steepest descent algorithm. In the pre-balancing phase, the system is gradually heated to 300 K at 1 atm using a 200 ps NVT set and 200 ps NPT set. A 50 ns production run followed to collect the balanced configuration for each 2 ps interval. If the system appeared unbalanced, the simulation was extended to 300 ns for further observation. A speed scaling thermostat 49 with a time constant equal to 0.1 ps was used to keep the temperature constant throughout the simulation. To maintain the pressure, the Berendsen pressure coupler was used in the NPT pre-balance run and the Parrinello-Rahman pressure coupler was used in the production run, with the pressure time constant and isothermal compression rate set at 2 ps and 4.5×10^-5 bar^-1, respectively. In the whole simulation process, the integral time step of the motion equation is 2 fs and the cut-off value for the non-bonding interaction was 12 Å. The particle grid Ewald algorithm 53 was used to calculate long range electrostatic interactions.

1. RMSD of NEP-GLP-1 Complex

The NEP-GLP-1 complex demonstrated remarkable structural stability during the 300 ns molecular dynamics simulation, providing a structural basis for its effective cleavage by NEP. The simulation trajectory revealed that the backbone RMSD of the complex remained stable within a range of 0.2–0.3 nm (average 0.24 ± 0.019 nm) throughout the process, with fluctuations of less than 0.5 Å after 100 ns, indicating no significant conformational shift in the binding pose.

Figure 6. Interaction analysis between NEP and GLP-1. A. RMSD plot of the NEP-GLP-1 complex during the 0-300 ns molecular dynamics simulation. B. Three-dimensional structure of the molecular docking results between NEP and GLP-1.

2. RMSD of NEP-GLP-1Z Complex

The simulation results indicated that GLP-1Z is expected to exhibit resistance to NEP-mediated proteolysis. After 300 ns of molecular dynamics simulation, the NEP-GLP-1Z complex showed a higher RMSD value (0.28±0.016 nm) compared to the NEP-GLP-1 complex. This reflects a decrease in the stability of its binding conformation, which is a key structural feature contributing to reduced NEP cleavage efficiency.

Figure 7. Interaction analysis between NEP and GLP-1Z. A. RMSD curve of the NEP-GLP-1Z complex during the 0–300 ns molecular dynamics simulation. B. Three-dimensional structure diagram of the molecular docking results between NEP and GLP-1Z.

3. Hydrogen Bonding Characteristics of the NEP-GLP-1 and NEP-GLP-1Z Complexes

The hydrogen bond network is also crucial for maintaining the conformation of the zinc-binding motif, the hydrophobic pocket, and mediating the specific catalytic mechanism of the NEP active site. To characterize changes in the binding mode between GLP-1(Z) and NEP, we analyzed a 300 ns MD simulation trajectory. A distance of ≤ 3.5 Å between donor and acceptor atoms was used as the criterion for hydrogen bond formation. The results indicated that H7, K26, and K34 of GLP-1 are the primary residues involved in hydrogen bonding with NEP. Given that position 7 is a well-established key residue for GLP-1 receptor interaction, our analysis suggests K26 and K34 as promising modifiable sites. In the engineered GLP-1Z, residues H8, E21, and K26 still exhibit a high propensity to form hydrogen bonds with NEP, facilitating binding and subsequent proteolysis. Therefore, further modifications at these sites could be explored to reduce affinity for NEP.

Figure 8. Hydrogen Bond Features of the NEP-GLP-1 and NEP-GLP-1Z Complexes. A. Critical hydrogen bonds at the NEP-GLP-1 binding interface. B. Critical hydrogen bonds at the NEP-GLP-1Z binding interface.

IV. Non-Canonical Amino Acid Library Screening and ΔΔG Validation

Based on the superior incorporation efficiency of 5-HTP compared to sTyr in EcN (refer to the Results section), we selected 5-HTP as the candidate molecule for computational screening using Rosetta. After obtaining its stable conformation through structural optimization and quantum mechanical (QM) calculations, we automatically generated Rosetta-compatible parameter files and a rotamer library via scripting, completing the necessary configurations to achieve successful incorporation of 5-HTP in GLP-1(Z).

To systematically evaluate the effects of 5-HTP incorporation on both GLP-1(Z) proteins' stability and their binding activity with the GLP-1R, we performed computational screening using both Cartesian ΔΔG and Flex ΔΔG protocols. The wild model of GLP-1Z (7-35) was predicted using AlphaFold3, while the crystal structures of wild GLP-1 (7-36) and GLP-1R were obtained from the PDB database (IDs: 3C5T and 6X18[4], respectively). We generated the initial GLP-1Z/GLP-1R complex by docking with ClusPro 2.0, followed by structural preprocessing to repair missing atoms and remove redundant information, ensuring structural standardization. Subsequently, The complex conformation was further optimized using MD simulations to resolve steric clashes and achieve thermodynamic stability. Finally, we calculated binding free energy changes (ΔΔG) for all possible 5-HTP scanning mutations to quantify their effects on protein stability and receptor binding.

1.Cartesian ΔΔG

Cartesian ΔΔG is a computational method for calculating changes in protein stability (ΔΔG) resulting from protein mutations, commonly implemented in the Rosetta software suite. The algorithm was developed by Phil Bradley and Yuan Liu, with related research findings published in journals such as JCTC. Compared to the traditional ddg_monomer protocols, Cartesian ΔΔG accommodates backbone flexibility through Cartesian-space optimization, eliminating the need for extensive backbone restraints and thereby allowing slight backbone movements.

The algorithm employs a two-step relaxation protocol: initial Fastrelax optimization identifies optimal side-chain conformations, followed by Cartesian-space Fastrelax to thoroughly relax the region near the mutation site. The energy difference is then evaluated and the ΔΔG values are accordingly corrected. This methodology offers faster computation and higher predictive accuracy, with improved Pearson correlation between predicted ΔΔG and experimental ΔΔG compared to ddg_monomer.

The preprocessed docking files of GLP-1_GLP-1R and GLP-1Z_GLP-1R, after removing GLP-1R, were utilized for Cartesian ddG calculations. Changes in protein stability was calculated as ΔΔG = ΔG_mut - ΔG_wt, where ΔG_mut and ΔG_wt represent mutant and wild-type free energies, respectively. A negative ΔΔG value indicates mutations that enhance protein stability. Final results for 5-HTP substitutions (denoted as X) are summarized in the table below.

Table 1. Cartesian ΔΔG Results of 5-HTP Scanning in GLP-1

Sample	ΔG_wt	ΔG_mut	ΔΔG
A24X	-84.8245	-84.8506	-0.0261
K26X	-84.757	-84.7596	-0.0026
S18X	-84.7043	-84.7056	-0.0013
I29X	-84.5713	-84.5721	-0.0008
E21X	-84.8288	-84.8293	-0.0005
V33X	-84.6135	-84.614	-0.0005
A30X	-84.7493	-84.7495	-0.0001
H7X	-84.3435	-84.3435	0
A8X	-84.4511	-84.4511	0
E9X	-84.3108	-84.3108	0
G10X	-84.449	-84.449	0
T11X	-84.7435	-84.7435	0
F12X	-84.7531	-84.7531	0
T13X	-84.9555	-84.9555	0
S14X	-84.6961	-84.6961	0
D15X	-84.7453	-84.7453	0
V16X	-84.7403	-84.7403	0
S17X	-84.6666	-84.6666	0
F28X	-84.7311	-84.7311	0
L32X	-84.566	-84.566	0
K34X	-84.7391	-84.7391	0
R36X	-85.873	-85.873	0
L20X	-84.8571	-84.857	0.0001
G22X	-84.8856	-84.8853	0.0003
Y19X	-84.8436	-84.8421	0.0015
A25X	-84.6625	-84.6608	0.0016
Q23X	-84.8486	-84.8468	0.0018
E27X	-84.9411	-84.9385	0.0026
G35X	-84.7036	-84.6961	0.0075
W31X	-84.7928	-84.6141	0.1786

Seven sites with ΔΔG < 0 were identified: A24, K26, S18, I29, E21, V33, and A30.

Table 2. Cartesian ΔΔG Results of 5-HTP Scanning in GLP-1Z

Sample	ΔG_wt	ΔG_mut	ΔΔG
E21X	-61.314	-61.4345	-0.1205
F28X	-61.0155	-61.112	-0.0965
L20X	-61.3185	-61.3703	-0.0518
R34X	-61.0486	-61.0863	-0.0376
A30X	-60.9166	-60.953	-0.0363
K27X	-61.0788	-61.1105	-0.0316
K26X	-61.0341	-61.06	-0.0258
G22X	-61.2405	-61.2535	-0.013
A24X	-61.099	-61.104	-0.005
F12X	-61.0661	-61.07	-0.0038
T13X	-61.3933	-61.3963	-0.003
T11X	-61.2211	-61.222	-0.0008
S14X	-61.374	-61.3746	-0.0006
G35X	-61.0816	-61.0823	-0.0006
H7X	-61.7971	-61.7971	0
H8X	-61.9836	-61.9836	0
E9X	-61.3193	-61.3193	0
Q23X	-61.2055	-61.2055	0
W31X	-60.9095	-60.9095	0
V33X	-60.8015	-60.8015	0
L32X	-60.8223	-60.8221	0.0001
V16X	-61.3806	-61.3798	0.0008
D15X	-61.4256	-61.424	0.0016
S17X	-61.342	-61.3403	0.0016
A25X	-61.0651	-61.063	0.0021
G10X	-61.2198	-61.2141	0.0056
S18X	-61.22	-61.1986	0.0213
I29X	-60.9678	-60.9293	0.0385
Y19X	-61.3928	-61.3441	0.0486

Fourteen sites with ΔΔG < 0 were identified: E21, F28, L20, R34, A30, K27, K26, G22, A24, F12, T13, T11, S14, and G35.

2. Flex ΔΔG

Flex ΔΔG is a computational method within the Rosetta software suite designed to evaluate the effects of protein mutations on stability or binding affinity. Its distinguishing feature is the allowance of flexible adjustments in residues surrounding the mutation site (including side chains and partial backbone regions) during energy optimization. Relevant methodological studies have been reported in journals such as Protein Science.

Unlike rigid ΔΔG calculations that rely on a fixed backbone conformation, Flex ΔΔG employs energy minimization and Monte Carlo sampling to explore low-energy conformational states near the mutation site. The calculation begins by constructing the initial mutant structure, followed by flexible relaxation of the local environment to search for the optimal conformation. The ΔΔG value is then derived from the optimized mutant and wild-type structures.

This approach more accurately reflects the dynamic characteristics of biomolecules and demonstrates superior predictive accuracy for scenarios involving local structural changes induced by mutations, particularly outperforming rigid models in the analysis of mutations at protein-protein interaction interfaces.

The preprocessed docking files of the GLP-1_GLP-1R and GLP-1Z_GLP-1R complexes were used as input for Flex ΔΔG calculations.

Table 3. Flex ΔΔG Results for 5-HTP Scanning in GLP-1

Sample	ΔΔG
E27X	-0.334663
A30X	-0.170607
A8X	-0.05427
H7X	0.200502
G10X	0.227868
K26X	0.272857
K34X	0.350991
R36X	0.35444
E9X	0.455663
V33X	0.489714
L32X	0.538231
W31X	0.5449
T11X	0.57746
G35X	0.590981
I29X	0.695001
F28X	0.753165
A25X	0.762696
E21X	0.867378
T13X	1.296619
G22X	1.501862
Q23X	1.548774
L20X	1.588415
F12X	1.654764
S18X	1.788207
S14X	1.799825
Y19X	1.84515
S17X	1.979547
D15X	2.118907
V16X	2.212899
A24X	2.294355

Three sites with ΔΔG < 0 were identified: E27, A30, and A8.

Table 4. Flex ΔΔG Results for 5-HTP Scanning in GLP-1Z

Sample	ΔΔG
H7X	-0.195403
R34X	-0.10771
G10X	-0.094491
L20X	-0.088267
H8X	-0.081194
W31X	-0.05782
I29X	-0.054197
V33X	-0.019801
Q23X	-0.004053
G35X	-0.0011
E9X	0.013912
K27X	0.025279
L32X	0.053631
F28X	0.055388
A30X	0.069835
K26X	0.241896
A24X	0.333949
V16X	0.385698
S17X	0.445652
E21X	0.500679
D15X	0.531144
A25X	0.61103
G22X	0.87842
Y19X	0.895563
F12X	1.057014
S14X	1.06598
S18X	1.781751
T11X	2.135362
T13X	2.245099

Ten sites with ΔΔG < 0 were identified:H7、R34、G10、L20、H8、W31、I29、V33、Q23、G35.

Based on the results from both Cartesian ΔΔG and Flex ΔΔG calculations, we propose that the L20 and R34 sites in GLP-1Z hold promise for optimization via 5-HTP mutation. Both sites exhibited negative ΔΔG values in both computational methods, suggesting that the mutations could enhance both the intrinsic stability of the peptide and its binding affinity to GLP-1R.

Figure 9. Heatmap of Rosetta ΔΔG analysis of 5-HTP incorporation in GLP-1(Z) with Cartesian ΔΔG and Flex ΔΔG.

Based on the Rosetta software platform, we systematically evaluated the mutational effects of 5-HTP in both GLP-1 and GLP-1Z using Cartesian ΔΔG and Flex ΔΔG methods. This screening identified promising mutation sites (L20 and R34 in GLP-1Z) that are predicted to enhance both structural stability and receptor binding affinity. The negative ΔΔG values at these sites suggest that 5-HTP substitution is conformationally favorable. Future work should focus on validating the actual function of these mutants through wet-lab experiments and exploring possible synergistic effects of multi-site combinatorial mutations to further optimize the drug-like properties of GLP-1-based proteins.

V. Conclusions and Outlook

This study integrated molecular dynamics (MD) simulations with Rosetta free energy calculations to evaluate the conformational stability and binding properties of GLP-1(Z), and to identify potential sites for incorporating the non-canonical amino acid 5-HTP. MD simulations revealed that the GLP-1Z–NEP complex exhibited higher conformational flexibility, suggesting reduced binding stability that may contribute to resistance against NEP-mediated proteolysis. Notably, the MD results were not entirely consistent with the subsequent Rosetta free energy calculations, indicating methodological discrepancies that require integrated interpretation.

Rosetta-based Cartesian ΔΔG and Flex ΔΔG scanning screened for favorable mutations in GLP-1Z and yielded two high-potential sites, L20 and R34, both with negative ΔΔG values. These sites are predicted to enhance both protein stability and binding affinity to GLP-1R. Future work should include dry-lab testing of multi-site combinatorial mutations and wet-lab validation of selected mutants. Although experimental verification was not completed within the project timeline, our project provides a clear direction for the rational optimization of GLP-1(Z). Future work could employ rational engineering at these validated sites to develop GLP-1-based peptides with enhanced therapeutic potential.

Modeling of Controlled Release Kinetics of Hydrogel Microspheres

---

I.Overview

Hydrogel is a three-dimensional network structure formed by hydrophilic polymer chains through physical interaction or chemical bonding. It can produce different responses when stimulated by different environmental stimuli, so it is widely used in the field of medicine as a drug controlled release carrier . For intestinal microbial delivery, hydrogels show great advantages due to their mild preparation conditions, good biocompatibility and targeted release capability. Current research focuses on the physical embedding of probiotics in hydrogel networks to achieve the protection and Targeted release of probiotics[8].

Drug release in hydrogels is a complex process involving multiple mechanisms such as diffusion and swelling. Current research on hydrogel drug delivery models mainly focuses on small molecule drugs. Due to the small size of Escherichia Coli, some models based on controlled release of small molecules can also be applied to predict the release behavior of Escherichia Coli under certain conditions. We hope to make our work more concrete through the tools of mathematical modeling, and provide theoretical support and direction guidance in the preparation and characterization of microspheres.

II.Diffusion-controlled model of spherical hydrogels

In order to quantitatively describe the release behavior of probiotics from hydrogel microspheres, an idealized theoretical model was first established. This model is based on the following basic assumptions: 1. The hydrogel microspheres are homogeneous and isotropic perfect spheres, and do not change their geometry and size due to swelling and degradation during the release of probiotics. 2. E. coli enters the external environment by simple diffusion. 3. Probiotics were uniformly distributed in the hydrogel with an initial concentration of $C_0\left(C\right(r,0)=C_0)$, and maintained a good impregnation at the interface, that is, the release space of microorganisms was infinite, and the concentration of probiotics on the surface and outside the hydrogel was always 0.

Since E. coli is relatively small, its migration behavior in the gel network can be compared to the diffusion process of small molecules. Can be used in three-dimensional rectangular coordinate system based on the principle of molecular free diffusion Fick's second law partial $\frac{\partial C}{\partial t}=D(\frac{{\partial }^{2}C}{\partial x^2}+\frac{{\partial }^{2}C}{\partial y^2}+\frac{{\partial }^{2}C}{\partial z^2})$ ($D$ is the diffusion coefficient, $\frac{\partial C}{\partial t}$ is the rate of change in the concentration of the diffuser with time).

Considering the spherical symmetry of the carrier, the calculation can be simplified by converting this equation to the spherical coordinate system. In this case, the diffusion equation can be written as follows:

$$\frac{M_t}{M_{\infty }}=1-\frac{6}{{\pi }^2}\sum _{n=1}^{\infty }\frac{1}{n^2}{e}^{-\frac{D{n}^2{\pi }^{2}t}{R^2}}$$

Where $M_t$ is the cumulative amount of probiotics released, $M_{\infty }$ is the maximum amount of probiotics that can be released by the carrier, $R$ is the radius of the hydrogel microspheres (in mm), and t is the release time (in hours).

The diffusion coefficient $D$ needs to be measured experimentally. After the $D$ value was measured, the release kinetics curves of the same hydrogel material with different microsphere sizes could be predicted using the above analytical solution, thereby providing a theoretical basis for the carrier design and the regulation of release behavior.

III.Swelling-controlled model of spherical hydrogels

Drug release from hydrogels is often controlled by more than simple diffusion. Hydrogels tend to have good swelling properties. After water absorption and expansion, the internal pores of the hydrogel increased, and the embedded probiotics were more likely to be released. At this point, the diffusion coefficient $D$ and the radius $R$ of the hydrogel microspheres are both functions of time. The above diffusion model is difficult to obtain an analytical solution. The swelling model starts from an entropy perspective to better explain the swelling controlled process of hydrogel drug release.

Current microscopic modeling of hydrogel drug release during swelling can be based on three theories: rubberlike elasticity theory, equilibrium swelling theory and mesh size theory[9]. According to Rubberlike elasticity theory, the elasticity of rubber-like substances is mainly due to the change of entropy rather than the change of internal energy. The swelling of the hydrogel comes from the change of entropy, which is based on the movement of chain segments in its cross-linking network, and its properties determine the elasticity of the gel. The equilibrium swelling theory suggests that the equilibrium swelling degree of hydrogels is determined by the balance between the osmotic pressure-driven swelling affinity and the elastic shrinkage affinity of the polymer network. The swelling affinity is determined by the entropy increase of water molecules entering the hydrogel network spontaneously, while the shrinkage affinity is caused by the entropy decrease caused by chain segment motion, which can be explained by rubberlike elasticity theory. The mesh size theory simplifies the microstructure of hydrogels to a porous structure composed of grids. The mesh size of hydrogels can be estimated by comprehensive use of the theoretical model, so as to calculate the swelling performance and release potential of hydrogels.

The swelling process of hydrogels is complex and depends on factors such as its polymer network structure, cross-linking density, and environmental conditions. Simulations are generally carried out after experimental measurements combined with the above theories. The contraction force and mesh size of sodium alginate gel microspheres used in this experiment are greatly affected by the concentration of sodium alginate and cross-linking agent Ca²⁺, while the swelling affinity is mainly affected by the environmental solution[10]. There are two empirical formulas for the shrinkage/swelling ratio:

$$f_{sw}=e^{-\frac{a_{sw}}{t}}$$
$$f_{sh}=e^{-\frac{a_{sh}}{t}}$$

where $a_{sw}$ is swelling affinity (in minutes), $a_{sh}$ is shrinkage affinity (in minutes), and $t$ is the time spent by the sample in water (in minutes).

By measuring the swelling curve of the hydrogel, the swelling affinity of the microspheres and the shrinkage affinity in the corresponding liquid can be measured.

IV.Test of experimental data

The Peppas model[11] is one of the most widely used models in the field of drug controlled release. Before the majority of the drug is released completely, this model describes the cumulative release behavior of the drug from the polymer carrier with a simple power function relationship:

$$\frac{M_t}{M_{\infty }}=k{t}^n$$

where $n$ is the release exponent, $k$ is the kinetic constant, and $t$ is the release time in hours. This equation can be linearized as follows:

$$\text{lg}\left(\frac{M_t}{M_{\infty }}\right)=\text{lg}k+n\cdot \text{lg}t$$

In general, Peppas model is quite accurate if $\frac{M_t}{M_{\infty }}<\mathrm{60}\%$. The Peppas model was used to linearly fit the experimental data for alginate microspheres release (see our Results page for details), which yielded n=4.16 and an R² of 0.988. Since the release index is significantly greater than 1, it can be speculated that the release process of the hydrogel prepared in this experiment is mainly controlled by swelling.

Figure 10. Peppas model curve of cross-linked hydrogel with 1% SA and 0.1 mol/L CaCl₂

In addition, the swelling rate of the hydrogel crosslinked with 1% SA and 0.1mol/L CaCl₂ was determined by the drying method. After swelling the hydrogel microspheres fully in sterile PBS for 1h, the microspheres were removed and weighed. The mass of the microspheres was measured by drying them in an oven at 55 ° C for 48h until the mass of the microspheres no longer changed. The formula $\mathrm{SW}=\frac{W_p}{W_d}$ (where $W_d$ is the weight of the dry hydrogel, and $W_p$ is the weight of the swollen hydrogel after immersion.) The swelling rate $\mathrm{SW}$ of the hydrogel was calculated.

The dry weight of 4.20g fully swelled hydrogel microspheres was only 0.04g, and the $\mathrm{SW}$ exceeded 10000%. This indicates that the microspheres have a strong swelling ability, which is more conducive to the release of probiotics. Due to the limitation of experimental conditions, we were not able to test the swelling curve of the lyophilized gel, but it can be speculated that the lyophilized microspheres should have a higher release index, mainly swelling release.

V.Conclusion

We concluded that the release of probiotics from the hydrogel microspheres in our project was primarily governed by a swelling-controlled mechanism. However, we also observed a clear tendency for bacterial proliferation within the hydrogel matrix. This indicates that a comprehensive model must incorporate the dynamics of bacterial growth on the porous material instead of relying solely on either a diffusion-controlled or swelling-controlled model.

Due to the lack of instruments for material characterization in the laboratory, we were unable to characterize the physical properties of the hydrogel microspheres, which resulted in the inability to calculate data such as mesh size and conduct in-depth swelling control model calculations. However, the model indicates that the shrinkage affinity and the cross-link density can be adjusted by adjusting the concentration of reactants and cross-linking agents to control the release curve. In addition, we can also use electrostatic generation method, phase separation method, multi-phase emulsion method and other methods and related experimental instruments to make microspheres with more uniform geometry and larger dynamic radius, so as to conveniently and accurately control the targeted release performance of hydrogels. We look forward to further enriching the experimental and theoretical modeling content when subsequent conditions permit, and constructing high-performance hydrogel microspheres suitable for our strain.

Calculations Based on the Logistic Growth Model and Application of Turbidostats in GlucoXpert

---

Overview

To evaluate the growth performance of our engineered E. coli for its application in a turbidostat-based in vitro system, we characterized its growth curve and performed material balance analyses. These modeling experiments provided key parameters, including the maximum specific growth rate ${\mu }_{\max }$ and maximum biomass $X\mathrm{_max}$ under simulated experimental conditions. The resulting data offer a critical foundation for scaling up the strain for industrial production and, ultimately, for achieving its long-term colonization and drug-induced synthesis in vivo.

Our methodology included determining the growth curve of the strain (exemplified by the AND-gate circuit strain), performing material balance calculations under batch cultivation, and simulating and computing its behavior under turbidostatic conditions. By integrating the characterization of genetic circuits with the implementation of corresponding hardware, we validated combinational data of dry-lab and wet-lab experiments, thereby guiding the practical application of our project.

I. Core Foundation: Logistic Cell Growth Model

The Logistic model assumes that cell growth is limited by environmental factors. The specific growth rate decreases as cell concentration increases, eventually reaching the maximum cell concentration ($X\mathrm{_max}$).

1. Integral Form of Cell Concentration vs. Time (Core of Growth Curve)

The essence of the Logistic model is to describe the change in cell concentration $X\left(t\right)$ with time $t$[12]. The integral form (directly applicable for fitting growth curves) is:

$$X\left(t\right)=\frac{X_0\cdot X_{\max }}{X_0+(X_{\max }-X_0)\cdot e^{-{\mu }_{\max }\cdot t}}$$

Parameter Definitions:

$X_0$: Initial cell concentration (cell quantity at $\mathrm{t\; =\; 0}$, determined experimentally or predefined);

$X_{\max }$: Maximum cell concentration (cell quantity in the stationary phase, a key parameter for fitting growth curves);

${\mu }_{\max }$: Maximum specific growth rate (units: $\mathrm{1/h}$, reflecting the potential growth capacity of cells, obtained by fitting);

$t$: Cultivation time (units: $h$).

Application:

To determine key growth parameters in batch culture, we prepared media with varying concentrations of the limiting substrate (for LB medium, this refers to the combined mass concentration of tryptone and yeast extract) and monitored the change in biomass (measured as OD₆₀₀) over time for each initial concentration. The data were fitted using the built-in logistic growth model in GraphPad Prism to derive the maximum biomass ($X_{\max }$) and maximum specific growth rate (${\mu }_{\max }$) , providing essential parameters for subsequent material balance calculations.

Figure 11. Growth curve fitting of E. coli carrying the AND-gate plasmids under a limiting substrate concentration of 15 g/L.

Cultures were inoculated at an initial OD₆₀₀ of 0.02 and incubated at 37 °C to generate the growth curve shown in Figure 11. Using the 15 g/L substrate condition as an example, the fitting yielded a high correlation (R² = 0.9961) with ${\mu }_{\max }$ = 1.394 h^-1 and $X_{\max }$ = 1.011.

2.Specific Growth Rate (μ)

The specific growth rate $\mu$ (units: $\mathrm{1/h}$) represents the "growth rate per unit cell" and reflects the current growth state of cells. It is derived by differentiating the Logistic model:

$$\mu =\frac{1}{X}\cdot \frac{dX}{dt}={\mu }_{\max }\cdot (1-\frac{X}{X_{\max }})$$

Physical Meaning:

When $X\ll X_{\max }$, $\mu \approx {\mu }_{\max }$ (corresponding to the exponential growth phase); when $X\to X_{\max }$, $\mu \to 0$ (corresponding to the stationary phase), which is consistent with actual growth patterns.

3. Expression for Cell Growth Rate (dX/dt)

This directly describes the rate of change in total cell quantity, derived by rearranging the above differential equation:

$$\frac{dX}{dt}={\mu }_{\max }\cdot X\cdot (1-\frac{X}{X_{\max }})$$

This equation serves as the core driving term for subsequent cell material balance calculations.

II. Material Balance in Batch Cultivation (Basic Scenario)

In batch cultivation, the reactor volume $V$ is approximately constant (evaporation and sampling losses are neglectable). The material balance follows the fundamental principle: Rate of Accumulation = Rate of Input – Rate of Output + Rate of Generation – Rate of Consumption. For a batch system, where no material enters or leaves during cultivation, this simplifies to: Rate of Accumulation = Rate of Generation – Rate of Consumption. The material balances for biomass, substrate, and product are derived as follows:

1. Cell Material Balance

In batch cultivation, cellular mass is generated solely through growth and is not subject to consumption (cell death is neglected. If significant, a death term a death term $k_{d}X$ can be incorporated). The material balance equation is:

$$\frac{d(X\cdot V)}{dt}=V\cdot \frac{dX}{dt}$$

Given the assumptions of constant volume $V$ and substrate-limited growth, the expression for the specific growth rate $\mu$ is given by the logistic model, thus simplifying to:

$$\frac{dX}{dt}={\mu }_{\max }\cdot X\cdot (1-\frac{X}{X_{\max }})$$

This is identical to the Logistic growth rate equation, verifying the self-consistency of the model.

If cell death is considered , a specific death rate $k_d$ (units: $\mathrm{1/h}$) is introduced, modifying the equation to:

$$\frac{dX}{dt}={\mu }_{\max }\cdot X\cdot (1-\frac{X}{X_{\max }})-k_d\cdot X$$

2. Substrate Material Balance

Substrates are mainly consumed for cell growth, with an additional portion potentially utilized for maintenance metabolism (which can be incorporated via a maintenance term). The core of the balance is: Substrate consumption rate = Substrate demand rate for cell growth.

First, the yield coefficient yield coefficient $Y_{X/S}$ (units: $g$ cells/$g$ substrate)is defined as the mass of cells produced per mass of substrate consumed, and is a key experimentally determined parameter (like, $Y_{X/S}\approx 0.5 g/g$ for E. coli utilizing glucose).

Neglecting maintenance metabolism, the substrate balance is:

$$\frac{d(S\cdot V)}{dt}=-V\cdot \frac{1}{Y_{X/S}}\cdot \frac{dX}{dt}$$

Assuming constant volume $V$, this simplifies to:

$$\frac{dS}{dt}=-\frac{1}{Y_{X/S}}\cdot {\mu }_{\max }\cdot X\cdot (1-\frac{X}{X_{\max }})$$

When maintenance metabolism is considered, representing substrate consumed for cellular maintenance independent of growth, a maintenance coefficient $m_S$ (units: $g$ substrate/($g$ cells·$h$)) is introduced. The substrate consumption rate then becomes:

$$\frac{dS}{dt}=-(\frac{\mu }{Y_{X/S,true}}+m_S)\cdot X$$

Where:

$Y_{X/S,true}$: True yield coefficient (yield considering only growth-related consumption, slightly higher than $Y_{X/S}$);

$\mu ={\mu }_{\max }(1-X/X_{\max })$: Logistic-specific growth rate.

Initial-Final Relationship:

Integrating the substrate balance from $\mathrm{t\; =\; 0}$ to $t$) relates the initial substrate concentration $S_0$, residual substrate concentration $S\left(t\right)$, and cell concentration $X\left(t\right)$:

$$S_0-S\left(t\right)=\frac{X}{\left(}+m_S\cdot {\int }_0^{t}X\left(\tau \right)d\tau$$

The integral term ${\int }_0^{t}X\left(\tau \right)d\tau$ represents the cumulative "cell concentration-time integral," accounting for total maintenance-associated substrate consumption.

3. Product Material Balance

Product formation (like ethanol, antibiotics and GLP-1) has two distinct relationships with cell growth, requiring different formulas for description:

(1) Growth-Associated Products (like Ethanol)

Product synthesis occurs concurrently with cell growth. This type of formation is modeled using the product yield coefficient $Y_{P/X}$ (units: $g$ product / $g$ cells):

$$\frac{dP}{dt}=Y_{P/X}\cdot \frac{dX}{dt}$$

Substituting the Logistic growth rate gives:

$$\frac{dP}{dt}=Y_{P/X}\cdot {\mu }_{\max }\cdot X\cdot (1-\frac{X}{X_{\max }})$$

Integrating with initial product concentration $P_0=0$ gives:

$$P\left(t\right)=Y_{P/X}\cdot \left[X\right(t)-X_0]$$

(2) Non-Growth-Associated Products (like Penicillin)

Product formation is decoupled from active growth and may continue during the stationary phase. It is described by a specific product formation rate $q_P$ (units: $g$ product / $g$ cells·$h$), often considered constant:

$$\frac{dP}{dt}=q_P\cdot X$$

Substituting the Logistic expression for $X\left(t\right)$ and integrating gives:

$$P\left(t\right)=q_P\cdot {\int }_0^{t}X\left(\tau \right)d\tau$$

The integral term must be evaluated either analytically or numerically based on the logistic expression for $X\left(\tau \right)$.

III.Extended Scenarios

In practical bioprocessing, cultivation volume often varies, necessitating modifications to the classical balance equations. Furthermore, the standard Logistic model can be extended with correction factors to better represent complex physiological conditions.

1. Fed-Batch Cultivation (Volume Variation)

In fed-batch cultivation, a substrate solution is fed into the reactor at a time-dependent feed rate $F\left(t\right)$ (units: $\mathrm{L/h}$). The system volume evolves as:

$$V\left(t\right)=V_0+{\int }_0^{t}F\left(\tau \right)d\tau$$

where $V_0$ is the initial volume.

The biomass balance must account for potential cell input from the feed stream:

$$\frac{d(X\cdot V)}{dt}=\mu \cdot X\cdot V-F\left(t\right)\cdot X_{feed}$$

Here, $X_{feed}$ denotes the cell concentration in the feed (typically zero).

The corresponding substrate balance becomes:

$$\frac{d(S\cdot V)}{dt}=F\left(t\right)\cdot S_{feed}-(\frac{\mu }{Y_{X/S,true}}+m_S)\cdot X\cdot V$$

where $S_{feed}$ represents the substrate concentration in the feed.

2. Extensions of the Logistic Model (Correction)

When the basic Logistic model fails to fit experimental data, correction terms can be introduced:

Product inhibition correction:

$$\mu ={\mu }_{\max }\cdot (1-\frac{X}{X_{\max }})\cdot (1-\frac{P}{P_{\max }})$$

where $P_{\max }$ is the maximum tolerated product concentration.

Substrate limitation correction:

$$\mu ={\mu }_{\max }\cdot \frac{S}{K_S+S}\cdot (1-\frac{X}{X_{\max }})$$

Here, combined with the Monod model, $K_S$ denotes the substrate half-saturation constant.

3.Turbidostat Calculations(Monod Model)

A turbidostat is a type of continuous cultivation system that maintains a constant cell concentration ($X=\mathrm{constant}$) dynamically adjusting the dilution rate $D$ (units: $\mathrm{1/h}$). Unlike chemostats that fix $D$ to control substrate levels, turbidostats employ optical density feedback: when $X$ exceeds the set value, the dilution rate increases to dilute cells back to the set levels; when $X$ is too low, $D$ decreases to permit growth. The core equations and operating principles are outlined below[13].

(1) Basic Assumptions for Turbidostats

Steady-state operation: $\frac{dX}{dt}=0$ and $\frac{dS}{dt}=0$ (cell and substrate concentrations are constant over time);

Constant volume: $V=\mathrm{constant}$ (feed rate $F$ equals effluent rate $F_{eff}$);

No cell death: $k_d=0$ (negligible for steady-state growth);

Feed contains no cells: $X_{feed}=0$.

(2) Core Definitions for Turbidostats

Dilution rate ($D$): Ratio of feed rate to reactor volume, $D=\frac{F}{V}$ (units:$\mathrm{1/h}$). This is the key adjustable parameter in turbidostats;

Set cell concentration ($X_s$): Target cell density maintained by the turbidostat (predefined based on process needs);

Specific growth rate adaptation: In turbidostats, $\mu$ is no longer described by the basic Logistic model (which assumes $X$ changes with time). Instead, $\mu$ is linked to substrate concentration via the Monod model (more accurate for continuous substrate-limited growth):

$$\mu ={\mu }_{\max }\cdot \frac{S}{K_S+S}$$

where $K_S$ (units: $\mathrm{g/L}$) is the Monod half-saturation constant (substrate concentration at which $\mu =\frac{1}{2}{\mu }_{\max }$)[14].

(3) Steady-State Material Balance for Turbidostats

Cell Balance (Key for Turbidostat Operation) At steady state, the rate of cell accumulation equals zero. The cell balance equation is derived as follows: Total cell input (feed) + Total cell generation (growth) = Total cell output (effluent) + Total cell loss (death). Substituting assumptions ($X_{feed}=0$, $k_d=0$, $F=F_{eff}$):

$$0+\mu \cdot X\cdot V=D\cdot X\cdot V+0$$

Simplifying as (divide both sides by $X\cdot V$, $X\ne 0$):

$$\mu =D$$

Critical Interpretation:

In a turbidostat, the specific growth rate of cells $\mu$ is exactly equal to the dilution rate $D$ at steady state. The system automatically adjusts $D$ to match $\mu$, ensuring $X$ remains at the set value $X_s$.

Substrate Balance (Linked to Cell Concentration) At steady state, substrate accumulation is zero. The substrate balance accounts for three terms: substrate input (feed), substrate consumption (growth + maintenance), and substrate output (effluent):

$$D\cdot S_{feed}=(\frac{\mu }{Y_{X/S,true}}+m_S)\cdot X+D\cdot S$$

Substitute $\mu =D$ (from cell balance) and rearrange to solve for the residual substrate concentration $S$ in the reactor:

$$S=S_{feed}-\frac{X}{Y_{X/S,true}}-\frac{m_S}{D}\cdot X$$

Special Case: No Maintenance Consumption as $m_S=0$, like rapid growth phase, the equation simplifies to:

$$S=S_{feed}-\frac{X}{Y_{X/S}}$$

This shows that residual substrate $S$ decreases as the set cell concentration $X$ increases, consistent with the trade-off between cell density and substrate utilization.

Product Balance (Growth-Associated and Non-Growth-Associated)

For turbidostats, product balance follows the same classification as batch cultivation, but steady-state conditions simplify the rate calculations.

1. Growth-Associated Products (like Ethanol)

Product formation is tied to cell growth. At steady state, the product formation rate equals the product effluent rate:

$$Y_{P/X}\cdot \mu \cdot X=D\cdot P$$

Substitute $\mu =D$ and solve for steady-state product concentration $P$:

$$P=Y_{P/X}\cdot X$$

2. Non-Growth-Associated Products (like Penicillin)

Product formation depends on total cell mass (not growth rate). At steady state:

$$q_P\cdot X=D\cdot P$$

Solve for $P$:

$$P=\frac{q_P}{D}\cdot X$$

Key Insight:

For non-growth-associated products, $P$ is inversely proportional to $D$, lowering the dilution rate to retain cells longer increases product concentration, which is critical for antibiotic fermentation.

(4) Practical Calculations for Turbidostats

HiZJU-China developed a turbidostat model based on the aforementioned theories using E. coli carrying two AND-gate plasmids cultivated in LB medium. Growth curves were obtained at 37 ℃ with initial OD₆₀₀ = 0.02 across LB concentrations gradients of 3–21 g/L.

Figure 12. Growth curves of E. coli carrying AND-gate plasmids in LB medium at varying substrate concentrations of 3–21$\mathrm{g/L}$.

After analyzing this growth curve, we generated a $\mathrm{\mu -S}$ saturation curve (A) and a Lineweaver-Burk double reciprocal plot (B) to determine the parameters[15] $K_s$ and ${\mu }_m$, as shown in Figure 13.

Figure 13. Kinetic analysis of bacterial growth using the Monod model. (A) $\mathrm{\mu -S}$ saturation curve. (B) Lineweaver-Burk double-reciprocal plot.

Through fitting, we obtained the values: $K_s$ = 3.696 ($\mathrm{g/L}$) and ${\mu }_{\max }$ = 1.588 ($\mathrm{1/h}$).

Known Parameters:

Feed substrate concentration: $S_{feed}=15 g/L$

Set cell concentration: $X_s=1.4 \left({\mathrm{the\; value\; of\; OD}}_{600}\right)$

Cell yield coefficient: $Y_{X/S}=0.3 g/g$

Estimated Product yield coefficient: $Y_{P/X}=0.01 g/g$

Maximum specific growth rate: ${\mu }_{\max }=1.588 1/h$

Monod half-saturation constant: $K_S=3.696 g/L$

eactor volume: $V=50 \mathrm{mL}$

Calculation Steps:

1. Calculate residual substrate $S$ (neglect maintenance):

$$S=S_{feed}-\frac{X_s}{Y_{X/S}}=15-\frac{1.4}{0.3}=10.33 g/L$$

2. Calculate specific growth rate $\mu$ (Monod model):

Since $S\gg K_S$ (even if $S=0$ in this case, we use $\mu \le {\mu }_{\max }$):

$$\mu ={\mu }_{\max }\cdot \frac{S}{K_S+S}\approx 1.274 1/h$$

3. Calculate dilution rate $D$ (steady-state condition $\mu =D$):

$$D=\mu =1.274 1/h$$

4. Calculate feed rate $F$ ($D=\frac{F}{V}$):

$$F=D\cdot V=1.674\cdot 0.05=0.0637 L/h$$

5. Calculate steady-state product concentration $P$ (growth-associated):

$$P=Y_{P/X}\cdot X_s=0.01\cdot 1.4=0.014 g/L$$

Result Summary:

Estimate based on the above data, the turbidostat operates at a feed rate of 0.0637 L/h, producing 0.014 $\mathrm{g/L}$ GLP-1.

(5) Comparison Between Turbidostats and Chemostats

To clarify the operational rationale behind our turbidostat-based approach, we compare it against the industrial-standard chemostat across several key parameters. This analysis helps elucidate our strategy for in vitro simulation of dynamic gut microbiota homeostasis and identifies the critical control parameters required.

Parameter	Turbidostat	Chemostat
Controlled Variable	Cell concentration $X$ (constant)	Dilution rate $D$ (constant)
Adjustable Parameter	Dilution rate $D$ (varies to maintain $X$)	None (fixed $D$; $X$ and $S$ adjust to steady state)
Steady-State $\mu$ vs. $D$	$\mu =D$ (always true)	$\mu =D$ (true, but $D$ is fixed)
Residual Substrate $S$	Varies with $X$set (higher $X$set → lower $S$)	Fixed (depends on $D$; higher $D$ → higher $S$)
Application Scenarios	High-cell-density cultivation Processes requiring constant $X$ (like, enzyme production) Studies on cell physiology at fixed density	Substrate-limited growth studies Steady-state product screening Large-scale fermentation with fixed $D$ (like ethanol)
Stability	Less stable (prone to fluctuations if $D$ adjusts too slowly)	More stable (fixed $D$ reduces feedback delays)

IV. Key Parameter Acquisition and Formula Application Process

The systematic acquisition of kinetic and stoichiometric parameters follows this workflow:

1. Experimental Data Collection: Measure time-course data during batch cultivation, including $X\left(t\right)$ (cell concentration, converted from OD₆₀₀ measurements), $S\left(t\right)$ (substrate concentration), and $P\left(t\right)$ (product concentration);




2. Logistic Growth Model Fitting: Fit the integrated logistic equation to the $X\left(t\right)$ data:


$$X\left(t\right)=\frac{X_0{X}_{\max }}{X_0+(X_{\max }-X_0)e^{-{\mu }_{\max }t}}$$


to determine the maximum biomass $X_{\max }$ and the maximum specific growth rate ${\mu }_{\max }$;




3. Stoichiometric Parameter Calculation: Use $S_0-S\left(t\right)$ and $X\left(t\right)-X_0$ to fit $Y_{X/S}=\frac{X\left(t\right)-X_0}{S_0-S\left(t\right)}$. Similarly, obtain $Y_{P/X}$ or $q_P$;




4. Model Validation and Application: Substitute parameters into material balance formulas to calculate $\frac{dX}{dt}$, $\frac{dS}{dt}$, or turbidostat $D$/$P$ and compare the predictions with experimental results to modify the model.

V. Summary Formula Table (Core Summary)

Type	Formula	Key Parameters
Logistic growth curve	$X\left(t\right)=\frac{X_0{X}_{\max }}{X_0+(X_{\max }-X_0)e^{-{\mu }_{\max }t}}$	$X_0,X_{\max },{\mu }_{\max }$
Specific growth rate (Logistic)	$\mu ={\mu }_{\max }(1-\frac{X}{X_{\max }})$	${\mu }_{\max },X,X_{\max }$
Cell growth rate	$\frac{dX}{dt}={\mu }_{\max }X(1-\frac{X}{X_{\max }})$	Same as above
Substrate consumption rate (no maintenance)	$\frac{dS}{dt}=-\frac{1}{Y_{X/S}}{\mu }_{\max }X(1-\frac{X}{X_{\max }})$	$Y_{X/S}$
Growth-associated product formation rate	$\frac{dP}{dt}=Y_{P/X}{\mu }_{\max }X(1-\frac{X}{X_{\max }})$	$Y_{P/X}$
Non-growth-associated product formation rate	$\frac{dP}{dt}=q_{P}X$	$q_P$
Turbidostat steady-state $\mu$ vs $D$	$\mu =D$	$\mu ,D$
Turbidostat residual substrate	$S=S_{feed}-\frac{X}{Y_{X/S,true}}-\frac{m_S}{D}\cdot X$	$S_{feed},X,Y_{X/S},m_S,D$
Turbidostat product (non-growth-associated)	$P=\frac{q_P}{D}\cdot X$	$q_P,D,X$

Through the above formula system, the biochemical reaction process based on the Logistic model can be fully described, including batch cultivation dynamics and turbidostat steady-state behavior. This provides a complete quantitative framework for analyzing cell growth, substrate consumption, and product formation, serving as core mathematical tools for reactor design, process optimization, and scale-up in biochemical engineering.

Conclusions

Through the simulation of bacterial growth kinetics, we have characterized the growth behavior of engineered E. coli carrying genetic circuits for GLP-1 synthesis and identified key parameters for both batch and turbidostat cultivation. These findings provide a critical foundation for optimizing turbidostat design and support the future goal of applying this system to intestinal microbial colonization.

Based on the determined maximum biomass, we can rationally design the hardware feeding strategy and predict the steady-state biomass under turbidostat operation. Furthermore, the model allows us to estimate the production levels of growth-associated products such as GLP-1 and GFP under induction, facilitating the identification of optimal induction timing and dosage.

By integrating dry and wet lab approaches, we have established a deeper understanding of the physiological properties of the engineered strain. This knowledge serves as a essential basis for subsequent in vitro co-culture simulations and will ultimately inform strategies for achieving effective in vivo colonization and function.

References

---

[1]Gare, C. L., White, A. M., & Malins, L. R. (2025). From lead to market: chemical approaches to transform peptides into therapeutics. Trends in biochemical sciences, 50(6), 467–480.

[2]Tiraboschi, G., Jullian, N., Thery, V., Antonczak, S., Fournie-Zaluski, M. C., & Roques, B. P. (1999). A three-dimensional construction of the active site (region 507-749) of human neutral endopeptidase (EC.3.4.24.11). Protein engineering, 12(2), 141–149.

[3]Ma, Y. P., Gu, J. W., & Ma, J. (2019). A GLP-1 mutant and its preparation method and use [Chinese Patent Publication No. CN 110256553 A]. National Intellectual Property Administration of the People's Republic of China (CNIPA).

[4]Zhang, X., Belousoff, M. J., Zhao, P., Kooistra, A. J., Truong, T. T., Ang, S. Y., Underwood, C. R., Egebjerg, T., Šenel, P., Stewart, G. D., Liang, Y. L., Glukhova, A., Venugopal, H., Christopoulos, A., Furness, S. G. B., Miller, L. J., Reedtz-Runge, S., Langmead, C. J., Gloriam, D. E., Danev, R., … Wootten, D. (2020). Differential GLP-1R Binding and Activation by Peptide and Non-peptide Agonists. Molecular cell, 80(3), 485–500.e7. https://doi.org/10.1016/j.molcel.2020.09.020

[5]Oefner, C., Pierau, S., Schulz, H., & Dale, G. E. (2007). Structural studies of a bifunctional inhibitor of neprilysin and DPP-IV. Acta crystallographica. Section D, Biological crystallography, 63(Pt 9), 975–981. https://doi.org/10.1107/S0907444907036281

[6]Runge, S., Thøgersen, H., Madsen, K., Lau, J., & Rudolph, R. (2008). Crystal structure of the ligand-bound glucagon-like peptide-1 receptor extracellular domain. The Journal of biological chemistry, 283(17), 11340–11347. https://doi.org/10.1074/jbc.M708740200

[7]van Zundert, G. C. P., Rodrigues, J. P. G. L. M., Trellet, M., Schmitz, C., Kastritis, P. L., Karaca, E., Melquiond, A. S. J., van Dijk, M., de Vries, S. J., & Bonvin, A. M. J. J. (2016). The HADDOCK2.2 Web Server: User-Friendly Integrative Modeling of Biomolecular Complexes. Journal of molecular biology, 428(4), 720–725. https://doi.org/10.1016/j.jmb.2015.09.014

[8]Toshifumi Udo,Zijin Qin,Yang Jiao,Rakesh K. Singh & Fanbin Kong.(2025).A comprehensive review of mathematical modeling in probiotic microencapsulation.Food Engineering Reviews,17(2),1-21.

[9]Nathan R. Richbourg & Nicholas A. Peppas.(2020).The swollen polymer network hypothesis: Quantitative models of hydrogel swelling, stiffness, and solute transport. Progress in Polymer Science,105.

[10]Aleš Ručigaj & Tilen Kopač.(2025).Mathematical modeling of swelling and shrinking dynamics in pH-sensitive hydrogels composed of alginate, anionic cellulose, and chitosan.Polymer,328,128452-128452.

[11]Peppas NA.(1985).Analysis of Fickian and non-Fickian drug release from polymers. Pharm Acta Helv,60(4),110-1.

[12]Bacaër, N. (2011). Verhulst and the logistic equation (1838). Retrieved from https://webpages.ciencias.ulisboa.pt/~mcgomes/aulas/dinpop/Mod13/Verhulst.pdf

[13]Bryson, V., & Szybalski, W. (1952). Microbial selection. Science, 116(115).

[14]Monod, J. (1949). The growth of bacterial cultures. Annual Review of Microbiology, 3(1), 371–394. https://doi.org/10.1146/annurev.mi.03.100149.0021

[15]Lineweaver, H., & Burk, D. (1934). The determination of enzyme dissociation constants. Journal of the American Chemical Society, 56(3), 658–666. https://doi.org/10.1021/ja01318a036