Model

General Description of Modeling

Our model serves two main purposes:

1. Quantitative Description of Project Design: Due to safety considerations, we were unable to conduct animal experiments to demonstrate the processes occurring during the operation of the project. Modeling can help in understanding therapeutic pathways, provide a quantitative perspective, and tell our story better.
2. Computational Methods for Project Engineering: If the project can be carried out, the model can help determine the parameters in the implementation process of the project, reduce the calculation amount in the experimental process, connect the wet experimental independent system , and make the design mathematically encapsulated as new components.

Our model can be divided into four interconnected parts, representing the inhalation of muscone, its binding to receptors, intracellular signal transduction and lactate secretion triggered by receptor activation, and the absorption of lactate. These models provide a comprehensive understanding of the project and yield valuable computational results.

Compartment Model for Muscone Inhalation

1.Model Description

The main focus of our project is the use of muscone as a signaling molecule to activate engineered yeast in the gut for therapeutic purposes. Therefore, it is crucial to provide a quantitative description and computational support for the diffusion of muscone in the body. This model describes the entire process from the inhalation of muscone to its increased concentration in the intestinal tract. We will establish a multi-compartment model that includes the following main processes:

1. Inhalation Process: Muscone is inhaled in the form of an aerosol into the lungs.
2. Pulmonary Process: Muscone distributes in the alveoli and may be exhaled, adhered to, or permeated into the microvessels.
3. Adhesion Process: A portion of muscone adheres to the respiratory mucosa and then diffuses into the systemic circulation.
4. Alveolar Microvessel Process: Muscone permeates into the alveolar microvessels and gradually enters the systemic circulation.
5. Systemic Circulation Process: Muscone distributes in the systemic circulation and is transported to various parts of the body through the bloodstream.
6. Intestinal Process: Muscone enters the target intestine through the mesenteric microvascular network, where its concentration begins to increase.

Fig 2 Processes in the Inhalation Model

Corresponding to the above processes, five compartments need to be established for simulation, where \(t\) represents the time variable:

2. Initial Settings and Assumptions

At \(t=0\), the amount of muscone in all compartments is \(0\).

Assuming that the total amount of inhaled muscone is \(Q_{\text{inhale}}\) (\(\text{mg}\)), which is assumed to be \(100\text{mg}\). Only \(0.5\%\) of muscone enters the systemic circulation through adhesion. In this model, since muscone only acts as a signaling molecule to activate yeast to synthesize lactate, we only consider the metabolism and excretion of muscone in the systemic circulation. We only focus on the short-term process of muscone appearing in the intestine from scratch, and the subsequent process of reaching a certain concentration can be ignored.

3. Model Equations

Inhalation Equation for Muscone

\[ V_{\text{inhale}}(t) =\frac{Q_{\text{inhale}}}{5}(u(t)-u(t-5)) \]

Explanation: This describes the rate equation for inhaling muscone over five seconds, where the total amount \( Q \) remains constant. The function \( u(t) \) is a step function, which takes the value of \( \frac{Q_{\text{inhale}}}{5} \) from \( t=0s \) to \( t=5s \), and is \( 0 \) otherwise, simulating the scenario of resting human respiration.

Compartment 0: \( Q_A(t) \)

\[ \frac{dQ_A(t)}{dt} = V_{\text{inhale}}(t) - \left( k_{\text{exhale}} + k_{\text{perm}} \right) Q_A(t) \]

Explanation: The amount of muscone in the alveoli increases through inhalation and decreases due to exhalation, adhesion to the respiratory mucosa, and permeation into the alveolar capillaries.

Parameters:

• \( k_{\text{exhale}} \): Since most of the muscone is rapidly exhaled, this value is relatively large, taken as \( 10 \ \text{min}^{-1} \)
• \( k_{\text{perm}} \): The rate of muscone permeation into the capillaries, affected by its physicochemical properties, is taken as \( 0.005 \ \text{min}^{-1} \)

Compartment 1: \( Q_M(t) \)

\[ \frac{dQ_M(t)}{dt} = 0.0005 \cdot k_{\text{adh}} V_{\text{inhale}}(t) - k_{\text{diffMC}} Q_M(t) \]

Explanation: The increase in muscone on the mucosa comes from adhesion in the alveoli, and the decrease is due to diffusion into the systemic circulation.

Parameters:

• \( k_{\text{adh}} \): The adhesion process is relatively slow, and only \( 0.5\% \) of muscone enters the systemic circulation through this pathway, taken as \( 0.001 \ \text{min}^{-1} \)
• \( k_{\text{diffMC}} \): Diffusion from the mucosa to the systemic circulation is slow, taken as \( 0.01 \ \text{min}^{-1} \)

Compartment 2: \( Q_L(t) \)

\[ \frac{dQ_L(t)}{dt} = k_{\text{perm}} Q_A(t) - k_{\text{diffLC}} Q_L(t) \]

Explanation: The increase in muscone in the alveolar capillaries comes from permeation in the alveoli, and the decrease is due to diffusion into the systemic circulation.

Parameters:

• \( k_{\text{perm}} \): Same as Compartment 0
• \( k_{\text{diffLC}} \): The diffusion rate from alveolar capillaries to the systemic circulation is relatively slow, taken as \( 0.05 \ \text{min}^{-1} \)

Compartment 3: \( Q_C(t) \)

\[ \frac{dQ_C(t)}{dt} = k_{\text{diffMC}} Q_M(t) + k_{\text{diffLC}} Q_L(t) - k_{\text{dist}} Q_C(t) - k_{\text{excrete}} Q_C(t) \]

Explanation: The increase in muscone in the systemic circulation comes from the input of mucosa and alveolar capillaries, and the decrease is due to distribution to the intestinal mesenteric microvascular network and excretion through various routes.

Parameters:

• \( k_{\text{diffMC}} \): Same as Compartment 1
• \( k_{\text{diffLC}} \): Same as Compartment 2
• \( k_{\text{dist}} \): The rate constant of muscone distribution from the systemic circulation to the intestinal mesenteric microvascular network, taken as \( 0.001 \ \text{min}^{-1} \)
• \( k_{\text{excrete}} \): Muscone is excreted from the systemic circulation through epidermal volatilization, urine, continuous respiration, etc., taken as \( 0.05 \ \text{min}^{-1} \)

Compartment 4: \( Q_I(t) \)

\[ \frac{dQ_I(t)}{dt} = k_{\text{dist}} Q_C(t) - k_{move}Q_I(t) \]

Explanation: The increase in muscone in the intestine comes from the distribution of the systemic circulation, and the decrease is due to metabolism and excretion through intestinal fluid and peristalsis.

\( k_{\text{dist}} \): Same as Compartment 3
\( k_{move} \): The metabolism and excretion of muscone in the intestine, taken as \( 0.02 \ \text{min}^{-1} \)

4. System of Equations:

In summary, we can write a system of ordinary differential equations and import it into MATLAB for simulation:

\[ \begin{align*} Q_{\text{inhale}}(t) & = 100(mg)(Assumption) \\ V_{\text{inhale}}(t) & =\frac{Q_{\text{inhale}}}{5}(u(t)-u(t-5)) \\ \frac{dQ_A(t)}{dt} & = V_{\text{inhale}}(t) -\left( k_{\text{exhale}} + k_{\text{perm}} \right) Q_A(t) \\ \frac{dQ_L(t)}{dt} & = k_{\text{perm}} Q_A(t) - k_{\text{diffLC}} Q_L(t) \\ \frac{dQ_M(t)}{dt} & = 0.0005\cdot k_{\text{adh}} V_{\text{inhale}}(t) - k_{\text{diffMC}} Q_M(t) \\ \frac{dQ_C(t)}{dt} & = k_{\text{diffMC}} Q_M(t) + k_{\text{diffLC}} Q_L(t) - k_{\text{dist}} Q_C(t) - k_{\text{excrete}} Q_C(t) \\ \frac{dQ_I(t)}{dt} & = k_{\text{dist}} Q_C(t)-k_{move}Q_I(t) \\ \end{align*} \]

Fig 3 Stimulation Result of the Muscone Inhalation Model

We simulated the distribution of muscone in the systemic circulation and obtained the concentration change curve of muscone in the systemic circulation. According to the model, after one breath, traces of muscone can spread into the intestine, similarly, the concentration change caused by continuous muscone is simulated by changing the inhalation equation, and the concentration of muscone in the intestine can be obtained in combination with experiment. Because there is no animal experimental support, the data are manually drafted, and the calculation method is more meaningful than the calculation results.

Molecular Dynamics Simulations of Muscone with its Receptor

One of the major contributions of our project was the development of a signaling pathway activated by gas molecules acting as switches, so we were interested in the behavior of the muscone when bound to its receptor. Molecular dynamics simulations can help us to visualize the binding process, track the parameters of structural changes, protein folding, and receptor-ligand interactions, and demonstrate their biological significance.

1. Molecular Model Construction and Preparation

• Before conducting molecular dynamics simulations, it is essential to accurately construct and prepare the three-dimensional molecular model of the research object. In this study, our goal is to simulate the interaction between muscone and the olfactory receptor Or5an6 (MOR215-1).

Constructing the Three-Dimensional Structures of Muscone and the Receptor:

• Muscone:
• Obtain the molecular structure information of muscone from the PubChem database. PubChem is a widely used chemical information database that provides data on molecular structures, chemical properties, and biological activities.
• Molecular structure obtained from PubChem link


Fig 4 The structural formula of the Muscone

• Next, use the chemical modeling tool Avogadro to convert this structure into a PDB format file (muscone.pdb) for further simulation and analysis. This step ensures that the three-dimensional geometric structure of muscone is accurate and suitable for subsequent molecular docking and dynamics simulations.
• Avogadro converts to pdb format muscone.pdb


Fig 5 The 3D geometry of Muscone

• Receptor Or5an6 (MOR215-1):
  • Retrieve the nucleic acid sequence of the MOR215-1 gene from the NCBI database and find its protein sequence information. With support from Uniprot, further obtain the amino acid sequence of this protein, which is an important component for subsequent docking simulations. In addition, obtain the three-dimensional structural model of the Or5an6 receptor (MOR215-1.pdb) from the AlphaFold database, whose model quality is widely recognized, especially in the field of protein structure prediction.
  • Nucleic acid sequence information in NCBI link
  DNA sequence of Or5an6
  >NC_000085.7:12371629-12372576 Mus musculus strain C57BL/6J chromosome 19, GRCm39 ATGCCTGGAGGGAGGAATAGCACAGTCATCACCAAGTTCATCCTTGTGGGATTCTCAGATTTTCCAAAGC TCAAGCTGGTTCTCTTTGTTATCTTCCTGGGAAGTTATCTCTCCACAGTGGTGTGGAACTTGGGCCTCAT CATCTTGATTAGGATTGACCCTTACCTACACACACCTATGTACTTCTTCCTCAGCAATTTGTCATTTTTA GATTTCTGTTACATTTCATCTACAACCCCTAAAATGCTCTCGGGATTCTTCCAGAAGTCTAAATCTATCT CCTTTGTTGGGTGCACCATGCAGTACTTCATCTTCTCAAGCCTGGGTCTGTCCGAATGCTGCCTTCTGGC AGCCATGGCTTATGACCGGTATGCTGCCATTTGTAATCCTCTTCTCTACACAGCCATCATGTCCCCGTCA CTCTGTGTGCACATGGTGGTTGGAGCCTATAGTACTGGTCTCTTGGGTTCATTGATTCAACTGTGTGCTA TACTTCAGCTCCATTTCTGTGGGCCAAATATTATAAACCATTTCTTTTGTGACCTGCCTCAGCTATTAGT TCTTTCCTGCTCTGAAACCTTTCCCCTGCAAGTCTTGAAATTTGTAATAGCAGTGATTTTTGGGGTGGCA TCTGTCATTGTTATCCTGATATCCTATGGTTATATCATTGGCACAATCCTGAATATCAGCTCAGTAGAAG GTAGGTCCAAGGCATTCAATACCTGTGCCTCTCACCTGACAGCAGTCACCCTCTTTTTTGGATCAGGACT CTTTGTCTATATGCGCCCCAGCTCCAACAGTTCCCAGGGTTATGACAAGATGGCTTCCGTGTTCTATACA GTGGTGATTCCCATGTTGAATCCTCTGATTTATAGTCTCAGGAACAAGGAAATAAAAGATGCTCTTCAGA GATGTAAAAATAAGTGCTTTTCTCAGTGCCACTGTTAG
  • Protein sequence information in Uniprot link
  Protein sequence of Or5an6
  >sp|Q8VFV4|O5AN6_MOUSE Olfactory receptor 5AN6 OS=Mus musculus OX=10090 GN=Or5an6 PE=2 SV=1 MPGGRNSTVITKFILVGFSDFPKLKLVLFVIFLGSYLSTVVWNLGLIILIRIDPYLHTPM YFFLSNLSFLDFCYISSTTPKMLSGFFQKSKSISFVGCTMQYFIFSSLGLSECCLLAAMA YDRYAAICNPLLYTAIMSPSLCVHMVVGAYSTGLLGSLIQLCAILQLHFCGPNIINHFFC DLPQLLVLSCSETFPLQVLKFVIAVIFGVASVIVILISYGYIIGTILNISSVEGRSKAFN TCASHLTAVTLFFGSGLFVYMRPSSNSSQGYDKMASVFYTVVIPMLNPLIYSLRNKEIKD ALQRCKNKCFSQCHC
  • Protein structure information in the AlphaFold database link MOR215-1.pdb

Fig 6 Protein structure of MOR215-1

System Preparation:

• To study how muscone binds to the receptor, molecular docking tools such as AutoDock and Vina are used to determine potential binding conformations and obtain docking data:
  • Mac software installation: AutoDock, Vina, XQuartz, MGLTools,
  • The 32-bit version of MGLTools cannot be used directly on Mac, but can be opened using open /Users/Shared/MGLTools/1.5.7/bin/adt
  • AutoDock generates PDBQD files: MOR215_1.pdbqt and muscone.pdbqt
    • In the same directory: adt, autodock4, autogrid4, vina
  • Vina configures the docking parameter file config.txt in the same directory as follows (based on actual coordinates), runs the docking simulation, and generates the file muscure.pdbqt. This file records the possible binding modes of muscone with the receptor at the defined coordinate position and dimensions (size_x, size_y, size_z). config.txt
    receptor = MOR215_1.pdbqt ligand = muscone.pdbqt center_x = -0.112 center_y = -0.219 center_z = -2.092 size_x = 40 size_y = 40 size_z = 40 exhaustiveness= 8 out=muscure.pdbqt

    ./vina --config config.txt

    Obtain muscure.pdbqt

• Visualization Structure Optimization and Inspection:
  • Use Pymol to analyze the generated docking structures, identify potential binding conformations, especially focusing on interactions at amino acids Arg-51 and Tyr-271, and determine stable low free energy conformations. Then extract this conformation as separate PDB files, isolating MOR.pdb and MUS.pdb.
    • Open muscure.pdbqt and MOR215-1.pdb
    • Select the appropriate conformation, display interactions Action - find - polar contacts - to any atoms
    • Add Label, export the image
  • Alternative Conformations
    • Arg-51
      

      Fig 7 Alternative Conformation with Arg-51

    • Tyr-271
      

      Fig 8 Alternative Conformation with Tyr-271

    • Select conformations with low free energy that are more likely, extract in pdb format, isolating MOR.pdb and MUS.pdb

2. Force field parameterization

Select Force Field:

• To perform molecular dynamics simulations, it is necessary to choose an appropriate molecular force field to describe the interactions between molecules within the system. CHARMM36 was selected as the force field for protein and small molecule (musk ketone) interactions. However, due to the absence of direct parameters for musk ketone in existing force fields, custom parameters need to be generated to supplement it.

Generate Force Field Parameters:

• Use Avogadro to convert to .mol2 format, adjust file information, and then use the software CGenFF to generate its CHARMM36 force field parameter file MUS_fix.mol2 and parameter file MUS.str. This step includes calculating the sorting of the chemical information file (sort_mol2_bonds.pl) to ensure file correctness.

perl sort_mol2_bonds.pl MUS.mol2 MUS_fix.mol2

3. Preprocessing

Build the system:

• Generate the topology file MOR_processed.gro for the receptor using GROMACS's pdb2gmx command.

gmx pdb2gmx -f MOR.pdb -o MOR_processed.gro -ter

• Convert the force field parameters of muscone to a format recognizable by GROMACS to generate its topology data using the CGenFF helper script.

python cgenff_charmm2gmx_py3_nx2.py MUS MUS_fix.mol2 MUS.str charmm36-jul2022.ff

Merge the system:

• Prepare the complete solvent system required for simulations using the editconf and solvate commands, merging the topology files of muscone mus.gro and receptor MOR_processed.gro into a single system complex.gro.

gmx editconf -f mus_ini.pdb -o mus.gro
gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0

• Solvate the system in solvent (such as water). From topol.top, it is known that the net charge is 9; add counterions to neutralize the system.

gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

[ molecules ]
; Compound        #mols
Protein_chain_A     1
MUS                 1
SOL         31227
CL               9

Energy minimization:

• Perform energy minimization on the overall system to eliminate unreasonable conflicts in the initial geometry. Achieve rapid convergence of energy through the gradient descent algorithm and ensure all atoms in the system are reasonably positioned within the force field. Analyze the energy minimization results to ensure both the maximum force and potential energy are within reasonable thresholds.

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

Steepest Descents converged to Fmax < 1000 in 1182 steps
Potential Energy  = -1.5345670e+06
Maximum force     =  8.9675085e+02 on atom 4987
Norm of force     =  1.3456365e+01

gmx energy -f em.edr -o potential.xvg
#11 0
xmgrace potential.xvg
dit xvg_show -f potential.xvg



Fig 9 Potential Energy Minimization

4. Molecular Dynamics Simulation

System Equilibration:

• To achieve thermal and mechanical equilibrium of the system, simulations are conducted in two stages: NVT (constant temperature) and NPT (constant pressure) equilibration. The system temperature is gradually increased to the target of 300K to reach stable conditions, while analyzing the curves of temperature, pressure, and density over time to ensure the stability of the system.

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -deffnm nvt
gmx energy -f nvt.edr -o temperature.xvg
#16 0
dit xvg_show -f temperature.xvg

gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr -maxwarn 1
gmx mdrun -deffnm npt
gmx energy -f npt.edr -o pressure.xvg
#17 0
gmx energy -f npt.edr -o density.xvg
#23 0
dit xvg_show -f temperature.xvg

Fig 10 Curve of the pressure over time

Fig 11 Curve of the density over time

Production Simulation:

• Under the conditions of equilibrium, a long-term production simulation is conducted. This simulation observes the time evolution characteristics of the dynamic interactions between muscone and the receptor, typically requiring simulation times ranging from hundreds of nanoseconds to several microseconds to ensure the reliability and reproducibility of the results. By sampling key frame data of the system during the dynamic process, it aids in further analyzing the intermolecular interactions and dynamic conformational changes.

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

Fig 12 Molecular dynamics simulation process

Fig 13 Results of the molecular dynamics simulations

5. Post-Processing and Analysis

• After the completion of molecular dynamics simulations, detailed post-processing and analysis of the data is a key step in understanding the behavior of molecular systems. This section will demonstrate the use of GROMACS tools to analyze trajectory data, and visualize it using VMD and PyMOL.

Trajectory Analysis:

To gain deeper insights into the interactions between muscone and the receptor, visualization tools are used to make the simulation process intuitive, identifying key interaction sites and structural changes.

• Using the gmx trjconv tool, the calculated trajectory data is centered and periodic boundary conditions are removed, generating the centered trajectory file md_0_10_center.xtc and its initial frame start.pdb.

gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact
#1 0
gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o start.pdb -dump 0
#0

• In PyMOL, key interactions of the protein-ligand system can be visualized through selection and rotation commands, identifying important residues and binding pockets that play significant roles during the simulation process.

# Pymol
select water, resn SOL
select ions, resn CL
select protein, not water and not ions
select ligand, resn MUS
deselect
cmd.rotate('x', 45)
cmd.rotate('y', 45)

Fig 14 Visualization of the protein-ligand system

• Finally, generate a fitted trajectory file suitable for analysis and animation production.

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans
#4 0
gmx trjconv -s md_0_10.tpr -f md_0_10_fit.xtc -o traj.pdb -dt 10
#0

Fig 15 Trajectory Analysis

RMSD (Root Mean Square Deviation) Analysis

• RMSD provides a fundamental metric for measuring structural deviation during the simulation process. By calculating the deviation of structures at different time points, it is possible to assess the stability of the system.

gmx distance -s md_0_10.tpr -f md_0_10_center.xtc -select 'resname "MUS" and name O plus resid 51 and name NH1' -oall
gmx make_ndx -f em.gro -o index.ndx
> 13 & a O  
> 1 & r 51 & a NH1
> 1 & r 51 & a CZ
> 20 | 21 | 22
> q
gmx angle -f md_0_10_center.xtc -n index.ndx -ov angle.xvg 
gmx make_ndx -f em.gro -n index.ndx
> 13 & ! a H*
> name MUS_Heavy
> q
gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_mus.xvg

xmgrace rmsd_mus.xvg
dit xvg_show -f rmsd_mus.xvg
xmgrace rmsd_mus.xvg

• From the RMSD curve, it can be observed that the system stabilizes after approximately 1 nanosecond, with a fluctuation range of around 1.3 Å. This indicates that the structures of the protein and ligand remained relatively stable throughout the simulation process, with no significant conformational changes occurring, supporting the potential biological significance of the results.

Fig 16 RMSD Analysis

Radius of Gyration (Rg) Calculation

• The radius of gyration (Rg) is used to assess the compactness of a protein and is an important indicator of protein folding or unfolding.
• If the protein folding is stable, its radius of gyration Rg will maintain a relatively stable value. If the protein unfolds, its Rg will change over time.

gmx gyrate -s md_0_10.tpr -f md_0_10_fit.xtc -o gyrate.xvg
#1
xmgrace gyrate.xvg

• In the simulation, the Rg value of the protein remained between 2.2 and 2.25 nanometers, indicating that the protein maintained a compact folded state during the simulation period at 300 K, with no significant unfolding or expansion occurring.

Fig 17 Radius of Gyration Calculation

Protein-Ligand Interaction Energy

• By calculating the Coulomb and Lennard-Jones interaction energies within the system, the binding energy between musk ketone and the receptor is quantified.

gmx make_ndx -f em.gro -o index.ndx
> 1 | 13
gmx grompp -f ie.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o ie.tpr
gmx mdrun -deffnm ie -rerun md_0_10.xtc -nb cpu

gmx energy -f ie.edr -o interaction_energy.xvg

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Coul-SR:Protein-MUS       0.0439222       0.38    2.88794    1.41131  (kJ/mol)
LJ-SR:Protein-MUS          -81.3724        1.6    9.75743    1.99715  (kJ/mol)
#21 | 22

xmgrace interaction_energy.xvg
dit xvg_show -f interaction_energy.xvg

Fig 18 Protein-Ligand Interaction Energy

• The calculation results show that the Lennard-Jones interaction dominates in stable binding, averaging -81.37 kJ/mol, indicating that hydrophobic interactions play a central role in the binding between the ligand and the receptor.

Ordinary Differential Equation of the signal transduction of the yeast MAPK pathway

1. Model Description

In our project, we express the muscone receptor (GPCR) on the yeast cell membrane. After a certain concentration of muscone diffuses into the intestine and binds to the receptor, it activates the receptor, which in turn activates the G protein. The G protein dissociates into α and βγ subunits, with the βγ subunit releasing and activating Ste20 and the scaffold protein Ste5. Ste5 can undergo oligomerization and other behaviors, recruiting Ste11, Ste7, and Fus3 near the plasma membrane. The cascade reaction is initiated by Ste20, and the signal is transmitted along the Ste11-Ste7-Fus3 cascade. Fus3 activates the transcription factor pFUS1, and the downstream gene is LahA, which expresses lactate dehydrogenase LDH, catalyzing the conversion of pyruvate to lactate. This model simulates the changes in the concentrations and phosphorylation states of molecules in the signaling transduction pathway by writing out chemical reactions and converting them into ordinary differential equations, in order to obtain the quantitative relationship between muscone activation and lactate secretion. The model includes the following main processes:

1. Activation of Muscone Receptor: The muscone receptor Ste2, derived from mouse olfactory epithelium, is a G protein-coupled receptor (GPCR) that is expressed on the cell membrane and receives signals. Its domains consist of α, β, and γ, where the Gα subunit is called Gpa1, and the Gα and Gγ subunits are Ste4 and Ste18, respectively, both anchored in the cell membrane, without discussing the scenario of their separation. After binding with muscone, Gpa1 will release Ste4-Ste18.
2. Formation of Scaffold: The released Ste4-Ste18 can bind to Ste5, and the Ste5 protein can undergo dimerization, oligomerization, and other behaviors, forming a scaffold near the cell membrane and recruiting proteins related to the cascade phosphorylation.
3. Cascade Reaction: The scaffold composed of Ste5 can recruit Ste11 (MAPKKK), Ste7 (MAPKK), and Fus3 (MAPK). Each of these three proteins has multiple phosphorylation modification sites, and the efficiency of catalyzing phosphorylation varies under different modification scenarios. Furthermore, the three proteins independently bind to Ste5, and a reaction can only occur when two adjacent proteins are simultaneously present on the scaffold, making this signaling pathway highly specific.
4. Activation of pFUS1: The transcription factor pFUS1 is activated by Fus3, and the downstream gene is LahA, which expresses lactate dehydrogenase to produce lactate.

2. Basic Assumptions

1. Since the model only simulates the signal transduction shortly after muscone activation, it does not consider protein synthesis and degradation, assuming that the concentrations of each protein remain stable during this time.
2. It is assumed that all proteins involved in the cascade reaction have the same dephosphorylation rate, denoted by \(k_{cat_{dephosph}}\).
3. The behavior of all molecules in the system is random and not influenced by environmental factors.

3. Model Equations

Activation of muscone Receptor

After Ste2 binds with muscone, it interacts with the G protein, causing the exchange of GDP bound to the G protein with GTP in the cytoplasm, releasing Ste4 and Ste18. After Gpa1 catalyzes the conversion of GTP to GDP, it can return to the cytoplasm and rebind, forming a G protein trimer. Since the original signaling pathway is the yeast pheromone signaling pathway, with the ligand being the pheromone, this section uses Pheromone to represent the molecules that activate the receptor.

Table 1: Variables of Receptor Activation Model

Table 2: Parameters of Receptor Activation Model

There are \(0.3{\mu M}\) of Pheromone and \(1{\mu M}\) of inactive G proteins. Known variables are entered, other variables are set to zero, and unknown parameters are defined. After starting the simulation, reactions occur according to the equations listed.

Fig 19 Receptor Activation

Formation of the Scaffold

Explanation: The binding of Ste4Ste18 with Ste5 and the oligomerization of Ste5 is a process that is not completely independent. Many equations can be derived through combinations, but here we only consider the dimerization process, and each reaction is reversible. Since Ste5 actually binds to Ste4, we abbreviate Ste5 as S5 and Ste4 as S4 in the equations.

Table 3: Variables of Scaffold Formation Model

Table 4: Parameters of Scaffold Formation Model

Assume that before signal transduction starts, there are only free Ste5 and just released Ste4Ste18 in the cell, with concentrations both equal to 1, and parameters are assumed. After starting the simulation, reactions occur according to the listed equations, and after a period of time, the concentrations reach equilibrium.

Fig 20 Scaffold Formation

Cascading Reactions

Reactions:

Explanation

Only the Ste5 bound to the scaffold has significance in recruiting Ste11, Ste7, and Fus3, and the binding to these three proteins is independent. Therefore, the Ste5 on the scaffold can be treated as three copies to calculate its binding with Ste11, Ste7, and Fus3 separately. The three proteins are activated through cascading phosphorylation initiated by Ste20, and the conditions for the reactions to occur are that the kinases are activated and bound to the scaffold. Each protein has different forms of phosphorylation modifications, which may have different catalytic reaction rates; thus, they need to be listed separately.

Ordinary Differential Equations

The forms of multiple reactions are similar; here, only a portion is selected for demonstration.

Taking Ste11 as an example to illustrate the binding of the kinase with Ste5:

Variables

Table 5: Variables of Ste11 Binding Model

Parameters

Table 6: Parameters of Ste11 Binding Model

Using Ste11 catalyzing the phosphorylation of Ste7 as an example to illustrate the phosphorylation process:

Variables

Table 7: Variables of Ste7 Phosphorylation Model

Parameters

\(kcat_{Ste11pS{Ste7_{pS}}}\): Represents the catalytic efficiency in this case.

Initial Conditions

The concentrations of the three kinases are known, assuming their initial state has not undergone phosphorylation. Some enzyme activity parameters are known, and other parameters are roughly estimated to the same order of magnitude.

Fig 21 Combination of kinases and scaffold

Fig 22 Ste11 phosphorylation

Fig 23 Ste7 phosphorylation

Fig 24 Fus3 phosphorylation

Lactate Absorption Model

1. Model Description

Our project alleviates IBD symptoms by secreting lactate in the intestine to weaken autoimmunity, but it may face two aspects of doubt: first, why can't lactate or lactate bacteria probiotics be taken directly; second, will the considerable secretion of lactate cause acidosis in the human body? We hope to model our project to describe how it has a better sustained release effect compared to direct lactate consumption, more precise control compared to probiotic intake, and to avoid adaptation of the immune system and gut microbiota. Additionally, we need to develop a computational method to achieve precise control over lactate secretion to regulate treatment time and prevent acidosis.

2. Basic Assumptions

1. Only the absorption process of lactate is described, without considering other effects of lactate on the human body.
2. It is assumed that the location where lactate acts on immune cells is separated from the intestinal environment.
3. It is assumed that the secretion rate of lactate is uniform, and activated yeast cells secrete a total amount of lactate \(a\) within time \(t_0\), secreting \(\frac{a}{n}\) of lactate in the time interval \(\frac{t_0}{n}\).

3. Model Equation

According to Fick's law :

\[ \frac{dQd}{dt} = -D \frac{dC}{dx} \]

\[ \frac{dQd}{dt} = K\times Qd \]

Direct Administration

In the case of direct lactate intake, the content of lactate in the intestine can be described by the following equation:

\[ Q_d = (Q_{d_0} + a)e^{-(k_1 + k_2)t} \]

Explanation: The absorption rate is proportional to the concentration of lactic acid, and the concentration of lactate declines in an exponential form.

Parameters:

• \( Q_d \): Remaining lactate content in the intestinal environment
• \( Q_{d_0} \): Initial lactate content in the intestinal environment
• \( a \): Total amount of lactate ingested
• \( k_1 \): Absorption rate of lactate
• \( k_2 \): Rate at which lactate is eliminated due to metabolism and excretion
• \( t \): Time

Induced Secretion

The remaining lactate content in the intestinal environment has a recursive relationship over time:

\[ Q_{d_i} = \left(Q_{d_{i-1}} + \frac{a}{n}\right)e^{-(k_1 + k_2)(t - (i-1)\frac{t_0}{n})} \]

We can obtain the expression:

\[ Q_{d_i} = \frac{a}{n} \sum_{m=1}^{i-1} e^{-(k_1 + k_2)\left(mt - \left(j \frac{(m+2)(m+1)}{2} \frac{t_0}{n}\right)\right)} \]

Fig 25 Lactate Absorption

By simulating the absorption process of lactate, we can conclude that in the case of direct administration, the concentration of lactate decreases exponentially over time, while in the case of induced secretion, the concentration of lactate slowly increases over time and reaches equilibrium after a certain period.

References

[1] Borghardt JM, Weber B, Staab A, Kloft C. Pharmacometric Models for Characterizing the Pharmacokinetics of Orally Inhaled Drugs. AAPS J. 2015 Jul;17(4):853-70. doi: 10.1208/s12248-015-9760-6. Epub 2015 Apr 7. PMID: 25845315; PMCID: PMC4477002.

[2] Thomson TM, Benjamin KR, Bush A, Love T, Pincus D, Resnekov O, Yu RC, Gordon A, Colman-Lerner A, Endy D, Brent R. Scaffold number in yeast signaling system sets tradeoff between system output and dynamic range. Proc Natl Acad Sci U S A. 2011 Dec 13;108(50):20265-70. doi: 10.1073/pnas.1004042108. Epub 2011 Nov 23. PMID: 22114196; PMCID: PMC3250143.

[3] Sheng Li. Global quantitative study on the non-equilibrium cell fate decision-making of yeast under pheromone perturbation[D].Jilin University,2023.DOI:10.27162/d.cnki.gjlin.2023.007553.

[4] Scott BM, Gutiérrez-Vázquez C, Sanmarco LM, da Silva Pereira JA, Li Z, Plasencia A, Hewson P, Cox LM, O'Brien M, Chen SK, Moraes-Vieira PM, Chang BSW, Peisajovich SG, Quintana FJ. Self-tunable engineered yeast probiotics for the treatment of inflammatory bowel disease. Nat Med. 2021 Jul;27(7):1212-1222. doi: 10.1038/s41591-021-01390-x. Epub 2021 Jun 28. PMID: 34183837.

[5] Tuncbag N, Braunstein A, Pagnani A, Huang SS, Chayes J, Borgs C, Zecchina R, Fraenkel E. Simultaneous reconstruction of multiple signaling pathways via the prize-collecting steiner forest problem. J Comput Biol. 2013 Feb;20(2):124-36. doi: 10.1089/cmb.2012.0092. PMID: 23383998; PMCID: PMC3576906.

[6] Kanehisa, M. (n.d.). Pathway: Sce04011. KEGG. Retrieved October 1, 2024, from https://www.kegg.jp/pathway/sce04011