Overview

The model simulates real-world states and processes using mathematical and computational techniques, offering deep insights into complex problems across various fields through simplification and abstraction. In the field of synthetic biology, modeling has become a crucial tool for studying complex phenomena such as protein structure, promoter sequences, gene expression, and system behavior. It not only reduces the costly and time-consuming trial-and-error stages in experimental processes but also helps scientists more efficiently understand biological systems, optimize project design, and thereby enhance the success rate and efficiency of research.

2024 JLU-NBBMS Model Section is neatly divided into the following three main operational segments:

1. Targeting Module Optimization

2. Gene Selection and Analysis

3. Safety Assessment

These three components work in close coordination, providing strong support for project design. Modeling and experimentation complement each other: experiments not only present the problems and requirements for modeling but also provide critical data for model construction and refinement. At the same time, modeling significantly enhances experimental efficiency and accuracy by optimizing experimental conditions, predicting outcomes, completing visualizations, and validating hypotheses. This close interaction systematizes and refines research progress, facilitating more innovative breakthroughs.

In the Targeting module, structure simulation, molecular docking, and molecular dynamics simulation were employed to enhance Salmonella's targeting ability. Through pan-cancer analysis, integrins were found to be highly expressed in various tumors, leading to the selection of RGD peptides as the targeting unit. Multiple RGD variants were analyzed for affinity differences with integrins. Using AlphaFold3, a fusion protein structure containing RGD, iRGD, and RGDFK was constructed and optimized. Molecular docking and dynamics simulations confirmed that Lpp'-OmpA-RGD exhibited the most stable binding with integrin, making it the final choice for the targeting unit.

In the Gene Selection module, we initially selected tumors from 14 different sites. For each tumor type, we conducted analysis and identified candidate genes suitable for RNAi. Specifically, we first obtained bulk transcriptome data from the tumors, followed by differential expression analysis. We then performed gene enrichment analysis to understand the relevant pathways for future use by iGEM participants. Finally, we conducted prognosis analysis and identified potential RNAi targets.

In the Safety part, we initially used mathematical modeling, including ordinary and partial differential equations (ODEs and PDEs), to explore immune system dynamics following Salmonella injection, focusing on macrophage-cytokine interactions. By integrating experimental data, we simulated bacterial diffusion, macrophage activation, and cytokine regulation, identifying key parameters via sensitivity analysis.

Subsequently, we simplified the PLGA nanoparticle model and applied six mathematical frameworks to investigate drug release mechanisms, such as diffusion, erosion, and swelling. Model fitting and prediction helped optimize the drug delivery system.

Finally, we developed an ODE-based model to quantify Salmonella growth dynamics and delayed lysis within host cells, revealing that delayed lysis significantly influences infection dynamics, thus supporting its potential for enhancing safety.

Targeting Design and Analysis Based on Structural Simulation, Molecular Docking, and Molecular Dynamics Simulations

1.Introduction

In wet lab experiments, we aim to use Salmonella-based bacteriotherapy for tumor treatment. Simultaneously, we seek to enhance the tumor-targeting specificity of Salmonella. Based on this objective, we designed the first part of our in dry lab experiment: “Targeting”.

Preliminary experimental results indicated that integrins, composed of alpha(α) and beta(β) chains, are highly expressed in the majority of tumor tissues across various cancers (Figure 1). Therefore, we intend to utilize integrins as the target for achieving targeted therapy with Salmonella.

Figure 1. Pan-Cancer Differential Expression of Integrin Genes.

The RGD (Arg-Gly-Asp) sequence is a short peptide composed of three amino acids: arginine (Arg), glycine (Gly), and aspartic acid (Asp). The RGD sequence is an essential functional fragment of extracellular matrix proteins, such as fibronectin, and it can recognize and bind to integrin receptors on the cell surface. This interaction regulates various biological processes, including cell adhesion, migration, proliferation, and differentiation.[1] Due to its specific binding to integrins, the RGD sequence is widely used in biomedical applications, including drug delivery, tissue engineering, and anti-tumor therapy.[2]

With ongoing development, several modified forms of RGD peptides have been generated, such as cyclo(-RGDfK) and iRGD, among others (Figure 2).

Figure 2. Standard RGD Peptide and Modified RGD Peptides.

It is important to clarify that these modified RGD peptides differ primarily in function from the standard RGD peptide, rather than in structural affinity. For instance, iRGD is a 9-amino acid cyclic peptide (sequence: CRGDKGPDC) and a molecular mimic.[3] This peptide demonstrates superior homing efficacy to tumor tissues compared to the standard RGD peptide and can more extensively penetrate extravascular tumor tissue. The underlying reason is its dual functionality: initially, it targets tumors through RGD-mediated homing, and subsequently, it binds to the neuropilin-1 receptor, activating the CendR pathway, which enhances drug penetration and transport within the tumor microenvironment.[4]

In this context, since we are more focused on the affinity between modified RGD peptides and integrins, we will not consider functional aspects here but instead concentrate on binding affinity. We can assume that the affinity of RGD peptides for integrins is positively correlated with the targeting efficacy of Salmonella.

2. Objective of the Experiment

The primary objective of this experiment is to evaluate the binding affinity of the standard RGD peptide and various modified RGD peptides to integrins through structural modeling, molecular docking analysis, and molecular dynamics simulations. The findings aim to provide a foundation for subsequent gene screening, drug release simulations, and therapeutic efficacy assessments.

3.Experimental Procedures and Result Analysis

3.1 Structural Simulation

3.1.1 Methods

The experiment initially involves structural modeling of the standard RGD peptide and various modified RGD peptides (including c(-RGDfK) and iRGD). The specific steps include:

Structural Simulation Cycle 1: We used AlphaFold3 (https://alphafoldserver.com/) to construct three fusion proteins: Lpp-OmpA-RGD, Lpp-OmpA-RGDFK, and Lpp-OmpA-iRGD. We then calculated their ipTM values (predicted accuracy of intra-unit placement).

Before proceeding to “Structural Modeling Cycle 2,” we reviewed the literature and found that the Lpp-OmpA outer membrane display system is a technology for displaying exogenous proteins or peptides on the outer membrane of Gram-negative bacteria. This system utilizes the N-terminal signal sequence of Lpp (lipoprotein) and the transmembrane domain of OmpA (outer membrane protein A) to anchor and stably display the target protein on the bacterial outer membrane. Lpp directs the protein to the outer membrane, while OmpA provides a stable anchoring point. This system is widely used in antibody screening, vaccine development, and enzyme engineering. By displaying antigens or enzymes on the cell surface, the Lpp-OmpA system enhances screening efficiency, stimulates immune responses, and improves biocatalytic efficiency, demonstrating significant potential in biotechnology and industrial applications. In practice, the system comprises the signal sequence of Lpp lipoprotein (20 amino acids), the first 9 amino acids of Lpp lipoprotein, and the N-terminal transmembrane domain sequence of OmpA.[5] The modified fusion proteins were named Lpp'-OmpA-RGD, Lpp'-OmpA-RGDFK, and Lpp'-OmpA-iRGD, and based on this, we proceeded with “Structural Modeling Cycle 2.”

Structural Simulation Cycle 2: We used AlphaFold2 (https://neurosnap.ai/) to re-simulate and analyze the structures of Lpp'-OmpA-RGD, Lpp'-OmpA-RGDFK, and Lpp'-OmpA-iRGD.

3.1.2 Results

The results of Structural Simulation Cycle 1 (Figure 3) indicate that all configurations have ipTM values of 0.41, which is below 0.5. This suggests that the predicted models may be inaccurate and carry a significant risk of failure.

Figure 3. Protein Structure Simulation Cycle 1.

The results of the structural simulation cycle 2 (Figure 4) show that the ipTM values are 0.58, 0.59, and 0.54, which are higher than 0.5 but lower than 0.8, indicating that the predictions remain in an unreliable range.

Figure 4. Protein Structure Simulation Cycle 2.

3.1.2 Discussion

Both simulations indicate that the low ipTM values suggest poor stability and accuracy of these structural prediction models, necessitating further optimization and validation. However, the stability of monomers does not fully determine the stability of the oligomers; factors such as surrounding environment pH, temperature, and ionic concentration also play a role in protein oligomer stability. Therefore, we decided to proceed with molecular docking studies to investigate the binding affinity.

3.2 Molecular docking

The docking analysis was primarily conducted through molecular dynamics simulations with integrins, using the structure composed of the alpha and beta chains from the molecular structure with ID 1L5G (https://www.rcsb.org/structure/1L5G) as the receptor molecule. Preparatory work for the docking experiments included constructing the following ligands (see table below):

FTable 1. Ligand protein information.

The ligands were then docked individually with the integrin receptor.

3.2.1 Methods

The detailed process of using AlphaFold3 for protein docking, PyMOL for visualization, and PDBePISA for protein binding energy and interface amino acid analysis is as follows:

1. Using AlphaFold3 for Protein Docking

AlphaFold3 is a deep learning-based protein structure prediction tool that can predict the three-dimensional structure of individual proteins as well as protein-protein complexes (docking predictions) [6].

Preparing Input Sequences:

Collect the amino acid sequences of two or more target proteins and prepare these sequences as input for AlphaFold3.

Running AlphaFold3:

Input the prepared amino acid sequences into AlphaFold3 and select the appropriate settings for protein docking simulations. AlphaFold3 will generate predicted three-dimensional structures of the resulting complex.

Reviewing Results:

Upon completion of the run, AlphaFold3 will produce multiple models. Select the most likely model (typically based on pLDDT scores or other model evaluation criteria) as the basis for subsequent analysis.

2. Using PyMOL for Visualization

PyMOL is a powerful molecular visualization tool that can be used to display and analyze the protein complex structures generated by AlphaFold3 [7].

Loading Structure Files:

Load the PDB files generated by AlphaFold3 into PyMOL using the 'load' command or through the graphical interface.

Analysis and Visualization:

Utilize various PyMOL features, such as coloring and displaying different rendering modes (e.g., wireframe, surface, stick models), to analyze the protein complex structure. Measure distances between key residues within the complex using the 'distance' command, or use the 'show' and 'hide' commands to highlight important regions.

Generating Images:

Adjust the viewpoint and rendering settings to meet visualization requirements. Use the 'ray' command to produce high-resolution rendered images, and save the images using the 'png' command.

3. Using PDBePISA for Protein Binding Energy and Interface Amino Acid Analysis

PDBePISA (Proteins, Interfaces, Structures, and Assemblies) is a tool for analyzing protein complexes, particularly for calculating binding energies and identifying interface amino acids [8].

Submitting Structures to PDBePISA:

Upload the PDB files generated by AlphaFold3 to the PDBePISA server. Select the functionality for analyzing the complex and initiate the calculation.

Interface Analysis:

PDBePISA will identify all interfaces formed within the complex and list the amino acid residues at these interfaces. Key amino acids at the interface can be determined based on changes in solvent-accessible surface area (ΔASA).

Calculating Binding Energy:

PDBePISA will calculate and provide the binding energies of the various interfaces within the complex, typically including ΔG (change in free energy), which helps in assessing the stability and binding affinity of the complex.

Interpreting Results:

Based on the analysis results from PDBePISA, key interface amino acids in the complex can be identified and the binding energies assessed, which is crucial for understanding protein-protein interactions.

3.2.2 Results

Molecular docking reveals differences in binding modes and affinities between various RGD-modified peptides and integrins.

Figure 5illustrates the binding interactions of Protein1, Protein2, and Protein3 with integrin. Visual inspection indicates that Protein1 (Lpp'-OmpA-RGD) exhibits a more intimate binding with integrin.

Figure 5. Molecular Docking Visualization.

Additionally, the contact area (interference area) and binding free energy provide further data-based evidence. Figure 6A clearly shows that the contact areas of the three binding interfaces between Protein1 and integrin are larger than those between Protein2 or Protein3 and integrin. Moreover, Figure 6A indicates that the absolute values of the binding free energy ΔG for Protein1 with integrin are greater compared to those for Protein2 and Protein3 with integrin.

Figure 6. Molecular Docking Parameters.

As shown in Figure 7, the stability of the binding between Protein1 and integrin can also be inferred from the number of hydrogen bonds and the distances between them. The figure demonstrates that Protein1 forms more hydrogen bonds with integrin, and these hydrogen bonds are closer in distance compared to those formed by Protein2 and Protein3, suggesting a more stable interaction.

Figure 7. Molecular Docking Hydrogen Bond Diagram.

3.2.3 Discussion

Simply considering the contact area (Interference area), binding free energy (ΔG), and the number/length of hydrogen bonds, it can be concluded that Lpp'-OmpA-RGD is the most effective. However, it can also be observed from the above figure that Lpp'-OmpA-RGDFK exhibits relatively good values for these parameters. Since proteins are dynamic during actual molecular binding (i.e., they are in a fluctuating state), using molecular dynamics simulations to model the dynamic process would be more convincing and facilitate a more accurate assessment of affinity.

3.3 Molecular dynamics simulation

Molecular dynamics simulations further analyzed the interactions of the above ligands (Lpp'-OmpA-RGD, Lpp'-OmpA-RGDFK, Lpp'-OmpA-iRGD) with integrins in a dynamic environment to observe their stability and binding behavior in the simulated environment.

3.3.1 Methods

1.Molecular Simulation

Molecular simulations were performed using Amber 22 (San Francisco, CA, USA) [7]. The ff19SB force field [8] was used to calculate system force field parameters. The TIP3P water model was employed for solvation, and counterions were added to neutralize the system. Once the system energy was minimized, the system was heated from 0 K to 310.15 K within 500 ps. The system was then equilibrated under the NVT ensemble, followed by a pre-equilibration at 310.15 K. Finally, a 100 ns molecular simulation was conducted under the NPT ensemble with periodic boundary conditions. All covalent bonds involving hydrogen were constrained using the SHAKE method. Dynamics results were analyzed using AmberTools23[9][10].

Note: In this simulation, no multiple rounds of repetition were performed. Therefore, the results have a certain degree of randomness.

Note:

1. The NVT ensemble (Canonical Ensemble) maintains a constant number of particles (N), volume (V), and temperature (T) during the simulation. NVT ensemble is used for preliminary processes in the production run, such as system energy minimization, heating, and pre-equilibration.

2. The production run simulation uses the NPT ensemble (Isothermal-Isobaric Ensemble), which maintains constant number of particles (N), pressure (P), and temperature (T). This setup is appropriate for simulating biomolecules under realistic conditions, as most experimental conditions are conducted at constant temperature and pressure.

3. In the simulation, the formal production step has a timestep of 0.002 fs, and the simulation runs for 250,000,000 steps. The time conversions are: 1000 fs = 1 ps, 1000 ps = 1 ns, 1000 ns = 1 μs, 1000 μs = 1 ms, 1000 ms = 1 s.

4. The pressure is set to 1.0 bar = 100,000 Pascals (Pa) = 100 kilopascals (kPa), which is considered to be 1 standard atmosphere.

5. Counterions Na⁺ and Cl^- were added using the Monte Carlo method in the solvent, with a concentration of 0.15 M.

2.Root Mean Square Deviation (RMSD) Analysis and Root Mean Square Fluctuation (RMSF) Analysis

RMSD measures the root-mean-square deviation of the target structure relative to the reference frame. An RMSD value of 0 indicates perfect superposition. It is defined as follows (Equation 1):

RMSF is a measure of the deviation between the position of particle \(i\) and its reference position, averaged over each residue. An RMSF value of 0 indicates no atomic fluctuation. It is defined as follows (Equation 2):

3. Binding Free Energy Analysis

In Amber 22, the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) calculations were performed using MMPBSA.py. The binding free energy of the interaction is derived as follows (Equation 3):

In Equation (3), the free energy G is calculated according to Equation (4):

Here, ΔE_Gas represents the gas-phase molecular mechanical energy, ΔG_Sovation represents the solvation free energy, and TΔS represents the entropy term.

ΔE_Gas can be further divided into three components: van der Waals energy (ΔE_vdW), electrostatic energy (ΔE_ELE) and gas-phase internal energy (ΔE_INT).

The solvation free energy G_Solvation is calculated using a continuum solvent model and consists of polar contributions G_PB and nonpolar contributions G_nonpolar.

The conformational entropy -TΔS is estimated through quasi-harmonic analysis of the simulation trajectory. Since its contribution to the results is minimal, it can be neglected.

Binding Free Energy Analysis Definitions:

Internal: Internal energy term

VDWAALS: van der Waals energy term

EEL: Electrostatic energy term

EPB: Electrostatic component of the solvation energy (calculated by the Poisson-Boltzmann (PB) equation)

ENPOLAR: Nonpolar component of the solvation energy

G_gas: VDWAALS + EEL, the energy in the gas phase, which can be understood as the force field energy

G_solv: EPB + ENPOLAR, solvation energy

TOTAL: Sum of all energy terms

DELTA_TOTAL: Binding free energy

4. Energy decomposition contribution analysis

By decomposing the solvent-accessible surface area of the protein onto each residue, the overall binding energy is divided into the binding energies of individual residues.

5. Hydrogen bond analysis

A hydrogen atom covalently bonded to an electronegative atom X can form a hydrogen bond with another highly electronegative and small-radius atom Y when X and Y are in close proximity, resulting in the X-H...Y hydrogen bond interaction.

The criteria for calculating and determining hydrogen bonds are as follows:

(1) The distance between atom X and atom Y is less than 3.0 Å (or 3.5 Å);

(2) The angle between X-H in the hydrogen bond donor and Y in the hydrogen bond acceptor is less than 30 degrees.

6. Radius of gyration (Rg) analysis

The radius of gyration (Rg) can represent the compactness of the protein structure and is also used to characterize changes in the degree of looseness of the protein's peptide chain during simulations.

In molecular dynamics (MD) simulations, the compactness and rigidity of the protein backbone can be evaluated using Rg, and it can also indicate the strength of protein stability under different temperatures. A smaller Rg value corresponds to a more compact and rigid protein structure. It is important to note that Rg is not an indicator of protein flexibility, which is better reflected by RMSF (Root Mean Square Fluctuation).

The calculation formula is as follows:

3.3.2 Results

1.Molecular Simulation

Figure 8. Trajectory visualization.

As shown in the figure, darker colors indicate stronger structural fluctuations. From the trajectory visualization animation, Lpp'-OmpA-RGD and Lpp'-OmpA-iRGD exhibit higher stability.

2. Root Mean Square Deviation (RMSD) Analysis and Root Mean Square Fluctuation (RMSF) Analysis Results

Figure 9. Root Mean Square Deviation (RMSD).

RMSD represents the distance between the same atom in different structures. The RMSD of a protein can reveal the positional changes between the protein's conformations during the simulation and its initial conformation. The trend in the RMSD of both the protein and the ligand is an important indicator of whether the simulation has reached stability. If the ligand may have a destabilizing effect on the protein, the protein's RMSD can not only reflect the stability of the protein structure but also further indicate whether the ligand has a destabilizing effect on the protein.

A lower RMSD value indicates high stability of the protein, while a higher RMSD suggests that the backbone has undergone conformational changes during the simulation.

As shown in Figure 9, the RMSD of Lpp'-OmpA-RGD and Lpp'-OmpA-RGDFK reached a stable state over time, whereas the RMSD of Lpp'-OmpA-iRGD did not stabilize, indicating that the structure post-binding did not reach equilibrium.

Figure 10. Root Mean Square Fluctuation (RMSF).

RMSF represents the average atomic positional changes over time and can characterize the flexibility and degree of motion of individual amino acids throughout the simulation.

The RMSF of a protein is used to determine the fluctuation of each residue during the simulation. Higher RMSF values are associated with regions of greater flexibility, such as loop regions (and intrinsically disordered regions). Additionally, amino acid domains distant from the binding site may also exhibit higher RMSF values.

In this simulation, residues 1-1592 are integrins, while residues 1593 and beyond correspond to different proteins (Lpp'-OmpA-RGD, Lpp'-OmpA-RGDFK, and Lpp'-OmpA-iRGD). Figure 10 shows that the flexibility or degree of motion in the iRGD region is higher, indicating greater instability.

3. Binding free energy analysis results (Differences (Complex - Receptor - Ligand))

Figure 11. Absolute value of binding free energy (ΔG).

According to the binding free energy analysis results, the absolute values of binding free energy for Lpp'-OmpA-RGD, Lpp'-OmpA-RGDFK, and Lpp'-OmpA-iRGD decrease progressively. From an energy perspective, this suggests that the structure of Lpp'-OmpA-RGD is more likely to reach equilibrium. (For detailed data, see "Supplementary Materials.")

4. Results of energy decomposition contribution analysis

Figure 12. Energy decomposition contribution diagram.

As shown in Figure 12, the RGD peptide exhibits a relatively higher contribution to the binding free energy from different amino acid residues. This suggests that, from an energy perspective, the structure of Lpp'-OmpA-RGD is more likely to reach equilibrium.

5. Hydrogen bond analysis

Figure 13. Hydrogen bond analysis diagram.

Hydrogen bonds and hydrophobic interactions play crucial roles in maintaining protein conformations. Additionally, hydrogen bonds between amino acids can alter the internal geometry of the protein. The following diagram shows the variation in the number of hydrogen bonds in the ligand-protein complex. As shown in Figure 13, Protein1 forms more hydrogen bonds with integrin, indicating a more stable binding.

6. Results of radius of gyration (Rg) analysis

Figure 14. Distribution of the radius of gyration (Rg).

The radius of gyration (Rg) characterizes the compactness of the protein structure. A smaller range of Rg (narrower peak) indicates greater compactness or rigidity of the protein. From the figure, it can be observed that the Rg distribution of Lpp'-OmpA-RGD is more concentrated, indicating that the complex of Lpp'-OmpA-RGD is the most compact and stable.

3.3.3 Discussion

Based on the molecular dynamics results, Lpp'-OmpA-RGD shows the highest stability upon binding with integrin, which indirectly indicates that the two have the best affinity.

4.Conclusion

Based on the results from structural simulations, molecular docking, and molecular dynamics simulations, the modified Lpp-OmpA-RGD (i.e., Lpp'-OmpA-RGD) exhibits better binding effectiveness and targeting ability.

Reference:

[1]Plow EF, Haas TA, Zhang L, Loftus J, Smith JW. Ligand binding to integrins. J Biol Chem. 2000 Jul 21;275(29):21785-8. doi: 10.1074/jbc.R000003200. PMID: 10801897.

[2]Ley K, Rivera-Nieves J, Sandborn WJ, Shattil S. Integrin-based therapeutics: biological basis, clinical use and new drugs. Nat Rev Drug Discov. 2016 Mar;15(3):173-83. doi: 10.1038/nrd.2015.10. Epub 2016 Jan 29. PMID: 26822833; PMCID: PMC4890615.

[3]Sugahara KN, Teesalu T, Karmali PP, Kotamraju VR, Agemy L, Greenwald DR, Ruoslahti E. Coadministration of a tumor-penetrating peptide enhances the efficacy of cancer drugs. Science. 2010 May 21;328(5981):1031-5. doi: 10.1126/science.1183057. Epub 2010 Apr 8. PMID: 20378772; PMCID: PMC2881692.

[4]Ruoslahti E. Tumor penetrating peptides for improved drug delivery. Adv Drug Deliv Rev. 2017 Feb;110-111:3-12. doi: 10.1016/j.addr.2016.03.008. Epub 2016 Apr 1. PMID: 27040947; PMCID: PMC5045823.

[5]Francisco JA, Campbell R, Iverson BL, Georgiou G. Production and fluorescence-activated cell sorting of Escherichia coli expressing a functional antibody fragment on the external surface. Proc Natl Acad Sci U S A. 1993 Nov 15;90(22):10444-8. doi: 10.1073/pnas.90.22.10444. PMID: 8248129; PMCID: PMC47793.

[6]Abramson J, Adler J, Dunger J, et al. Accurate structure prediction of biomolecular interactions with AlphaFold 3[J]. Nature, 2024: 1-3.

[7]"Open-Source PyMOL". Schrodinger, Inc. 5 November 2021. Retrieved 7 November 2021.

[8]Paxman JJ, Heras B. Bioinformatics Tools and Resources for Analyzing Protein Structures. Methods Mol Biol. 2017;1549:209-220. doi: 10.1007/978-1-4939-6740-7_16. PMID: 27975294.

[9]Case, David A., et al. "AMBER 22 Reference Manual." (2022).

[10]Tian, Chuan, et al. "ff19SB: Amino-acid-specific protein backbone parameters trained against quantum mechanics energy surfaces in solution." Journal of chemical theory and computation 16.1 (2019): 528-552.

[11]Roe, Daniel R., and Thomas E. Cheatham III. "PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data." Journal of chemical theory and computation 9.7 (2013): 3084-3095.

[12]Case, David A., et al. "AmberTools." Journal of chemical information and modeling 63.20 (2023): 6183-6191.

Supplementary Materials

Overview

In current cancer research, open databases like TCGA^[1] and GEO^[2] provide large-scale transcriptomic data from various tumor patients, greatly advancing cancer studies. By analyzing this data, potential therapeutic targets can be identified, enabling RNA interference to inhibit oncogenes.

We selected 14 tumor types for analysis based on WHO’s tumor mortality ranking^[3], focusing on their low survival rates and high treatment difficulty.

As a user, you can select a tumor of interest. For example, lung cancer includes four treatment options.

1.Carboplatin+PTX(paclitaxel)

2.Cisplatin+Etoposide100

3.Cisplatin+Etoposide137

4.Cisplatin+Navelbine

On our website, users can view the main chemotherapy options for each tumor, with each treatment displayed as images A, B, C, and D.

A.DEG: This section identifies genes with significant expression changes between tumor and normal tissues under different treatments, presented as volcano and heatmaps.

B.Gene Enrichment Analysis: This explores whether the identified differentially expressed genes are enriched in specific biological processes or pathways by linking them to known functions and diseases.

C.Prognosis Analysis: Combines differential gene expression with clinical data to assess if gene expression correlates with patient survival or disease progression.

D.Gene Candidates: Based on differential expression, enrichment, and prognosis analysis, potential drug resistance genes are ranked by HR value. These genes may serve as RNAi targets for future experimental validation and treatment strategies.

Our analysis lays the foundation for RNAi-based tumor target screening, identifying candidate genes and providing a data platform for users to explore tumor targets, supporting personalized treatment strategies.

Lung

Carboplatin+Paclitaxel

Cisplatin+Etoposide 1

Cisplatin+Etoposide 2

Cisplatin+Vinorelbine

Fluorouracil+Leucovorin+Oxaliplatin

Capecitabine

2-Fluoropyrimidine

Cisplatin

AC+lxabepilone

FEC+Docetaxel

TFAC

TFEC

Cyclophosphamide+Doxorubicin 1

Cyclophosphamide+Doxorubicin 2

Fluorouracil

Gemcitabine

Capecitabine

Bicalutamide+Leuprolide

Leuprolide

Cisplatin

Cytarabine+ldarubicin

Doxorubicin+Mercaptopurine+
Methotrexate+Prednisone+Vincristine

Lomustine+Procarbazine+Vincristine

Temozolomide

Carboplatin+Gemcitabine

Cisplatin+Doxorubicin+Methotrexate+Vinblastin

DEG:

Cisplatin+Gemcitabine

Carboplatin+Paclitaxel 1

CarboplatinPaclitaxel 2

Cisplatin+Pemetrexed

Reference

[1] https://www.cancer.gov/ccg/research/genome-sequencing/tcga

[2] https://www.ncbi.nlm.nih.gov/geo/

[3] https://gco.iarc.fr/en

Safety Module

Immune Reaction

1. Introduction

The immune system plays a pivotal role in defending the host against pathogens such as Salmonella, which can cause severe infections leading to significant morbidity and mortality. The complexity of the immune response involves various cellular and molecular interactions, including the activation and regulation of macrophages, the secretion of cytokines, and the overall modulation of the inflammatory response. Understanding these interactions at a quantitative level is essential for developing effective therapeutic strategies and enhancing our knowledge of host-pathogen dynamics.

In recent years, mathematical modeling has emerged as a powerful tool to simulate and predict the behavior of immune responses under various conditions. By leveraging mathematical equations to represent biological processes, models can provide insights into the temporal and spatial dynamics of immune cell populations and cytokine levels. These models serve as a bridge between experimental observations and theoretical predictions, enabling researchers to explore scenarios that may be challenging or impossible to replicate in laboratory settings.

This study focuses on the mathematical modeling of the immune regulatory system's response to Salmonella infection, specifically examining the roles of resting and activated macrophages, as well as the dynamics of key cytokines such as TNF-⍺, IL-6, IL-8, and IL-10. Using a combination of ordinary differential equations (ODEs) and partial differential equations (PDEs), we aim to simulate the spatial and temporal behaviors of these immune components, providing a comprehensive understanding of their interactions during an immune response.

Our model incorporates well-established biological mechanisms, including the logistic growth of bacterial populations, macrophage activation, and cytokine-mediated regulation. By parameterizing the model based on experimental data and literature values, we ensure that the simulations reflect realistic physiological conditions. The results from these simulations offer valuable insights into how the immune system orchestrates its response to Salmonella, highlighting the critical factors that influence the outcome of the infection.

The following sections will outline the assumptions underlying our model, the specific parameters used, and the detailed mathematical formulations of the ODE and PDE models. We will then present the results of our simulations, analyzing the spatial diffusion of bacteria and macrophages, the dynamic interactions between resting and activated macrophages, and the cytokine response under different conditions. Finally, we will discuss the stability of the model parameters and their impact on the overall immune response.

2. Research Objectives

The primary objectives of this study are to investigate the dynamic behaviors and interactions within the immune regulatory system during a Salmonella infection. The study focuses on multiple aspects of the immune response, utilizing both partial differential equations (PDEs) and ordinary differential equations (ODEs) to model and simulate these processes. The specific objectives are outlined as follows:

A. Examine the Spatial Diffusion of Bacteria and Macrophages:

To simulate and analyze the spatial diffusion of Salmonella bacteria, resting macrophages, and activated macrophages within the tissue over time.

B. Investigate the Dynamic Interactions Between Resting and Activated Macrophages:

To model and explore the time-dependent conversion of resting macrophages into activated macrophages and their respective roles in the immune response. Apply the ODE model to simulate the dynamic changes in macrophage populations, focusing on their activation and deactivation processes.

C. Analyze the Impact of Immune Response on Salmonella Bacterial Concentration:

To compare the concentration of Salmonella bacteria in the tissue with and without the presence of an immune response.

D. Simulate and Compare Cytokine Dynamics Against Experimental Data:

To simulate the production and regulation of key cytokines (TNF-⍺, IL-6, IL-8, and IL-10) during an immune response and compare the results with experimental observations.mAssess the accuracy of the model in replicating cytokine dynamics and its potential for predicting immune response outcomes.

E. Conduct a Parameter Sensitivity Analysis:

To evaluate the stability and sensitivity of the model parameters, particularly focusing on the most influential parameters that govern the immune response dynamics.

3. Model Assumptions

3.1 Biological Assumptions

3.1.1 Constant Parameter Stability

It is assumed that all biological parameters (e.g., rates of macrophage activation, cytokine production, bacterial growth, and decay rates) remain constant throughout the simulation. This implies that external factors such as temperature, pH levels, and nutrient availability do not fluctuate significantly during the infection process.

3.1.2 Isolated System

The immune response is modeled as an isolated system where no external infections or pathogens are introduced during the simulation period. This means that the primary focus remains solely on the interaction between the host's immune system and the Salmonella bacteria, without considering additional infections or external interventions.

3.1.3 Uniform Tissue Environment

The tissue in which the immune response occurs is assumed to be homogenous and isotropic, meaning that the properties of the tissue (e.g., permeability, diffusivity) are uniform throughout the spatial domain. This simplifies the diffusion equations by assuming that the movement of cells and molecules is consistent across the tissue.

3.1.4 Neglect of Adaptive Immune Response

The model focuses primarily on the innate immune response, particularly the roles of macrophages and cytokines. The adaptive immune response, including the activation of T-cells and the production of antibodies, is not explicitly modeled. This assumption is made to concentrate on the initial phase of infection and the early immune response.

3.1.5 No Cell Death Beyond Modeled Decay

The model assumes that macrophage death or other cell loss occurs only through the decay terms included in the equations. Other forms of cell death, such as apoptosis or necrosis triggered by factors not modeled (e.g., toxins or other pathogens), are not considered.

3.2 Mathematical Assumptions

3.2.1 Linear Superposition of Diffusion

The diffusion of bacteria and immune cells (macrophages) is modeled using linear partial differential equations, assuming that the diffusion process follows Fick's law and can be linearly superimposed. This means that the effects of diffusion from multiple sources are additive and do not interact non-linearly.

3.2.2 Time-Independent Coefficients

The diffusion coefficients, decay rates, and other parameters are assumed to be time-independent. This simplifies the equations and allows for easier integration and analysis of the model over time.

3.2.3 Boundary Conditions

The model assumes zero-flux boundary conditions at the edges of the spatial domain, meaning that no bacteria or immune cells cross the boundaries of the tissue. This assumption is made to focus on the internal dynamics of the infection site without considering interactions with adjacent tissues.

3.2.4 Spatial and Temporal Discretization

In solving the PDEs, the spatial and temporal domains are discretized into finite grids. It is assumed that the grid resolution is sufficiently fine to capture the critical dynamics of the system without introducing significant numerical errors.

4. Parameter Settings

We extracted and obtained the following experimental parameters from specific experimental data and paper data for numerical fitting and data analysis.

Table 1. Parameters settings

5. Model Formulation

In this section, we present the detailed mathematical formulation of the immune system response to Salmonella infection. The model is divided into two main parts: the partial differential equation (PDE) diffusion model, which captures the spatial dynamics of bacteria and immune cells, and the system of ordinary differential equations (ODEs), which describe the temporal dynamics of cytokine production and macrophage activation.

5.1 PDE Diffusion Model

The PDE diffusion model is employed to simulate the spatial distribution of Salmonella bacteria, resting macrophages, and activated macrophages within the tissue. This model captures the diffusion process, which is governed by the following equation:

Where:

u(x,t) represents the concentration of bacteria or macrophages at position x and time t. D is the diffusion coefficient (e.g. , D_A for bacteria, D_MR for resting macrophages, and D_MA for activated macrophages). μ is the natural decay rate of the respective species.∇²u(x,t) denotes the Laplacian, representing the diffusion process in space.The diffusion coefficient D_A, D_MR, and D_MA dictate how rapidly each species spreads throughout the tissue. The decay rate μ represents the loss of concentration due to natural decay or immune system activity.

5.2 Hill Function Overview

To capture the regulatory effects of cytokines on macrophage activation and other immune processes, we employ Hill functions. These functions describe the saturation effects observed in biological systems, where a response to a stimulus (e.g., cytokine concentration) increases gradually and eventually plateaus.

A typical Hill function is given by:

Where:

x is the concentration of the cytokine or other regulatory molecule. K is the concentration at which the response is half-maximal. n is the Hill coefficient, which determines the steepness of the response curve.

In our model, Hill functions are used to describe the activation and suppression effects of cytokines like TNF-α, IL-6, IL-8, and IL-10 on macrophages and other immune components.

5.3 ODE Model Formulation

The temporal dynamics of the immune response are captured using a system of seven ordinary differential equations (ODEs). Each equation represents the rate of change of a specific biological component over time, accounting for both the activation and decay processes.

5.3.1 Bacterial Growth and Decay

The growth and decay of Salmonella bacteria are modeled by the following equation:

Where:

A represents the bacterial concentration. β_A is the bacterial replication rate. k_A is the carrying capacity of the bacteria. μ_A is the natural decay rate of bacteria. λ_MR and λ_MA are the phagocytosis rates by resting and activated macrophages, respectively. This equation models the logistic growth of bacteria, limited by the carrying capacity k_A, and the reduction in bacterial concentration due to immune system activity.

5.3.2 Resting Macrophage Dynamics

The dynamics of resting macrophages (MRM_RMR) are governed by:

Where:

μ_MR is the resting macrophages' natural decay rate. M_{R_max} is the maximum capacity of resting macrophages. γ_MA and k_MTNF represent the activation rates of macrophages by TNF-α and other factors. This equation captures the recruitment and activation of resting macrophages in response to bacterial presence.

5.3.3 Activated Macrophage Dynamics

The dynamics of activated macrophages (MRM_RMR) are described by:

Where:

μ_MA is the natural decay rate of activated macrophages. The other terms describe the activation of macrophages by TNF-α and their subsequent reduction. M_A(0) = mean(M_R) ̅is an initial condition. This equation reflects the increase in activated macrophages due to the immune response and their eventual decay.

5.3.4 TNF-α Dynamics

The production and decay of TNF-α are governed by:

Where:

k_TNFM is the TNF-α upregulation factor. k_TNF is the TNF-α decay rate. q_TNF represents the baseline concentration of TNF-α. This equation models the cytokine's production in response to immune stimuli and its subsequent decay.

5.3.5 IL-6 Dynamics

The production and decay of IL-6 are described by:

Where:

k_IL6M and k_IL6TNF are the rates of IL-6 production. k_IL6 is the decay rate of IL-6. q_IL6 represents the baseline concentration of IL-6.This equation captures the regulation of IL-6 in response to TNF-α and other cytokines.

5.3.6 IL-8 Dynamics

The production and decay of IL-8 are modeled by:

Where:

k_IL8M and k_IL8TNF are the rates of IL-8 production. k_IL8 is the decay rate of IL-8. k_IL8 represents the baseline concentration of IL-8. This equation models IL-8 production and its regulation by TNF-α and IL-10.

5.3.7 IL-10 Dynamics

The production and decay of IL-10 are described by:

Where:

k_IL10M and k_IL10IL6 are the rates of IL-10 production. k_IL10 is the decay rate of IL-10. k_IL10 represents the baseline concentration of IL-10.This equation captures the regulatory role of IL-10 in modulating the immune response.

This comprehensive set of equations forms the mathematical model that describes the immune response to Salmonella infection. The PDE model captures the spatial dynamics of bacterial and immune cell diffusion, while the ODE model captures the temporal dynamics of cytokine production and macrophage activation. Together, these models provide a detailed framework for understanding the complex interactions within the immune system during infection.

6. Results

In this section, we present the results derived from the mathematical models described earlier. The analysis is divided into several subsections, each focusing on different aspects of the immune response to Salmonella infection. These subsections cover the results from the PDE diffusion model, the ODE model for bacterial proliferation, the dynamics of macrophages, and the cytokine activity.

6.1 Results from the PDE Diffusion Model

The PDE diffusion model simulates the spatial distribution of Salmonella bacteria, resting macrophages, and activated macrophages within the tissue. The model captures the diffusion processes, revealing how these entities spread and interact over time. The results were obtained using the finite difference method to approximate the solutions of the PDEs.

Bacterial Diffusion: The results show that Salmonella bacteria, starting from a localized infection site, rapidly diffuse throughout the tissue. The rate of diffusion is governed by the bacterial diffusion coefficient D_A. The spatial spread of the bacteria is initially fast but slows down as the bacterial concentration approaches the carrying capacity of the environment.

Macrophage Diffusion: Resting macrophages (MR) exhibit slower diffusion compared to bacteria, as reflected by their lower diffusion coefficient D_MR. However, once activated (MA), macrophages spread more rapidly due to a higher diffusion coefficient (D_MA). The model illustrates the dynamic response of macrophages to the presence of bacteria, where activated macrophages rapidly move towards the infection site.

Spatial Distribution: The results indicate that activated macrophages tend to accumulate around areas of high bacterial concentration, creating a focused immune response. This spatial localization is critical for effective pathogen clearance, as it ensures that immune resources are concentrated where they are most needed.

Figure 1. Spatial distribution for activated macrophages

Here, we select this representative image of resting macrophages diffusion simulation when D_MR=4.32 × 10⁻² mm³/day and μ_MR=0.033 day⁻¹ as a typical example to show the trend.

6.2 Bacterial Proliferation (ODE 1)

The first ODE model describes the temporal dynamics of Salmonella bacterial proliferation in the tissue. This model incorporates bacterial growth, natural decay, and reduction due to immune system activity.

In the absence of an immune response, the bacteria exhibit logistic growth, characterized by an initial exponential increase in concentration, followed by a plateau as the population reaches the carrying capacity (k_A). This behavior is typical of bacterial populations in a resource-limited environment.

When the immune response is active, the model shows a significant reduction in bacterial growth. The phagocytosis by resting and activated macrophages(λ_MR and λ_MA) introduces a decay term that counteracts the bacterial replication, leading to a lower steady-state concentration or complete clearance of bacteria depending on the parameters.

Mathematical analysis of the ODE reveals that the system can have multiple steady states, depending on the initial conditions and parameter values. For instance, if the immune response is strong enough (high λ_MA and λ_MR), the bacteria may be completely eradicated. Otherwise, a non-zero steady-state bacterial concentration may persist.

Figure 2. Bacteria concentration-Time curve

6.3 Dynamics of Macrophages (ODE 2/3)

The second and third ODE models describe the dynamics of resting and activated macrophages during the immune response to Salmonella infection.

The results show that resting macrophages (M_R) decrease over time as they become activated in response to bacterial presence. The activation rate is influenced by factors such as M_R is accompanied by a corresponding increase in M_A.

Activated macrophages (M_A) initially increase rapidly as resting macrophages are recruited and activated. Over time, M_A reaches a peak and then gradually declines due to natural decay and the reduction in bacterial concentration (which serves as the activation stimulus).

The figures depict the time course of both M_R and M_A, showing how the balance between these two populations shifts during the immune response. The transition from a resting to an activated state is a crucial aspect of effective immune function, as activated macrophages are more potent in pathogen clearance.

Figure 3. The Status of macrophages-Time curve

6.4 Cytokine Activity and Experimental Comparison (ODE 4/5/6/7)

The fourth to seventh ODE models focus on the production and decay of key cytokines: TNF-α, IL-6, IL-8, and IL-10. These cytokines play critical roles in modulating the immune response.

TNF-α Dynamics: TNF-α is produced rapidly in response to bacterial infection and macrophage activation. The model shows that TNF-α levels rise quickly, peaking early in the infection, and then decline as the bacterial load decreases and the immune response stabilizes.

IL-6 and IL-8 Dynamics: Both IL-6 and IL-8 follow a similar pattern, with rapid production followed by a decline. These cytokines are involved in promoting inflammation and recruiting additional immune cells to the infection site. The model captures the transient nature of their activity, which is crucial for avoiding chronic inflammation.

IL-10 Dynamics: IL-10, a regulatory cytokine, exhibits a delayed response compared to TNF-α, IL-6, and IL-8. The model shows that IL-10 levels increase as the inflammatory response peaks, helping to modulate and eventually suppress excessive immune activity.

The simulated cytokine levels are compared with experimental data to validate the model. The figures show good agreement between the model predictions and observed cytokine dynamics, demonstrating the model's accuracy in capturing the key aspects of the immune response.

Further mathematical analysis of the ODEs reveals how the interaction between cytokines influences the overall immune response. Stability analysis shows that the cytokine network can have stable or oscillatory behavior depending on the parameter values, which may correspond to different infection outcomes. The results underscore the complexity of cytokine interactions and their critical role in orchestrating an effective immune response while avoiding overactivation that could lead to tissue damage.

Figure 4. The Cytokine concentration-Time curve

Figure 5. The Cytokine dynamics-Time curve

7. Sensitivity Analysis Process and Conclusions

7.1 Process Overview

The sensitivity analysis was conducted on the cytokine models for TNF-α, IL-6, IL-8, and IL-10. The parameter chosen for adjustment was the upregulation factor k_TNFM, which directly influences the production of TNF-α by activated macrophages.The selected parameter, k_TNFM, was adjusted by ±10%. Specifically:

·Up 10%: The value of k_TNFM was increased by 10%.

·Down 10%: The value of k_TNFMwas decreased by 10%.

The ODEs describing the dynamics of TNF-α, IL-6, IL-8, and IL-10 were solved numerically using MATLAB’s ode45 solver for the base parameter values as well as the ±10% adjusted values.

7.2 Mathematical Analysis

7.2.1 Impact on TNF-α Dynamics

The sensitivity analysis shows that a ±10% change in 〖k﷩TNFM〗 leads to a visible shift in the concentration of TNF-α over time. Specifically, an increase in k_TNFM results in a higher peak and a slower decay in TNF-α concentration, while a decrease results in a lower peak and faster decay. This suggests that the production of TNF-α is moderately sensitive to changes in the upregulation factor.

7.2.2 Impact on IL-6 Dynamics

IL-6 dynamics show a similar trend, but the sensitivity to k_TNFM is less pronounced compared to TNF-α. This indicates that while TNF-α influences IL-6 production, the overall effect on IL-6 is buffered by other regulatory mechanisms in the model.

7.2.3 Impact on IL-8 Dynamics

The IL-8 concentration appears to be highly responsive to changes in k_TNFM, with the upregulation having a more substantial impact on both the peak concentration and the steady-state level. This is consistent with the role of IL-8 as a key pro-inflammatory cytokine, which is tightly regulated by TNF-α.

7.2.4 Impact on IL-10 Dynamics

IL-10 dynamics are slightly influenced by the changes in 〖k﷩TNFM〗, with a minor shift in the concentration over time. The regulation of IL-10, being a primarily anti-inflammatory cytokine, is less directly tied to TNF-α, which explains the lower sensitivity observed.

7.3 Conclusion

The parameter k_TNFM shows varying degrees of influence on the dynamics of the cytokines studied. TNF-α and IL-8 show the highest sensitivity, reflecting their direct roles in the pro-inflammatory response. IL-6 and IL-10, while affected, exhibit more stable behavior, indicating a more complex regulation that dampens the effect of single parameter changes.

The sensitivity analysis reveals that the cytokine network is robust but still responsive to key regulatory factors. Small changes in these parameters can lead to significant shifts in cytokine levels, which could impact the overall immune response. This highlights the importance of accurate parameter estimation in the model, as well as the potential for targeting these pathways therapeutically.

Figure 6. The Sensitivity analysis results

8. Model Strengths, Limitations, and Potential Extensions

8.1 Strengths

8.1.1 Comprehensive Modeling

The model integrates both spatial (PDE) and temporal (ODE) dynamics, providing a detailed representation of how the immune system responds to Salmonella infection. This dual approach allows for the exploration of both localized immune responses and overall cytokine dynamics.

8.1.2 Modularity

The model's structure is modular, allowing for the easy addition or modification of components, such as introducing new cytokines or immune cells. This flexibility makes it suitable for studying a wide range of infections and immune responses.

8.1.3 Predictive Capability

The model successfully predicts the dynamics of key cytokines (TNF-α, IL-6, IL-8, and IL-10) and macrophages, which are critical for understanding the immune response. The sensitivity analysis further demonstrates the model’s ability to identify influential parameters, providing insights into potential therapeutic targets.

8.2 Limitations

Like all mathematical models, certain simplifications and assumptions are necessary to make the model tractable. These include the use of constant parameters and the exclusion of certain biological processes, which might not fully capture the complexity of the immune response.The model primarily focuses on macrophages and a limited set of cytokines. While these components are crucial, the exclusion of other immune cells and signaling molecules may limit the model’s applicability to more complex or systemic infections.

8.3 Potential Extensions

8.3.1 Incorporation of Additional Immune Cells

The model can be extended to include other important immune cells, such as T cells, dendritic cells, or neutrophils. This would allow for a more comprehensive simulation of the immune response and provide insights into interactions between different cell types.

8.3.2 Adaptive Immune Response

Expanding the model to incorporate the adaptive immune response, including T-cell activation and antibody production, would offer a more complete picture of how the immune system combats Salmonella over a longer period.

8.3.3 Therapeutic Simulations

The model could be further developed to simulate the effects of various therapeutic interventions, such as cytokine inhibitors or macrophage-targeted therapies. This would provide a valuable tool for predicting treatment outcomes and optimizing therapeutic strategies.

9. Conclusion

This study presents a comprehensive mathematical model to simulate the immune system's response to Salmonella infection, integrating both spatial (PDE) and temporal (ODE) dynamics. The model successfully captures the intricate interactions between bacteria, macrophages, and key cytokines such as TNF-α, IL-6, IL-8, and IL-10, offering valuable insights into the mechanisms underlying immune regulation and pathogen clearance.

Through a series of numerical simulations, the model illustrates how Salmonella bacteria proliferate within the tissue, how macrophages become activated and respond to the infection, and how cytokines orchestrate the immune response. The sensitivity analysis further highlights the critical parameters that influence these dynamics, demonstrating the model’s predictive power and identifying potential targets for therapeutic intervention.

While the model incorporates certain simplifications and focuses on specific immune components, its strengths in biological relevance, modularity, and validation against experimental data underscore its utility as a research tool. The model also offers a flexible framework that can be extended to explore other pathogens, immune cells, and therapeutic strategies.

In conclusion, this model provides a powerful platform for advancing our understanding of the immune response to Salmonella and serves as a foundation for future studies aimed at optimizing treatment approaches and exploring the broader implications of immune regulation in infectious diseases.

9. Conclusion

·Talaei, K., Garan, S.A., Quintela, B.M., Olufsen, M.S., Cho, J., Jahansooz, J.R., Bhullar, P.K., Suen, E.K., Piszker, W.J., Martins, N.R.B., Moreira de Paula, M.A., dos Santos, R.W., & Lobosco, M. (2021). A Mathematical Model of the Dynamics of Cytokine Expression and Human Immune Cell Activation in Response to the Pathogen Staphylococcus aureus. Frontiers in Cellular and Infection Microbiology, 11, 711153.

·León-Triana, O., Sabir, S., Calvo, G.F., Belmonte-Beitia, J., Chulián, S., Martínez-Rubio, Á., Rosa, M., Pérez-Martínez, A., Ramirez-Orellana, M., & Pérez-García, V.M. (2021). CAR T cell therapy in B-cell acute lymphoblastic leukaemia: Insights from mathematical models. Communications in Nonlinear Science and Numerical Simulation, 94, 105570.

Drug Release

1. Introduction

In recent years, the development of drug delivery systems (DDS) has become a focal point in pharmaceutical research, particularly for improving the bioavailability and therapeutic efficacy of drugs. Polymeric nanoparticles, especially those based on Poly(lactic-co-glycolic acid) (PLGA), have garnered significant attention due to their biocompatibility, biodegradability, and versatility in drug encapsulation. PLGA has been extensively used as a carrier in drug delivery systems, allowing for controlled and sustained drug release, which is crucial for maintaining optimal therapeutic levels over extended periods.

One of the critical aspects of designing effective drug delivery systems is understanding and predicting the kinetics of drug release from the polymeric matrix. This involves the use of mathematical models that can accurately describe the drug release mechanisms and kinetics. Various models have been proposed to elucidate the drug release behavior from PLGA-based systems, each offering insights into different aspects of the release process.

In this study, we focus on six well-established mathematical models: the Higuchi model, Hixson-Crowell model, Square-root of Mass model, Three-second Root of Mass model, Weibull model, and Korsmeyer-Peppas model. These models were chosen based on their frequent application in previous studies and their ability to represent different mechanisms of drug release, such as diffusion, erosion, and swelling.

The Higuchi model assumes a drug release mechanism dominated by diffusion, where the amount of drug released is proportional to the square root of time. The Hixson-Crowell model describes drug release from systems where the surface area and radius of the particles change over time, often used when the drug release is controlled by dissolution or erosion processes. The Square-root of Mass and Three-second Root of Mass models, derived from the Hixson-Crowell model, offer additional perspectives on the relationship between drug release and time, particularly in systems with changing surface areas.

The Weibull model is an empirical model that is highly flexible and can describe various release profiles, including immediate and delayed release. This model is often used due to its versatility in fitting a wide range of drug release data. Finally, the Korsmeyer-Peppas model is a semi-empirical model that describes drug release behavior involving both diffusion and swelling mechanisms, making it applicable to a broad range of drug delivery systems.

This study aims to apply these models to the drug release kinetics from PLGA-based nanoparticles, utilizing parameters and datasets similar to those employed in previous studies. By comparing the fit and predictive capabilities of these models, we seek to identify the most appropriate mathematical framework for describing drug release from PLGA nanoparticles, thereby contributing to the optimization of drug delivery systems.

2. Model Assumptions

In the context of modeling drug release kinetics from PLGA-based nanoparticles, several assumptions are made to simplify the complex interactions and mechanisms involved. These assumptions are critical for applying the selected mathematical models and are outlined as follows:

A. Homogeneous Polymer Matrix: The PLGA matrix is assumed to be homogeneous, meaning that the distribution of the drug within the polymer is uniform.

B. Constant Diffusion Coefficient: The diffusion coefficient of the drug within the PLGA matrix is assumed to remain constant over time. This implies that the drug molecules move through the polymer at a steady rate, unaffected by changes in the polymer's physical or chemical properties during degradation.

C. No Initial Burst Release: It is assumed that there is no significant initial burst release of the drug. This means that the drug release begins gradually and follows the kinetics described by the models without a large initial release phase, which can occur due to surface-bound drug molecules.

D. Negligible Polymer Swelling: The PLGA nanoparticles are assumed not to swell significantly in the release medium. This assumption is important for models like the Higuchi and Korsmeyer-Peppas models, where swelling could otherwise alter the drug release kinetics by affecting the diffusion path.

E. Controlled Release Mechanism: The drug release is assumed to be controlled primarily by diffusion and/or erosion mechanisms. For the Hixson-Crowell and related mass-root models, it is assumed that the particle size decreases uniformly over time due to polymer erosion.

F. Isothermal Conditions: The release experiments are assumed to be conducted under constant temperature conditions, typically at physiological temperature (37°C). This assumption ensures that the diffusion coefficient and degradation rate of PLGA remain constant, simplifying the modeling process.

3. Parameters Extraction

The following table outlines the key parameters used in the six selected mathematical models for analyzing drug release kinetics from PLGA-based nanoparticles. The parameters, their symbols, units, and reference values are based on the original research and adapted for this study.

Table 1. Parameters settings

4. Mathematical Models

In this section, we present the mathematical models used to describe the drug release kinetics from PLGA-based nanoparticles. These models offer different perspectives on the release mechanisms, such as diffusion, erosion, and swelling.

4.1 Higuchi Model

Formula:

Description:

The Higuchi model is one of the most widely used models for describing drug release from matrix systems. It assumes that the drug release is governed by Fickian diffusion and that the matrix is homogeneous. The model indicates that the amount of drug released is proportional to the square root of time, reflecting a diffusion-controlled process. This model is particularly useful for systems where the drug is dispersed in an insoluble matrix.

·f(t): The cumulative fraction of drug released at time t.

·K_H: The Higuchi release constant, which reflects the drug’s diffusion rate in the matrix.

The Higuchi model suggests that drug release is driven by the concentration gradient, and as time progresses, the release rate decreases. This model is particularly suitable for understanding the initial phase of drug release from solid matrices.

4.2 Hixson-Crowell Model

Formula:

Description:

The Hixson-Crowell model accounts for changes in surface area and particle size as the drug dissolves. This model is often applied when the drug release is controlled by the erosion or dissolution of the matrix, rather than purely by diffusion. It assumes that the cube root of the remaining drug fraction decreases linearly over time.

·f(t): The cumulative fraction of drug released at time t.

·K_HC: The Hixson-Crowell constant, related to the dissolution or erosion rate.

This model is particularly relevant for systems where the drug particles undergo a reduction in size and surface area, such as tablets or spherical particles that erode over time.

4.3 Square-root of Mass Model

Formula:

Description:

The Square-root of Mass model is a variant of the Hixson-Crowell model, emphasizing the relationship between the square root of the remaining mass and time. This model is particularly useful for systems where drug release follows a square-root time dependence, often due to the combined effects of diffusion and erosion.

·f(t): The cumulative fraction of drug released at time t.

·K_SM: The release constant for the square-root of mass model, indicative of the rate at which the remaining drug mass decreases.

This model is beneficial for systems where drug release slows down over time as the mass and surface area of the drug decrease.

4.4 Three-second Root of Mass Model

Formula:

Description:

This model, like the Square-root of Mass model, is derived from the Hixson-Crowell model but focuses on the three-second root of the remaining drug mass. It provides a different perspective on the erosion or dissolution-controlled drug release.

·f(t): The cumulative fraction of drug released at time t.

·K_3SM: The release constant for the three-second root of mass model, which indicates the rate of drug release based on the three-second root of the remaining drug mass.

This model is particularly relevant for more complex erosion processes where the drug release is not linearly related to time but follows a cubic relationship.

4.5 Weibull Model

Formula:

Description:

The Weibull model is an empirical model that is highly flexible and can describe a variety of drug release profiles, including immediate, delayed, and sigmoidal release patterns. The model’s parameters allow for a detailed description of the shape and timing of the release curve.

·f(t): The cumulative fraction of drug released at time t.

·t_d: The lag time before the onset of drug release.

·τ: A scale parameter related to the characteristic time of release.

·β: A shape parameter that indicates the nature of the release curve (exponential, sigmoidal, etc.).

The Weibull model’s flexibility makes it applicable to a wide range of drug delivery systems, making it particularly useful for fitting complex release data where other models may not provide an adequate description.

4.6 Korsmeyer-Peppas Model

Formula:

Description:

The Korsmeyer-Peppas model is a semi-empirical model that describes drug release from polymeric systems, where the release mechanism may involve both diffusion and swelling. The release exponent nnn provides insight into the type of release mechanism: n=0.5n = 0.5n=0.5 indicates Fickian diffusion, while 0.5 < n < 10.5 < n < 10.5 < n < 1 suggests anomalous transport (combination of diffusion and swelling).

·f(t): The cumulative fraction of drug released at time t.

·K_KP: The release rate constant for the Korsmeyer-Peppas model.

·n: The release exponent, indicating the mechanism of drug release.

This model is particularly useful for characterizing the release from swellable matrices where the release mechanism is not purely diffusion but also involves matrix relaxation or erosion.

5. Results

The simulation results provide a comprehensive overview of the drug release kinetics from PLGA-based nanoparticles, as predicted by the six selected mathematical models. The plotted curves illustrate the differences in release profiles, allowing for a detailed comparison and analysis of the models.

The Higuchi model shows a characteristic curve where the drug release is proportional to the square root of time. This model suggests that the release rate decreases over time, which is consistent with the diffusion-controlled release from a homogeneous matrix. The curve starts with a steep slope, indicating a higher release rate at the beginning, which gradually tapers off as time progresses.

The Hixson-Crowell model produces a release curve that assumes a linear reduction in particle size and surface area over time. The model predicts a somewhat linear relationship between the cube root of the remaining drug and time. This model’s release profile is less steep initially compared to the Higuchi model, indicating a slower initial release rate that becomes more consistent over time.

Similar to the Hixson-Crowell model, the Square-root of Mass model suggests a square root dependency on time. The curve is very close to that of the Hixson-Crowell model but slightly deviates due to the different power dependency. The initial release rate is moderate, with a steady release observed over time, indicating that this model is also suited for describing systems with combined diffusion and erosion mechanisms.

The Three-second Root of Mass model is closely related to the Square-root of Mass model but with a cubic root dependency. The curve generated by this model follows a similar trajectory to the Square-root of Mass model but shows an even slower release rate over time, especially at the later stages. This model’s slight variation from the previous two emphasizes the effect of different exponent values on the release kinetics, demonstrating that small changes in the model's structure can lead to noticeable differences in predictions.

The Weibull model’s flexibility is evident in its ability to describe a wide range of release profiles. The curve in this scenario shows an initial lag phase (due to t_d) followed by a sigmoidal release pattern. The shape parameter β influences the curvature, making this model particularly versatile for fitting various experimental data. The model’s ability to capture both immediate and delayed release phases makes it a powerful tool for characterizing complex drug release behaviors.

The Korsmeyer-Peppas model, with its power-law form, displays a release curve that depends on the exponent n. In this case, the curve suggests a non-linear release profile indicative of anomalous transport, which could involve a combination of diffusion and swelling mechanisms. The release rate in this model is higher initially and slows down over time, similar to the Higuchi model but with greater flexibility to describe non-Fickian transport.

When comparing the release profiles, it is clear that the Higuchi and Korsmeyer-Peppas models predict a relatively rapid initial release followed by a slower, sustained release. This is consistent with diffusion-dominated processes in polymer matrices. The Hixson-Crowell, Square-root of Mass, and Three-second Root of Mass models, on the other hand, predict more consistent and moderate release rates, which are characteristic of systems where erosion or surface reduction plays a significant role.

The Weibull model stands out for its flexibility, effectively capturing the initial delay in drug release and providing a good fit for both rapid and gradual release phases. This makes the Weibull model particularly useful when the release mechanism is not well-defined or when multiple mechanisms contribute to the release kinetics.

The comparison between the Square-root of Mass model and the Three-second Root of Mass model shows that the change in the exponent from 1/2 to 2/3 results in a subtle but significant difference in the release profiles. The latter model predicts a slower release, emphasizing the sensitivity of the release kinetics to the specific mathematical form used.

The simulation results highlight the strengths and limitations of each model in describing the drug release from PLGA-based nanoparticles. The Higuchi and Korsmeyer-Peppas models are well-suited for systems where diffusion is the primary release mechanism, while the Hixson-Crowell and mass-root models are better for describing erosion-controlled release. The Weibull model’s versatility makes it a strong candidate for fitting experimental data across various conditions.

Ultimately, the choice of model depends on the specific characteristics of the drug delivery system being studied. Understanding the nuances of each model allows researchers to select the most appropriate tool for predicting drug release behavior, which is crucial for the design and optimization of effective drug delivery systems.

Figure 1. Release profiles evaluation for varieties of model

6. Discussion and Future Perspectives

The mathematical modeling of drug release kinetics from PLGA-based nanoparticles, as explored through the six models in this study, provides a solid foundation for understanding and predicting the behavior of drug delivery systems. However, there are several avenues for extending and improving these models to enhance their applicability and accuracy in real-world scenarios.

6.1 Incorporating Multiphase Release Mechanisms

The current models typically focus on single-phase release mechanisms, such as diffusion or erosion. However, in many real-world drug delivery systems, multiple release phases are present, including initial burst release, sustained release, and final tail release. To improve the accuracy of predictions, future models could incorporate multiphase mechanisms, combining elements from different models to better reflect the complexity of drug release in heterogeneous systems.

6.2 Nanoparticle Interaction Models

As drug delivery systems become increasingly sophisticated, with multifunctional nanoparticles that can target specific cells, respond to external stimuli, and release multiple drugs, the corresponding models must evolve to capture these complexities. Future models should consider not only the drug release kinetics but also the interactions between nanoparticles and their biological environment, including cell uptake, endosomal escape, and intracellular trafficking.

6.3 Integration with Pharmacokinetic-Pharmacodynamic (PK-PD) Models

For a comprehensive understanding of drug efficacy, it is essential to link drug release models with pharmacokinetic-pharmacodynamic (PK-PD) models. This integration would allow researchers to predict not only how the drug is released but also how it is absorbed, distributed, metabolized, and excreted in the body, and ultimately, its therapeutic effect. Developing such integrated models would provide a powerful tool for drug development and clinical decision-making.

7. References

· Pourtalebi Jahromi, L., Ghazali, M., Ashrafi, H., & Azadi, A. (2020). A comparison of models for the analysis of the kinetics of drug release from PLGA-based nanoparticles. Heliyon, 6(3), e03451.

· Gutiérrez-Valenzuela, C.A., Guerrero-Germán, P., Tejeda-Mansir, A., Esquivel, R., Guzmán-Z, R., & Lucero-Acuña, A. (2016). Folate Functionalized PLGA Nanoparticles Loaded with Plasmid pVAX1-NH36: Mathematical Analysis of Release. Applied Sciences, 6(12), 364.

· D.N. Kapoor, A. Bhatia, R. Kaur, R. Sharma, G. Kaur, S. Dhawan, PLGA: a unique polymer for drug delivery, Ther. Deliv. 6 (1) (2015) 41–58.

Delayed Lysis

1. Introduction

Salmonella enterica, a significant pathogen, is responsible for a wide range of infections in humans and animals, often leading to severe gastrointestinal diseases. Understanding the intracellular dynamics of Salmonella is crucial for developing more effective treatment strategies, particularly in the context of its ability to evade the host immune response and establish persistent infections. One of the critical processes in the intracellular lifecycle of Salmonella is the delayed lysis of host cells, a mechanism that allows the bacteria to replicate within the host cell before causing cell death and subsequently spreading to new cells. This delayed lysis is a complex process influenced by various factors, including nutrient availability and the bacterium's growth rate within the host cell.

The delayed lysis model for Salmonella within host cells is particularly relevant for understanding the pathogen's replication dynamics and its interaction with the host environment. This model allows us to quantify the growth of Salmonella inside host cells and to predict when lysis is likely to occur, thereby offering insights into how the timing of this event influences infection outcomes. By incorporating key biological parameters, such as the growth rate of Salmonella and the availability of nutrients within the host cell, this model provides a comprehensive framework for analyzing the delayed lysis process.

2. Model Assumptions

To accurately model the intracellular behavior of Salmonella enterica and its delayed lysis within host cells, several key assumptions are made:

A. Neglect of External Factors: The model assumes that external factors, such as the host immune response, do not significantly influence the intracellular growth of Salmonella enterica or the timing of lysis. This simplification allows the focus to remain on the intrinsic dynamics of bacterial growth and lysis.

B. Fixed Delay Time for Lysis: A fixed delay time (τ) is assumed before lysis is triggered. During this delay, Salmonella enterica continues to grow until a critical density is reached, after which the host cell lysis occurs. This delay represents the time required for the bacteria to reach sufficient numbers to initiate the lysis process.

C. Constant Lysis Rate: Once lysis is triggered, the rate at which host cells undergo lysis (δ) is constant. This means that, after the delay period, the lysis of host cells proceeds at a steady rate, regardless of the bacterial population size at that moment.

These assumptions provide a simplified yet powerful framework for understanding the key processes governing the delayed lysis of host cells by Salmonella enterica.

3. Parameter Settings

The parameters have been selected based on experimental data and relevant literature, providing a realistic foundation for the model. The values are adjustable depending on specific experimental conditions or host cell types

Table 1. Parameters settings

4. Mathematical Model

We develop a mathematical model to describe the intracellular growth dynamics of Salmonella enterica and the delayed lysis of the host cell. This model is based on ordinary differential equations (ODEs) that capture the essential features of bacterial replication and the timing of host cell lysis. The model is built on two key differential equations that describe the growth of Salmonella enterica within the host cell and the delayed lysis of the host cell.

4.1 Bacterial Growth Equation

·S: The population of Salmonella enterica within the host cell.

·μ: The maximum growth rate of Salmonella enterica .

·C: The concentration of nutrients available to Salmonella enterica within the host cell .

·K: The half-saturation constant, representing the nutrient concentration at which the bacterial growth rate is half of its maximum value .

·δ: The lysis rate constant , representing the rate at which host cells begin to undergo lysis as the bacterial population grows.

The term S∙μ∙C/(C+K) describes the logistic growth of Salmonella enterica within the host cell, where bacterial growth is constrained by the availability of nutrients. The term δS counts for the eventual lysis of the host cell, which reduces the bacterial population once lysis begins. The equation captures the dynamics of bacterial growth within the host cell, which initially follows a logistic pattern as nutrients are consumed. As Salmonella enterica replicates, the available nutrients within the host cell are gradually depleted, leading to a slowing of the growth rate.

4.2 Delayed Lysis Equation

This equation models the delayed lysis of the host cell.

·L: Represents the process of lysis of the host cell.

·δ: The same lysis rate constant as in the bacterial growth equation.

·τ: The delay time, representing the period during which Salmonella enterica continues to grow within the host cell before initiating lysis.

In this equation, S(t-τ) represents the bacterial population at a time τ hours earlier. The delayed lysis occurs after this fixed time delay τ, reflecting the period required for the bacterial population to reach a critical density capable of triggering the lysis of the host cell.

The overall aim of the ODE model is to provide a quantitative framework for understanding the intracellular lifecycle of Salmonella enterica. By simulating the delayed lysis process, we can predict how changes in key parameters, such as the growth rate or the delay time, influence the timing and extent of host cell destruction.

5. Results and Conclusion

The model developed to describe the intracellular dynamics of Salmonella enterica and the delayed lysis of host cells was implemented and solved. The results of this simulation are presented in two plots: the growth of Salmonella enterica within the host cell over time, and the delayed lysis process that occurs after a specified delay period.

The first plot illustrates the growth of Salmonella enterica within the host cell over a 20-hour period. Initially, the bacterial population increases rapidly, following a logistic growth pattern. This exponential phase is driven by the high availability of intracellular nutrients and the intrinsic replication rate of the bacteria. As time progresses, the growth rate begins to decelerate as the nutrients within the host cell are depleted, causing the population to approach a saturation point.

This deceleration is reflected in the gradual flattening of the growth curve. The model assumes that nutrient availability is the primary limiting factor for bacterial growth, which aligns with biological observations of Salmonella replication within host cells. The impact of the nutrient limitation is evident as the curve begins to plateau, indicating that the bacterial population is reaching its carrying capacity within the host cell.

The second plot depicts the delayed lysis process, which is triggered after a fixed delay time (τ) of 3 hours. During this delay period, the bacterial population continues to grow undisturbed, accumulating within the host cell. Once the delay time is reached, lysis begins, and the rate of lysis is determined by the lysis rate constant (δ).

The lysis process is characterized by a sharp decline in the bacterial population, as the host cell ruptures and releases the bacteria. This sudden drop in population size reflects the catastrophic impact of lysis on the intracellular bacteria, as well as the host cell. The timing of the lysis event is critical; if the delay time is too short, the bacterial population may not reach a sufficiently high density to ensure effective dissemination. Conversely, a longer delay allows the population to grow larger, potentially leading to more severe infection outcomes as more bacteria are released during lysis.

The simulation results demonstrate the interplay between bacterial growth and delayed lysis within the host cell. The logistic growth model effectively captures the initial expansion of the Salmonella enterica population, while the delayed lysis model introduces a critical time-dependent factor that influences the ultimate fate of the host cell.

One of the key findings from this analysis is the importance of the delay time (τ) in determining the infection dynamics. The delay allows the bacterial population to build up to a level that maximizes the impact of lysis, facilitating the spread of the infection to neighboring cells. This mechanism is likely an evolutionary adaptation that Salmonella enterica uses to optimize its survival and transmission within the host. Moreover, the lysis rate constant (δ) plays a significant role in determining the rate at which the host cell is destroyed once lysis is initiated. The results of this mathematical modeling study provide valuable insights into the intracellular dynamics of Salmonella enterica and the delayed lysis of host cells.

Figure 1. Growth-lysis comparison curves

6. References

•Fazzino, L., Anisman, J., Chacón, J. M., Heineman, R. H., & Harcombe, W. R. (2020). Lytic bacteriophage have diverse indirect effects in a synthetic cross-feeding community. The ISME Journal, 14(1), 123-134.

•Deatherage DE, Barrick JE. Identification of mutations in laboratory evolved microbes from next-generation sequencing data using breseq. Methods Mol Biol 2014;1151:165–88.