Helicobacter pylori infection affects over half of the population. This seemingly tiny bacterium can cause chronic gastritis, gastric ulcers, and even increase the risk of gastric cancer. Once infected, it often remains dormant in the body for a long time, severely impacting quality of life.

Traditionally, we have relied primarily on combination therapy, primarily antibiotics, which initially achieved eradication rates of approximately 90%. However, with the rise of drug resistance, efficacy has steadily declined. Furthermore, antibiotics indiscriminately damage the intestinal flora, leading many patients to experience side effects such as diarrhea and bloating. Later, probiotics were tried as adjunctive therapy or even as monotherapy, but the results have been less than ideal, with eradication rates often below 30%. The core problem is that probiotics, as large particles, have difficulty penetrating gastric mucus and reaching the deeper layers where H. pylori accumulates. To this end, we designed a highly efficient delivery system. The stomach’s highly acidic environment and thick mucus layer have always been two major obstacles for drugs to reach the site of H. pylori infection. We encapsulated the modified yeast into microspheres using a pH-responsive gel—a gel that remains stable in gastric acid (pH < 4) and only dissolves and releases upon reaching the gastric mucosa (pH 6-7). Crucially, we loaded the microspheres with calcium carbonate-based micromotors, which provide propulsion and help the microspheres penetrate the viscous mucus layer and reach the “lesions” where H. pylori accumulates.

Thus, in the dry lab, we aimed to model the drug’s journey from entry into the body and gradual diffusion into the stomach wall, its penetration of the gastric mucus layer under the propulsion of the micromotors, its adhesion to the H. pylori, activation of the GPCR response pathway, and secretion of AiiA (killer protein), as well as the impact of drug intervention on the spread of H. pylori in the human population. In the following content, we will follow this logic to introduce our drylab content.

Dissolution Kinetics Modeling of pH-Sensitive Microspheres

First, we needed to ensure that our living drug would not lose its activity in the highly acidic environment of gastric fluid. Therefore, we devised a gel coating method. This prevented the Saccharomyces boulardii from dying in gastric fluid. Furthermore, to ensure that it could sense and respond properly deep within the gastric mucus layer, we also needed to ensure that the outer shell would readily degrade in a neutral environment. Therefore, we first needed to simulate the degradation of the entire gel under different pH conditions.

Mathematical Modeling

We model the dissolution of alginate-based drug-loaded microspheres using a first-order kinetic model. The main assumption is that the dissolution rate is proportional to the amount of undissolved material remaining at any given time. Mathematically, this can be expressed as:

where:

• $M(t)$: mass (or fraction) of undissolved microspheres at time tt
• $k$: dissolution rate constant, which depends on pH

Solving the differential equation by separation of variables gives:

Assuming $M_0 = 1$ ( 100% undissolved at $t=0$ ), it further simplifies to:

In real scenarios, the rate constant k varies with pH, so we further model:

This allows the prediction of drug release profiles at arbitrary pH levels.

Computational Method with Simulation

Workflow:

1. Collect experimental data of undissolved mass vs. time at different pH values.
2. Fit each data to the exponential decay model to extract the rate constant $k$.
3. Fit the $k$ pH relationship using a suitable nonlinear function (see Figure 1).
4. Use the fitted $k(pH)$ to simulate release curves under any target pH. Here we consider k as a function of ph.

Test data:

Tabulated experimental data (The values in the table represent the percentage of undepolymerization.) Here, due to reasons such as the amount of experiments, we were unable to conduct large-scale experiments. Therefore, in order to illustrate the reliability of our method, we used simulated data as a demonstration.

Result

We used the above data to calculate and obtained the following results.

Figure 1: Dissolution rate constant k as a function of pH.

After obtaining the functional relationship of $k$ with respect to pH, we can then predict the change in the percentage of undegraded gel microspheres over time at any given pH.

Figure 2: Model-predicted undissolved fraction curve at pH = 2.0.

Figure 3: Model-predicted undissolved fraction curve at pH = 6.0.

In real-world experiments, we found that in a strongly acidic environment, the microspheres remained undegraded for up to an hour. In a neutral environment, however, the gel microspheres completely degraded in about 90 seconds, which was close to our prediction.

Propulsion Model of Micropheres in Acidic Solution

Background

Next, we simulated the velocity of the microspheres after entering the human body, driven by the calcium carbonate reaction. We simulated the propulsive force of calcium carbonate-coated microspheres in gastric fluid and reviewed the classical Stokes solution for low Reynolds number flow around a sphere. We derived the equations of motion and simulated the velocity-time behavior to obtain the velocity of the microspheres in the liquid environment.

Parameters and Physical Assumptions

We consider a spherical $CaCO_3$ microsphere with:

• Mass: $m = 1\times 10^{-7}\,\mathrm{kg}$ (100 $\mu g$)
• Density: $\rho = 3.000\,\mathrm{g/cm^3}$

At the same time, by consulting the literature, we can also obtain the parameters of the gastric mucus layer：

Its radius is determined by:

Assume the CaCO₃ reacts with excess acid at pH 2, producing CO₂ by

Assuming the entire CaCO₃ sphere reacts over time $T_{\mathrm{react}}$, the maximum CO₂ generation rate is:

where $M_{\mathrm{CO}_2}$ and $M_{\mathrm{CaCO}_3}$ are the molar masses of CO₂ (44.01 g/mol) and CaCO₃ (100.09 g/mol), respectively.

The thrust is generated by the liberation of CO₂, imparting momentum to the particle (neglecting losses):

where $v_{\mathrm{ex}}$ denotes the effective exhaust velocity of CO₂ as it escapes (estimation: $v_{\mathrm{ex}}$ can be approximated as gas flow velocity at interface, typically 10–30 m/s per microfluidics literature).

Modeling the Drag Force in Fluid

The drag force is described using Stokes law.

Stream Function in Cylindrical Coordinates

For an incompressible, axisymmetric, low-Reynolds-number flow past a fixed sphere of radius $R$ in a uniform stream $u$, introduce the Stokes stream function $\psi(r,z)$ satisfying:

which ensures $\nabla\cdot\mathbf{u}=0$ by construction.

Vorticity and Governing Equations

The sole nonzero vorticity is:

and it satisfies:

Exact Solutions for Stream Function and Velocities

The well-known solution is:

leading to velocity components $u_r(r,z)$ and $u_z(r,z)$.

Vorticity, Pressure, and Drag Force

The azimuthal vorticity and pressure fields are:

In spherical coordinates:

The total drag force—evaluated via surface integration of the stress tensor—is:

where $\mu$ is the dynamic viscosity of the surrounding solution (for water at 25 °C, $\mu=5\times10^{-3}$ Pa·s⁴; acidic media similar). This classical Stokes drag formula is valid for $Re \leq 1$ ¹.

Equation of Motion

Applying Newton’s second law:

with $F_{\mathrm{thrust}}$ as above. If $\dot{m}_{\mathrm{CO}_2}$ is constant, the analytical solution for velocity is straightforward. If time-varying, a numerical solution is required. So we will get the numerical solution through numerical simulation.

Result

In summary, we wrote a program in Python to solve the ODE equations above. We then plotted the velocity-time curve.

Figure 4: Predicting the speed of a microsphere.

The results showed that the microspheres’ speed in the stomach was relatively stable, maintaining a constant speed of around 20 micrometers per second when only the micromotor (calcium carbonate) and acid reaction were considered. To further ensure the accuracy and reliability of our model, we consulted relevant literature. Existing articles have shown that similar methods have also achieved microsphere movement speeds of 20-150 $\mu m/s$ micrometers per second²,³, demonstrating the validity of our model.

At the same time, take into account the possible adverse effects of gas release，we can simply calculate the amount of gas released during the reaction of each microsphere.

Figure 5: CO₂ release dynamics predition.

One microsphere will only produce about 0.01 ml of gas, which means that taking 2g of the drug will only produce 80 ml of gas, which is safe enough for the human body.

Mathematical Modeling of Microsphere Diffusion in the Stomach

Background

To study the movement, diffusion, and contact of gel microspheres in the gastric cavity in drug delivery systems, we establish a dynamic model combining Brownian motion, gastric fluid flow, and boundary conditions. The goal is to predict the spatial distribution of microspheres, adhesion ratio, and excretion over a given timescale.

Model Assumptions

1. Geometry: The gastric cavity is defined by two eighth-order polynomial curves (lesser curvature and greater curvature), forming a closed region $\Omega$.
2. Movement: Microsphere movement is driven by Brownian diffusion and gastric fluid flow.
3. Boundary conditions: Microspheres near the gastric wall have a probability of adhesion; if they reach the outlet, they are considered excreted.

Mathematical Model

1. Gastric Cavity Boundary

Here $s$ denotes the lesser curvature and $l$ denotes the greater curvature. The gastric cavity region is:

2. Microsphere Dynamics

Each microsphere position $(x,y)$ satisfies the update equations:

• Diffusion term (Brownian motion):

• Flow term:

where $v$ is the average flow velocity and $\hat{\mathbf{v}}(x,y)$ is the locally normalized velocity direction.

3. State Transitions

Microspheres exist in three states:

• Free state ($S=0$): moving freely;
• Attached state ($S=1$):

• Exited state ($S=2$):

Simulation Results

We perform simulation according to the above formula. We also observed the effects of the drug 60 seconds and 10 minutes after it entered the body. Through numerical simulations (initial microsphere number $N=1000$, total simulation time 600s), the following conclusions are drawn:

Figure 6: The diffusion of microspheres in the stomach was observed in 60 seconds.

Figure 7: The diffusion of microspheres in the stomach.

1. Microspheres gradually deposit in the gastric fundus region under diffusion and flow effects, continually contacting the gastric wall. At 60 seconds, a large number of microspheres have diffused to the surface of the stomach wall.
2. Over time, the number of free microspheres decreases, the proportion of attached microspheres increases, and some microspheres are excreted through the outlet.
3. Long-term simulation results show that about 80% of gel microspheres eventually contact and adhere to the gastric wall, while the remaining microspheres mainly exist in free or excreted forms.²⁸

Physical Modeling of the Penetration Process of Microspheres through Gastric Mucus Layer

Background

The objective of this section is to establish a physical model to describe and analyze the penetration process of a microsphere drug delivery system through the gastric mucus layer. In our design, the drug is encapsulated in a pH-responsive hydrogel microsphere, with calcium carbonate-based micromotors loaded on its surface. Upon entering the stomach, the microspheres remain stable in the acidic environment and release the drug when approaching the gastric mucosa. A key challenge lies in effectively penetrating the thick gastric mucus layer to reach the H. pylori aggregation region.

We aim to answer the following questions through our physical model:

1. What are the conditions for microsphere penetration of the mucus layer?
2. How does the incident angle $\theta$ affect penetration efficiency? Is there a maximum angle $\theta_{max}$ such that penetration is possible when $\theta \leq \theta_{max}$? What is the corresponding penetration probability?
3. How do mucus thickness $H$, viscosity $\eta$, microsphere radius $R$, and propulsion speed $v$ affect the penetration conditions?

Figure 8: Illustration of microspheres through gastric mucus layer.

Model Assumptions

To simplify the complex real-world situation, we make the following assumptions:

• Microsphere assumption:
  The microsphere is a rigid sphere with radius $R$ and mass $m$. It moves at a constant velocity $v$ with a fixed incident angle $\theta$. The propulsion force always aligns with the direction of motion.

• Medium assumption:
  The mucus layer has thickness $H$, modeled as a uniform Newtonian fluid with constant viscosity $\eta$, independent of velocity. The mucus layer boundaries are clearly defined.

• Mechanical assumption:
  Upon entering the mucus layer, the microsphere is subject to:

  • Propulsion force $F_p$;
  • Viscous resistance $F_d$, described by Stokes’ law: $$\begin{align} F_d = 6\pi \eta R v \end{align}$$
  • The relationship between propulsion force and velocity has been simulated in previous work, showing that in the mucus layer region, velocity $v \approx 50\ \mu\text{m/s}$.

  Gravity and buoyancy are neglected.

• Energy criterion assumption:
  The criterion for penetration is that the microsphere’s kinetic energy in the vertical direction must be greater than or equal to the energy needed to overcome resistance and interfacial effects.

Mathematical Modeling Process

Time and Angle Relationship from the Perspective of Uniform Motion

Arrival Time and Incident Angle

A microsphere moves with velocity $v$ at constant speed into a mucus layer of thickness $H$, with incident angle $\theta$ (relative to the perpendicular direction).

Its travel path length is:

The time to reach the gastric wall is:

Rearranging:

where $t \geq \tfrac{H}{v}$ is required for physical meaning.

Reachable Angle Condition

If the microsphere has a maximum sustainable travel time $T$, then only when

can it reach the gastric wall. That is:

Therefore, the critical incident angle is:

Boundary cases:

• If $T < H/v$, no angles can reach;
• If $T = H/v$, only vertical incidence $\theta=0$ can reach;
• If $T > H/v$, all $\theta \in [0,\theta_{\max}]$ are reachable.

Penetration Probability (Uniform Angle Distribution Assumption)

Assuming the incident angle is uniformly distributed over $[0,\pi/2]$, the reachable probability is:

Numerical Calculation

Experiments and literature show that $H=150\ \mu\text{m}, T=90\ \text{s}$, with $v \approx 50\ \mu\text{m/s}$ in the mucus layer region:

penetration probability:

Results and Physical Interpretation

• Increasing viscosity $\eta$, radius $R$, or thickness $H$ increases required energy, reducing the critical angle $\theta_{\max}$, requiring incidence closer to vertical.
• Increasing velocity $v$ increases kinetic energy, increasing $\theta_c$, making penetration easier at larger angles.
• Under current assumptions, most microspheres can penetrate the gastric mucus layer and reach the gastric wall.

Overall, we establishes an energy-based physical model describing the microsphere penetration process through the gastric mucus layer, deriving an analytical expression for the critical incident angle $\theta_{\max}$ and estimating the penetration range. The model reveals how velocity, radius, viscosity, and thickness affect penetration ability, providing a theoretical basis for optimizing microsphere design and guiding subsequent experiments.

Molecular Dynamics Simulation Analysis of the HopQ-C1ND Protein Complex

Background

HopQ, a type I adhesin of H. pylori, mediates bacterial colonization and regulates host immune responses by specifically binding to the N-terminal domain (C1ND) of CEACAM1 on human host cells. CEACAM1 is a core molecular pair in H. pylori pathogenesis⁵. We achieved targeted adhesion to H. pylori by displaying C1ND on the surface of yeast.

Molecular dynamics (MD) simulations can track the motion of atoms under physiological conditions and dynamically analyze the structural stability of protein complexes. We used MD simulations to study the dynamic behavior of the HopQ-C1ND complex under physiological conditions, providing theoretical validation for our project’s adhesion system.

Indicators of the Molecular Dynamics Simulation

1. The protein complex is recognized to maintain stable overall conformation under simulated physiological conditions (300 K, 1 atm, aqueous environment), with backbone root-mean-square deviation (RMSD) stabilizing within a physiologically reasonable range of ≤0.35 nm in the later simulation stages.

2. The binding interface region of the complex is considered to exhibit higher structural flexibility compared to core domains, with significantly increased root-mean-square fluctuation (RMSF) values, and this flexibility directly correlates with binding specificity and dynamic regulatory functions.

3. Complex stability is supposed to be maintained through synergistic effects of multiple non-covalent interactions, with van der Waals forces and electrostatic forces as primary driving forces, and the total binding free energy showing a significantly negative value, indicating thermodynamically spontaneous and stable binding.

4. After energy minimization and NVT/NPT equilibration, physical parameters including temperature, pressure, and density are required to converge to physiologically steady-state ranges (temperature fluctuation ≤1 K, pressure approaching 1 atm, density ≈1000 kg/m³), to unsure biological relevance of simulation results.

5. The solvent-accessible surface area (SASA) and radius of gyration (Rg) of the complex are required to exhibit small fluctuation amplitudes during simulation, with hydrogen bond counts maintained at high levels, which collectively reflect conformational stability.

Methods

This study conducted all-atom molecular dynamics simulations using GROMACS 2018.8 software¹⁰, with the following workflow:

Target Protein Structure Acquisition and Screening

Target Protein Structure Acquisition and Screening High-resolution experimental structures were retrieved and screened from the Protein Data Bank (PDB)⁷,⁸ :

• HopQ monomer structure: PDB ID 5LP2 (X-ray diffraction, 2.60 Å resolution, no mutations, expressed in E. coli);
• C1ND (CEACAM1 N-domain) structure: PDB ID 4WHD (X-ray diffraction, 2.50 Å resolution, homodimer, human origin);
• HopQ-CEACAM1 complex reference structure: PDB ID 6AW2 (X-ray diffraction, 2.68 Å resolution, heterodimer, containing complete binding interface).

All structures were validated using the wwPDB Validation tool, with Ramachandran plot favored regions ≥95% and no significant structural defects.

Protein Structure Preprocessing

1. Structure splitting and integration: Binding interface coordinates of HopQ and CEACAM1 were extracted from 6AW2, with 5LP2 and 4WHD as templates to generate initial HopQ-C1ND complex structure using PyMOL;

2. Residue repair: Missing residues 120-125 in 5LP2 and broken peptide bonds in 4WHD were completed using the SWISS-MODEL server⁶,⁹, with modeling templates selected from homologous structures with ≥85% sequence identity (GMQE=0.85, QMEAN=-1.2);

3. Structure optimization: Atomic overlaps in the repaired structure were eliminated through energy minimization to ensure bond lengths and angles conform to standard force field parameters.

Figure 9: The structure of HopQ-CEACAM1 complex.

Topology File Construction and Force Field Selection

1. Force field selection: The OPLS-AA/L all-atom force field was chosen for its high accuracy in describing protein-ligand interactions, widely applied in protein complex simulations;

2. Topology file generation: GROMACS’ gmx pdb2gmx module was used to assign force field parameters to the complex, generating topology (top) and coordinate (gro) files;

3. Ion and solvent models: The SPC water model was used to describe water molecules, with Cl⁻ ions added to neutralize total system charge (0.15 mol/L ion concentration, simulating physiological salt environment).

   gmx pdb2gmx -f 6AW2.pdb -o 6AW2_processed.gro -water tip3p `

Simulation Environment Construction

1. Boundary conditions: Cubic periodic boundary conditions were employed, with minimum distance between box edges and protein surface set to 10 Å to ensure complete solvent encapsulation;

   gmx editconf -f 6AW@_processed.gro -o 6AW2_newbox.gro -c -d 1.2 -bt cubic

2. Solvent filling: Approximately 25,000 SPC water molecules were added using the gmx solvate module;

   gmx solvate -cp 6AW2_newbox.gro -cs spc216.gro -o 6AW2_solv.gro -p topol.top

   

   Figure 10: HopQ-CEACAM1 complex after solvent filling.

3. Ion addition: Water molecules in the solvent box were replaced with Cl⁻ ions using gmx genion to neutralize system charge (total charge +4, 4 Cl⁻ ions added).

   gmx grompp -f inputs/ions.mdp -c 6AW2_solv.gro -p topol.top -o ions.tpr
   
   gmx genion -s ions.tpr -o 6AW2_solv_ions.gro -p topol.top -pname NA -nname CL -neutral

Figure 11: HopQ-CEACAM1 complex after ion addition.

Pre-Simulation Equilibration

1. Energy minimization: Combined “steepest descent (1000 steps) + conjugate gradient (5000 steps)” methods were used for energy minimization, with termination condition of maximum force ≤1000 kJ/(mol·nm), resulting in final system potential energy below -4.4×10⁶ kJ/mol;

   gmx grompp -f inputs/minim.mdp -c 6AW2_solv_ions.gro -p topol.top -o em.tpr
   
   gmx mdrun -v -deffnm em
   
   gmx energy -f em.edr -o potential.xvg

Figure 12: Energy change during energy minimization.

2. NVT equilibration: 200 ps simulation under constant volume (V) and temperature (T=300 K) conditions, using V-rescale thermostat for temperature control, 2 fs integration step, with protein backbone atoms constrained (LINCS algorithm) to stabilize system temperature at 300±1 K;

   gmx grompp -f inputs/nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
   
   gmx mdrun -deffnm nvt
   
   gmx energy -f nvt.edr -o temperature.xvg

   

   Figure 13: Temperature change during NVT equilibration.

3. NPT equilibration: 200 ps simulation under constant pressure (P=1 atm) and temperature (T=300 K) conditions, using Parrinello-Rahman barostat with compressibility 4.5×10⁻⁵ bar⁻¹ to stabilize system density at 1000±5 kg/m³.

   gmx grompp -f inputs/npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
   
   gmx mdrun -deffnm npt
   
   gmx energy -f npt.edr -o pressure.xvg

   

   Figure 14: Density change during NPT equilibration.

   

   Figure 15: Pressure change during NPT equilibration.

Production Simulation and Result Analysis

1. Production run: 100 ns NPT ensemble simulation with parameters consistent with equilibration, trajectory saved every 10 ps;

   gmx grompp -f inputs/md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_10.tpr
   
   gmx mdrun -deffnm md_0_10

2. Trajectory analysis: Core indicators calculated using GROMACS modules:

   gmx rms: Backbone Cα atom RMSD calculation;

   gmx rmsf: Per-atom RMSF calculation;

   gmx gyrate: Total Rg and axial Rg calculation;

   gmx sasa: Total SASA and binding interface SASA calculation;

   gmx hbond: Hydrogen bond counting (criteria: donor-acceptor distance ≤0.35 nm, angle ≥150°);

   gmx mmpbsa: Binding free energy calculation using MM-GBSA method;

3. Visualization and plotting: Protein structures visualized with PyMOL, and indicator curves plotted with QtGrace.

Parameter Table

Results

Validation of Simulation System Equilibrium

1. Energy convergence characteristics: During energy minimization, total system potential energy decreased from initial -2.0×10⁶ kJ/mol to -4.4×10⁶ kJ/mol, with bond length, bond angle, and dihedral energies converging synchronously (fluctuation amplitude ≤5%), indicating complete elimination of unreasonable interatomic interactions. Potential energy remained stable during NVT and NPT equilibration with no significant jumps, confirming establishment of energy steady state.

2. Temperature and pressure dynamics: Temperature rapidly increased to 300 K during initial NVT equilibration (0-50 ps), then stabilized at 300±0.5 K with Gaussian distribution (mean 300 K, standard deviation 0.3 K); pressure gradually converged to 1±50 atm from initial fluctuations (-400~200 atm) during NPT equilibration, with amplitude ≤20 atm by 200 ps; density stabilized at 1000±2 kg/m³ after 100 ps NPT equilibration, consistent with liquid water density. All parameters met physiological steady-state standards, validating the rationality of the simulation environment.

Overall Structural Stability of Protein Complex

1. RMSD time evolution: Backbone Cα atom RMSD showed “three-stage” behavior:

   Initial stage (0-20 ns): Rapid increase (0.10→0.28 nm), corresponding to initial conformational adjustment in aqueous environment;

   Transition stage (20-80 ns): Slow increase with minor fluctuations (0.28→0.30 nm), representing local structural optimization;

   Stable stage (80-100 ns): Stabilized at 0.29±0.01 nm with fluctuation amplitude ≤0.02 nm and no significant drift, indicating achievement of stable overall complex structure.

   

   Figure 16: RMSD analysis during simulation.

2. Radius of gyration (Rg) characteristics: Total Rg remained stable at 3.02±0.03 nm throughout 100 ns simulation with only 0.99% coefficient of variation, demonstrating good overall compactness without significant expansion or contraction. Axial Rg analysis showed:

   X-axis Rg: 1.5-2.8 nm, mean 2.1±0.17 nm, coefficient of variation 8.2%;

   Y-axis Rg: 2.0-2.9 nm, mean 2.4±0.16 nm, coefficient of variation 6.5%;

   Z-axis Rg: 1.5-2.8 nm, mean 2.2±0.17 nm, coefficient of variation 7.9%.

   Maximum X-axis Rg fluctuation, combined with secondary structure analysis, is hypothesized to be related to flexible movement of HopQ’s α-helix terminal segment (5LP2: residues 250-260).

   

   Figure 17: Radius of gyration during simulation.

3. Solvent-accessible surface area (SASA): Total SASA averaged 195±15 nm² with range 180-210 nm² and 7.7% coefficient of variation. Binding interface SASA decreased to 45±5 nm² in later simulation stages, approximately 18% lower than initial stage (55 nm²), indicating increased burial of interface residues within the molecule after binding, reducing solvent exposure and further enhancing complex stability.

   

   Figure 18: Solvent-accessible surface area during simulation.

Local Structural Flexibility Characteristics (RMSF Analysis)

RMSF profiles revealed significant flexibility differences across complex regions, categorized into three types:

1. Low flexibility regions (RMSF≤0.2 nm): Mainly distributed in complex core domains, including 4WHD’s β-sheet region (residues 10-50) and 5LP2’s β-barrel region (residues 80-120). Restricted atomic movement in these regions forms the structural scaffold, providing stable support for binding.

2. Moderate flexibility regions (0.2-0.4 nm): Located at junctions between α-helices and β-sheets (e.g., 5LP2 residues 150-170). Moderate flexibility in these regions enables conformational fine-tuning of core domains, avoiding binding energy loss from rigid structures.

3. High flexibility regions (RMSF≥0.4 nm): Precisely localized at the HopQ-C1ND binding interface, including 5LP2 residues 210-230 (RMSF=0.55±0.03 nm) and 4WHD residues 80-100 (RMSF=0.52±0.04 nm). Sequence analysis showed these regions are enriched in polar amino acids (Asp, Lys, Asn), whose dynamic movement modulates hydrogen bond and electrostatic interaction formation, enhancing binding specificity.

Figure 19: RMSF analysis during simulation.

Intermolecular Interaction Characteristics

1. Hydrogen bond network: Hydrogen bond count remained stable at 2000-2500 pairs during simulation, averaging 2250±150 pairs. Trajectory tracking identified 12 high-frequency stable hydrogen bonds (duration ≥80 ns), with core pairs including:

   • 5LP2: Arg185 - 4WHD: Glu92 (92 ns duration);
   • 5LP2: Asn201 - 4WHD: Lys88 (88 ns duration);
   • 5LP2: Tyr215 - 4WHD: Asp76 (85 ns duration).

   These hydrogen bonds are located at the binding interface center, representing key interactions maintaining complex stability.

   

   Figure 20: Hydrogen bond analysis during simulation.

2. Binding free energy analysis: MM-GBSA calculations showed total binding free energy (TOTAL) of -34864.90±500 Kcal/mol, with significantly negative value indicating thermodynamically spontaneous and stable binding. Energy component contributions were:

   • Van der Waals energy (VDWAALS): -11490.15 Kcal/mol, accounting for 32.9% of total binding energy as primary driving force;

   • Electrostatic energy (EEL): -6965.76 Kcal/mol, 19.9% contribution, synergizing with van der Waals forces to enhance binding;

   • Generalized Born solvation energy (EGB): -4679.48 Kcal/mol, with polar solvation effects favoring complex stability;

   • Surface energy (ESURF): 155.09 Kcal/mol, with positive nonpolar solvation contribution offset by EGB, resulting in total solvation free energy (GSOLV) of -4524.39 Kcal/mol;

   • Gas-phase free energy (GGAS): -3574.99 Kcal/mol, reflecting stability of intramolecular interactions.

   

   Figure 21: Free energy analysis during simulation.

Data Sources

1. Protein structure data: Protein Data Bank (PDB) database, including:

   • HopQ: PDB ID 5LP2 Opens in a new tab , released 2016;
   • C1ND: PDB ID 4WHD Opens in a new tab , released 2015;
   • HopQ-CEACAM1 complex: PDB ID 6AW2 Opens in a new tab , released 2018.

2. Structure repair tool: SWISS-MODEL server Opens in a new tab , accessed August 2025, modeling quality assessment: GMQE=0.85, QMEAN=-1.2.

3. Simulation and analysis software: GROMACS 2018.8 Opens in a new tab , force field parameters from OPLS-AA/L, original literature 2001; QtGrace 0.2.6 (for data visualization); PyMOL 2.5 (for structure display).

4. Force field and model parameters: SPC water model and OPLS-AA/L force field parameters from official GROMACS force field library (D:/gromacs/gmx2018.8/share/gromacs/top).

Modeling and Analysis of GPCR Signaling Pathway

Background

In the microscopic biological world, G protein-coupled receptors (GPCRs) function as cellular “signal antennas” that receive extracellular signals and initiate intracellular cascade reactions. This process plays a core role in numerous life activities, including yeast mating and mammalian physiological regulation. In our project, GPCRs serve as the central component of the sensing system. However, the activation of GPCRs and the transmission of downstream pathways involve highly dynamic and complex conformational changes, which cannot be fully captured by static experimental structures alone.

To address this challenge, we screened for the most sensitive and effective engineered GPCRs using molecular docking and molecular dynamics (MD) simulations. Furthermore, we deciphered the complete signal transduction mechanism from perception to response via ordinary differential equation (ODE) mathematical modeling, enabling quantitative analysis of the signaling process. Ultimately, we innovatively reengineered the yeast mating signal system into a H. pylori recognition system. Our GPCR signaling pathway modeling, integrated with wet-lab validation, provides a promising universal strategy for engineering yeast GPCR sensing systems to achieve specific functions, offering valuable insights for future iGEM teams focusing on GPCR modification.

Molecular Dynamics (MD) Simulation Analysis of GPCRs

Sources of GPCR Variants

To engineer the endogenous yeast GPCR (Ste2) for the recognition of α-methylhistamine, we adopted three parallel design strategies to screen for the most sensitive and efficient GPCR sensors:

Strategy 1: Chimeric GPCR

Principle: This strategy leverages the modular nature of GPCRs. We fused the ligand-recognition domains—extracellular domains (ECDs) and transmembrane domains (TMDs)—of human histamine receptors (e.g., hH1R, hH2R, hH3R, hH4R) with the intracellular domains (ICDs) of the yeast Ste2 receptor, which mediate interactions with the downstream G protein Gpa1, thereby constructing novel chimeric proteins.

Strategy 2: Exogenous GPCR + Chimeric G-protein

Principle: G proteins also exhibit prominent modularity, with their GPCR selectivity primarily determined by several key amino acids at the C-terminus. We retained intact, unmodified human GPCRs (e.g., hH3R) while engineering the endogenous yeast Gα protein (Gpa1). Specifically, we replaced the 5 C-terminal amino acids of Gpa1 with the corresponding sequences from three human Gα subtypes (Gαi, Gαs, Gαq), generating three chimeric G proteins: Gpai, Gpas, and Gpaq.²⁶

Strategy 3: Literature-Validated Positive Controls

To verify the reliability of our simulation system, we incorporated two systems previously validated in the literature²⁷ :

mhH3R-Gpai: An optimized modified histamine H3 receptor with high ligand sensitivity, paired with the chimeric Gpai protein.

P2Y2-Gpap: The human P2Y2 receptor, an extracellular ATP (eATP) sensor. Literature has confirmed its synergistic function with the chimeric G protein Gpap in yeast, where it successfully activates downstream pathways upon ATP binding. This system served as the “positive control” for our simulation workflow to ensure our analytical methods could accurately capture known activation events.

Methods

Preparation

1. GPCRs and G proteins: The 3D structures of target GPCRs (hH1R, hH2R, hH3R, hH4R, mhH3R, P2Y2) and G proteins (Gpai, Gpaq, Gpap, Gpas) were prioritized from the Protein Data Bank (PDB)⁷,⁸. For proteins lacking experimental structures or containing missing loops (e.g., chimeric proteins), homology modeling tools SWISS-MODEL were used for structure completion and construction. All structures underwent rigorous inspection prior to use.

2. Ligand Preparation: The mol2 structure file of the target ligand Nα-methylhistamine (NAMH) was retrieved from the ZINC database or CHEMICAL BOOK website. The Antechamber tool and General Amber Force Field (GAFF) were used to generate ligand topology files¹⁴, which were then converted to GROMACS-compatible itp force field files using the ACPYPE script¹⁵,¹⁶.

3. Construction of Protein-Membrane Complexes: The “Membrane Builder” module of the CHARMM-GUI online server was employed for system construction. The GPCR structure was embedded in a pre-equilibrated phospholipid bilayer mimicking the yeast membrane composition (POPE: POPI: POPC: POPS: ERG = 6: 2: 2: 1: 6; see the figure below for details)¹⁷ . Subsequently, the corresponding Gα protein (e.g., Gpai) was positioned on the intracellular side of the GPCR, with reference to resolved GPCR-G protein complex structures to ensure a physiologically relevant initial orientation.

   

   Figure 22: The composition table of membrane.

Then we can get the results：

Figure 23: The generation report of the membrane.

Figure 24: Sketch map of membrane-protein complex.

Molecular Docking Process

Semi-flexible docking was performed using AutoDock Vina to predict the optimal binding conformation and affinity of NAMH within the active pockets of different GPCR variants.

Workflow:

1. The GPCR was treated as a rigid receptor, and NAMH as a flexible ligand.

2. A sufficiently sized grid box was defined in the region enclosed by transmembrane helices, based on either the binding sites of known homologous ligands or cavity recognition algorithms.

3. Docking calculations were executed to generate multiple binding conformations.

4. The conformation with the lowest binding affinity (kcal/mol) and the most clustered, physiologically reasonable structure was selected as the initial structure for subsequent MD simulations.

   

   Figure 25: Molecular docking site.

Molecular Dynamics Simulation Process

Software and Force Fields: All MD simulations were performed using the GROMACS 2020 package¹⁰. The CHARMM36m force field was applied to proteins, lipids, and ions, and the TIP3P water model was used for solvation.

   grep -v "HSM" input.pdb | grep "ATOM" > protein.pdb

   grep "HSM" input.pdb > ligand.pdb
   
   gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top -ff    amber99sb-ildn -water tip3p -ignh

   gmx editconf -f ligand.pdb -o ligand.gro

   sed -i '/forcefield.itp/a #include "ligand.itp"' topol.top

   echo "" >> topol.top 

   echo "; Ligand" >> topol.top

   echo "HSM         1" >> topol.top

   PROTEIN_ATOMS=$(tail -n 2 protein.gro | head -n 1 | awk '{print $1}')

   LIGAND_ATOMS=$(tail -n 2 ligand.gro | head -n 1 | awk '{print $1}')

   BOX_VECTOR=$(tail -n 1 protein.gro)

   TOTAL_ATOMS=$((PROTEIN_ATOMS + LIGAND_ATOMS

   echo "Complex of Protein and Ligand" > complex.gro

   echo " ${TOTAL_ATOMS}" >> complex.gro

   sed -n '3,$(($PROTEIN_ATOMS+2))p' protein.gro >> complex.gro

   sed -n '3,$(($LIGAND_ATOMS+2))p' ligand.gro >> complex.gro

   echo "${BOX_VECTOR}" >> complex.gro

1. System Preparation:

• Solvation: The constructed GPCR-membrane-G protein complex was placed in a sufficiently large simulation box and solvated with TIP3P water molecules.

     gmx editconf -f complex.gro -o complex_boxed.gro -c -d 1.0 -bt cubic
  
     gmx solvate -cp complex_boxed.gro -cs spc216.gro -o system_solvated.gro -p topol.top

• Neutralization: Na⁺ or Cl⁻ ions were added to neutralize the total charge of the system and achieve a physiological salt concentration (~0.15 M).

  gmx grompp -f ions.mdp -c system_solvated.gro -p topol.top -o ions.tpr 
  
  echo "SOL" | gmx genion -s ions.tpr -o system_ions.gro -p topol.top -neutral

• Energy Minimization: The steepest descent algorithm was used to minimize the system energy, eliminating unfavorable steric clashes and atomic interactions.

  gmx grompp -f em.mdp -c system_ions.gro -p topol.top -o em.tpr 
  
  gmx mdrun -deffnm em

1. Equilibration:

• NVT Equilibration: A 1-nanosecond (ns) simulation was conducted under constant volume and temperature (310 K). Position restraints were applied to the heavy atoms of proteins and lipids to allow water molecules and ions to equilibrate around the solute.

  gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
  
  gmx mdrun -deffnm nvt

• NPT Equilibration: A 10-ns simulation was performed under constant pressure (1 bar) and temperature (310 K). Position restraints were gradually weakened and ultimately removed to stabilize the system density and pressure.

  gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -p topol.top -o npt.tpr 
  
  gmx mdrun -deffnm npt

3. Simulation:

• A 40,000-ps (40 ns) production simulation was performed for each docked system under the NPT ensemble. The Parrinello-Rahman pressure coupling and Nosé-Hoover temperature coupling methods were utilized. Long-range electrostatic interactions were calculated via the Particle Mesh Ewald (PME) algorithm, and a cutoff radius was set for van der Waals interactions and short-range Coulomb forces. Trajectory coordinates were saved every 10 picoseconds (ps).

  gmx grompp -f md.mdp -c npt.gro -p topol.top -o md.tpr 
  
  gmx mdrun -deffnm md

4. Result Analysis and Output:

• GROMACS built-in tools (e.g., gmx rms, gmx rmsf, gmx covar, gmx anaeig, gmx hbond) were used to analyze trajectory files, extract key parameters, and export data in xvg format for subsequent plotting and analysis.

  # --- 1.preparation---
  
  echo -e "resname ${LIGAND_ITP_MOLNAME}\nname 19 MyLigand\n\"Protein\" | \"MyLigand\"\nname 20 Complex\nq" | gmx make_ndx -f em.gro -o index.ndx
  
  check_file_exists "index.ndx"
  
  echo "Complex System" | gmx trjconv -s md.tpr -f md.xtc -o md_whole.xtc -n index.ndx -pbc whole -center -ur compact
  
  check_file_exists "md_whole.xtc"
  
  echo "Backbone System" | gmx trjconv -s md.tpr -f md_whole.xtc -o md_final.xtc -fit rot+trans
  
  check_file_exists "md_final.xtc"
  
  # --- 2. RMSD analysis ---
  
  echo "Backbone Backbone" | gmx rms -s em.tpr -f md_final.xtc -o rmsd.xvg -tu ns
  
  # --- 3. RMSF analysis ---
  
  echo "C-alpha" | gmx rmsf -s md.tpr -f md_final.xtc -o rmsf.xvg -res
  
  # --- 4. Rg analysis---
  
  echo "Complex" | gmx gyrate -s md.tpr -f md_final.xtc -o gyrate.xvg
  
  # --- 5. TM2-TM6 Distance analysis ---
  
  gmx make_ndx -f em.gro -o index_tm2-tm6.ndx <<EOF
  
  a CA & r ${TM2_RES}
  
  name 21 TM2_CA
  
  a CA & r ${TM6_RES}
  
  name 22 TM6_CA
  
  q
  
  EOF
  
  check_file_exists "index_tm2-tm6.ndx"
  
  echo -e "TM2_CA\nTM6_CA" | gmx distance -s md.tpr -f md_final.xtc -n index_tm2-tm6.ndx -oav distance_TM2-TM6.xvg

Results

Molecular Docking

AutoDock enabled the determination of the spatial binding modes of α-methylhistamine with GPCRs, and the binding affinities were quantified using the built-in scoring function of AutoDock Vina, as summarized in the table below.

Figure 26: Binding energy without G protein (P2Y2 to ATP, others to NAMH).

Figure 27:Binding energy with G Protein (P2Y2 to ATP, others to NAMH).

Three hH3R combinations (hH3R-Gpai, hH3R-Gpas, hH3R-Gpaq), along with mhH3R-Gpai and P2Y2-Gpap, were selected as candidates based on their favorable binding affinities. However, binding affinity alone is insufficient; the key criterion is whether ligand binding induces conformational changes to initiate downstream signal transduction. Thus, we further performed MD simulations for these five candidate combinations.

Molecular Dynamics Simulation Results

Quality Control:

1. Root Mean Square Deviation (RMSD): RMSD of the protein backbone was calculated to assess structural stability during simulations. All five systems reached equilibrium within ~10–20 ns, with RMSD values plateauing, indicating stable system structures without extreme events such as unfolding.

   

   Figure 28: RMSD analysis during simulation.

2. Radius of Gyration (Rg): Rg reflects the compactness of the protein structure. In all simulations, the total Rg values (black line in figures) remained stable without significant fluctuations, confirming that proteins maintained their native folded conformation and compactness throughout the 40-ns simulation.

   

   Figure 29: Radius of gyration during simulation.

3. Root Mean Square Fluctuation (RMSF): RMSF analysis revealed residue-specific flexibility. All systems exhibited typical GPCR dynamic features: low fluctuations (structural rigidity) in transmembrane helices, and high fluctuations (structural flexibility) in intracellular/extracellular loops and N/C-termini.

   

   Figure 30: RMSF analysis during simulation.

4. Principal Component Analysis (PCA):

   Eigenvalues of the Covariance Matrix: For all systems, the first ~10 eigenvalues accounted for the majority of total system motion, followed by a rapid decay to near zero. This indicates that protein motion is dominated by a few major collective modes (principal components), justifying the use of the first two components (PC-1, PC-2) to describe conformational changes.

   

   Figure 31: Eigenvalues of the covariance matrix.

   Projection onto Eigenvectors: Trajectory projection onto PC-1 and PC-2 showed that systems did not fully converge to a single conformational cluster within 40 ns, but rather exhibited directional drift along PC-1. This drift suggests a slow, global conformational transition occurring on a timescale longer than the simulation. Specifically,CH3R, hH3R-Gpai, hH3R-Gpaq, hH3R-Gpas, mhH3R-Gpai, and the positive control P2Y2-Gpap all showed a continuous increasing trend along PC-1, indicating shared major motion modes. Behavior along PC-2 was more complex, with oscillations or gradual changes (e.g., hH3R-Gpai showed decreasing PC-2 values with oscillation, while P2Y2-Gpap showed an initial increase followed by a decrease), reflecting differences in secondary motion modes. PCA results confirmed high dynamics of these GPCR-G protein complexes, and the shared drift along PC-1 (consistent with the positive control) suggests conformational changes associated with receptor activation.

   

   Figure 32: Projection onto eigenvectors.

Activation State Analysis

GPCR activation is characterized by conformational rearrangement of transmembrane (TM) helices. According to a Nature study²⁵, the class D fungal GPCR Ste2 has a unique activation mechanism summarized as “obstruction removal and inward compaction”:

Inactive State: In the absence of agonist, the intracellular end of TM helix H7 is disordered and folded inward, physically blocking the G protein binding site.

Active State: Agonist binding triggers conformational changes, including a ~20 Å outward swing of the disordered H7 region (removing the obstruction) and a ~12 Å inward movement of the intracellular end of H6 (forming the G protein binding site).

We used the “12 Å inward movement of the H6 intracellular end” as the activation criterion, monitoring the H2-H6 distance during simulations:

Figure 33: The distance between H2 and H6.

P2Y2-Gpap System (Figure 33.g): As the positive control, this system displays a clear and rapid activation event. The H2-H6 distance sharply decreases from an initial ~3.5 nm (35 Å) to a stable conformation at ~1.5-2.0 nm (15-20 Å). This represents a significant reduction of approximately 1.5-2.0 nm (15-20 Å), confirming the validity of the simulation criteria for activation.

hH3R-Gpaq System (Figure 33.e): This system demonstrates a strong activation profile similar to the positive control. The distance undergoes a substantial reduction from ~8.0 nm (80 Å) to a stable range of ~5.0 nm (50 Å). The total decrease of ~3.0 nm (30 Å) clearly indicates a significant conformational change consistent with receptor activation.

CH3R System (Figure 33.b): This simulation also shows a notable distance reduction. The H2-H6 distance decreases from approximately 9.5 nm (95 Å) to a more compact state around 7.5 nm (75 Å), resulting in a total reduction of ~2.0 nm (20 Å), which suggests an activation-like conformational change.

hH3R-Gpai System (Figure 33.c): This system also meets the quantitative requirements for activation. The H2-H6 distance decreases from an initial value of approximately 10.5 nm (105 Å) to a range around 8.0 nm (80 Å). This results in a total reduction of ~2.5 nm (25 Å), which exceeds the threshold for an activation event. However, it’s worth noting that the final state exhibits more fluctuation compared to the stable conformations seen in the P2Y2-Gpap and hH3R-Gpaq systems.

CH1R (Figure 33.a), hH3R-Gpas (Figure 33.d), and mhH3R-pai (Figure 33.f) Systems: These three systems do not show clear evidence of sustained activation.

**Figure 33.a:** The distance fluctuates dynamically between 4.0 nm and 7.5 nm without settling into a stable, reduced state

**Figure 33.c:** Although there is a slight downward trend from ~10.5 nm to ~8.0 nm, the high degree of fluctuation and modest overall change suggest an unstable or non-activated state.

**Figure 33.f:** The distance remains in a relatively stable but dynamic range between 2.0 nm and 3.5 nm, showing no significant reduction.

Conclusion

The analysis reveals that the hH3R-Gpaq system (e) exhibits a clear and stable activation event, a behavior that is strongly mirrored by the positive control, P2Y2-Gpaq (g). Additionally, the hH3R-Gpai (c) and CH3R (b) systems show significant H2-H6 distance reductions (~2.5 nm and ~2.0 nm respectively) that meet the quantitative criteria for activation, although their final states are more dynamic. Conversely, the other hH3R-Gpas simulation (d) shows a distinct increase in distance, indicating a failure to activate. The remaining systems did not display conformational changes consistent with activation.

In summary, we established an efficient screening and analysis workflow using a dry-lab strategy combining molecular docking and MD simulations: First, molecular docking and binding affinity scoring enabled initial screening of numerous GPCR-G protein candidates, identifying high-potential combinations. Next, MD simulations probed the dynamic conformational evolution of engineered GPCRs upon NAMH binding. Finally, comprehensive evaluation of binding stability and conformational activity identified CH3R，hH3R-Gpai ,hH3R-Gpaq as the optimal candidate. This dry-lab strategy not only narrowed the scope of subsequent wet-lab validation but also revealed the molecular mechanisms of receptor-protein-ligand interactions, laying a solid theoretical foundation for functional validation of the sensing system.

ODE Modeling of G Protein Pathway Signal Transduction

When the H. pylori-associated signaling molecule hH3R diffuses to the gut and binds to the seven-transmembrane receptor Ste2, it activates downstream G proteins, recruits the scaffold protein Ste5 to form a signal complex, and further initiates the MAPK cascade reaction composed of Ste11-Ste7-Fus3. Ultimately, it activates the transcription factor Ste12 to regulate the expression of target genes¹². To clearly reveal the dynamic change rules of each molecular concentration, the rate-limiting steps of key reactions, and the synergistic mechanism of signal transmission in this process, this section adopts ordinary differential equations (ODEs) to construct a mathematical model, realizing the quantitative simulation and analysis of the GPCR signaling pathway, and laying a foundation for studying the sensitivity and response time of the sensing system in our project.

Figure 34: The GPCR signaling pathway in modified yeast.

Assumptions

The model is based on the following core assumptions that conform to biological common sense and physical laws:

1. Short-Term Steady-State Assumption: The simulation focuses on the short-term signal transmission process (within 300 seconds) after hH3R activation. De novo protein synthesis is not considered for the time being, and it is assumed that the total concentrations of core molecules such as Ste2, Ste5, and the MAPK family remain stable.

2. Unified Dephosphorylation Assumption: In the MAPK cascade reaction, the dephosphorylation of molecules such as Ste11, Ste7, and Fus3 is mediated by phosphatases. To simplify the model, it is assumed that they have the same basic dephosphorylation rate constant.

3. Molecular Behavior Assumption: The reactions of all molecules in the system follow the random collision theory and are not disturbed by external environmental factors such as temperature and pH. The reaction process is only determined by concentration and reaction rate constants.

4. Gβγ Concentration Correlation Assumption: The concentration of activated G protein βγ subunits (Gβγ) is directly linked to the concentration of the hH3R-Ste2 receptor complex, that is, the degree of receptor activation directly reflects the release amount of Gβγ.

5. Hill Coefficient Assumption: The binding of the transcription factor Ste12 to DNA does not consider cooperative effects for the time being, and the Hill coefficient is set to n=1, meaning that each binding site is independent of each other.

Core Formulas

The model is divided into five functional modules according to the signal transmission sequence. The core reaction equations, ODE equations, and parameter meanings of each module are as follows:

1. Receptor Activation Module

The binding and dissociation kinetic process of hH3R and the seven-transmembrane receptor Ste2, which is the initial step of signal transmission.

1. Reaction Equations

   • Binding reaction: $hH3R + Ste2 \stackrel{K_{hs}}{\to} hH3R-Ste2$

   • Dissociation reaction: $hH3R-Ste2 \stackrel{K_{sh}}{\to} hH3R + Ste2$

2. ODE Equations

   $$\begin{align} \frac{d[Ste2]}{dt} &= -K_{hs}[hH3R][Ste2] + K_{sh}[hH3R-Ste2] \\ \frac{d[hH3R-Ste2]}{dt} &= K_{hs}[hH3R][Ste2] - K_{sh}[hH3R-Ste2] \\ \frac{d[hH3R]}{dt} &= -K_{hs}[hH3R][Ste2] + K_{sh}[hH3R-Ste2] \end{align}$$

3. Parameter and Variable Meanings

   Symbol	Meaning
   $K_{hs}$	Reaction rate constant: Rate constant for hH3R and Ste2 to form the hH3R-Ste2 complex
   $K_{sh}$	Reaction rate constant: Rate constant for the dissociation of the hH3R-Ste2 complex into free hH3R and Ste2
   $[x]$	Concentration of the molecule in brackets (e.g., $[Ste2]$ represents the concentration of free Ste2 receptors)
   $hH3R-Ste2$	Receptor-ligand complex formed by the binding of hH3R and Ste2

2. Scaffold Formation Module

Activated Gβγ subunits recruit the scaffold protein Ste5 to form a multi-level assembly process of monomolecular binding, dimers, and advanced complexes, providing anchor sites for the MAPK cascade reaction.

1. Reaction Equations

   • Ste5 dimerization: $Ste5 + Ste5 \stackrel{K_{SS}}{\to} Ste5-Ste5$

   • Binding of Gβγ and Ste5: $G\beta\gamma + Ste5 \stackrel{K_{SG}}{\to} G\beta\gamma-Ste5$

   • Binding of Gβγ and double Ste5: $G\beta\gamma + Ste5 + Ste5 \stackrel{K_{GSS1}}{\to} G\beta\gamma-Ste5-Ste5$

   • Binding of Gβγ-Ste5 and Ste5: $G\beta\gamma-Ste5 + Ste5 \stackrel{K_{GSS2}}{\to} G\beta\gamma-Ste5-Ste5$

   • Binding of double Gβγ-Ste5: $G\beta\gamma-Ste5 + G\beta\gamma-Ste5 \stackrel{K_{GSSG1}}{\to} G\beta\gamma-Ste5-Ste5-G\beta\gamma$

   • Binding of Gβγ and Gβγ-Ste5-Ste5: $G\beta\gamma-Ste5-Ste5 + G\beta\gamma \stackrel{K_{GSSG2}}{\to} G\beta\gamma-Ste5-Ste5-G\beta\gamma$

2. ODE Equations

   $$\begin{align} \frac{d[Ste5-Ste5]}{dt} &= K_{SS}[Ste5]^2 \\ \frac{d[G\beta\gamma-Ste5]}{dt} &= K_{SG}[G\beta\gamma][Ste5] \\ \frac{d[G\beta\gamma-Ste5-Ste5]}{dt} &= K_{GSS1}[G\beta\gamma][Ste5]^2 + K_{GSS2}[G\beta\gamma-Ste5][Ste5] \\ \frac{d[G\beta\gamma-Ste5-Ste5-G\beta\gamma]}{dt} &= K_{GSSG1}[G\beta\gamma-Ste5]^2 + K_{GSSG2}[G\beta\gamma-Ste5-Ste5][G\beta\gamma] \\ \frac{d[Ste5]}{dt} &= -2K_{SS}[Ste5]^2 - K_{SG}[G\beta\gamma][Ste5] \\ &\quad -2K_{GSS1}[G\beta\gamma][Ste5]^2 - K_{GSS2}[G\beta\gamma-Ste5][Ste5] \\ [G\beta\gamma] &= (hH3R-Ste2)_{stable} \end{align}$$

3. Parameter and Variable Meanings

   Symbol	Meaning
   $K_{SS}$	Reaction rate constant: Rate constant for two free Ste5 molecules to form a Ste5-Ste5 dimer
   $K_{SG}$	Reaction rate constant: Rate constant for Gβγ and Ste5 to form the Gβγ-Ste5 complex
   $K_{GSS1}$	Reaction rate constant: Rate constant for Gβγ and two free Ste5 molecules to form Gβγ-Ste5-Ste5
   $K_{GSS2}$	Reaction rate constant: Rate constant for Gβγ-Ste5 and free Ste5 to form Gβγ-Ste5-Ste5
   $K_{GSSG1}$	Reaction rate constant: Rate constant for two Gβγ-Ste5 molecules to form Gβγ-Ste5-Ste5-Gβγ
   $K_{GSSG2}$	Reaction rate constant: Rate constant for Gβγ-Ste5-Ste5 and Gβγ to form Gβγ-Ste5-Ste5-Gβγ
   $G\beta\gamma$	Activated G protein βγ subunits, whose steady-state concentration is equal to the steady-state concentration of the hH3R-Ste2 complex
   $Ste5-Ste5$	Homodimer formed by Ste5
   $G\beta\gamma-Ste5$	Gβγ-Ste5 scaffold complexes in different assembly forms

3. MAPK Cascade Reaction Module

The scaffold complex recruits Ste11 (MAPKKK), Ste7 (MAPKK), and Fus3 (MAPK), and realizes signal amplification and transmission through sequential phosphorylation.

1. Reaction Equations

   • Binding and dissociation of Ste5 and Ste11: $Ste5_{off-st1} + Ste11_{off} \underset{K_{offSte5-Ste11}}{\stackrel{K_{onSte5-Ste11}}{\rightleftharpoons}} Ste5-Ste11$

   • Phosphorylation of Ste11: $Ste11 \stackrel{ste20}{\to} Ste11pS \stackrel{ste20}{\to} Ste11pSpS \stackrel{ste20}{\to} Ste11pSpSpT$

   • Phosphorylation of Ste7: $Ste7 \xrightarrow{Ste11pS, Ste11pSpS, Ste11pSpSpT} Ste7pS \xrightarrow{Ste11pS, Ste11pSpS, Ste11pSpSpT} Ste7pSpT$

   • Phosphorylation of Fus3: $Fus3 \stackrel{Ste7pS, Ste7pSpT}{\to} Fus3pY \stackrel{Ste7pS, Ste7pSpT}{\to} Fus3pYpT$

2. ODE Equations

   • Kinetics of Binding and Dissociation Between Ste5 and Ste11

     $$\begin{align} \frac{d[Ste5_{off-st1}]}{dt} &= K_{offSte5-Ste11}[Ste5-Ste11] - K_{onSte5-Ste11}[Ste5_{off-st1}][Ste11_{off}] \\ \frac{d[Ste11_{off}]}{dt} &= K_{offSte5-Ste11}[Ste5-Ste11] - K_{onSte5-Ste11}[Ste5_{off-st1}][Ste11_{off}] \\ \frac{d[Ste5-Ste11]}{dt} &= -K_{offSte5-Ste11}[Ste5-Ste11] + K_{onSte5-Ste11}[Ste5_{off-st1}][Ste11_{off}] \end{align}$$

   • Phosphorylation Kinetics (Taking Key Molecules as Examples)

     $$\begin{align} \frac{d[Ste11pSpSpT]}{dt} &= K_{cat\_Ste11}[Ste11pSpS] - K_{dephos}[Ste11pSpSpT] \\ \frac{d[Ste7pS]}{dt} &= K_{cat\_Ste11pSSte7pS}[Ste11pS] \cdot \frac{[Ste11_{total}]}{[Ste5-Ste11]} \cdot \frac{[Ste5-Ste7]}{[Ste7_{total}]} \cdot \left(\frac{[Ste7pS]}{[Ste7_{total}]}\right)^2 + ... \\ \frac{d[Fus3pYpT]}{dt} &= K_{cat\_Fus3}[Fus3pY] - K_{dephos}[Fus3pYpT] \end{align}$$

3. Parameter and Variable Meanings

   Symbol	Meaning
   $K_{onSte5-Ste11}$	Binding rate constant: Rate constant for free Ste5 ($Ste5_{off-st1}$) and free Ste11 ($Ste11_{off}$) to bind
   $K_{offSte5-Ste11}$	Dissociation rate constant: Rate constant for the dissociation of the Ste5-Ste11 complex into free Ste5 and Ste11
   $K_{cat\_Ste11}$	Catalytic rate constant: Efficiency constant for Ste11pSpS to catalyze the formation of Ste11pSpSpT
   $K_{cat\_Ste11pSSte7pS}$	Catalytic rate constant: Efficiency constant for Ste11pS to catalyze the formation of Ste7pS from Ste7
   $K_{dephos}$	Dephosphorylation rate constant: Unified dephosphorylation rate constant for phosphorylated molecules such as Ste11, Ste7, and Fus3
   $Ste5_{off-st1}$	Free scaffold protein Ste5 not bound to Ste11
   $Ste11_{off}$	Free kinase Ste11 not bound to Ste5
   $Ste11pS/Ste11pSpS/Ste11pSpSpT$	Phosphorylated forms of Ste11 at the S302 site, double S sites, and double S+T sites, respectively
   $Ste7pS/Ste7pSpT$	Phosphorylated forms of Ste7 at the S359 site and S+T sites, respectively
   $Fus3pY/Fus3pYpT$	Phosphorylated forms of Fus3 at the Y site and Y+T sites, respectively
   $Ste11_{total}/Ste7_{total}$	Total concentrations of Ste11 and Ste7 in the system (assumed to be constant)

4. Ste12 Activation Module

After the nuclear translocation of phosphorylated Fus3, the Ste12-Dig complex dissociates, releasing free Ste12 to initiate downstream transcription.

1. Reaction Equations

   • Nuclear translocation of Fus3: $Fus3_{cyt} \stackrel{K_{in}}{\to} Fus3_{nuc}$

   • Nuclear export of Fus3: $Fus3_{nuc} \stackrel{K_{out}}{\to} Fus3_{cyt}$

   • Dissociation of Ste12-Dig: $Fus3p + Ste12-Dig1-Dig2 \stackrel{K_{phos1}/K_{phos2}}{\to} Fus3 + Dig1 + Dig2 + Ste12$

   • Rebinding of Ste12: $Ste12 + Dig1 + Dig2 \stackrel{K_{bind}}{\to} Ste12-Dig1-Dig2$

2. ODE Equations

   $$\begin{align} \frac{d[Fus3_{cyt}]}{dt} &= -K_{in}[Fus3_{cyt}] + K_{out}[Fus3_{nuc}] \\ \frac{d[Fus3_{nuc}]}{dt} &= K_{in}[Fus3_{cyt}] - K_{out}[Fus3_{nuc}] \\ \frac{d[C]}{dt} &= -(K_{phos1} + K_{phos2})[Fus3_{nuc}][C] \\ \frac{d[Dig1]}{dt} &= K_{phos1}[Fus3_{nuc}][C] - K_{dephos1}[Msg5][Dig1] \\ \frac{d[Dig2]}{dt} &= K_{phos2}[Fus3_{nuc}][C] - K_{dephos2}[Ptp2][Dig2] \\ \frac{d[Ste12]}{dt} &= (K_{phos1} + K_{phos2})[Fus3_{nuc}][C] - K_{bind}[Ste12][Dig1][Dig2] \end{align}$$

3. Parameter and Variable Meanings

   Symbol	Meaning
   $K_{in}$	Transport rate constant: Rate constant for phosphorylated Fus3 to enter the nucleus ($Fus3_{nuc}$) from the cytoplasm ($Fus3_{cyt}$)
   $K_{out}$	Transport rate constant: Rate constant for phosphorylated Fus3 to be exported from the nucleus to the cytoplasm
   $K_{phos1}/K_{phos2}$	Phosphorylation rate constants: Rate constants for Fus3 to catalyze the phosphorylation of Dig1 and Dig2, respectively
   $K_{dephos1}/K_{dephos2}$	Dephosphorylation rate constants: Rate constants for Msg5 to catalyze the dephosphorylation of Dig1 and Ptp2 to catalyze the dephosphorylation of Dig2, respectively
   $K_{bind}$	Binding rate constant: Rate constant for free Ste12 to rebind with phosphorylated Dig1 and Dig2
   $C$	Ste12-Dig1-Dig2 ternary complex (inactive form of Ste12 bound to Dig1 and Dig2)
   $Fus3_{cyt}/Fus3_{nuc}$	Concentrations of phosphorylated Fus3 in the cytoplasm and nucleus, respectively
   $Msg5/Ptp2$	Phosphatases, whose concentrations are assumed to be constant
   $Dig1/Dig2$	Inhibitory proteins bound to Ste12, which dissociate from Ste12 after phosphorylation

5. Transcription-Translation Module

Free Ste12 acts as a transcription factor to initiate the transcription of the target gene (GeneA). The transcribed mRNA is then translated into protein, accompanied by the degradation of nucleic acid and protein.

1. Reaction Equations

   • Transcription: $GeneA + Ste12 \xrightarrow{K_{trans}} mRNA_A + Ste12$

   • mRNA degradation: $mRNA_A \stackrel{K_{degmRNA}}{\to} Degradation\ Products$

   • Translation: $mRNA_A \xrightarrow{K_{transl}} ProteinA$

   • Protein degradation: $ProteinA \stackrel{K_{degprotein}}{\to} Degradation\ Products$

2. ODE Equations

   $$\begin{align} \frac{d[mRNA_A]}{dt} &= K_{trans} \cdot \frac{[Ste12]^n}{[Ste12]^n + K_d^n} - K_{degmRNA}[mRNA_A] \\ \frac{d[ProteinA]}{dt} &= K_{transl}[mRNA_A] - K_{degprotein}[ProteinA] \end{align}$$

3. Parameter and Variable Meanings

   Symbol	Meaning
   $K_{trans}$	Transcription rate constant: Rate constant for Ste12 to initiate the transcription of GeneA into $mRNA_A$
   $K_{transl}$	Translation rate constant: Rate constant for $mRNA_A$ to be translated into ProteinA
   $K_{degmRNA}$	Degradation rate constant: Degradation rate constant of $mRNA_A$
   $K_{degprotein}$	Degradation rate constant: Degradation rate constant of ProteinA
   $K_d$	Dissociation constant: Reflects the binding characteristics of Ste12 to the regulatory sequence of GeneA
   $n$	Hill coefficient: Reflects the cooperative effect of Ste12 binding, set to $n=1$ (no cooperative effect) in this model
   $mRNA_A$	Messenger RNA transcribed from GeneA
   $ProteinA$	Target protein translated from $mRNA_A$ (AiiA enzyme or adhesion molecule)

Parameter Table

The model parameters, integrating literature reports and reasonable assumptions, cover two core types: reaction rate constants and initial concentrations, as detailed below: