Modeling & Analysis

We compiled anticancer peptide drugs that are approved or in clinical development from public databases.

Data sources:

1. ClinicalTrials.gov: The world's largest registry of clinical trials, covering the vast majority of publicly developed drug studies.
2. THPdb: A curated database of peptide therapeutics that are approved or in clinical phases.

Procedure:

We retrieved records using the keywords "anticancer" AND "peptide drug," removed duplicates, and organized the resulting entries by drug class.

Figure 1.1. Proportions of mainstream classes of peptide-based anticancer therapeutics

From the figure, antibody modalities account for the largest share (red, blue, and purple segments), indicating that industry practice predominantly relies on antibodies as the principal therapeutic form. Although these peptide-related agents offer high target specificity, they typically entail very high R&D costs, narrow applicability, and limited tissue penetration. For example, trastuzumab is a monoclonal antibody that binds HER2 and suppresses tumor growth via pathway blockade; in colorectal cancer (CRC), only ~5% of patients are eligible. Other peptide-related modalities—such as colony-stimulating factors and antibody–drug conjugates (ADCs)—carry notable adverse-event profiles and potential misuse risks. Consequently, a delivery technology that helps peptide drugs reach targets efficiently while maintaining low toxicity during transit is crucial for unlocking their anticancer potential.

In response, our project aims to develop a delivery vehicle with high specificity, low toxicity, and high bioavailability to guide existing peptide drugs to intracellular targets.

1. markedly improve target specificity while reducing toxicity.
2. help peptide drugs reach their targets and increase bioavailability.
3. tailor to colorectal cancer and enable noninvasive delivery.

Beyond statistical summaries from clinical drug databases, we leveraged GEO transcriptomic datasets to perform differential expression and pathway enrichment between colorectal cancer and normal tissues. Results showed:

1. Significant up-regulation of MMP genes (MMP1, MMP3, MMP7, MMP9; p < 0.0001), implicating key roles in ECM remodeling and invasion;
2. Prominent enrichment of ECM-receptor interaction, cell adhesion molecules (CAMs), and focal adhesion pathways—highly concordant with MMP function.

Grounded in the literature, we executed the following workflow:

1. Data acquisition and cohort stratification: We obtained transcriptomes for colorectal cancer and normal intestinal tissues (training set GSE10950). Clustering analyses confirmed clear separation of tumor vs. normal expression profiles, providing phenotypic support for downstream analyses.
2. Differentially expressed genes (DEGs) and necroptosis-related DEGs (NRDEGs): Using statistical criteria, we identified DEGs between tumor and normal tissues, with specific focus on MMP2 and MMP9. Volcano and violin plots provided complementary views to highlight markedly overexpressed proto-oncogenic candidates.
3. Gene set enrichment analysis (GSEA): Moving beyond single genes, we profiled pathway-level enrichments for ECM-receptor interaction, CAMs, and focal adhesion in tumor tissues. We further evaluated upstream regulation—e.g., TGF-β and VEGF—on MMP expression, assembling a pro-oncogenic signaling network.
4. Protein-protein interaction (PPI) network and core-module mining: We built a PPI network from DEGs to pinpoint modules coordinating with MMPs upstream and downstream, revealing putative functional partners and core gene clusters implicated in tumorigenesis.
5. Validation cohort analysis: We validated core DEGs and pathway signatures in an independent dataset ("Transcriptomic and Cellular Content Analysis of Colorectal Cancer by Combining Multiple Independent Cohorts"). Diagnostic performance of shortlisted genes was assessed to confirm robustness and generalizability.

Objective:

The GEO repository is the world's largest public gene-expression database. Leveraging published transcriptomes avoids redundant experiments, reduces costs, and ensures reproducibility. We selected GSE10950 as the primary training cohort (owing to its larger sample size) and used the supplementary data from "Transcriptomic and Cellular Content Analysis of Colorectal Cancer by Combining Multiple Independent Cohorts" as an independent validation cohort to test robustness and mitigate single-dataset bias.

Procedure:

We retrieved GSE10950 (tumor vs. normal colorectal tissue expression profiles) from GEO as the training set and downloaded the supplementary files of the study above as the validation set. Data preprocessing and curation were conducted with UltraEdit and Excel.

Objective:

• Identify DEGs that differ significantly between colorectal cancer and normal tissues; such genes can serve as molecular markers of disease phenotypes and may drive tumorigenesis.
• Focus on the matrix metalloproteinase (MMP) family to narrow from genome-wide signals to genes directly tied to our core biology (ECM remodeling, aberrant adhesion), thereby filtering noise and sharpening the target space.
• Utilize visualization tools in Sangerbox (http://sangerbox.com/)—heatmaps, volcano plots, violin plots—to depict DEG counts and expression patterns, laying the groundwork for downstream GSEA and PPI network modeling.

Procedure:

We uploaded tumor and normal CRC datasets to Sangerbox and applied default criteria (FDR/adjusted p < 0.05 with |log2FC| thresholds) to call significant DEGs and annotate NRDEGs (necroptosis-related DEGs). Particular attention was given to MMP2 and MMP9 and other MMP family members. Using volcano plots, heatmaps, violin plots, and ROC curves, we characterized their tumor-biased overexpression and nominated putative pro-oncogenic candidates.

Figure 2.1. Differential expression heatmap

Caption: Blue = Tumor; Red = Normal.

Sample grouping	Left/top annotation bars denote sample classes (Tumor vs. Normal)	Hierarchical clustering shows clear separation
Gene expression	The heatmap color scale encodes expression levels: blue = low, red = high	A clear "tumor vs. normal" differential pattern is observed

Table 2.1. Notes on the differential-expression heatmap

Sample stratification is clear, with a marked tumor--normal separation. In the dendrogram, tumor (blue bar) and normal (red bar) samples are fully segregated, indicating a substantial divergence in their gene-expression programs.

Biological implication: Colorectal cancer exhibits expression signatures that are distinctly different from normal intestinal tissue, providing phenotypic support for downstream differential analyses.

The figure also reveals "tumor high-/low-expression gene clusters." The heatmap shows four expression modules (upper-left red, upper-right blue, lower-left blue, lower-right red): the upper-left represents a tumor high-expression cluster (red blocks enriched), whose genes may be pro-oncogenic (e.g., the matrix metalloproteinase (MMP) family and cell-proliferation-related genes). The lower-left represents a tumor low-expression cluster (blue blocks enriched) → these genes may be tumor suppressors or normal tissue-specific genes.

Figure 2.2. Volcano plot of differential expression

Axis/Feature	Description	Observation
X-axis (−log10 (p-value))	Statistical significance of differential expression: larger values (further right) indicate stronger significance (smaller p-values)	Range 0–26, with most genes at p < 0.05
Y-axis (log2 fold change)	Magnitude of expression change: larger values (higher up) denote greater tumor-vs-normal differences	Range −8 to 8, with red (upregulated) points clustering in the upper region
Color (direction of change)	Red = tumor high expression (up-regulated); green = tumor low expression (down-regulated). Red points outnumber green ones	Many matrix metalloproteinase (MMP) family genes appear in red

Table 2.2. Notes on the volcano plot

Key findings:

• In the volcano plot, MMP-9 appears as a significant red point, indicating marked up-regulation in tumors with strong statistical support (high −log10(p-value)).
• By contrast, MMP-2 appears in black, not passing the preset significance thresholds (e.g., combined FDR/adjusted p and |log2FC| criteria), indicating that MMP-2 does not reach statistical significance in this training cohort.

Figure 2.3. Violin plots for MMP9 vs. MMP2 enrichment/expression

Figure 2.4. Box plots for MMP9 vs. MMP2 enrichment/expression (left: MMP9; right: MMP2)

Plot Type	Feature Description	Observation
Violin shape	Distribution density of gene expression: the wider the violin, the greater the sample density at that expression level	Here, the tumor group (blue) shows generally wider violins, indicating a larger fraction of samples at higher expression ranges and greater heterogeneity in tumors
Box plot	Statistical features of expression: the box denotes the interquartile range (IQR, 25%–75%), the horizontal line inside marks the median, and whiskers span 1.5×IQR (points beyond are outliers)	In this figure, tumor box plots sit higher than normal

Table 2.3. Notes on violin and box plots

Figure 2.5. Gene enrichment analysis plot for matrix metalloproteinase 9 (MMP9) vs matrix metalloproteinase 2 (MMP2) (left: MMP9; right: MMP2)

The AUC of matrix metalloproteinase 9 (MMP9) is around 0.82, indicating potential as a CRC molecular biomarker and as the targeting moiety for the TRACER-sensitive peptide.

Objective:

Unlike single-gene differential analysis, GSEA evaluates predefined gene sets (e.g., pathways, biological processes) to assess overall expression trends in tumor vs. normal groups and determine whether pathways are activated or suppressed.

The core goal here is to verify whether pathways associated with high expression of matrix metalloproteinases (MMPs)—such as ECM remodeling and aberrant cell adhesion—are significantly enriched in colorectal cancer. Results show these pathways are markedly activated in tumor samples, supporting the hypothesis that the matrix metalloproteinase (MMP) family participates in ECM remodeling, adhesion abnormalities, and their upstream regulatory signals.

We also observed enrichment of regulatory pathways such as TGF-β and VEGF, providing clues for subsequent mechanistic analyses (e.g., links between ECM remodeling and signaling pathways).

Procedure:

Using the Sangerbox web platform, we performed GSEA on the uploaded differentially expressed genes. We selected the c2 (curated pathway) collection from MSigDB, ran the analysis, and obtained enriched pathways and biological processes. We focused on core MMP-related pathways—ECM-receptor interaction, cell adhesion molecules (CAMs), and focal adhesion—and concurrently evaluated enrichment of upstream signaling pathways such as TGF-β and VEGF.

Figure 2.6. KEGG pathway enrichment bubble chart

Caption: Pathway names—e.g., Human T-cell leukemia virus 1 infection, Axon guidance—are on the y-axis (enriched KEGG pathways).

• Metrics: X-axis GeneRatio = the proportion of differentially expressed genes (DEGs) in a pathway relative to all DEGs (larger values indicate a stronger association between the pathway and the DEGs).
• Bubble size (Count): Number of DEGs contained in the pathway (larger bubbles indicate broader pathway impact).
• Color (−log10(p-value)): Statistical significance of pathway enrichment (darker colors indicate smaller p-values and more significant enrichment).

Below is the analysis of specific KEGG pathways:

Figure 2.7. Detailed signaling pathway analysis

Pathway name	ES (enrichment score)	NP (nominal P value)	Interpretation	Mechanistic link
ECM RECEPTOR INTERACTION	ES = 0.4557	NP = 0.0978	Pathway enriched in tumors (positive ES, right-shifted curve), near significance	Matrix metalloproteinase (MMP) activity degrades ECM components (e.g., collagen, fibronectin); as this pathway underpins cell--ECM adhesion, MMP function directly impacts its activity
CELL ADHESION MOLECULES	ES = 0.4368	NP = 0.0432	Pathway significantly enriched in tumors (positive ES, p < 0.05)	Focal adhesions are key for cell--ECM linkage; ECM degradation by matrix metalloproteinases (MMPs) disrupts focal adhesions, and pathway activation can in turn regulate MMP expression
FOCAL ADHESION	ES = 0.4168	NP = 0.0039	Pathway significantly enriched in tumors (positive ES, p < 0.01)	Cell adhesion molecules (CAMs; e.g., integrins) act via the ECM; ECM degradation by matrix metalloproteinases (MMPs) alters CAM function, with the pathway and MMPs co-regulating cell migration
VEGF SIGNALING PATHWAY	ES = 0.2971	NP = 0.1381	Pathway enriched in tumors (positive ES, right-shifted curve), near significance	The VEGF pathway promotes angiogenesis; matrix metalloproteinases (MMPs) degrade vascular basement membrane (facilitating neovascularization), and VEGF can upregulate MMP gene expression
PATHWAYS IN CANCER	ES = −0.2050	NP = 0.3496	Pathway enriched in normals (negative ES, left-shifted curve), not significant	"Pathways in cancer" provides an overview; high matrix metalloproteinase (MMP) expression is one hallmark of pathway activation, contributing to ECM degradation and angiogenesis
WNT SIGNALING PATHWAY	ES = 0.2393	NP = 0.2373	Pathway enriched in tumors (positive ES, right-shifted curve), not significant	The TGF-β pathway regulates the balance of ECM synthesis and degradation; matrix metalloproteinases (MMPs) act as downstream effectors (e.g., TGF-β upregulating MMP2)
TGF-β SIGNALING PATHWAY	ES = 0.3194	NP = 0.1357	Pathway enriched on the tumor side but not strictly significant (p ≈ 0.1; suggests association with tumor activation, pending validation)	Activation of the Wnt pathway promotes tumor cell proliferation and migration; matrix metalloproteinases (MMPs) degrade the ECM to aid invasion, and the pathway can directly regulate MMP gene expression (e.g., Wnt upregulating MMP7)

Table 2.4. Notes on specific signaling pathways

Conclusion:

• ECM-receptor interaction, cell adhesion molecules (CAMs), and focal adhesion pathways are significantly enriched in tumors, representing the biological context of MMP functions.
• TGF-β and VEGF pathways show trend-level enrichment on the tumor side and can directly upregulate MMP gene expression (e.g., TGF-β promoting transcription of MMP2 and MMP9 via SMAD signaling).

Value: Extending from "single genes" to "pathway networks" reveals the upstream regulatory routes underlying high expression of MMPs, corroborating that MMP levels are markedly elevated in colorectal cancer tissues compared with normal tissues.

Objective:

Proteins form functional networks through their interactions, and a PPI network can reveal cooperative relationships among differentially expressed genes. The network constructed in this study shows that matrix metalloproteinase 9 (MMP9), MMP1, MMP3, MMP7, MMP10, MMP2 and other matrix metalloproteinase (MMP) family genes are highly interconnected with collagen family members (e.g., COL1A1, COL4A1, COL6A3) and integrins (the ITGA series), suggesting cooperative roles in ECM remodeling and aberrant cell adhesion.

By screening hub genes (e.g., matrix metalloproteinase 9 (MMP9), COL1A1, ITGA5), key nodes in the network can be pinpointed; dysfunction of these genes may be decisive for stromal remodeling, tumor invasion, and metastasis in colorectal cancer, providing a prioritized list for subsequent experimental validation and target studies.

Module analysis (e.g., modules enriched with multiple collagens and integrins) helps focus on the core biological unit of ECM structural remodeling and extracellular matrix signaling, implying that these genes may act in concert to regulate high matrix metalloproteinase (MMP) expression and ECM abnormalities, thereby guiding deeper mechanistic investigation.

Procedure:

Construct a PPI network of differentially expressed genes—focusing on the matrix metalloproteinase (MMP) family and their upstream/downstream ECM-related genes—using the STRING database.

Figure 2.8. PPI analysis graph

Legend: Green = gene neighborhood/co-occurrence; Red = gene fusion; Blue = co-expression; Purple = experimental evidence; Light blue = curated database associations; Yellow = text mining; Black = co-expression or mutual regulation; Pink = protein homology.

Nodes with more edge colors and denser connections indicate proteins playing more central roles in the network.

The PPI network shows that the matrix metalloproteinase (MMP) family (e.g., MMP9, MMP1, MMP3) forms a highly interactive network with collagen family members (COL1A1, COL6A3, COL4A1) and integrins (ITGA series), suggesting cooperative regulation in extracellular matrix (ECM) remodeling, cell adhesion, and tumor invasion/metastasis.

Among them, MMP9, COL1A1, and ITGA5 may serve as potential hub genes, offering candidate directions for further mechanistic studies and targeted interventions. However, this is not the main focus of our project, so we did not pursue it further.

Bioinformatics findings may be influenced by dataset-specific factors; consistency in an independent validation cohort can markedly improve the reliability of conclusions, demonstrating that the matrix metalloproteinase 9 (MMP9) gene is significantly differentially expressed in colorectal cancer, rather than being an artifact of a single dataset.

We used the independent validation dataset "Transcriptomic and Cellular Content Analysis of Colorectal Cancer by Combining Multiple Independent Cohorts" to externally validate the core differentially expressed genes and pathways identified in the training set.

In this validation cohort, a Venn diagram was used to intersect differentially expressed genes (DEGs) from the training set with those from the validation set; matrix metalloproteinase (MMP) family genes—specifically MMP9—consistently appeared in the intersection, suggesting concordant expression patterns across cohorts and potential biological significance. A volcano plot displayed the overall distribution of gene expression between tumor and normal tissues, visually showing the numbers and significance of up- and down-regulated genes, thereby aiding assessment of genome-wide differential trends within the validation set.

This "discovery in the training set—replication in the validation set" strategy, via multidimensional visualization with Venn and volcano plots, not only intuitively presents the expression and significance of core genes in the validation cohort, but also objectively evaluates their diagnostic potential and clinical applicability, further strengthening the robustness and generalizability of the conclusions.

Figure 2.9. Volcano plot of the validation cohort

Figure 2.10. Venn diagram of differentially expressed genes in the training vs. validation cohorts

Caption: The matrix metalloproteinase 9 (MMP9) gene appears in the intersection of the two sets.

These results confirm that MMP9 is a core gene markedly activated in colorectal cancer and reveal its associated pathways and immune microenvironment features. We will translate these findings into modeling parameters: using MMP9 expression level as the activation probability and the ECM degradation rate as an input to the tissue model, thereby making subsequent structural modeling and pharmacokinetic simulations more biologically realistic.

After validating TRACER's feasibility at the molecular level, we face a practical engineering question: how to produce this therapeutic efficiently and robustly? Our strategy employs a quorum-sensing triggered synchronized lysis circuit (SLC) in engineered bacteria. This is a complex dynamical system whose behavior (e.g., lysis period and efficiency) depends strongly on strain traits and environment.

Accordingly, the goal here shifts from "validation" to "prediction and guidance." We construct an ordinary differential equation (ODE) model to address key questions and provide a clear "operating manual" for the wet lab:

• How will the system behave? Under different strain backgrounds, will synchronized lysis be sustained or damp quickly?
• How should experiments be run? What sampling frequency captures the critical dynamics?
• How should results be interpreted? How do we decide from data whether the circuit behaves as intended?

To capture the core dynamics of SLC, we formulate an ODE system with four key variables: cell density (N), AHL concentration (H), lysis protein concentration (L), and LuxI synthase concentration (I). Based on literature review, the relevant processes can be described by ODEs as follows:

Meaning:

• N: Cell number (approximated by OD or CFU).
• H: Concentration of AHL.
• L: Lysis protein (readouts via SDS-PAGE/immunoassay or fluorescent tags).
• I: LuxI synthase.
• α_H: Promoter activity activated by AHL (Hill model with coefficient 4 to capture cooperativity).
• μ_L: Cell death rate mediated by the lysis protein (Hill model).

Because circuit behavior may vary across chassis strains, we avoid a single parameter set and instead use two literature-informed representative scenarios to probe behavioral bounds.

We identified two parameter sets from the literature:

Set A (literature E. coli-like) ^[16]

Set A (literature E. coli-like) models efficient AHL synthesis and clearance. It is expected to yield sustained/stable oscillations—an ideal, periodically "harvestable" production mode.

Parameter	Value	Parameter	Value
μG	0.2	N0	20
k	10	L0	1
n	2	b	15
μ	12	CL	0.5
CI	1	α0	1
αH	20	H0	4.5
γL	2	γI	2
γC	12

Expected phenomenon: quorum-synchronized lysis with oscillations persisting over multiple cycles.

Parameter set B (S. Typhimurium-like in the literature) ^[17]

Parameter set B (S. Typhimurium-like in the literature) models a weaker AHL signaling system. It is expected to yield damped oscillations—i.e., oscillations that decay to a steady state after several cycles—representing a plausible but less efficient mode.

Parameter	Value	Parameter	Value
μG	0.2	N0	20
k	4	L0	1
n	4	b	15
μ	1	CL	0.5
CI	1	α0	0.4
αH	8	H0	2
γL	2	γI	1
γC	12

Expected behavior: With weaker synthesis/activation and a steeper lysis threshold (n↑, k↓), oscillations are strongly damped and approach steady state.

(1) Population rhythms of cells and AHL

Figure 6.1-A (Set A): N shows clear periodic rise-fall; each surge of H increases promoter activity, then L rises and the death rate increases, triggering synchronized lysis; after N drops, H is cleared and the next cycle begins.

Interpretation: A classic self-excited oscillator; OD/CFU should show >2 cycles of drop peaks, phase-led by AHL peaks in LC-MS (AHL peak precedes lysis peak).

Figure 6.1-B (Set B): After the first one or two peaks, amplitudes successively shrink and converge to steady state.

Interpretation: Damped oscillations, suggesting weaker AHL induction or a higher lysis threshold in this chassis.

These simulations translate into three practical guidelines:

• Set expectations: Damped oscillations need not mean failure; they may reflect chassis traits (à la Set B), guiding host choice/engineering.
• Optimize sampling: With a ~3-4 h period, measure OD600 every 15-20 min; densify sampling near peaks/valleys for AHL and protein to capture dynamics.
• Define quantitative calls:
  • If over 3 cycles ΔN/N_peak > 30% and amplitude decay < 20% → "sustained oscillation" (Set A-like).
  • If amplitude decays > 50% within 3 cycles → "damped" (Set B-like).

(2) Phase relation of L and LuxI/II

Figure 6.2-A (Set A): I rises before L (driven by promoter activity), L peaks ~half a cycle later; each L peak coincides with increased death rate and a sharp N drop.

Interpretation: Clear phase ordering.

Figure 6.2-B (Set B): The peak value decreases gradually and the peak-to-peak interval slightly lengthens.

Unlike the continuous oscillation in Figure A, the oscillation in this scenario exhibits significant damping attenuation, with the cracking peak decreasing gradually and eventually stabilizing. This indicates that the strength of the positive feedback loop driving cracking is insufficient to maintain sustained synchronous behavior.

According to the analysis of model parameters, the two main inhibitory factors that may cause this phenomenon are:

1. The degradation of key proteins is too fast: the ClpXP system responsible for degrading LuxI proteins may be too efficient (equivalent to a "brake" that is too heavy), resulting in AHL signaling molecules being unable to effectively accumulate to the threshold that triggers cleavage.
2. The synthesis rate of key proteins is too low: The basal synthesis rate of LuxI or lytic proteins (equivalent to the upper limit of the "throttle") itself is low, and even with strong activation of AHL signals, the system's response ability is limited and insufficient to support the next round of strong lysis.

This simulation result warns us that when selecting or designing host strains, we must pay attention to their endogenous protein degradation and synthesis efficiency, ensuring that the positive feedback loop can be effectively triggered and amplified, otherwise the expected periodic lysis effect may not be achieved.

In Sum,

Parameter set	Ref.	Traits	Simulation	Experimental implication
Set 1 (E. coli)	Nat. Commun. 2020	High AHL synthesis/clearance; strong LuxI	Stable oscillations of N/H/L/I (~3-4 h period)	Predicts periodic autolysis; offers quantitative guidance
Set 2 (S. Typhimurium)	Nat. Microbiol. 2017	Low synthesis/clearance; weak LuxI	Damped oscillations toward steady state	Oscillations may vanish in some chassis; choose hosts carefully

Table 6.1. Two chassis, parameterized behaviors, and guidance

Suggested workflow: predicted period T_pred ≈ 120-180 min (placeholder). Measure OD600 every 10-15 min and photo-record clarity.

Record sheet (to backfill):

Period T (experiment): ________ min; amplitude (cycles 1/2): ________ / ________;

TBP = to be

Item	Symbol	Value
Average Period	T_measure	Pending
Amplitude of cycle	Amp_measure	Pending
Damping	Yes/No	Pending

Parameter set A (E. coli-like): predicts sustained synchronized-lysis oscillations;

Parameter set B (S. Typhimurium-like): predicts rapidly damped oscillations approaching steady state;

Different chassis exhibit markedly different dynamics in the synchronized lysis circuit (SLC). The model suggests that sustained oscillations are maintained only with specific hosts (e.g., E. coli-like parameters as in set A). Therefore, achieving the desired oscillatory behavior requires targeted screening and optimization of the host strain.

Due to time constraints, we were unable to obtain the results of the MTT experiment. However, we will conduct the experiment and perform pharmacokinetic modeling after Wiki Freeze. We will use PH-Sim to simulate different dosing regimens for individuals in East Asia and analyze the concentration-time spectra of TRACER-IL24 in plasma, distal colorectum and kidneys. We will calculate pharmacokinetic parameters, including AUC, t_max and c_max, to characterize in vivo PK features and provide data support and theoretical basis for future animal studies and potential therapeutic applications.

MTT results analysis

Process cell viability according to the formula below:

By fitting the dose-response curve to the raw data:

The fitting result yields the effect concentration parameter.

Convert the in vitro concentration effect relationship into effective drug dose mg/kg through in vivo pharmacokinetic PK:

• C: Target plasma concentration mg/L
• V_d: Distribution volume, drug dose/in vivo drug concentration=1.12L
• F: bioavailability, Value ≈ 0.2
• Weight: weight, modeled as East Asian, Default Weight = 60kg

Obtain the effective concentration parameter effective drug dose from above and conduct pharmacokinetic modeling.

1. Individual & Compound Create

Parameters	Type/Value	Source
Species	Human	Default
Population	East Asian	Default
Gender	Male	Default
Age	30.00 years	Default
Weight	60.03 kg	Default
Height	169.96 cm	Default
BMI	20.78 kg/m²	Default
Lipophilicity	-2.50 Log Units	Estimation
Fraction Unbound	0.25%	Estimation
Molecular Weight	23.80kDa	Database
Solubility at pH 7.00	0.21 mg/L	Estimation
Advanced Compound Parameters	/	Default

2. Administration Protocol Create

According to the calculated effective drug dose, four rounds of administration were set, and the administration method was set as Intravenous Infusion to simulate an intelligent microneedle delivery device. The infusion time for each round was set to 60.00 minutes, and the dose matched the effective drug dose.

3. Simulation Run

Curve Selection: 1. Distal Colon → Interstitial Concentration 2. Kidney → Interstitial Concentration 3. Peripheral Venous Blood → Plasma Concentration. Used to evaluate whether the drug efficacy target concentration can be reached, while assessing plasma and kidney exposure, and determining the efficacy vs. toxicity balance.

To explore the effects of different inoculation amounts on the population dynamics and lysis process of engineering bacteria, we used the Agent-Based Modeling (ABM) approach and simulated with NetLogo.

Figure 8.1. NetLogo simulation workflow

Initialization: Generate Num bacterial agents on a 2D grid as the initial inoculum (2000-10000 gradient). Assign each agent base attributes—lifespan, division cycle, and maximum health—with randomized variation to capture heterogeneity.

Nutrient-driven growth: At each timestep, bacteria consume nutrients on their current patch and decide to divide based on accumulated energy; biomass increase triggers the quorum-sensing module.

Quorum sensing: Viable cells secrete AHL, LuxR, and the AHL-LuxR complex at specified rates; molecules diffuse to neighboring patches and undergo degradation.

Lysis module: When local AHL-LuxR concentration remains above threshold, cells enter lysis state and synthesize lysis protein. Once this exceeds its threshold, the cell dies and instantly releases intracellular lysis protein as a burst, accelerating lysis in neighboring cells.

Outputs: At each step, record total cell count, biomass, and concentrations of signaling molecules and lysis protein to analyze synchronous lysis patterns.

(1) Relationship between inoculum size and lysis

Figure 8.2: The effect of different initial inoculation amounts on biomass accumulation and fragmentation cycles

[Caption]: Simulate setting Num=2000, 5000, 8000, 10000, 12000, 15000, 20000. The vertical axis represents biomass, and the horizontal axis represents time step (tick). As the vaccination amount increases, the lysis cycle is significantly shortened and the synchronization of the population is enhanced.

Figure 8.3: Inoculation volume fragmentation cycle

As the inoculation amount Num increases (2000 → 10000), the first lysis cycle is significantly shortened.

Figure 8.4: When the inoculation amount is Num=20000, the growth quantity of engineering bacteria.

[Caption]: When the inoculation amount is too high, too many lytic proteins are produced in the lysed area, and the rate of inactivation of lytic proteins is not sufficient to rapidly reduce the effective amount of lytic proteins, resulting in a low bacterial count after the first lysis and the system losing the observable second cycle.

The results indicate that the initial inoculation amount determines the accumulation speed of quorum sensing signals, which in turn affects the lysis cycle. In the future, the wet experimental section will introduce inoculation volume for validation.

(2) Visualization of spatiotemporal distribution:

Figure 8.5: Spatiotemporal visualization from NetLogo simulations

Caption: Each panel shows the lysis process under different values of Num. a.Num=2000; b.Num=2000; c.Num=2000; d.Num=2000. The upper-left subfigure plots total biomass and product over time; the upper-right shows spatial distributions (red = high lysis-protein concentration, green = viable cells). Results indicate that at higher inoculum levels, lysis events are more clustered and exhibit more stable periodicity.

• Spatial clustering of signaling molecules and lysis protein is visible in local patches, forming "hotspot" regions.
• Lysis events erupt locally and then spread to neighboring areas, driving synchronized lysis across the population.
• At higher inoculum, this spatial clustering effect is amplified, accelerating the onset of synchronized lysis.

The NetLogo simulations reveal, at the tissue scale, a quantitative relationship between inoculum size and lysis period:

• A higher initial inoculum accelerates quorum-signal build-up and shortens the lysis period.
• When the initial inoculum is too low, lysis initiation is delayed.

These findings provide a theoretical basis for large-scale industrial culture of engineered bacteria: by tuning the initial inoculum, the lysis rhythm can be precisely controlled to manage cultivation cycles and achieve optimal yield.

Our current work has built an unprecedented multi-scale modeling framework that not only provides multidimensional validation of the TRACER-IL24 delivery system but also forms a complete logical closed loop from "why do it (Why)" to "can it work (How)" and then to "how to use it (What's next)." The excellence of this framework lies in the fact that each layer of the model takes the conclusions of the previous layer as input. However, we must also soberly recognize that every model is a simplification of complex biological reality, with underlying assumptions and limitations that must be faced.

Our modeling story begins with the most macro-level clinical need and advances layer by layer, ultimately returning to the question of clinical application. The core strength of the work lies in deeply realizing the synergy between model and experiment—the model is no longer mere "validation," but rather a "guidance system" that runs throughout.

Although our modeling framework is logically coherent, each layer is built on certain idealized assumptions. These constitute the limitations of the current work and are precisely the starting points for future optimization.

Molecular-level limitations—challenges of static predictions and simplified environments:

Inherent uncertainty of structure prediction: AlphaFold provides a "static" most probable conformation. For the highly flexible peptides in TRACER, prediction confidence is inherently lower. The real activation process is dynamic, whereas our simulations can capture only two snapshots—"before activation" and "after activation."

Limitations of simulation duration: Our MD simulations are at the nanosecond (ns) scale, sufficient to verify stability and initial trends. However, observing a complete catalytic event or membrane traversal requires microseconds (μs) or longer, which is beyond our current computational resources.

Over-simplified environment: Our transmembrane simulations use an idealized pure lipid bilayer. Real cancer cell membranes are covered by a glycocalyx and embed numerous proteins, which can greatly affect TRACER's behavior.

System-level limitations—universality of parameters and neglect of individual variability:

Parameter dependence of ODE and ABM models: The ODE and NetLogo models of the population lysis circuit rely heavily on kinetic parameters from the literature (e.g., protein synthesis/degradation rates). These parameters come from specific strains and experimental conditions and may differ from the true parameters of our engineered strains. Before calibration with wet-lab data, the model's predictions serve more as qualitative trend guidance than precise quantitative forecasts.

"Average person" assumption in the PBPK model: Our pharmacokinetic model is based on a standard physiological database for an "East Asian healthy male." This overlooks large inter-patient variability (e.g., age, sex, body weight, hepatic/renal function) and, more importantly, does not account for the impact of disease states on physiological parameters (for example, tumors themselves may alter local blood flow and tissue permeability).

Clinical translation-level limitations—bridging the gap from organ concentration to tumor targeting:

Core challenge: uncertainty of target-site concentration: The greatest limitation of the PBPK model is that it predicts the average interstitial concentration of the drug in the "distal colorectum" organ region rather than the actual concentration within the tumor microenvironment. Tumor tissue has unique physiological features (e.g., high interstitial pressure, abnormal vasculature) that make permeation from blood into the tumor core much more difficult than in normal tissue. Therefore, while our current conclusion that "the target-site concentration may be insufficient" is an important warning, its precision still needs improvement.

Facing the above limitations is key to driving the project toward maturity. Future work will focus on calibrating models with experimental data and using more refined models to answer more complex questions.

Molecular level: toward more realistic dynamic simulations

Future plan: We will perform long-timescale simulations using enhanced sampling methods (such as umbrella sampling) to compute the transmembrane potential of mean force (PMF), thereby elevating the assessment of transmembrane potential from the qualitative "TRF framework" to a quantitative comparison of energy barriers. At the same time, we will verify the actual biological activity of IL-24 through wet-lab experiments and use the results to guide optimization of the linker design.

Production level: building an accurate "digital twin"

Future plan: The current parameter dependence of the models precisely underscores the importance of data back-filling. We will use wet-lab measurements such as growth curves and AHL concentrations to perform parameter fitting and calibration for the ODE and NetLogo models, turning them into "digital twin" models that can accurately predict the behavior of our specific engineered strains, thereby enabling precise control of the production process.

Clinical application level: building an integrated PBPK model for multiple administration routes

Core task—bridging the "tumor target" gap: Addressing the greatest limitation of the PBPK model, our future core effort is to establish an integrated multi-route PBPK model combining microneedles (systemic) and oral engineered bacteria (local), and on this basis introduce an independent tumor physiological compartment. This model will simulate how the drug overcomes physiological barriers of tumor tissue and predict whether the drug concentration in the tumor core can reach the EC90 level under the synergy of the two administration routes. This will be the most critical step in moving from "proof of concept" to "preclinical research."

Our current work has built an unprecedented multi-scale modeling framework; while its logical coherence is an advantage, its true methodological core and innovative highlight lies in the Transmembrane Readiness Framework (TRF) we proposed to solve key scientific questions. The entire modeling system can be viewed as a complete narrative constructed around the "input side" and "output side" of TRF.

In our project, the most critical design hypothesis is: "Can cleavage by MMP9 (matrix metalloproteinase MMP9) effectively 'activate' the transmembrane ability of the BRⅡ peptide?" This is a typical complex question that is difficult to answer with a single metric. The existing literature often uses single indicators in isolation—such as insertion depth, energy profiles, or tilt angle—lacking a "gold standard" for comprehensive evaluation.

Innovativeness and necessity of TRF:

To fill this gap, we proposed TRF. Its innovation lies in deconstructing the transmembrane process into three logically progressive layers for the first time: A—geometry and interactions (structural premise), B—thermo-dynamic (energy feasibility), and C—kinetics (process sustainability). This three-layer, eight-metric framework turns the vague problem of "transmembrane potential" into a set of specific, quantifiable, and comparable indicators. It is no longer a "yes/no" judgment but a scoring of "readiness." In our work, TRF clearly "adjudicated" that the post-cleavage Cut state is significantly superior to the Total state across all three layers. This conclusion provides the firmest microscopic theoretical support for the feasibility of the entire TRACER design.

TRF as the engine driving the project forward:

TRF is not merely a static evaluation tool; it in fact drives the logic of our preceding and subsequent work. Bioinformatics analysis (Chapter 2) and structural modeling (Chapter 3) provided validated input models for TRF (i.e., the Total and Cut conformations); and TRF's conclusion—"the activated state can indeed traverse the membrane"—directly forms the fundamental premise for the subsequent pharmacokinetic (PBPK) model (Chapter 7). Without the strong evidence from TRF, the assumptions in PBPK simulations that the drug can enter cells and exert effects would be castles in the air. Therefore, TRF is the bridge and engine that carries us from molecular design to clinical application prediction.

Although the TRF framework itself is systematic, the precision (clarity) of its conclusions is inevitably constrained by the "resolution" of the simulation data we feed into it. This constitutes a major limitation of the current work.

Qualitative thermodynamic evaluation: In TRF's B-layer (thermodynamics), we used steered MD (SMD) to qualitatively compare energetic resistance to transmembrane passage. This is computationally efficient but relatively coarse; it can only provide a ranking of "relative barrier heights," not precise free-energy barrier values (ΔG). This leaves our conclusion at the qualitative level of "Cut is easier than Total."

Idealization of the simulation environment: The entire TRF analysis was conducted in a simplified pure lipid bilayer without other membrane proteins. This ignores complex structures on real cancer cell membranes, such as the glycocalyx and receptors. It means our current "high TRF score" represents potential under ideal conditions, while behavior on crowded, complex real cell surfaces remains unknown.

Acknowledging these limitations provides a clear path for the future development and application of TRF. Our goal is to upgrade it from a project-specific solution to a more precise and general platform-style design tool.

Core upgrades: achieving quantitative and high-resolution TRF

Upgrade the B-layer—thermodynamic evaluation: Our top priority is to use enhanced sampling methods with higher computational cost but greater accuracy, such as umbrella sampling, to compute the potential of mean force (PMF) for peptide translocation. This will upgrade TRF from a qualitative comparison tool to a precise analytical framework that provides quantitative barrier values (units: kJ/mol). We will then not only know which is better but also "by how much," and we can directly compare with other thermodynamic data such as molecular binding energies.

Building more realistic simulation environments: In future simulations, we will construct cell membrane models closer to physiological reality—for example, introducing varying proportions of sphingolipids and cholesterol, and even embedding representative glycoproteins on the membrane surface—to test whether TRF evaluations remain robust under more complex conditions.

Application expansion: Turning TRF into a general-purpose design platform for the iGEM community

Evaluating other payloads: TRACER is a delivery platform. In the future, we can use a standardized TRF workflow to rapidly assess TRACER's potential to deliver other types of drugs. By simply swapping the payload, rerunning simulations, and applying TRF analysis, we can efficiently predict the transmembrane efficiency of new fusion proteins and greatly accelerate preliminary drug screening.

Guiding rational design and evolution of TRACER peptides: The most exciting application prospect of TRF is to use it as a computer-aided design (CAD) tool. We can perform in silico multi-site mutations on the BRⅡ peptide segment to generate dozens of virtual candidate designs, then rapidly score and rank them with TRF, and submit only the top-scoring designs for wet-lab validation. This will fundamentally change traditional trial-and-error-dependent peptide optimization and enable truly rational design.

In summary, the TRF paradigm is the core innovation of our work. It has now successfully validated the rationality of our design; in the future, with continuous quantitative upgrades and expansion of application scenarios, it has the potential to become a powerful and efficient computational simulation platform serving a broader scope of peptide drug design—perfectly embodying the iGEM project's core value of "contributing to the community."