The Modeling of Reactive Oxygen Species (ROS) Promoters

---

1.description

---

In our project, we conducted an in-depth analysis through modeling to study the concentration of reactive oxygen species (ROS) during enteritis and its impact on the expression levels of specific promoters. In this model, we simulated the dynamic changes in ROS concentration during enteritis based on both exogenous and endogenous factors. Building upon this model, we then considered three key influencing factors—ROS concentration, temperature, and pH—to simulate the expression levels of the three promoters we selected during the course of enteritis. Finally, we plotted the promoter activity throughout the entire enteritis process.

Specifically, we modeled the process of immune cell activation and the suppression of antioxidant enzyme expression by pro-inflammatory cytokines, using modified exponential and logistic functions to simulate the changes in ROS concentration. Additionally, we modeled the fluctuations in temperature and pH to better capture the influence of these external conditions on promoter expression. Ultimately, the model predicted response curves for the three promoters in relation to changes in ROS, temperature, and pH during the course of enteritis, providing a strong basis for further experimental validation.

---

2.Hypothesis

---

(1)During the process of inflammatory bowel disease (IBD), ROS concentration is influenced by two factors: endogenous and exogenous factors. For endogenous factors, we considered two categories: the activation of intestinal immune cells, which produce large amounts of ROS, and the suppression of antioxidant enzyme expression by pro-inflammatory cytokines such as tumor necrosis factor-α (TNF-α), interleukin-6 (IL-6), and interleukin-1β (IL-1β), leading to an increase in ROS concentration. As for exogenous factors, we considered the impact of temperature and pH on ROS concentration.

(2)The activation of immune cells, such as neutrophils and macrophages, directly increases ROS production, and this process can be represented by a linear function\( P_{\mathrm{immune}}(t) \).

(3)Pro-inflammatory cytokines such as TNF-α, IL-6, and IL-1β increase ROS concentration by inhibiting the expression of antioxidant enzymes, thereby reducing ROS clearance. This process can also be represented by the function \( P_{\mathrm{cytokine}}(t) \).

---

3.Parameter of model1

---

Parameter	Value	Description	Reference
\( \alpha_{\mathrm{immune}} \)	5	The rate of ROS production induced by immune cell activation	[1]
\( \beta_{\mathrm{immune}} \)	1	The rate of increase in ROS production during immune cell activation	assumption
\( P_{\mathrm{max}} \)	20	The maximum rate of ROS production during immune cell activation	assumption
\( \alpha_{\mathrm{cytokine}} \)	2	The proportional constant of the rate of ROS production induced by pro-inflammatory cytokines	[2]
\( E_{a} \)	5000	The activation energy for the effect of temperature on ROS production	[3]
\( T_{0} \)	310	Initial temperature	assumption
\( R_{0} \)	0	Initial ROS concentration	assumption
\( \Delta T_{\max} \)	2	The maximum amplitude of temperature variation	assumption
\( \gamma \)	0.1	The regulatory constant for the effect of pH on ROS production	[4]
\( \mathrm{pH}_0 \)	7.4	Initial pH value	assumption
\( \Delta\mathrm{pH}_{\max} \)	0.5	The maximum amplitude of pH variation	assumption
\( \omega_{\mathrm{pH}} \)	0.2	The frequency of pH fluctuations	assumption
\( \kappa_{\mathrm{pH}} \)	0.1	The decay rate of pH variation	assumption
\( K_{2} \)	0.1	The rate constant of natural ROS clearance	assumption
\( \kappa_{T} \)	0.3	The rate constant of temperature variation	assumption

---

4.Abbreviations of the model1

---

---

5.Equation

---

By integrating the effects of both endogenous and exogenous factors on ROS concentration, the rate of change in ROS concentration can be described by the following ordinary differential equation:

\[ \frac{dR(t)}{dt}=(P_{\mathrm{immune}}(t)+P_{\mathrm{cytokine}}(t))\cdot f_T(T(t))\cdot f_{\mathrm{pH}}(\mathrm{pH}(t))-k_2\cdot R(t) \]

Modeling of endogenous influencing factors:

Among them,\( P_{\mathrm{immune}}(t) \) represents the portion of the increase in ROS concentration due to immune cell activation. This typically rises with the inflammatory response, as activated immune cells (such as neutrophils and macrophages) produce large amounts of ROS to combat infections and damaged tissues. However, high ROS concentrations can inhibit the activity of immune cells, preventing an unlimited increase. In conclusion, the rate of ROS production due to immune cell activation initially shows exponential growth before stabilizing.

We used modified exponential and logistic functions to simulate this process:

\[ P_{\mathrm{immune}}(t)=\frac{\alpha_{\mathrm{immune}}\cdot e^{\beta_{\mathrm{imunune}}t}}{1+\frac{\alpha_{\mathrm{immune}}}{P_{\mathrm{max}}}\cdot(e^{\beta_{\mathrm{immune}}t}-1)} \]

The modified exponential function

\[ P_{\mathrm{immune}}(t)=\frac{P_{\max}}{1+e^{-\gamma(t-t_0)}} \]

Logistic function

And we plotted the two function curves respectively:

We found that the modified exponential function better fits the actual situation: the growth rate is faster, and it reaches a stable state in a shorter period of time.

Next, assuming that the concentration of pro-inflammatory cytokines changes over time, the inhibitory effect on antioxidant enzymes is represented by a linear relationship:

\[ P_{\mathrm{cytokine}}(t)=\alpha_{\mathrm{cytokine}}\cdot C(t) \]

Modeling of exogenous influencing factors:

1.The function representing the effect of temperature\( f_T\left(T\left(t\right)\right) \)

The effect of temperature on ROS production can still be represented using the Arrhenius equation:

\[ f_T(T(t))=e^{-\frac{E_a}{RT(t)}} \]

Where \( \text{Ea} \) is the activation energy, and R is the gas constant.

2.The function representing the effect of pH on ROS production\( f_{\mathrm{pH}}(\mathrm{pH}(t)) \)

Assuming the effect of pH on ROS concentration can be represented by a quadratic function, it simulates the maximum impact within an optimal pH range:

\[ f_{\mathrm{pH}}(\mathrm{pH}(t))=1-\gamma\cdot(\mathrm{pH}(t)-\mathrm{pH}_{opt})^2 \]

Where \( \gamma \) is the adjustment coefficient, and \( \mathrm{pH}_{opt} \) represents the optimal pH value.

The final model equation:

\[ \frac{dR(t)}{dt}=\left(\frac{\alpha_{\mathrm{immma}}\cdot e^{\beta_{\mathrm{mmma}}t}}{1+\frac{\beta_{\mathrm{mmma}}}{P_{\mathrm{max}}}\cdot(e^{\beta_{\mathrm{mmmax}}t}-1)}+\alpha_{\mathrm{cytoline}}\cdot C(t)\right)\cdot e^{-\frac{E_a}{RT(t)}}\cdot\left(1-\gamma\cdot\left(\mathrm{pH}(t)-\mathrm{pH}_{\mathrm{opt}}\right)^2\right)-k_2\cdot R(t) \]

Equation(2): Simulation of the expression of promoters dps, katG, and gorA during the course of enteritis:

During the course of enteritis, the main factors referenced are the dynamic changes in ROS concentration, temperature, and pH. The ROS concentration variation has already been modeled in Equation (1). Next, we only need to model temperature and pH.

Temperature Modeling:

Temperature during the IBD process may be affected by fever or localized inflammation. Let's assume the temperature variation over time follows the model below:

1. The initial temperature \( T_0 \) is the normal human body temperature (approximately 310K).

2. As inflammation progresses, the temperature will rise and gradually stabilize.

The temperature can be represented by a model of asymptotic variation:

\[ T(t)=T_0+\Delta T_{\max}\cdot\left(1-e^{-\kappa_Tt}\right) \]

pH Modeling: The pH of the intestine fluctuates during the course of IBD.

\[ \mathrm{pH}(t)=\mathrm{pH}_0+\Delta\mathrm{pH}_{\max}\cdot\sin(\omega_\text{pH}t)\cdot e^{-\kappa_\text{pH}t} \]

Promoter Activity and Selection Modeling:

Assume that the sensitivity of each promoter to ROS concentration, temperature, and pH can be represented by a response function.

The response function for promoter iii is defined as follows: \( A_i(R(t),T(t),\mathrm{pH}(t)) \) ,Its form is: \( A_{i}(t)=\alpha_{i}\cdot f_{R,i}(R(t))\cdot f_{T,i}(T(t))\cdot f_{\mathrm{pH},i}(\mathrm{pH}(t)) \)

The specific form of the response function:

The response of each promoter to different parameters can be modeled using a logistic function.

\[ f_{R,i}(R(t))=\exp\left(-\frac{(R(t)-R_{\mathrm{opt},i})^2}{2\sigma_{R,i}^2}\right) \]

\[ f_{T,i}(T(t))=\exp\left(-\frac{(T(t)-T_{\mathrm{opt},i})^2}{2\sigma_{T,i}^2}\right) \]

\[ f_{\mathrm{pH},i}(\mathrm{pH}(t))=\exp\left(-\frac{(\mathrm{pH}(t)-\mathrm{pH}_{\mathrm{opt},i})^2}{2\sigma_{\mathrm{pH},i}^2}\right) \]

Where \( R_{\mathrm{opt},i_,} T_{\mathrm{opt},i_,} \mathrm{pH}_{\mathrm{opt},i} \) are the optimal activation conditions for promoter i , and \( \sigma_{R,i},\sigma_{T,i},\sigma_{\mathrm{pH},i} \) are the response width parameters.

---

6.Results and Analysis

---

Finally, we obtained the response patterns of the three promoters during the process of enteritis, and the results are shown in the figure below:

In this figure, we can observe:

In summary, we found that katG exhibited the most significant expression during the course of enteritis, reaching its peak in the mid-stage of inflammation, indicating a strong oxidative stress response. In contrast, the expression of the dps and gorA promoters was more moderate, with gorA showing particularly weak response, suggesting its role in the stress response during enteritis may be limited. Based on the significant expression pattern of katG during enteritis and its key regulatory role in oxidative stress response, we ultimately chose katG as the experimental promoter to better study and regulate stress response mechanisms in enteritis.

---

Screening of ROS-Sensitive Promoters

---

1.description

---

Escherichia coli possesses its own unique set of operon systems that respond to high oxidative stress, among which OxyR has garnered our attention as the ROS sensor and transcription activator in E. coli. In the presence of H2O2 or ROS, two cysteine residues (Cys199 and Cys208) of the OxyR protein are oxidized to form a disulfide bond, leading to a conformational change in the protein. This results in binding to promoter sequences of some ROS-responsive genes, subsequently activating the expression of downstream antioxidant genes [5]. The OxyR transcription factor belongs to the LysR-type transcriptional regulators (LTTRs) and is a homologous tetramer, with each subunit comprising an N-terminal DNA-binding domain [6]. Here, we identify the optimal promoter sequence that ensures the highest binding efficiency and specificity with the OxyR transcription factor through transcription factor binding site prediction, molecular docking, and molecular dynamics simulation.

---

2. Motif-Based Prediction of OxyR Binding Sites

---

Since we use Escherichia coli as the engineering strain, we first retrieved the target genes and binding motifs of OxyR from the RegulonDB database, which is a transcriptional regulatory network database specific to E. coli (http://regulondb.ccg.unam.mx). As shown in the following diagram obtained from the database, OxyR potentially activates the expression of 32 target genes, represses the expression of 12 genes, has bidirectional interactions with 4 genes, and is associated with 3 operon-related genes. The three candidate promoters of interest in the wet lab also originate from the genes activated by OxyR.

Some studies have shown that OxyR can recognize nucleotide motifs composed of ATAGnt elements arranged at 10 bp intervals, but the specific motif remains unknown [7]. We retrieved the upstream sequences (including promoter sequences) of the aforementioned target genes' transcription start sites (TSS) from the National Center for Biotechnology Information (NCBI) and stored them in .fasta file format. Subsequently, all .fasta files were merged into a single large .fasta file. The MEME Suite software allows biologists to discover new motifs in unaligned sets of nucleotide or protein sequences and perform various other motif-based analyses. This suite provides motif discovery algorithms using probabilistic models (MEME) and discrete models (STREME) [8]. After deploying meme (v.5.7.7) on the server and setting the motif length parameter to 10bp~20bp, we performed the prediction and obtained five motifs as shown in the figure below. These motifs are arranged based on their E-value from low to high, with higher E-values indicating motifs that are less likely to have biological significance. Lower scores suggest that the motif was found in more sequences and matches well. The motifs identified from the total of 21 promoter sequences of OxyR target genes are considered potential binding sites where OxyR facilitates transcription activation. We conducted further analysis on these identified sites.

For the selected motif, we used FIMO to match the motif's position and matching degree on the promoter sequences of all OxyR target genes, scoring and calculating the p-value. The method involves calculating the score for a sequence's position with the motif by summing the appropriate entries in each column of the position-specific scoring matrix representing the motif. The p-value for motif occurrence is defined as the probability that a random sequence of the same length as the motif would match the sequence at that position with an equal or better score [9]. As shown in the figure below, we can see that this motif scores the highest in the promoter sequence of KatG.

---

3.Molecular Docking

---

After obtaining the motif of the targeted sequence for OxyR, we were able to identify the binding site of OxyR in the promoter sequence, which facilitates the selection of docking models. The protein motif of its DNA binding domain obtained from RegulonDB is of the H-T-H type, with the sequence as follows: FRRAADSCHVSQPTLSGQIR. Due to the absence of a complete structural file for OxyR in the PDB database, we turned to the AlphaFold Protein Structure Database to download the high-confidence predicted structure of OxyR (AFDB accession: AF-P0ACQ4-F1-v4). Subsequently, we employed the molecular docking software HDOCK developed by Huang Lab to assess the binding affinity between DNA and protein [10].We first input six candidate DNA sequence files along with the PDB file of the OxyR protein. In setting the binding sites, we entered the motif of OxyR's DNA binding domain. We then submitted the task for molecular docking on the online platform (http://hdock.phys.hust.edu.cn/). Each result provided 100 or more docking models to choose from. Our criterion for model selection was that the binding site of OxyR in the model must correspond to the location of the motif in the promoter. Based on this, the results of docking between the six promoter sequences and the OxyR protein were output as PDB files.We utilized Python (v.3.12.6) and the open-source version of PyMOL (v.3.0.0) for visualization of the PDB files, as shown in the figure below. Each docking model is associated with a docking score and confidence level. A lower docking score indicates tighter binding, while a higher confidence level denotes greater reliability of the docking model. We used bar charts to present the quality of the docking results. From the figure, it is evident that the binding score between the KatG gene's promoter sequence and OxyR is the lowest, with the highest confidence level, demonstrating that the binding affinity is the best. In summary, selecting the promoter sequence of the KatG gene as an element responding to ROS is the optimal choice.

---

Reference

[1]Morris G, Gevezova M, Sarafian V, Maes M. Redox regulation of the immune response. Cell Mol Immunol. 2022 Oct;19(10):1079-1101. doi: 10.1038/s41423-022-00902-0. Epub 2022 Sep 2. PMID: 36056148; PMCID: PMC9508259.

[2]Bhatt S, Nagappa AN, Patil CR. Role of oxidative stress in depression. Drug Discov Today. 2020 Jul;25(7):1270-1276. doi: 10.1016/j.drudis.2020.05.001. Epub 2020 May 8. PMID: 32404275.

[3]Vergnolle N. TRPV4: new therapeutic target for inflammatory bowel diseases. Biochem Pharmacol. 2014 May 15;89(2):157-61. doi: 10.1016/j.bcp.2014.01.005. Epub 2014 Jan 16. PMID: 24440740.

[4]Lee S, Shanti A. Effect of Exogenous pH on Cell Growth of Breast Cancer Cells. Int J Mol Sci. 2021 Sep 14;22(18):9910. doi: 10.3390/ijms22189910. PMID: 34576073; PMCID: PMC8464873.

[5]Pomposiello PJ, Demple B. Redox-operated genetic switches: the SoxR and OxyR transcription factors. Trends Biotechnol. 2001 Mar;19(3):109-14. doi: 10.1016/s0167-7799(00)01542-0. PMID: 11179804.

[6]Wei Q, Minh PN, Dötsch A, Hildebrand F, Panmanee W, Elfarash A, Schulz S, Plaisance S, Charlier D, Hassett D, Häussler S, Cornelis P. Global regulation of gene expression by OxyR in an important human opportunistic pathogen. Nucleic Acids Res. 2012 May;40(10):4320-33. doi: 10.1093/nar/gks017. Epub 2012 Jan 24. PMID: 22275523; PMCID: PMC3378865.

[7]Toledano MB, Kullik I, Trinh F, Baird PT, Schneider TD, Storz G. Redox-dependent shift of OxyR-DNA contacts along an extended DNA-binding site: a mechanism for differential promoter selection. Cell. 1994 Sep 9;78(5):897-909. doi: 10.1016/s0092-8674(94)90702-1. PMID: 8087856.

[8]Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, Ren J, Li WW, Noble WS. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res. 2009 Jul;37(Web Server issue):W202-8. doi: 10.1093/nar/gkp335. Epub 2009 May 20. PMID: 19458158; PMCID: PMC2703892.

[9]Grant CE, Bailey TL, Noble WS. FIMO: scanning for occurrences of a given motif. Bioinformatics. 2011 Apr 1;27(7):1017-8. doi: 10.1093/bioinformatics/btr064. Epub 2011 Feb 16. PMID: 21330290; PMCID: PMC3065696.

[10]Yan Y, Tao H, He J, Huang S-Y.* The HDOCK server for integrated protein-protein docking. Nature Protocols, 2020; doi: https://doi.org/10.1038/s41596-020-0312-x.

---

Adhesion Factor Modeling Based on Cellular Automata

---

1.description

---

In the wet experiment section, we designed three adhesion factors—CM, CMC, and CMCC—to allow engineering bacteria to persist longer in the intestine. However, we currently do not know which of these three adhesion factors has the best adhesion effect. Therefore, this model is primarily used to simulate the changes in cells within the intestine after adding these three adhesion factors, comparing their adhesion effects, and ultimately selecting an adhesion factor suitable for experimentation.

---

2.Cellular Automata Model Framework

---

Cellular Automaton is a discrete mathematical model that simulates the dynamic evolution of complex systems through simple local rules[1]. Below, I will introduce the basic components of the cellular automaton in this model:

(1)Cells: We assume that the intestine is an ideal cylindrical structure, with intestinal cells uniformly distributed along the side of the cylinder. By unfolding the side of the cylinder, we obtain a rectangular prism. This rectangular prism is divided into many small grids, where each small grid represents the minimum unit of activity for the dry experiment, and not every small grid corresponds to a single intestinal cell.

(2)Grid: The model uses a two-dimensional grid, where each cell can interact with the surrounding cells.

(3)State: Each cell can be in one of several finite states, which represent different physical or biological conditions in the model. In this model, the states may indicate the following situations: empty, probiotics dominant, pathogens dominant, and a mixture of probiotics and pathogens.

(4)Neighborhood: The neighborhood of a cell defines which surrounding cells it interacts with. In this experiment, the neighborhood includes the cell itself and all cells within a distance of 2 units up, down, left, and right (including diagonal directions).

---

3.Rules of the model

---

(1)Intestinal cells can attach to either probiotics or pathogens; however, if the intestinal cells absorb engineered bacteria that contain adhesive factors, their attachment ability will be significantly enhanced. The greater the number of engineered bacteria within a small grid, the stronger the attachment ability of the intestinal cells in that grid. There is an upper limit (saturation effect) to this enhancement of attachment ability by engineered bacteria; when the quantity of engineered bacteria exceeds a certain threshold, the attachment ability no longer increases.

(2)The initial state of the intestinal cells contains both probiotics and pathogens, as the simulation represents the intestinal environment of an individual with inflammatory bowel disease. Therefore, we assume that the initial ratio of probiotics to pathogens is 2:8[2], and they are randomly distributed across the intestinal cells.

(3)The probiotics and pathogens attached to the intestinal cells will proliferate over time. In a given small grid, if the number of probiotics exceeds the number of pathogens, the proliferation rate of the pathogens will slow down. Conversely, if the number of pathogens exceeds the number of probiotics, the proliferation rate of the probiotics will slow down[3].

(4)Aside from the intestinal cells, all other cells will die after a certain period of time.

(5)For the engineered bacteria that are attached to the intestinal cells, their proliferation rate is constant. However, their death does not affect the attachment ability of the intestinal cells, as it is the adhesive factors secreted by the engineered bacteria that play a crucial role. These adhesive factors do not disappear even if the engineered bacteria die.

(6)When probiotics dominate within a cell, they will diffuse into the neighboring cells. Engineered bacteria will also diffuse, but their diffusion is not influenced by the number of probiotics or pathogens in the neighborhood. Once the concentration of engineered bacteria within a cell reaches saturation, they will begin to diffuse into the neighboring cells.

(7)Engineered bacteria exert their effects by releasing adhesive factors, and different adhesive factors exhibit varying preferences for attaching to probiotics and pathogens. For each adhesive factor, specific parameters are set that dictate its adhesion capacity for probiotics and pathogens, with the parameter values derived from empirical wet lab data.

(8)Each cell can only accommodate a limited number of probiotics or pathogens. However, the introduction of engineered bacteria increases this upper limit. The more engineered bacteria present, the higher the upper limit for the capacity of probiotics and pathogens in that cell. Nonetheless, there is also a maximum number of engineered bacteria that each cell can hold.

The above are the main rules of this model. For the definitions of the various values in the rules, please refer to the code.

---

4.Results and Analysis

---

Based on the above model, we simulated the changes in the number of probiotics and pathogens on intestinal cells after introducing different engineered bacteria into the intestines of individuals with inflammatory bowel disease. We used a health status index to represent this, calculated as the number of probiotics minus the number of pathogens within each cell. Negative values are shown in blue, while positive values are shown in red, with deeper colors indicating a higher quantity of corresponding probiotics or pathogens. In the result graphs, we found that the effect of CM was the poorest, with pathogens ultimately dominating. The differences between CMC and CMCC were not significant at first glance; however, upon closer inspection, it is evident that CMC reached a healthy state—where probiotics dominated—more quickly than CMCC. Therefore, we chose CMC as the adhesive factor for this experiment.

In the construction of Model 2, we considered various biologically realistic processes, including the lifespan and proliferation mechanisms of bacteria, the enhancement effects of engineered bacteria, dynamic changes in adhesion capabilities, saturation effects, and resource limitations. This approach enhances the biological plausibility of the model's results and aids in understanding the complex dynamics within the intestinal microbial ecosystem. Additionally, compared to many continuous mathematical models (such as partial differential equation models), our model employs discrete states and rules, thereby avoiding the complexities associated with solving differential equations. Furthermore, since the update of each cell at every time step depends only on the states of neighboring cells, the computational process of the model can be parallelized across multiple processors, resulting in faster simulation speeds.

In summary, this model not only has advantages such as simulation capability, efficiency, and clarity, but it also simulates the environmental changes in the intestines of individuals with inflammatory bowel disease after the introduction of engineered bacteria with different adhesive factors. Among the three adhesive factors established in wet experiments, CMC was selected as the most suitable adhesive factor, guiding the next steps in the wet experiments.

---

Reference

[1] G. D. R. Codd (1968). Cellular Automata. Scientific American, 218(6), 119-132.

[2] Lloyd-Price J, Arze C, Ananthakrishnan AN, Schirmer M, Avila-Pacheco J, Poon TW, Andrews E, Ajami NJ, Bonham KS, Brislawn CJ, Casero D, Courtney H, Gonzalez A, Graeber TG, Hall AB, Lake K, Landers CJ, Mallick H, Plichta DR, Prasad M, Rahnavard G, Sauk J, Shungin D, Vázquez-Baeza Y, White RA 3rd; IBDMDB Investigators; Braun J, Denson LA, Jansson JK, Knight R, Kugathasan S, McGovern DPB, Petrosino JF, Stappenbeck TS, Winter HS, Clish CB, Franzosa EA, Vlamakis H, Xavier RJ, Huttenhower C. Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases. Nature. 2019 May;569(7758):655-662. doi: 10.1038/s41586-019-1237-9. Epub 2019 May 29. PMID: 31142855; PMCID: PMC6650278.

[3] Martin AJM, Serebrinsky-Duek K, Riquelme E, Saa PA, Garrido D. Microbial interactions and the homeostasis of the gut microbiome: the role of Bifidobacterium. Microbiome Res Rep. 2023 May 10;2(3):17. doi: 10.20517/mrr.2023.10. PMID: 38046822; PMCID: PMC10688804.

---

Modeling the Effects of Anti-inflammatory Peptides on Immune Response Processes

---

1.description

---

In this section, we used an agent-based model to simulate the effects of melittin (Mel) and indole-3-carbinol (Indo) on immune cells, cytokine levels, and colon length in mice. The simulation ran for 7 days with iterations every half day. We can view the immune system of the mice as a complex multi-agent system, where different types of immune cells and inflammatory factors interact to regulate the inflammatory state. Our main objective is to observe the recovery and damage of colon tissue by simulating the immune system's response to different drugs.

---

2.Definition of Model Body (Agents)

---

The agent-based model is a simulation model that studies the behavior of complex systems by simulating the interactions between individual agents. In the immune system, each "agent" can represent different immune cells (such as macrophages, T cells, B cells, etc.), and the interactions among these cells determine the overall immune response process. Below are the model agents in our modeling process, with each agent representing an important component of the intestinal immune system. We have simplified the intestine into the following six parts:

Colon Tissue: Used to record the health status of the intestine (e.g., colon length).

Neutrophils: A type of immune cell responsible for initial defense, primarily functioning to engulf and destroy pathogens. However, excessive activation can lead to tissue damage and inflammation.

Macrophages: Primarily classified into M1 type (pro-inflammatory) and M2 type (anti-inflammatory), they play roles at different stages of the immune response.

T Cells: Can be divided into effector T cells (which promote inflammation) and regulatory T cells (which suppress inflammation).

Pro-inflammatory Factors: Such as TNF-α, IL-6, etc., which can exacerbate intestinal inflammation.

Anti-inflammatory Factors: Such as IL-10, which can suppress the inflammatory response and promote tissue repair.

---

3.Agent Behavior Rules

---

Next, we need to define the behaviors of each agent and how they interact with one another. Below are the key rules established:

Neutrophils, Macrophages, and T Cells:

(1)Under conditions of intestinal inflammation, the number of these cells increases, leading to an exacerbation of inflammation.

(2)They will migrate to the site of inflammation, release pro-inflammatory or anti-inflammatory factors, and interact with other cells.

(3)The melittin (Mel) significantly reduces their numbers, while indole-3-carbinol (Indo) also reduces the number of these cells, but to a lesser extent.

Pro-inflammatory factors and anti-inflammatory factors:

(1)Pro-inflammatory factors are secreted by neutrophils, M1 macrophages, and effector T cells, further attracting more immune cells to the site of inflammation. This accelerates the division of immune cells, exacerbates the degree of inflammation, and leads to a reduction in colon length.

(2)Anti-inflammatory factors are secreted by neutrophils, M2 macrophages, and regulatory T cells. These factors inhibit pro-inflammatory signals and promote tissue repair, slowing down the proliferation of immune cells and reducing the destructiveness of the inflammatory response, thereby increasing colon length.

Colon Tissue:

(1)In the state of colitis, the length of the colon is typically shortened, indicating damage to the colon tissue due to inflammation.

(2)When pro-inflammatory factors decrease and anti-inflammatory factors increase, inflammation is reduced, and the length of the colon is restored.

---

4.Mechanism of Action of Melittin (Mel) and Indole-3-carbinol (Indo)

---

Melittin :

1.Inhibit the expression of pro-inflammatory factors, such as TNF-α and IL-6.

2.Significantly reduce the number of neutrophils, macrophages, and T cells.

3.Their effects are relatively strong, leading to a notable reduction in immune cells and pro-inflammatory factors, potentially accelerating intestinal recovery.

Indole-3-carbinol（Indo）:

1.Primarily acts by promoting the expression of anti-inflammatory factors (such as IL-10), balancing the immune response.

2.Slightly reduces the number of immune cells.

3.The effect of indole-3-carbinol (Indo) is milder in terms of immune modulation, and it can promote intestinal repair through anti-inflammatory factors.

---

5.Main Equations and Parameters of the Model

---

During the entire experiment simulation, the mice experienced inflammation for the first three days. At halfway through the third day, we introduced the medication to treat colitis and observe its effects. To quantify and simulate the process in the model, we need to set various parameters. Below are some key formulas and variables(Please refer to the code for specific parameter values):

The formula for updating colon length: The colon length starts from a normal value of 8.7 cm. As inflammation progresses, the colon length shortens, but after the administration of medication, the colon length begins to recover.

\[ L_{t+1}=\begin{cases}L_t-2\cdot r_d\cdot I_t,&\mathrm{if~}t<7\\L_t-2\cdot r_d\cdot I_t+0.55\cdot r_r\cdot A_t,&\mathrm{if~}t\geq7&\end{cases} \]

Parameter	Description
\( L_{t} \)	The colon length on day t
\( I_{t} \)	The level of pro-inflammatory factors on day t
\( A_{t} \)	The level of anti-inflammatory factors on day t
\( r_{d} \)	Rate of colon damage
\( r_{r} \)	Rate of colon recovery due to the medication
\( t_{0} \)	Time Point of Drug Efficacy

Immune Cell Count Update Formula:The dynamic changes in immune cell numbers are influenced by the natural cell generation rate, the effects of the medication, and the presence of both pro-inflammatory and anti-inflammatory factors. Pro-inflammatory factors accelerate immune cell proliferation, while anti-inflammatory factors inhibit their growth.

\[ C_{t+1}=\begin{cases}C_t+r_r\cdot R_1\cdot(1+\gamma_\text{pro}I_{\text{pro},t}-\gamma_\text{anti}I_{\text{anti},t})&\text{if }t <7\\C_t+r_r\cdot R_1\cdot(1+\gamma_\text{pro}I_{\text{pro},t}-\gamma_\text{anti}I_{\text{anti},t})-d_\text{eff}\cdot R_2&\text{if }t\geq7\end{cases} \]

Parameter	Description
\( C_{t} \)	Number of immune cells on day t
\( r_{r} \)	Natural rate of cell generation
\( R_{1} \)	Random fluctuations in cell generation
\( \gamma_{\mathrm{pro}} \)	Influence coefficient of pro-inflammatory factors
\( I_{\mathrm{pro},t} \)	Level of pro-inflammatory factors
\( \gamma_{\mathrm{anti}} \)	Influence coefficient of anti-inflammatory factors
\( I_{\mathrm{anti},t} \)	Level of anti-inflammatory factors
\( d_{\mathrm{eff}} \)	Attenuation coefficient of drug effects
\( R_{2} \)	Random fluctuations in drug effects

Cytokine Update Formula: The rules for updating cytokines (either pro-inflammatory or anti-inflammatory factors) include the activity of immune cells, the effects of the medication, and the weighted influence of immune cells on cytokine levels.

\[ I_{t+1}=I_t+A_\text{cell}+(w_\text{neutrophil}\cdot N_t+w_\text{macrophage}\cdot M_t+w_\text{T-cell}\cdot T_t)+\Delta I_\text{pro}+\Delta I_\text{anti} \]

Parameter	Description
\( I_{t} \)	Factor level on day t
\( A_\mathrm{cell} \)	Baseline influence of cell activity
\( N_{t} \)	Number of neutrophils on day t
\( w_{\text{neutrophil}} \)	Weight of neutrophils
\( M_{t} \)	Number of macrophages on day t
\( w_{\text{macrophage}} \)	Weight of macrophages
\( T_{t} \)	Number of T cells on day t
\( w_{\mathrm{T-cell}} \)	Weight of T cells
\( \Delta I_{\mathrm{pro}} \)	Drug effect on pro-inflammatory factors (effective from day 7 onwards)
\( \Delta I_{\mathrm{anti}} \)	Drug effect on anti-inflammatory factors (effective from day 7 onwards)

---

6.Results and Analysis

---

Figure: The diagram shows the changes in multiple physiological indicators over time after the addition of melittin and indole-3-carbinol during colitis. t = 6 represents the starting point of drug treatment. From the result graph, it can be observed that:

In conclusion, melittin (Mel) demonstrates stronger anti-inflammatory effects compared to indole-3-carbinol (Indo) in the treatment of colitis. Melittin can more rapidly reduce the numbers of neutrophils, macrophages, and T cells, and effectively suppress the levels of pro-inflammatory cytokines, thereby helping the intestines recover to a normal state more quickly.

During the construction of Model 3, we simulated the treatment process of colitis with melittin (Mel) and indole-3-carbinol (Indo), respectively. This model meticulously reflects the changes in various indicators throughout the treatment process. By comparing the therapeutic effects of different drugs on colitis through parameter values and equations, we selected the most appropriate peptide, which is melittin, used in our wet experiments. Therefore, Model 3 is instructive, as it validates the results of our wet experiments.

---

Evolution of Anti-Inflammatory Peptide Functions

---

1.description

---

Our project utilizes anti-inflammatory peptide that is a modified form of melittin, which is a peptide composed of 26 amino acid residues. We have combined two monomeric melittin peptides using a linker peptide sequence to create a composite melittin (D-Melittin) containing 75 amino acids. The anti-inflammatory mechanism of melittin is quite complex; in our modeling, we hypothesize that it inhibits inflammation pathways by binding to TLR4, thereby blocking the downstream activation of NF-kB [1]. Additionally, melittin exhibits certain cytotoxic abilities by forming pores in microbial membranes through interactions with phospholipids. However, due to the lack of species conservation in the targeting of phospholipids by melittin, it can also damage human cell membranes, which has limited its clinical applications [2]. In this project, the primary role of melittin is to attenuate the excessive immune response of immune cells rather than to kill microorganisms. We intend to conduct virtual screening through protein mutagenesis to test whether the mutations can enhance stability, improve binding affinity to TLR4, and reduce cytotoxic effects on cell membranes.

---

2. Virtual Screening for Enhanced Anti-Inflammatory Effects

---

Due to being a synthetically produced protein, there is no structural domain for this protein in the PDB database. Therefore, we utilized AlphaFold 3 for structural prediction, with the results shown in the figure below. We chose TLR4, which has been studied most extensively, as a representative of TLRs. TLR4 consists of 839 amino acids; however, there is no complete structure of TLR4 in the PDB database. We obtained its high-confidence structure from the AlphaFold predicted structure database (AFDB accession: AF-O00206-F1). To determine the binding site of the receptor, we retrieved its amino acid sequence from the UniProt database (ID: O00206) and then searched for different domains of the amino acid sequence in the InterPro database (version 101.0), as illustrated in the figure below. Among these, the sequence recognized by TLR4 for PAMPs corresponds to the LRR domain (leucine-rich repeat), which can be seen at positions 57-114 and 447-508 aa of TLR4. We used the HADDOCK 2.4 [3] online platform for molecular docking. The docking binding site was specified as the LRR domain of the receptor, while it was assumed that every amino acid of the melittin monomer could participate in the binding process. Our docking results will be applied in the next step of random mutagenesis.

We used Rosetta software (version 3.10) to perform saturation mutagenesis at 10 amino acid sites of the Melittin protein. These 10 amino acids include six charged residues at positions 22-27, specifically KRKRQQ; mutating these may affect binding due to changes in the surface charge of the protein. The four amino acids at positions 12-15, located at the turn of the peptide segment, are TGLP; mutating these may also impact binding by altering the spatial structure of the protein. The figure below shows the changes in the free energy of D-Melittin protein before and after mutation, as well as the overall change in free energy when D-Melittin binds to TLR4. From the three heatmap results, we can see that the effect of the mutation 27Q➡Y is the most favorable. The mutated Melittin was similarly subjected to structural prediction using AlphaFold 3, followed by molecular docking, with the docking three-dimensional structure shown in the figure below.

---

3. Verification of Reduced Cytotoxicity

---

Cytotoxicity is generally caused by the hemolytic symptoms resulting from the binding of bee venom peptides to the cell membrane of red blood cells, with the targets of the bee venom peptides being lipid molecules in the cell membrane. To simulate real human conditions, we used CHARMM-GUI to generate a phospholipid bilayer model of human cell membranes. The distribution and composition of various lipid molecules in the inner and outer membranes are shown in the table below, and the generated three-dimensional structural diagram is presented in the figure below.After screening a series of mutations for stability and the most stable binding, we performed molecular docking of the identified mutants with the phospholipid bilayer, while also comparing them with D-Melittin and Melittin proteins. Due to the random nature of binding with the phospholipid bilayer, no specific binding sites were designated; the docking results are shown in the figure below. From the results, it can be observed that the binding strength of the engineered bee venom peptides to the lipid molecules of the cell membrane is significantly lower than that of the unmodified bee venom peptides. Moreover, specific amino acid mutations in the modified bee venom peptides further decrease their binding strength to the cell membrane, leading to reduced cytotoxicity.In summary, the mutation of 27Q➡Y enhances the inhibitory capacity of the peptides' binding to TLR4 while also alleviating cytotoxicity arising from interactions with cell membranes to a certain extent. This provides a potential strategy for optimizing anti-inflammatory peptide components in engineered bacteria in the future.

---

Reference

[1]Ahmedy OA, Ibrahim SM, Salem HH, Kandil EA. Antiulcerogenic effect of melittin via mitigating TLR4/TRAF6 mediated NF-κB and p38MAPK pathways in acetic acid-induced ulcerative colitis in mice. Chem Biol Interact. 2020 Nov 1;331:109276. doi: 10.1016/j.cbi.2020.109276. Epub 2020 Sep 28. PMID: 33002459.

[2]Zhang HQ, Sun C, Xu N, Liu W. The current landscape of the antimicrobial peptide melittin and its therapeutic potential. Front Immunol. 2024 Jan 22;15:1326033. doi: 10.3389/fimmu.2024.1326033. PMID: 38318188; PMCID: PMC10838977.

[3]R.V. Honorato, M.E. Trellet, B. Jiménez-García1, J.J. Schaarschmidt, M. Giulini, V. Reys, P.I. Koukos, J.P.G.L.M. Rodrigues, E. Karaca, G.C.P. van Zundert, J. Roel-Touris, C.W. van Noort, Z. Jandová, A.S.J. Melquiond and A.M.J.J. Bonvin. The HADDOCK2.4 web server: A leap forward in integrative modelling of biomolecular complexes. Nature Prot., In Press (2024).