{% extends "layout.html" %} {% block title %}Model{% endblock %} {% block page_content %}
0%

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.

introduction_figure

fig 1 General Description of Model

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.
Inhalation-1

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:

  • Compartment 0 (Alveolar Space, \(A\)): \(Q_A(t)\) represents the amount of muscone in the alveoli (\(\text {mg}\)).
  • Compartment 1 (Respiratory Mucosa, \(M\)): \(Q_M(t)\) represents the amount of muscone adhered to the respiratory mucosa (\(\text {mg}\)).
  • Compartment 2 (Alveolar Capillaries, \(L\)): \(Q_L(t)\) represents the amount of muscone in the alveolar capillaries (\(\text{mg}\)).
  • Compartment 3 (Systemic Circulation, \(C\)): \(Q_C(t)\) represents the amount of muscone in the systemic circulation(\(\text{mg}\)).
  • Compartment 4 (Target Intestine, \(I\)): \(Q_I(t)\) represents the amount of muscone in the intestine(\(\text{mg}\)).
  • 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:

    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:

    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:

    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:

    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*} \]

    Inhalation-2

    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.

    % Define parameters Q_inhale = 100; % mg k_exhale = 10; k_perm = 0.005; k_adh = 0.001; k_diffMC = 0.01; k_diffLC = 0.05; k_dist = 0.001; k_excrete = 0.05; k_move = 0.02; % Define the time range tspan = [0 300]; % From 0 to 5 minutes initial_conditions = [0 0 0 0 0]; % The initial condition is 0 % solve ODE [t, y] = ode45(@(t,y) odefun(t, y, Q_inhale, k_exhale, k_perm, k_adh, k_diffMC, k_diffLC, k_dist, k_excrete, k_move), tspan, initial_conditions); % calculate V_inhale V_inhale = Q_inhale / 5 * (heaviside(t) - heaviside(t-5)); figure('Position', [100, 100, 1200, 1000]); % V_inhale(t) subplot(3,2,1) plot(t, V_inhale) title('V_{inhale}(t)') xlabel('Time (s)') ylabel('V_{inhale}') % Q_A(t) subplot(3,2,2) plot(t, y(:,1)) title('Q_A(t)') xlabel('Time (s)') ylabel('Q_A') % Q_L(t) subplot(3,2,3) plot(t, y(:,2)) title('Q_L(t)') xlabel('Time (s)') ylabel('Q_L') % Q_M(t) subplot(3,2,4) plot(t, y(:,3)) title('Q_M(t)') xlabel('Time (s)') ylabel('Q_M') % Q_C(t) subplot(3,2,5) plot(t, y(:,4)) title('Q_C(t)') xlabel('Time (s)') ylabel('Q_C') % Q_I(t) subplot(3,2,6) plot(t, y(:,5)) title('Q_I(t)') xlabel('Time (s)') ylabel('Q_I') sgtitle('Simulation Results') % ODE function dydt = odefun(t, y, Q_inhale, k_exhale, k_perm, k_adh, k_diffMC, k_diffLC, k_dist, k_excrete, k_move) V_inhale = Q_inhale / 5 * (heaviside(t) - heaviside(t-5)); dydt = zeros(5,1); dydt(1) = V_inhale - (k_exhale + k_perm) * y(1); % dQ_A/dt dydt(2) = k_perm * y(1) - k_diffLC * y(2); % dQ_L/dt dydt(3) = 0.0005 * k_adh * V_inhale - k_diffMC * y(3); % dQ_M/dt dydt(4) = k_diffMC * y(3) + k_diffLC * y(2) - k_dist * y(4) - k_excrete * y(4); % dQ_C/dt dydt(5) = k_dist * y(4) - k_move * y(5); % dQ_I/dt end

    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

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

    AF-Q8VFV4-F1-model_v4.pdb

    fig 6 Protein structure of MOR215-1

    System Preparation:

    2. Force field parameterization

    Select Force Field:

    Generate Force Field Parameters:

    perl sort_mol2_bonds.pl MUS.mol2 MUS_fix.mol2

    3. Preprocessing

    Build the system:

    Merge the system:

    Energy minimization:

    4. Molecular Dynamics Simulation

    System Equilibration:

    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
    pressure

    fig 10 Curve of the pressure over time

    density

    fig 11 Curve of the density over time

    Production Simulation:

    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
    MDS

    fig 12 Molecular dynamics simulation process

    density

    fig 13 Results of the molecular dynamics simulations

    5. Post-Processing and Analysis

    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.

    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
    # 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)
    start.pdb

    fig 14 Visualization of the protein-ligand system

    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
    Trajectory Analysis

    fig 15 Trajectory Analysis

    RMSD (Root Mean Square Deviation) Analysis

    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
    RMSD Plot

    fig 16 RMSD Analysis

    Radius of Gyration (Rg) Calculation

    gmx gyrate -s md_0_10.tpr -f md_0_10_fit.xtc -o gyrate.xvg
    #1
    xmgrace gyrate.xvg
    Additional Radius of Gyration Plot

    fig 17 Radius of Gyration Calculation

    Protein-Ligand Interaction Energy

    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
    Interaction Energy Plot

    fig 18 Protein-Ligand Interaction Energy

    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.
    MAPK Pathway

    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

    Reactions:

    \[ \begin{align*} \text{Pheromone} + \text{Ste2} & \rightarrow \text{PheromoneSte2} \\ \text{PheromoneSte2} & \rightarrow \text{Pheromone} + \text{Ste2} \\ \text{PheromoneSte2} + \text{Gpa1Ste4Ste18} & \rightarrow \text{PheromoneSte2Gpa1Ste4Ste18} \\ \text{PheromoneSte2Gpa1Ste4Ste18} & \rightarrow \text{PheromoneSte2Gpa1} + \text{Ste4Ste18} \\ \text{PheromoneSte2Gpa1} & \rightarrow \text{PheromoneSte2} + \text{Gpa1} \\ \text{Gpa1} + \text{Ste4Ste18} & \rightarrow \text{Gpa1Ste4Ste18} \end{align*} \]

    Explanation

    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.

    Ordinary Differential Equations

    \[ \begin{align*} \frac{d{P}}{dt} & = k_{off_{PS}}{PS} - k_{on_{PS}}{P}*{S} \\ \frac{d{S}}{dt} & = k_{off_{PS}}{PS} - k_{on_{PS}}{P}*{S} \\ \frac{d{PS}}{dt} & = k_{on_{PS}}{P}*{S} + k_{off_{SG}} {PSG} \\ & \quad - k_{off_{PS}}{PS} - k_{on_{SG}}{PS} * {GSS} \\ \frac{d{GSS}}{dt} & = k_{on_{GS}}{SS} * {G} - k_{on_{SG}}{PS} * {GSS} \\ \frac{d{PSGSS}}{dt} & = k_{on_{SG}}{PS} * {GSS} - k_{on_{GS}}{PSGSS} \\ \frac{d{PSG}}{dt} & = k_{on_{GS}}{PSGSS} - k_{off_{SG}} {PSG} \\ \frac{d{SS}}{dt} & = k_{on_{GS}}{PSGSS} - k_{on_{GS}}{SS} * {G} \\ \frac{d{G}}{dt} & = k_{off_{SG}} {PSG} - k_{on_{GS}}{SS} * {G} \\ \end{align*} \]

    Variables

    Table 1: Variables of Receptor Activation Model

    Variable Represents Molecule Concentration (\(\mu M\))
    \(P\) Pheromone -
    \(S\) Ste2 \(0.287\)
    \(PS\) PheromoneSte2 -
    \(GSS\) Gpa1Ste4Ste18 -
    \(PSGSS\) PheromoneSte2Gpa1Ste4Ste18 -
    \(PSG\) PheromoneSte2Gpa1 -
    \(SS\) Ste4Ste18 \(2\times 10^{-4}\)
    \(G\) Gpa1 \(2\times 10^{-4}\)
    Parameters

    Table 2: Parameters of Receptor Activation Model

    Parameter Meaning Value Unit
    \(k_{on_{PS}}\) Binding rate of Pheromone to Ste2 \(0.185\) \({\mu M}^{-1} \cdot s^{-1}\)
    \(k_{off_{PS}}\) Dissociation rate of PheromoneSte2 \(1 \times 10^{-3}\) \(s^{-1}\)
    \(k_{on_{SG}}\) Binding rate of PheromoneSte2 to Gpa1Ste4Ste18 - \({\mu M}^{-1} \cdot s^{-1}\)
    \(k_{off_{SG}}\) Dissociation rate of PheromoneSte2Gpa1 - \(s^{-1}\)
    \(k_{on_{GS}}\) Binding rate of Gpa1 to Ste4Ste18 - \({\mu M}^{-1} \cdot s^{-1}\)
    \(k_{off_{GS}}\) Dissociation rate of PheromoneGpa1Ste4Ste18 - \(s^{-1}\)
    Initial Conditions

    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.

    Receptor Activation

    fig 19 Receptor Activation

    % parameter setting kon_PS = 0.185; koff_PS = 1e-3; kon_SG = 0.1; koff_GS = 0.1; koff_SG = 0.1; kon_GS = 0.01; % initial value setting Pheromone = 0.3; Ste2 = 0.287; Pheromone_Ste2 = 0.0; Gpa1_Ste4Ste18 = 1.0; Pheromone_Ste2_Gpa1_Ste4Ste18 = 0.0; Pheromone_Ste2_Gpa1 = 0.0; Ste4Ste18 = 2e-4; Gpa1 = 2e-4; % time setting tspan = [0 500]; % Define the systems of ordinary differential equations odefun = @(t, y) [ koff_PS*y(3)-kon_PS*y(1)*y(2); koff_PS*y(3)-kon_PS*y(1)*y(2); kon_PS*y(1)*y(2)+koff_SG*y(6)-koff_PS*y(3)-kon_SG*y(3)*y(4); kon_GS*y(7)*y(8)-kon_SG*y(3)*y(4); kon_SG*y(3)*y(4)-koff_GS*y(5); koff_GS*y(5)-koff_SG*y(6); koff_GS*y(5)-kon_GS*y(7)*y(8); koff_SG*y(6)-kon_GS*y(7)*y(8); ]; % Find the numerical solution of the system of ordinary differential equations [t, Y] = ode45(odefun, tspan, [Pheromone,Ste2,Pheromone_Ste2,Gpa1_Ste4Ste18,Pheromone_Ste2_Gpa1_Ste4Ste18,Pheromone_Ste2_Gpa1,Ste4Ste18,Gpa1,]); % plot the results figure; plot(t, Y(:,1), 'DisplayName', 'Pheromone', 'LineWidth', 2); hold on; plot(t, Y(:,2), 'DisplayName', 'Ste2', 'LineWidth', 2); plot(t, Y(:,3), 'DisplayName', 'PS', 'LineWidth', 2); plot(t, Y(:,4), 'DisplayName', 'GS', 'LineWidth', 2); plot(t, Y(:,5), 'DisplayName', 'PSGS', 'LineWidth', 2); plot(t, Y(:,6), 'DisplayName', 'PSG', 'LineWidth', 2); plot(t, Y(:,7), 'DisplayName', 'Ste4Ste18', 'LineWidth', 2); plot(t, Y(:,8), 'DisplayName', 'Gpa1', 'LineWidth', 2); xlabel('Time'); ylabel('Concentration'); title('Dynamics of Variables Over Time'); legend('show'); grid on;

    Formation of the Scaffold

    Reactions:
    \[ \begin{align*} Ste5 + Ste5 & \leftrightarrows Ste5Ste5 \\ Ste4Ste18Ste5 + Ste5 & \leftrightarrows Ste4Ste18Ste5Ste5 \\ Ste4Ste18Ste5 + Ste4Ste18Ste5 & \leftrightarrows Ste4Ste18Ste5Ste5Ste4Ste18 \\ Ste4Ste18 + Ste5 & \leftrightarrows Ste4Ste18Ste5 \\ Ste4Ste18 + Ste5Ste5 & \leftrightarrows Ste4Ste18Ste5Ste5 \\ Ste4Ste18 + Ste4Ste18Ste5Ste5 & \leftrightarrows Ste4Ste18Ste5Ste5Ste4Ste18 \\ \end{align*} \]

    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.

    Ordinary Differential Equations:

    \[ \begin{align*} \frac{d{S5}}{dt} & = -2 k_{on_{S5:S5}}{S5}^2 + 2 k_{off_{S5:S5}}{S55} \\ & \quad -k_{on_{S4:S5}}{S5}*{S4} + k_{off_{S4:S5}}{S45} \\ & \quad -k_{on_{S4S5:S5}}{S5}*{S45}+k_{off_{S4S5:S5}}{S5}*{S455}\\ \frac{d{S55}}{dt} & = k_{on_{S5:S5}} {S5}^2- k_{off_{S5:S5}}{S55} \\ & \quad - k_{on_{S4:S5S5}} {S4}* {S55} + k_{off_{S4:S5S5}} {S455} \\ \frac{d{S45}}{dt} & = k_{on_{S4:S5}}{S5}*{S4}- k_{off_{S4:S5}} {S45} \\ & \quad -k_{on_{S4S5:S5}}{S5}*{S45}+k_{off_{S4S5:S5}}{S5}*{S455}\\ & \quad -2 k_{on_{S4S5:S5S4}}{S45}^2 + 2 k_{off_{S4S5:S5S4}}{S4554} \\ \frac{d{S455}}{dt} & = k_{on_{S4:S5S5}} {S4}* {S55} - k_{off_{S4:S5S5}} {S455} \\ & \quad +k_{on_{S4S5:S5}}{S5}*{S45}-k_{off_{S4S5:S5}}{S5}*{S455}\\ & \quad -k_{on_{S4:S5S5S4}}{S455}*{S4}+k_{off_{S4:S5S5S4}}{S4554}\\ \frac{d{S4554}}{dt} & = k_{on_{S4:S5S5S4}}{S455}*{S4}-k_{off_{S4:S5S5S4}}{S4554}\\ & \quad +k_{on_{S4S5:S5S4}}{S45}^2 - k_{off_{S4S5:S5S4}}{S4554} \\ \frac{d{S4}}{dt} & = -k_{on_{S4:S5}}{S5}*{S4}+ k_{off_{S4:S5}} {S45} \\ & \quad - k_{on_{S4:S5S5}} {S4}* {S55} + k_{off_{S4:S5S5}} {S455} \\ & \quad -k_{on_{S4:S5S5S4}}{S455}*{S4}+k_{off_{S4:S5S5S4}}{S4554}\\ \end{align*} \]
    Variables

    Table 3: Variables of Scaffold Formation Model

    Variable Represents Molecule
    \(S5\) Ste5
    \(S55\) Ste5Ste5
    \(S45\) Ste4Ste18Ste5
    \(S455\) Ste4Ste18Ste5Ste5
    \(S4554\) Ste4Ste18Ste5Ste5Ste4Ste18
    \(S4\) Ste4Ste18
    Parameters

    Table 4: Parameters of Scaffold Formation Model

    Parameter Meaning
    \(k_{on_{S5:S5}}\) Binding rate of Ste5 and Ste5
    \(k_{off_{S5:S5}}\) Dissociation rate of Ste5:Ste5
    \(k_{on_{S4:S5}}\) Binding rate of Ste4Ste18 and Ste5
    \(k_{off_{S4:S5}}\) Dissociation rate of Ste4Ste18:Ste5
    \(k_{on_{S4S5:S5}}\) Binding rate of Ste4Ste18Ste5 and Ste5
    \(k_{off_{S4S5:S5}}\) Dissociation rate of Ste4Ste18Ste5:Ste5
    \(k_{on_{S4:S5S5}}\) Binding rate of Ste4Ste18 and Ste5Ste5
    \(k_{off_{S4:S5S5}}\) Dissociation rate of Ste4Ste18:Ste5Ste5
    \(k_{on_{S4:S5S5S4}}\) Binding rate of Ste4Ste18Ste5Ste5 and Ste4Ste18
    \(k_{off_{S4:S5S5S4}}\) Dissociation rate of Ste4Ste18Ste5Ste5:Ste4Ste18
    \(k_{on_{S4S5:S5S4}}\) Binding rate of Ste4Ste18Ste5 and Ste4Ste18Ste5
    \(k_{off_{S4S5:S5S4}}\) Dissociation rate of Ste4Ste18Ste5:Ste5Ste4Ste18
    Initial conditions

    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.

    Scaffold Formation

    fig 20 Scaffold Formation

    % parameter setting kon_S5S5 = 0.1; koff_S5S5 = 0.01; kon_S4S5_S5 = 0.05; koff_S4S5_S5 = 0.01; kon_S4S5_S5S4 = 0.03; koff_S4S5_S5S4 = 0.01; kon_S4S5 = 0.2; koff_S4S5 = 0.1; kon_S4_S5S5 = 0.1; koff_S4_S5S5 = 0.05; kon_S4_S5S5S4 = 0.05; koff_S4_S5S5S4 = 0.02; % initial value setting Ste5 = 1.0; Ste5_Ste5 = 0.0; Ste4Ste18_Ste5 = 0.0; Ste4Ste18_Ste5_Ste5 = 0.0; Ste4Ste18_Ste5_Ste5_Ste4Ste18 = 0.0; Ste4Ste18 = 1.0; % time setting tspan = [0 100]; % Define the systems of ordinary differential equations odefun = @(t, y) [ -2*kon_S5S5*y(1)^2+2*koff_S5S5*y(2)-kon_S4S5*y(1)*y(6)+koff_S4S5*y(3)-kon_S4S5_S5*y(1)*y(3)+koff_S4S5_S5*y(4); kon_S5S5*y(1)^2-koff_S5S5*y(2)-kon_S4_S5S5*y(2)*y(6)+koff_S4_S5S5*y(4); kon_S4S5*y(1)*y(6)-koff_S4S5*y(3)-kon_S4S5_S5*y(1)*y(3)+koff_S4S5_S5*y(4)-2*kon_S4S5_S5S4*y(3)^2+2*koff_S4S5_S5S4*y(5); kon_S4_S5S5*y(2)*y(6)-koff_S4_S5S5*y(4)+kon_S4S5_S5*y(1)*y(3)-koff_S4S5_S5*y(4)-kon_S4_S5S5S4*y(4)*y(6)+koff_S4_S5S5S4*y(5); kon_S4_S5S5S4*y(4)*y(6)-koff_S4_S5S5S4*y(5)+kon_S4S5_S5S4*y(3)^2-koff_S4S5_S5S4*y(5); -kon_S4S5*y(1)*y(6)+koff_S4S5*y(3)-kon_S4_S5S5*y(2)*y(6)+koff_S4_S5S5*y(4)-kon_S4_S5S5S4*y(4)*y(6)+koff_S4_S5S5S4*y(5); ]; % Find the numerical solution of the system of ordinary differential equations [t, Y] = ode45(odefun, tspan, [Ste5, Ste5_Ste5, Ste4Ste18_Ste5, Ste4Ste18_Ste5_Ste5, Ste4Ste18_Ste5_Ste5_Ste4Ste18, Ste4Ste18]); % plot the results figure; plot(t, Y(:,1), 'DisplayName', 'Ste5', 'LineWidth', 2); hold on; plot(t, Y(:,2), 'DisplayName', 'Ste5Ste5', 'LineWidth', 2); plot(t, Y(:,3), 'DisplayName', 'Ste4Ste18Ste5', 'LineWidth', 2); plot(t, Y(:,4), 'DisplayName', 'Ste4Ste18Ste5Ste5', 'LineWidth', 2); plot(t, Y(:,5), 'DisplayName', 'Ste4Ste18Ste5Ste5Ste4Ste18', 'LineWidth', 2); plot(t, Y(:,6), 'DisplayName', 'Ste4Ste18', 'LineWidth', 2); xlabel('Time'); ylabel('Concentration'); title('Dynamics of Variables Over Time'); legend('show'); grid on

    Cascading Reactions

    Reactions:

    \[ \begin{align*} Ste5_{off_{Ste11}} + Ste11_{off} & \leftrightarrows Ste5Ste11 \\ Ste5_{off_{Ste7}} + Ste7_{off} & \leftrightarrows Ste5Ste7 \\ Ste5_{off_{Fus3}} + Fus3_{off} & \leftrightarrows Ste5Fus3 \\ \end{align*} \]

    \[ \begin{align*} Ste11 & \xrightarrow {Ste20} Ste11_{pS} \\ Ste11_{pS} & \xrightarrow {Ste20} Ste11_{pSpS} \\ Ste11_{pSpS} & \xrightarrow {Ste20} Ste11_{pSpSpT} \\ \end{align*} \]

    \[ \begin{align*} Ste7 & \xrightarrow {Ste11_{pS},Ste11_{pSpS},Ste11_{pSpSpT}} Ste7_{pS} \\ Ste7_{pS} & \xrightarrow {Ste11_{pS},Ste11_{pSpS},Ste11_{pSpSpT}} Ste7_{pSpT}\\ \end{align*} \]

    \[ \begin{align*} Fus3 & \xrightarrow {Ste7_{pS},Ste7_{pSpT}} Fus3_{pY} \\ Fus3 & \xrightarrow {Ste7_{pS},Ste7_{pSpT}} Fus3_{pT} \\ Fus3_{pY} & \xrightarrow {Ste7_{pS},Ste7_{pSpT}} Fus3_{pYpT} \\ Fus3_{pT} & \xrightarrow {Ste7_{pS},Ste7_{pSpT}} Fus3_{pYpT} \\ \end{align*} \]

    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:

    \[ \begin{align*} \frac{dSte5_{off_{Ste11}}}{dt} & = k_{off_{Ste5Ste11}}Ste5Ste11 - k_{on_{Ste5Ste11}}Ste5_{off_{Ste11}} * Ste11_{off} \\ \frac{dSte11_{off}}{dt} & = k_{off_{Ste5Ste11}}Ste5Ste11 - k_{on_{Ste5Ste11}}Ste5_{off_{Ste11}} * Ste11_{off} \\ \frac{dSte5Ste11}{dt} & = - k_{off_{Ste5Ste11}}Ste5Ste11 + k_{on_{Ste5Ste11}}Ste5_{off_{Ste11}} * Ste11_{off} \\ \end{align*} \]

    Variables

    Table 5: Variables of Ste11 Binding Model

    Variable Represents Molecule
    \(Ste5_{off_{Ste11}}\) Unbound kinase Ste5
    \(Ste11_{off}\) Unbound scaffold Ste11
    \(Ste5Ste11\) Bound Ste5 and Ste11

    Parameters

    Table 6: Parameters of Ste11 Binding Model

    Parameter Meaning Units
    \(k_{off_{Ste5Ste11}}\) Dissociation rate of Ste5Ste11 \({s}^{-1}\)
    \(k_{on_{Ste5Ste11}}\) Association rate of Ste5 and Ste11 \({\mu M}^{-1}·s^{-1}\)

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

    \[ \frac{dSte7_{pS}}{dt} = kcat_{Ste11pS{Ste7_{pS}}}Ste11_{pS}*\frac{Ste5Ste11}{Ste11_{total}}*\frac{Ste5Ste7}{Ste7_{total}}*\frac{Ste7_{pS}}{Ste7_{total}}+\ldots \]

    Variables

    Table 7: Variables of Ste7 Phosphorylation Model

    Variable Represents Molecule
    \(Ste7_{pS}\) Phosphorylated Ste7 at S359
    \(Ste11_{pS}\) Phosphorylated Ste11 at S302
    \(Ste5Ste11\) Ste11 bound to Ste5
    \(Ste5Ste7\) Ste7 bound to Ste5
    \(Ste7_{total}\) Total amount of Ste7

    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.

    Combination of kinases and scaffold

    fig 21 Combination of kinases and scaffold

    Ste11 phosphorylation

    fig 22 Ste11 phosphorylation

    Ste11 phosphorylation

    fig 23 Ste7 phosphorylation

    Fus3 phosphorylation

    fig 24 Fus3 phosphorylation

    % initial value setting scaffold_Ste5 = 0.5; Ste5_off_Ste11 = scaffold_Ste5; Ste11_off = 0.2; Ste5_Ste11 = 0.0; Ste5_off_Ste7 = scaffold_Ste5; Ste7_off = 0.05; Ste5_Ste7 = 0.0; Ste5_off_Fus3 = scaffold_Ste5; Fus3_off = 0.088; Ste5_Fus3 = 0.0; Ste11 = Ste11_off; Ste11_pS = 0; Ste11_pSpS = 0; Ste11_pSpSpT = 0; Ste7 = Ste7_off; Ste7_pS = 0; Ste7_pSpT = 0; Fus3 = Fus3_off; Fus3_pY = 0; Fus3_pT = 0; Fus3_pYpT = 0; % parameter setting kon_Ste5_Ste11 = 0.1; % Ste5_off_Ste11 + Ste11_off -> Ste5_Ste11 koff_Ste5_Ste11 = 0.2; % Ste5_Ste11 -> Ste5_off_Ste11 + Ste11_off kon_Ste5_Ste7 = 0.3; % Ste5_off_Ste7 + Ste7_off -> Ste5_Ste7 koff_Ste5_Ste7 = 0.4; % Ste5_Ste7 -> Ste5_off_Ste7 + Ste7_off kon_Ste5_Fus3 = 0.5; % Ste5_off_Fus3 + Fus3_off -> Ste5_Fus3 koff_Ste5_Fus3 = 0.6; % Ste5_Fus3 -> Ste5_off_Fus3 + Fus3_off kcat_Ste20Ste11_pS = 1.8438; % Ste11 -> Ste11_pS kcat_Ste20Ste11pS_pS = 1.8438; % Ste11_pS -> Ste11_pSpS kcat_Ste20Ste11pSpS_pT = 1.8438; % Ste11_pSpS -> Ste11_pSpSpT kcat_Ste11pSSte7_pS = 1; % Ste11_pS + Ste7 -> Ste7_pS kcat_Ste11pSSte7pS_pT = 0.5; % Ste11_pS + Ste7_pS -> Ste7_pSpT kcat_Ste11pSpSSte7_pS = 2; % Ste11_pSpS + Ste7 -> Ste7_pS kcat_Ste11pSpSSte7pS_pT = 1.5; % Ste11_pSpS + Ste7_pS -> Ste7_pSpT kcat_Ste11pSpSpTSte7_pS = 3; % Ste11_pSpSpT + Ste7 -> Ste7_pS kcat_Ste11pSpSpTSte7pS_pT = 2.5; % Ste11_pSpSpT + Ste7_pS -> Ste7_pSpT kcat_Ste5Ste7pSFus3_pY = 2.6079; % Ste7_pS + Fus3 -> Fus3_pY kcat_Ste5Ste7pSFus3_pT = 2.6079; % Ste7_pS + Fus3 -> Fus3_pT kcat_Ste5Ste7pSFus3pY_pT = 2.6079; % Ste7_pS + Fus3_pY -> Fus3_pYpT kcat_Ste5Ste7pSFus3pT_pY = 2.6079; % Ste7_pS + Fus3_pT -> Fus3_pYpT kcat_Ste5Ste7pSpTFus3_pY = 0.86812; % Ste7_pSpT + Fus3 -> Fus3_pY kcat_Ste5Ste7pSpTFus3_pT = 6.7879; % Ste7_pSpT + Fus3 -> Fus3_pT kcat_Ste5Ste7pSpTFus3pY_pT = 3.5252; % Ste7_pSpT + Fus3_pY -> Fus3_pYpT kcat_Ste5Ste7pSpTFus3pT_pY = 1.5; % Ste7_pSpT + Fus3_pT -> Fus3_pYpT kcat_dephosph = 0.1; % time setting tspan = [0 40]; % Define the systems of ordinary differential equations odefun = @(t, y) [ koff_Ste5_Ste11*y(3)-kon_Ste5_Ste11*y(1)*y(2); % Ste5_off_Ste11 y(1) koff_Ste5_Ste11*y(3)-kon_Ste5_Ste11*y(1)*y(2); % Ste11_off y(2) -koff_Ste5_Ste11*y(3)+kon_Ste5_Ste11*y(1)*y(2); % Ste5_Ste11 y(3) koff_Ste5_Ste7*y(6)-kon_Ste5_Ste7*y(4)*y(5); % Ste5_off_Ste7 y(4) koff_Ste5_Ste7*y(6)-kon_Ste5_Ste7*y(4)*y(5); % Ste7_off y(5) -koff_Ste5_Ste7*y(6)+kon_Ste5_Ste7*y(4)*y(5); % Ste5_Ste7 y(6) koff_Ste5_Fus3*y(9)-kon_Ste5_Fus3*y(7)*y(8); % Ste5_off_Fus3 y(7) koff_Ste5_Fus3*y(9)-kon_Ste5_Fus3*y(7)*y(8); % Fus3_off y(8) -koff_Ste5_Fus3*y(9)+kon_Ste5_Fus3*y(7)*y(8); % Ste5_Fus3 y(9) -kcat_Ste20Ste11_pS*scaffold_Ste5*(y(3)*y(10)/(y(2)+y(3))^2)+kcat_dephosph*y(11); % Ste11 y(10) kcat_Ste20Ste11_pS*scaffold_Ste5*(y(3)*y(10)/(y(2)+y(3))^2)-kcat_Ste20Ste11pS_pS*scaffold_Ste5*(y(3)*y(11)/(y(2)+y(3))^2)-kcat_dephosph*y(11)+kcat_dephosph*y(12); % Ste11_pS y(11) kcat_Ste20Ste11pS_pS*scaffold_Ste5*(y(3)*y(11)/(y(2)+y(3))^2)-kcat_Ste20Ste11pSpS_pT*scaffold_Ste5*(y(3)*y(12)/(y(2)+y(3))^2)-kcat_dephosph*y(12)+kcat_dephosph*y(13); % Ste11_pSpS y(12) kcat_Ste20Ste11pSpS_pT*scaffold_Ste5*(y(3)*y(12)/(y(2)+y(3))^2)-kcat_dephosph*y(13); % Ste11_pSpSpT y(13) -kcat_Ste11pSSte7_pS*y(11)*(y(3)/(y(2)+y(3)))*(y(6)*y(14)/(y(5)+y(6))^2)-kcat_Ste11pSpSSte7_pS*y(12)*(y(3)/(y(2)+y(3)))*(y(6)*y(14)/(y(5)+y(6))^2)-kcat_Ste11pSpSpTSte7_pS*y(13)*(y(3)/(y(2)+y(3)))*(y(6)*y(14)/(y(5)+y(6))^2)+kcat_dephosph*y(15); % Ste7 y(14) kcat_Ste11pSSte7_pS*y(11)*(y(3)/(y(2)+y(3)))*(y(6)*y(14)/(y(5)+y(6))^2)-kcat_Ste11pSSte7pS_pT*y(11)*(y(3)/(y(2)+y(3)))*(y(6)*y(15)/(y(5)+y(6))^2)+kcat_Ste11pSpSSte7_pS*y(12)*(y(3)/(y(2)+y(3)))*(y(6)*y(14)/(y(5)+y(6))^2)-kcat_Ste11pSpSSte7pS_pT*y(12)*(y(3)/(y(2)+y(3)))*(y(6)*y(15)/(y(5)+y(6))^2)+kcat_Ste11pSpSpTSte7_pS*y(13)*(y(3)/(y(2)+y(3)))*(y(6)*y(14)/(y(5)+y(6))^2)-kcat_Ste11pSpSpTSte7pS_pT*y(13)*(y(3)/(y(2)+y(3)))*(y(6)*y(15)/(y(5)+y(6))^2)-kcat_dephosph*y(15)+kcat_dephosph*y(16); % Ste7_pS y(15) kcat_Ste11pSSte7pS_pT*y(11)*(y(3)/(y(2)+y(3)))*(y(6)*y(15)/(y(5)+y(6))^2)+kcat_Ste11pSpSSte7pS_pT*y(12)*(y(3)/(y(2)+y(3)))*(y(6)*y(15)/(y(5)+y(6))^2)+kcat_Ste11pSpSpTSte7pS_pT*y(13)*(y(3)/(y(2)+y(3)))*(y(6)*y(15)/(y(5)+y(6))^2)-kcat_dephosph*y(16); % Ste7_pSpT y(16) -kcat_Ste5Ste7pSFus3_pY*y(15)*(y(6)/(y(5)+y(6)))*(y(9)*y(17)/(y(8)+y(9))^2)-kcat_Ste5Ste7pSFus3_pT*y(15)*(y(6)/(y(5)+y(6)))*(y(9)*y(17)/(y(8)+y(9))^2)-kcat_Ste5Ste7pSpTFus3_pY*y(16)*(y(6)/(y(5)+y(6)))*(y(9)*y(17)/(y(8)+y(9))^2)-kcat_Ste5Ste7pSpTFus3_pT*y(16)*(y(6)/(y(5)+y(6)))*(y(9)*y(17)/(y(8)+y(9))^2)+kcat_dephosph*y(18)+kcat_dephosph*y(19); % Fus3 y(17) kcat_Ste5Ste7pSFus3_pY*y(15)*(y(6)/(y(5)+y(6)))*(y(9)*y(17)/(y(8)+y(9))^2)-kcat_Ste5Ste7pSFus3pY_pT*y(15)*(y(6)/(y(5)+y(6)))*(y(9)*y(18)/(y(8)+y(9))^2)+kcat_Ste5Ste7pSpTFus3_pY*y(16)*(y(6)/(y(5)+y(6)))*(y(9)*y(17)/(y(8)+y(9))^2)-kcat_Ste5Ste7pSpTFus3pY_pT*y(16)*(y(6)/(y(5)+y(6)))*(y(9)*y(18)/(y(8)+y(9))^2)-kcat_dephosph*y(18)+kcat_dephosph*y(20); % Fus3_pY y(18) kcat_Ste5Ste7pSFus3_pT*y(15)*(y(6)/(y(5)+y(6)))*(y(9)*y(17)/(y(8)+y(9))^2)-kcat_Ste5Ste7pSFus3pT_pY*y(15)*(y(6)/(y(5)+y(6)))*(y(9)*y(19)/(y(8)+y(9))^2)+kcat_Ste5Ste7pSpTFus3_pT*y(16)*(y(6)/(y(5)+y(6)))*(y(9)*y(17)/(y(8)+y(9))^2)-kcat_Ste5Ste7pSpTFus3pT_pY*y(16)*(y(6)/(y(5)+y(6)))*(y(9)*y(19)/(y(8)+y(9))^2)-kcat_dephosph*y(19)+kcat_dephosph*y(20); % Fus3_pT y(19) kcat_Ste5Ste7pSFus3pY_pT*y(15)*(y(6)/(y(5)+y(6)))*(y(9)*y(18)/(y(8)+y(9))^2)+kcat_Ste5Ste7pSFus3pT_pY*y(15)*(y(6)/(y(5)+y(6)))*(y(9)*y(19)/(y(8)+y(9))^2)+kcat_Ste5Ste7pSpTFus3pY_pT*y(16)*(y(6)/(y(5)+y(6)))*(y(9)*y(18)/(y(8)+y(9))^2)+kcat_Ste5Ste7pSpTFus3pT_pY*y(16)*(y(6)/(y(5)+y(6)))*(y(9)*y(19)/(y(8)+y(9))^2)-2*kcat_dephosph*y(20); % Fus3_pYpT y(20) ]; % Find the numerical solution of the system of ordinary differential equations [t, Y] = ode45(odefun, tspan, [Ste5_off_Ste11; Ste11_off; Ste5_Ste11; Ste5_off_Ste7; Ste7_off; Ste5_Ste7; Ste5_off_Fus3; Fus3_off; Ste5_Fus3; Ste11; Ste11_pS; Ste11_pSpS; Ste11_pSpSpT; Ste7; Ste7_pS; Ste7_pSpT; Fus3; Fus3_pY; Fus3_pT; Fus3_pYpT]); figure; hold on; plot(t, Y(:, 1), 'LineWidth', 2, 'DisplayName', 'Ste5_{off_{Ste11}}'); plot(t, Y(:, 2), 'LineWidth', 2, 'DisplayName', 'Ste11_{off}'); plot(t, Y(:, 3), 'LineWidth', 2, 'DisplayName', 'Ste5Ste11'); plot(t, Y(:, 4), 'LineWidth', 2, 'DisplayName', 'Ste5_{off_{Ste7}}'); plot(t, Y(:, 5), 'LineWidth', 2, 'DisplayName', 'Ste7_{off}'); plot(t, Y(:, 6), 'LineWidth', 2, 'DisplayName', 'Ste5Ste7'); plot(t, Y(:, 7), 'LineWidth', 2, 'DisplayName', 'Ste5_{off_{Fus3}}'); plot(t, Y(:, 8), 'LineWidth', 2, 'DisplayName', 'Fus3_{off}'); plot(t, Y(:, 9), 'LineWidth', 2, 'DisplayName', 'Ste5Fus3'); xlabel('Time (t)'); ylabel('Concentration'); title('Variable Change Graph'); legend show; grid on; hold off; figure; hold on; plot(t, Y(:, 10), 'LineWidth', 2, 'DisplayName', 'Ste11'); plot(t, Y(:, 11), 'LineWidth', 2, 'DisplayName', 'Ste11_{pS}'); plot(t, Y(:, 12), 'LineWidth', 2, 'DisplayName', 'Ste11_{pSpS}'); plot(t, Y(:, 13), 'LineWidth', 2, 'DisplayName', 'Ste11_{pSpSpT}'); xlabel('Time (t)'); ylabel('Concentration'); title('Ste11 Variable Change Graph'); legend show; grid on; hold off; figure; hold on; plot(t, Y(:, 14), 'LineWidth', 2, 'DisplayName', 'Ste7'); plot(t, Y(:, 15), 'LineWidth', 2, 'DisplayName', 'Ste7_{pS}'); plot(t, Y(:, 16), 'LineWidth', 2, 'DisplayName', 'Ste7_{pSpT}'); xlabel('Time (t)'); ylabel('Concentration'); title('Ste7 Variable Change Graph'); legend show; grid on; hold off; figure; hold on; plot(t, Y(:, 17), 'LineWidth', 2, 'DisplayName', 'Fus3'); plot(t, Y(:, 18), 'LineWidth', 2, 'DisplayName', 'Fus3_{pY}'); plot(t, Y(:, 19), 'LineWidth', 2, 'DisplayName', 'Fus3_{pT}'); plot(t, Y(:, 20), 'LineWidth', 2, 'DisplayName', 'Fus3_{pYpT}'); xlabel('Time (t)'); ylabel('Concentration'); title('Fus3 Variable Change Graph'); legend show; grid on; hold off;

    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} \]

    Because the distance between diffusion is very small, the concentration difference between the two sides of the system replaces the concentration gradient, so this formula can be simplified to:

    \[ \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:

    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)} \]

    Lactate Absorption

    fig 25 Lactate Absorption

    % Parameter settings Qd0 = 0; % Initial lactic acid level a = 50; % Total amount of lactic acid administered k1 = 0.1; % Absorption rate of lactic acid k2 = 0.05; % Metabolism and excretion rate of lactic acid t0 = 20; % Total time n = 10; % Number of induced secretions time_step = 0.1; % Time step % Time array t = 0:time_step:t0; % Direct administration simulation Qd_direct = (Qd0 + a) * exp(-(k1 + k2) * t); % Induced secretion simulation Qd_induced = zeros(1, length(t)); for i = 1:n % Calculate the time point for each induced secretion current_time = (i-1) * (t0 / n); if current_time <= t0 % Update concentration for j = 1:length(t) if t(j) >= current_time Qd_induced(j) = Qd_induced(j) + a/n * exp(-(k1 + k2) * (t(j) - current_time)); end end end end % Plot results figure; hold on; plot(t, Qd_direct, 'r-', 'LineWidth', 2, 'DisplayName', 'Direct Administration'); plot(t, Qd_induced, 'b-', 'LineWidth', 2, 'DisplayName', 'Induced Secretion'); xlabel('Time (t)'); ylabel('Intestinal Lactate Level (Qd)'); title('Lactate Concentration Change Curve'); legend; grid on; hold off;

    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


    {% endblock %}