A new mRNA vaccine targeting multiple cancer markers for colorectal cancer: computer design and immune response testing

The UniProt database (https://www.uniprot.org/↗) was utilized to acquire amino acid sequences of the six target proteins using their respective accession numbers (1): THBS2 (P35442), (2) FSTL3 (P95633), (3) TNNT1 (P13805), (4) BGN (P21810), (5) CTHRC1 (Q96CG8), and (6) NOX4 (Q9NPH5).

Following the prediction of linear B-cell epitopes (LBL), cytotoxic T lymphocyte (CTL), and helper T lymphocyte (HTL) epitopes, we systematically evaluated the immunobiological properties of candidate epitopes. Antigenicity was assessed using the VaxiJen v3.0 web server (http://www.ddg-pharmfac.net/vaxijen/↗), configured with the tumor protein antigen model at a threshold score of 0.5. This machine learning platform employs alignment-free algorithms trained on bacterial, viral, and tumor antigen datasets to predict whole-protein immunogenicity. Allergenicity screening was performed through AllerTop v2.1 (https://www.ddg-pharmfac.net/allertop_test/↗), which utilizes amino acid propensity scales and auto-cross covariance transformation for epitope safety evaluation. Toxicity predictions were conducted via ToxinPred (https://webs.iiitd.edu.in/raghava/toxinpred/↗), employing support vector machine models trained on toxic/non-toxic peptide datasets. Epitopes demonstrating strong antigenicity (VaxiJen score ≥0.5), non-allergenicity, and non-toxicity were retained for downstream analyses.

Geographical and environmental factors drive substantial diversity in human leukocyte antigen (HLA) alleles across global populations. To evaluate the epidemiological relevance of our vaccine candidate, we conducted global population coverage analysis using the Immune Epitope Database (IEDB) Population Coverage Calculation Tool (http://tools.iedb.org/population/↗). This platform calculates cumulative coverage probabilities for selected HLA class I (CTL) and class II (HTL) epitopes based on their binding affinities to region-specific HLA allele distributions across 16 geographical regions. The analysis integrated frequency data of HLA-A, -B, -C, and -DRB1 alleles from ethnically diverse populations to estimate potential vaccine efficacy thresholds.

To evaluate the binding affinity between selected T lymphocyte epitopes and their corresponding major histocompatibility complex (MHC) class I alleles, we performed systematic molecular docking simulations. MHC-I crystal structures were retrieved from the RCSB Protein Data Bank (PDB). Structural preprocessing involved removal of water molecules and non-essential ligands using PyMOL v2.5.7 (28). Epitope sequences underwent 3D conformational modeling via PEP-FOLD 3.5 server (https://bioserv.rpbs.univ-paris-diderot.fr/services/PEP-FOLD3/↗), followed by complex assembly with MHC-I proteins. Energy minimization was conducted using Rosetta's Relax protocol (v3.13) to optimize molecular geometries (29). Finally, epitope-MHC docking simulations were executed in local refinement mode through RosettaDock, employing Monte Carlo minimization algorithms to sample low-energy binding conformations.

The 5'UTR sequence in particular is a major determinant of translation efficiency and thus an intriguing target for engineering (30). The ribosomal 43 S pre-initiation complex (PIC) scans the 5'UTR in the 5'-to-3' direction until a start codon is found. Therefore, 5'UTRs can affect translation by capturing PICs prematurely via upstream start codons (uAUGs) and ORFs (uORFs) (31). To identify 5′-UTR variants with enhanced translational efficiency, we employed RNA-FM (32), a pretrained large language model, to generate embeddings for 83,919 synthetic human 5′-UTRs of 75 distinct lengths from the combinatorial library generated by Sample et al. (33). The convolutional neural network (CNN) architecture was constructed based on Sample et al. to screen candidate 5′-UTRs. The CNN model was constructed as three 1D convolutional layers with 120 filters and a ReLU activation for each layer. The third convolution layer output one channel, which was fed into two fully-connected layers with one output node as the final prediction. The inputs of the model were RNA-FM embeddings (640 dim). In this study, the data was partitioned into a training set and a validation set with a ratio of 8:2. This partitioning strategy aimed to ensure sufficient data for model training while enabling effective validation of the model's generalization ability.

Our 5′-UTR candidate pool integrates four distinct sources: 1) 212 variable-length sequences from Chu et al. (34), 2) 8 Kozak sequence-containing variants from Li et al. (35), 3) The 5′-UTR of human α-globin genes, and 4) the 5′-UTR of BNT162b2(SARS-CoV-2 vaccine, from Pfizer/BioNTech) (36).

A potent mRNA vaccine construct requires systematic integration of five critical components (1): an open reading frame (ORF) encoding antigenic elements (2), the Untranslated Regions (UTRs) flanking the coding regions, (3) a Kozak sequence incorporating the start codon (GCCACCAUGG) to enhance translational initiation (37), (4) functionally optimized linkers, and (5) regulatory termination signals. The proposed construct features a 5'→3' architecture comprising a modified 5′m⁷GCap structure, followed by an optimized 5' untranslated region (5'UTR) and Kozak sequence to maximize ribosomal engagement. The ORF initiates with a signal peptide for intracellular trafficking, connected via an EAAAK rigid linker to helper T lymphocyte (HTL) epitopes interconnected through GPGPG spacers that maintain domain autonomy. This HTL cluster transitions via KK linkers to linear B-cell epitopes (LBL), followed by AAY-linked cytotoxic T lymphocyte (CTL) epitopes strategically positioned upstream of an MHC class I trafficking domain (MITD). The MITD sequence enhances immunogenicity through dual mechanisms: 1) optimizing antigen presentation efficiency via endosomal targeting motifs, and 2) directing vaccine components to antigen-processing compartments. Epitope segregation through GPGPG, KK, and AAY linkers ensures proper conformational folding while preserving immunological functionality. The construct terminates with a UAA stop codon, a stabilizing 3' untranslated region (3'UTR), and a 120-nucleotide poly(A) tail to ensure mRNA stability and translational fidelity. This multi-layered design leverages structural bioengineering principles to balance epitope accessibility, intracellular trafficking efficiency, and immune activation capacity.

Following vaccine sequence assembly, systematic bioinformatic validation was performed to assess four critical properties: (1) antigenicity, (2) allergenicity, (3) toxicity, and (4) physicochemical stability. Antigenicity prediction employed dual machine learning platforms—VaxiJen v3.0 (threshold: 0.5; tumor antigen model) and ANTIGENpro (SCRATCH Protein Predictor suite)—using the translated amino acid sequence excluding trafficking/processing elements (tPA and MITD domains). Allergenicity screening was conducted through AllerTop v2.1 (https://www.ddg-pharmfac.net/allertop_test/↗) using default parameters. Toxicity profiling utilized ToxinPred (https://webs.iiitd.edu.in/raghava/toxinpred/↗), implementing SVM models trained on experimentally validated toxic peptides. Physicochemical characterization was performed via the ProtParam tool (https://web.expasy.org/protparam/↗), quantifying six essential parameters: amino acid composition, molecular weight, theoretical isoelectric point (pI), instability index (II), aliphatic index (AI), and grand average of hydropathicity (GRAVY). This multi-platform validation strategy ensures structural integrity while confirming immunological safety and biophysical stability of the vaccine candidate.

To predict the immunogenic potential of the mRNA vaccine construct in humans, we conducted in silico immune simulations using the C-IMMSIM immunoinformatics platform (https://kraken.iac.rm.cnr.it/C-IMMSIM/↗). The vaccination regimen was configured to follow the recommended dosage schedule of current vaccines: three 1,000-unit doses administered at time-steps 1, 84 (3 weeks), and 168 (6 weeks). All parameters were maintained at default settings. Dynamic simulation outputs quantified key immunological metrics: 1) Antigen-specific lymphocyte proliferation rates, 2) Cytokine production profiles, 3) Memory cell differentiation kinetics, and 4) Antibody titer trajectories. This computational framework enables systematic evaluation of the vaccine's capacity to elicit coordinated humoral and cellular immune responses while predicting long-term immunological memory formation.

Leveraging the degeneracy of the genetic code, we implemented a multi-algorithm optimization strategy to enhance translational efficiency of the mRNA vaccine construct. Five codon optimization methods, including LinearDesign (38), were systematically applied to resolve synonymous codon conflicts while maximizing expression potential. Post-optimization sequences underwent rigorous biophysical characterization through three key metrics: 1) Minimum free energy (MFE) of mRNA secondary structures predicted by RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi↗), 2) GC content analysis (optimal range: 45-55%), and 3) Codon Adaptation Index (CAI) calculations relative to human codon usage tables. This combinatorial approach achieves dual objectives: (a) elimination of cryptic splice sites and ribosomal drop-off sequences through codon bias correction, and (b) stabilization of mRNA architecture via thermodynamic optimization of folding patterns.

The RNAfold web server was employed to predict the secondary structure of the mRNA vaccine construct using McCaskill's partition function algorithm (39). This computational framework calculates the minimum free energy (MFE) conformation through dynamic programming-based thermodynamic modeling of RNA folding pathways. The analysis yielded two critical outputs: 1) the MFE-optimized secondary structure visualization, and 2) quantitative thermodynamic stability metrics (ΔG in kcal/mol). Subsequent evaluation of mRNA structural features focused on identifying persistent stem-loop formations and regions of high base-pairing potential that might impede ribosomal scanning. This characterization enabled rational design optimization to balance thermodynamic stability (ΔG < -300 kcal/mol) with translational efficiency through strategic codon rearrangement in unstable regions (MFE > -150 kcal/mol). The refined mRNA architecture demonstrates enhanced resistance to endonucleolytic degradation while maintaining optimal ribosomal accessibility-critical parameters for ensuring structural integrity and sustained antigen expression in human physiological conditions.

To investigate the vaccine construct's innate immune activation potential, we performed protein-protein docking between predicted vaccine epitopes and human Toll-like receptor 3 (TLR3; PDB ID 1ZIW) and Toll-like receptor 4 (TLR4; PDB ID 3FXI) using the ClusPro 2.0 server (https://cluspro.org/↗). This rigid-body docking platform employs a hierarchical protocol combining fast Fourier transform (FFT) correlation approaches with Monte Carlo minimization to sample >10⁹ possible binding conformations.

B-cell epitopes were predicted for the amino acid sequences of THBS2, FSTL3, TNNT1, BGN, CTHRC1, and NOX4 using the ABCPred web server. For each protein, the top five scoring epitopes were initially retained based on prediction rank. Subsequent refinement employed a tripartite filtering strategy: 1) Antigenicity validation via VaxiJen v3.0 (threshold: 0.5; tumor antigen model), 2) Allergenicity screening using AllerTop v2.1, and 3) Toxicity profiling through ToxinPred. A total of six epitopes demonstrating strong antigenicity (VaxiJen score ≥0.5), non-allergenicity, and non-toxicity were selected for final inclusion (). This stringent selection protocol ensures exclusive retention of epitopes with optimal immunogenic potential while mitigating risks of hypersensitivity or cytotoxic responses. 1

Potential cytotoxic T lymphocyte (CTL) epitopes across the six target proteins were predicted using the MHC-I Binding Predictions tool within the IEDB Analysis Resource (http://tools.iedb.org/mhci/↗). Predicted epitopes were subsequently filtered through a triaxial validation pipeline: 1) Antigenicity scoring via VaxiJen v3.0 (threshold ≥0.5), 2) Allergenicity assessment using AllerTop v2.1, and 3) Toxicity profiling through ToxinPred. Epitopes were prioritized based on combined metrics: percentile binding rank (<1.0), antigenicity score (>0.7), and evolutionary conservation index (>0.8 in ConSurf analysis). Twelve high-affinity epitopes residing in phylogenetically conserved regions were ultimately selected for vaccine incorporation (Supplementary Table S1), ensuring broad HLA coverage and variant-resistant immunogenicity.

Potential helper T lymphocyte (HTL) epitopes were identified through comprehensive analysis of the six target proteins using the IEDB MHC-II Binding Predictions tool (http://tools.iedb.org/mhcii/↗) with NetMHCIIPan 4.1 algorithm. Following rigorous triaxial screening—antigenicity assessment (VaxiJen v3.0 score ≥0.5), allergenicity profiling (AllerTop v2.1), and toxicity evaluation (ToxinPred), 18 epitopes residing in phylogenetically conserved domains (ConSurf conservation score >0.85) were prioritized for vaccine inclusion (Supplementary Table S1). After calculation and filtering, 18 selected HTL epitopes could induce Th1/Th2-polarized cytokine responses (IFN-γ, IL-4, IL-10) through TCR-MHC II complex stabilization, thereby amplifying both humoral and cellular arms of vaccine-induced immunity via cognate T-B lymphocyte collaboration.

The worldwide epidemiological relevance of the vaccine construct was evaluated through HLA allele population coverage analysis using the IEDB Population Coverage Calculation Tool (http://tools.iedb.org/population/↗). The 30 incorporated HLA alleles demonstrated differential coverage across 16 geographical regions, with peak efficacy observed in Europe (89.13%) and North America (84.69%) (Table 1). Comparatively reduced coverage rates were identified in Central America (15.29%), South America (52.23%), and Northeast Asia (52.38%), reflecting regional disparities in HLA allele distribution patterns (details in Supplementary Figure S1).

Through the previous prediction and screening steps, 30 lymphocyte epitopes and their corresponding HLA alleles were identified (Supplementary Table S2), and six of these epitope-allele pairs underwent molecular docking simulations using Rosetta's LocalDock module. Crystal structures of selected HLA class I molecules (PDB IDs: 8EMF, 7LG0, 7L1C, 4LNR, 5VUD, 6MPP) were energy-minimized to remove steric clashes (<0.3 Å RMSD deviation). Docking metrics (Supplementary Table S3) revealed two critical energy parameters: (1) Total_Score, predominantly reflecting monomeric folding energy, and (2) Interface_Score (I_sc), quantifying interaction energy across the binding interface. As emphasized in Rosetta documentation, I_sc provides superior predictive value for epitope-HLA binding stability. The MSDTEEQEY epitope demonstrated optimal binding with HLA-A*01:01, achieving the lowest I_sc value (-115.51 REU) and complete structural accommodation within the HLA binding groove (Figure 2F).

Our CNN model was trained and validated over 50-epoch using 83,919 human 5′-UTR sequences (spanning 75 distinct lengths) from Sample et al.'s library. Optimal regression performance was achieved at epoch 15, yielding an R² value of 0.844 with corresponding error metrics (mean squared error [MSE] = 0.307, mean absolute error [MAE] = 0.397, root mean squared error [RMSE] = 0.554). External prediction conducted on 222 literature-curated 5′-UTR sequences identified 10 top-performing candidates, with complete rankings detailed in Table 2. The highest-ranked variant (sequence: GGGATCTTATTCCACCTTCTGAAGCTTCTGTCGAACCAGTTGTAAGGAGA) was ultimately selected as the 5′-UTR component for our CRC mRNA vaccine construct.

The mRNA vaccine construct was proposed to be arranged from the N to C terminus in the following order:5′m⁷GCap, 5′UTR, Kozak sequence, Signal peptide (tPA)–GPGPG Linker–LLQVVYLHSNNITKV–GPGPG Linker–TTLLDLQNNDISELR–GPGPG Linker–KISKIHEKAFSPLRK–GPGPG Linker–RDGFKGEKGECLRES–GPGPG Linker–GRDGFKGEKGECLRE–GPGPG Linker–MNSTINIHRTSSVEG–GPGPG Linker–FVSSMGSGNPAPGGV–GPGPG Linker–GFVSSMGSGNPAPGG–GPGPG Linker–VSSMGSGNPAPGGVC–GPGPG Linker–IQKIIGEKYHALNSR–GPGPG Linker–KPAEFTQHKFVKICM–GPGPG Linker–RWKLLFDEIAKYNRG–GPGPG Linker–NCPYVHNPAQIDTDN–GPGPG Linker–DNCPYISNANQADHD–GPGPG Linker–TAQLKQDGKSRGTLL–GPGPG Linker–SEKFDLMAKLKQQKY–GPGPG Linker–FDLMAKLKQQKYEIN–GPGPG Linker–QQRFRTEKERERQAK–KK Linker–DEEASGADTSGVLDPD–KK Linker–ACCQRWYFTFNGAECS–KK Linker–PQSCVVDQTGSAHCVV–KK Linker–IGRPRWKLLFDEIAKY–KK Linker–SGTQQRGRSCDVTSNT–KK Linker–YNRISHAQKFRKGAGK–AAY Linker–DPDSVTPTY–AAY Linker–KLQKLYISK–AAY Linker–KQKAQLRQR–AAY Linker–IPKGKQKAQL–AAY Linker–MRPGAPGPLW–AAY Linker–VMYRGRCRK–AAY Linker–SPFEESLNY–AAY Linker–QANFPQTWLW–AAY Linker–SVDFSGTFY–AAY Linker–RVSNDNQFLW–AAY Linker–MSDTEEQEY–AAY Linker–KEEEELVAL–AAY Linker–MITD sequence–Stop codon–3′UTR–Poly(A) tail. The schematic diagram of the vaccine construct is shown in Figure 3, where key components (signal peptide, HTL/LBL/CTL epitope clusters, linkers, and regulatory elements) are color-coded to clarify the 5′→3′ architecture.

Following assembly of the complete mRNA vaccine construct, we systematically evaluated the translated polypeptide sequence through three critical safety parameters: (1) antigenicity prediction using VaxiJen v3.0 (threshold ≥0.5; tumor antigen model) and ANTIGENpro (SCRATCH suite), (2) allergenicity screening via AllerTop v2.1, and (3) toxicity profiling through ToxinPred. Physicochemical characterization was performed using the ExPASy ProtParam server (https://web.expasy.org/protparam/↗), quantifying six biostability metrics: molecular weight (54603.7), theoretical pI (9.27), instability index (47.94), aliphatic index (60.57), grand average of hydropathicity (GRAVY: -0.776), and thermal stability (in vitro half-life >30 hours at 50 °C). As summarized in Table 3, the construct demonstrated strong antigenicity (VaxiJen: 0.7284; ANTIGENpro: 0.93), non-allergenicity, and non-toxicity.

Three-dose immunization simulations (1,000 vaccine units per dose) conducted via the C-IMMSIM platform revealed coordinated adaptive immune responses (Figure 4). Primary analysis demonstrated predominant IgM over IgG titers post-initial vaccination, with significant immunoglobulin amplification following booster doses (Figure 4A). Sustained antibody elevation post-antigen clearance suggests established immunological memory, enabling rapid anamnestic responses upon antigen re-exposure. T lymphocyte profiling showed progressive expansion of both activated and resting helper T cell populations, stabilizing at elevated levels within 30 days post-immunization (Figures 4D, E). Peak T helper cell proliferation occurred at day 10, coinciding with macrophage activation maxima (days 5-10). Resting macrophage populations initially declined (days 0-2) before rebounding through monocyte differentiation (day 3 onward). These data indicate successful antigen presentation and lymphocyte priming. The simulation data collectively confirm the mRNA vaccine's capacity to orchestrate robust humoral and cellular immunity against CRC.

To maximize translational efficiency of the mRNA vaccine construct, we implemented a comparative codon optimization strategy employing five algorithms: LinearDesign, JCat, OPTIMIZER, Gensmart, and ExpOptimizer. The optimized nucleotide sequences and corresponding metrics—Codon Adaptation Index (CAI), GC content, and minimum free energy (MFE)—are detailed in Table 4.

All algorithms except LinearDesign achieved comparable CAI values of approximately 0.8 (range: 0.81-0.96), while maintaining optimized GC content within the preferred human codon usage optimum (52.25-67.43%). Among conventional approaches, JCat demonstrated superior performance in both CAI (0.96) and GC content (67.43%). However, LinearDesign exhibited exceptional thermodynamic stability with an MFE of -1193.70 kcal/mol, exceeding other algorithms by more than 2.1-fold. Considering the critical importance of mRNA structural stability in vivo, we selected LinearDesign-optimized sequence for further development, thereby ensuring sustained antigen presentation and robust immunogenicity in human hosts.

The mRNA vaccine construct's secondary structure was computationally validated through a dual-platform approach: 1) LinearDesign's constraint programming framework predicted MFE conformations, and 2) RNAfold's partition function algorithm (ViennaRNA Package 2.6) provided comparative thermodynamic profiling (Figure 4). LinearDesign yielded an MFE of -1193.70 kcal/mol, while RNAfold analysis revealed enhanced stability (MFE = -1233.87 kcal/mol). The subtle differences in the secondary structures of the sequences obtained by the two optimization methods can be seen in Figure 5. This convergence between distinct computational methodologies (<3.4% MFE variance) confirms exceptional structural stability. The robust thermodynamic profile correlates with enhanced vaccine efficacy through three mechanisms: 1) Reduced secondary structure interference with translational initiation complexes, 2) Increased nuclease resistance via compact folding motifs, and 3) Maintenance of epitope codon optimality under physiological temperature fluctuations (37 ± 2°C). These features collectively ensure sustained antigen expression critical for eliciting durable immune responses.

Protein-peptide docking between the vaccine-encoded epitope structure and human Toll-like receptor 3 (TLR3; PDB ID 1ZIW) and Toll-like receptor 4 (TLR4; PDB ID 3FXI)was performed using ClusPro 2.0 (https://cluspro.org/↗), the peptide sequences produced after the translation of the mRNA vaccine are predicted by AlphaFold2 (40). with the highest-ranked TLR3-vaccine complex demonstrating a binding energy of -1378.8 kcal/mol and the highest-ranked TLR4-vaccine complex's binding energy of -1581.2 kcal/mol (Supplementary Table S4, Figure 6). These results validate the vaccine construct's capacity to engage TLR3 and TLR4-mediated endosomal sensing pathways, initiating robust dendritic cell maturation and subsequent Th1-polarized adaptive immunity–critical features for eliciting durable antitumor responses in colorectal cancer.