Vaccine candidate designed against carcinoembryonic antigen-related cell adhesion molecules using immunoinformatics tools

Abstract

Carcinoembryonic antigen-related cell adhesion (CEACAM) molecules belong to a family of membrane glycoproteins that mediate intercellular interactions influencing cellular growth, immune cell activation, apoptosis, and tumor suppression. Several family members (CEACAM1, CEACAM5, and CEACAM6) are highly expressed in cancers, and they share a conserved N-terminal domain that serves as an attractive target for cancer immunotherapy. A multi-epitope vaccine candidate against this conserved domain has been developed using immunoinformatics tools. Specifically, several epitopes predicted to interact with MHC class I and II molecules were linked together with appropriate linkers. The tertiary structure of the vaccine is generated by homology and ab initio modeling. Molecular docking of epitopes to MHC structures have revealed that the lowest energy conformations are the epitopes bound to the antigen-binding groove of the MHC molecules. Subsequent molecular dynamics simulation has confirmed the stability of the binding conformations in solution. The predicted vaccine has relatively high antigenicity and low allergenicity, suggesting that it is an ideal candidate for further refinement and development.

Communicated by Ramaswamy H. Sarma

Keywords:

• Carcinoembryonic antigen-related cell adhesion molecules
• vaccine design
• immunoinformatics
• molecular modeling

Introduction

The carcinoembryonic antigen-related cell adhesion molecules (CEACAM) family includes 12 members that generally have one (sometimes two) Immunoglobulin (Ig)-like variable (V)-set domain, but they differ in the number of Ig-like constant C2-set domains, as well as the membrane anchorage (Beauchemin et al., 1999). Four members (CEACAM5-8) are associated with the membrane through a glycosylphosphatidylinositol (GPI) linkage, whereas seven members (CEACAM1, 3, 4, 18-21) are anchored to the cellular membrane via bona fide transmembrane domains. Only one member of the CEACAM family, CEACAM16, is a secreted protein with no membrane anchorage (Beauchemin & Arabzadeh, 2013).

Functionally, the expression of CEACAM molecules starts early in human embryonic and fetal development (weeks 9-14) (Eades-Perner et al., 1994; Nap et al., 1988) and is significantly elevated in colorectal (Jothy et al., 1993), gastric (Kodera et al., 1993), lung (Singer et al., 2010), pancreatic (Gebauer et al., 2014), and skin (Khatib et al., 2011) carcinoma, and thus they belong to a larger family of ‘carcinoembryonic antigens (CEA)’ (Gold & Freedman, 1965). In normal adult tissues, CEA is localized in the stomach, tongue, esophagus, cervix, sweat glands and prostate, as well as in columnar epithelial and goblet cells of the colon (Hammarstrom, 1999). Research in the past five decades has established that CEACAM molecules are involved in diverse functions in cell adhesion and signaling, and play important roles in cancer progression, inflammation, angiogenesis, and metastasis (Beauchemin et al., 2013; Gray-Owen & Blumberg, 2006; Muenzner et al., 2005; Sadarangani et al., 2011; Tchoupa et al., 2014). Three CEACAM proteins (CEACAM1, CEACAM5, and CEACAM6) are considered as valid clinical biomarkers and recently emerged as attractive therapeutic targets for cancer immunotherapy (Dankner et al., 2017; Horst & Wagener, 2004; Kuespert et al., 2006). Indeed, one member of the CEACAM family, CEACAM5, was ranked 13^th out of 75 representative cancer antigens based on a suite of pre-defined and pre-weighted criteria (Cheever et al., 2009). Since CEACAM molecules share the conserved Ig-like V domain at the N terminus, a cancer vaccine that targets this region potentially has a universal antitumor effect on all cancers that overexpress CEACAM proteins.

The conventional approach to vaccine development typically involves time-consuming and expensive experimental studies, along with ethical concerns. As such, there is a growing interest in utilizing bioinformatics and immunoinformatics tools for vaccine design (Atapour et al., 2020; Dorosti et al., 2019; Hajighahramani et al., 2019; Khalili et al., 2015; Mirza et al., 2016; Pandey et al., 2018; Sabetian et al., 2019), which is further empowered by the recent development of synthetic genomics (Bambini & Rappuoli, 2009). Immunoinformatics is an emerging field that interfaces computer science and experimental immunology, aiming to use computational methods and resources for the understanding of immunological information (Korber et al., 2006; Tong & Ren, 2009). One of its primary goals is to develop algorithms to predict potential B- and T-cell epitopes, which reduces the time and cost required for laboratory analysis of antigens. The application of these in silico techniques for epitope mapping has accelerated the development of vaccines (Li Pira et al., 2010; Zhao et al., 2013).

In this study, a multi-epitope vaccine candidate has been designed targeting the conserved V-set domain of CEACAM molecules (Figure 1). Several CD4+ and CD8+ epitopes are predicted by a set of immunoinformatics tools, which are then linked together by appropriate linkers for the enhancement of epitope presentation and separation. To the best of our knowledge, this is the first vaccine candidate targeting the conserved N terminal domain of CEACAM molecules that are overexpressed in a variety of cancers, which deserves further refinement and development.

Figure 1. Research strategy workflow. Flowchart of the methods for the prediction and validation of the CEACAM vaccine candidate. Boxes with sharp edges denote methods or tools, whereas boxes with rounded edges represent the data cumulated from tools. The validation for molecular docking/molecular dynamics simulation includes RMSD, RMSF, radius of gyration, snapshots during simulations and binding free energies. The validation for molecular modeling includes Ramachandran plots.

Materials and methods

Sequence and structure retrieval of CEACAM molecules

The protein sequences of CEACAM1 (P13688), CEACAM5 (P06731) and CEACAM6 (P40199) were retrieved from the National Center for Biotechnology Information (NCBI) protein databases. The crystal structure of CEACAM1 Ig-like V-set domain (PDB ID 5DZL) was retrieved from the RCSB PDB database (Berman et al., 2000).

Multiple sequence alignment

Multiple sequence alignment of the N-terminal domain of CEACAM1, CEACAM5, and CEACAM6 was performed using the T-Coffee (Notredame et al., 2000) multiple sequence aligner hosted by the EBI. T-Coffee uses its progressive alignment algorithm to perform the alignment and is set to produce an alignment output in the ClustalW format.

Prediction of cytotoxic T-lymphocyte (CTL) and helper T-lymphocyte (HTL) epitopes

Five epitope prediction methods including IEDB (Moutaftsi et al., 2006), NetMHC 4.0 (Andreatta & Nielsen, 2016; Nielsen et al., 2003), BIMAS server (Parker et al., 1994), SYFPEITHI server (Rammensee et al., 1999), and ProPred server (Singh & Raghava, 2001) were harnessed to determine potential CTL and HTL epitopes within the CEACAM1 Ig-like V-set domain (see parameters and thresholds of the servers in Supplementary Table S1). These methods predict the peptide epitopes from the antigen of interest based on different algorithms. Specifically, the IEDB server takes in the antigen sequence and runs a sequence alignment method based on artificial neural networks (ANNs) to predict 8-11 amino acid long peptide epitopes. The length was restricted to 9 peptides per epitope for MHC class I and up to 25mers for MHC class II molecules. The NetMHC 4.0 server uses a preset training of 81 MHC alleles to produce novel 9-mer epitopic representations for MHC class I molecules and up to 15mer sequences for MHC class II molecules. The BIMAS server ranks the 9-mer sequences based on independent binding of individual peptide side-chains. The ProPred server uses quantitative matrices that identify promiscuous binding regions useful in selecting vaccine candidates.

Table 1. Top MHC-I-bound epitopes identified from the antigen (* represents different alleles).

Epitope	MHC I Allele	Allfele PDB ID	IEDB	ProPred	NetMHC 4.0	BIMAS	SYPEITHI
IYPNASLLI	HLA-A*24:02(A24)	5XOV	YES	YES	YES	YES	YES
QLFGYSWYK	HLA-A*03:01(A3)	3RL1	YES	YES	YES	YES	YES
QLFGYSWYK	HLA-A*11:01(A11)	n/a	YES	YES	YES	YES	YES
VTQNDTGFY	HLA-A*01:01(A1)	3BO8	YES	YES	YES	YES	YES
QLFGYSWYK	HLA-A*68:01(A68.1)	4HX1	YES	YES	YES	YES	YES
EATGQFHVY	HLA-A*26:01	n/a	YES	n/a	YES	n/a	Yes
VTQNDTGFY	HLA-A*30:02	n/a	YES	n/a	YES	n/a	n/a

Vaccine design and modeling of the 3 D structure of the vaccine

The predicted epitopes were arranged in the order of the antigen, CEACAM1 Ig-like V-set domain, sequence. The epitopes were linked to form a single vaccine candidate with AAY and GPGPG motifs as the linkers to fill gapped sequences between MHC class I and MHC class II epitopes respectively to enhance epitope presentation and separation.

The 3 D structure of the vaccine was predicted by both homology modeling and ab initio modeling. Modeller within the UCSF’s Chimera tool (Pettersen et al., 2004) was used to find the three-dimensional homology structure of the vaccine candidate using the original antigen structure (PDB ID 5DZL) as the template. The vaccine structure was also predicted using PEP-FOLD 3.0 server (Shen et al., 2014), which predicts structure based on structural alphabet (SA) letters describing the conformations of groups of four consecutive residues. PEP-FOLD 3.0 couples the predicted series of the SA letters with a greedy algorithm and a coarse-grained force field to predict a final 3 D structure.

The Rampage server (Biasini et al., 2014) was utilized to assess the quality of the predicted structures by showing the number of residues falling in the favorable and unfavorable regions based on the phi and psi angles of rotation in a molecule.

Physiochemical property, antigenicity, and allergenicity analysis of the vaccine

The analysis of the vaccine construct was done on the ProtParam server (Wilkins et al., 1999) which predicts the pI, solubility, Molecular Weight, Half-Life (in-vitro) and Grand Average of Hydropathicity Values (GRAVY) for the input protein sequence.

The antigenicity of the predicted vaccine’s peptide sequence was verified using the ANTIGENpro tool from UC Irvine's Scratch Protein Prediction server (Magnan et al., 2010). The peptide sequence of the vaccine structure was entered into the tool and an antigen probability score was produced.

The allergenicity of the candidate vaccine was verified using the AllerTop v. 2.0 web server (Dimitrov et al., 2013) from the Medical University of Sofia's Drug Design Group. This tool accepts the vaccine sequence as input and runs a sequence structure mining algorithm on a dataset of known allergens and non-allergens. The tool then determines which member of its database most resembled the vaccine sequence and reports if that protein, as well as the input vaccine sequence, is an allergen or not.

Interferon γ epitopes have been used in the field of immunology to induce the innate as well as the adaptive immune system to elicit anti-tumor pathways (Castro et al., 2018). Using the IFNepitope server (Dhanda et al., 2013), possible IFN γ epitopes were predicted from the vaccine construct by overlapping regions of protein and predicting their potency and ability to induce IFN γ. Predictions were calculated using a Support Vector Machine (SVM) predictor based on the IEDB’s helper T cell database.

Immune simulation

C-ImmSim 10.1 Server (http://150.146.2.1/C-IMMSIM) was used to interpret the host immune response to the antigen (i.e. the vaccine construct). C-ImmSim is an agent-based computational immune response simulator that utilizes position-specific score matrix (PSSM) and machine learning methods for predicting epitope and immune interactions, respectively (Rapin et al., 2010). The parameters in the server were set based on the predominant HLA alleles of predicted epitopes (Tables 1 and 2). The host HLA selection parameter for MHC class I was set on A1010, A2402, and B0702 and for DR MHC class II was sat on DBR1_1101. In accordance with the 3-dose schedule of another cancer vaccine, HPV vaccine, recommended by CDC, the second dose should be given 1-2 months after the first dose, and the third dose should be given 6 months after the first dose. All simulation parameters were set at default with time steps set at 1, 126, 546 (each time step is 8 h and time step 1 is injection at time = 0). Therefore, the intervals between the first and the second dose, as well as between the second and the third dose were 6 and 20 weeks respectively.

Table 2. Top MHC-II-bound epitopes identified from the antigen (* represents different alleles).

Peptide	Allele	Allele PDB ID	IEDB	NetMHC 4.0	NetMHC Core
TQNDTGFYTLQVIKS	HLA-DRB1*11:01	6CPN	Yes	Yes	FYTLQVIKS
QQLFGYSWYKGERVD	HLA-DRB1*09:01	n/a	Yes	Yes	YSWYKGERV
VDGNRQIVGYAIGTQ	HLA-DRB1*08:01	n/a	Yes	n/a
AEGKEVLLLVHNLPQ	HLA-DRB1*04:04	6BIY	Yes	Yes	VLLLVHNLP
KEVLLLVHNLPQQLF	HLA-DRB4*01:03	n/a	Yes	Yes	LLLVHNLPQ
AEGKEVLLLVHNLPQ	HLA-DRB1*04:02	n/a	Yes	Yes	LLLVHNLPQ
GKEVLLLVHNLPQQL	HLA-DRB1*13:02	n/a	Yes	Yes	LVHNLPQQL

Epitope structure modeling

The epitope structures were modeled by the PEPFOLD 3.0 server (Shen et al., 2014), which runs a de novo prediction algorithm with the query sequence against a coarse-grained force field to predict the three-dimensional model of the query and validate using hidden Markov models for existing sequences.

Vaccine-HLA docking

The ClusPro server (Kozakov et al., 2013; 2017; Vajda et al., 2017) was employed to analyze the docked binding affinity scores between a given ligand and receptor molecules. The server predicts the energy of the structures based on five different stability forces, namely, Van der Waals forces, electrostatic forces, Decoys as Reference States (DARS) energy algorithm, attractive and repulsive forces. The predicted structures are Fast Fourier Transformed to produce the top optimal cluster with the lowest energy at their centers.

Using the ClusPro server, the epitopes were docked to the Human Leukocyte Antigen (HLA) allele structures available in PDB (Tables 1 and 2). The allele structures for docking were prepared by Chimera. Specifically, the heteroatoms (atoms other than carbon) from the HLA structures were removed and the beta-microglobulin supporting structure was retained due to their function of providing stability to the HLA alleles in the host. Hydrogen atoms were added to the structures using the Tools  AddH menu option to create a pre-docking structure. Then several clusters were analyzed to verify the location and binding of the epitope to each HLA allele structure.

The toll-like receptor 4 (TLR4) was used to dock with the epitopes. A 3-dimensional PDB model (PDB ID: 4G8A) of the TLR 4 was used to achieve this. The PDB structure was first cleaned by removing any solvent and non-standard molecules. Then, the variable extracellular domain of the TLR4 (chain B of 4G8A) was extracted and used to dock with the 3-D structures of epitopes.

Vaccine validation using ligand interaction study and molecular dynamics simulation

The Ligplot Tool (Wallace et al., 1995) was employed to visualize the interactions and assess the stability of a protein-ligand structure based on H-bonds and hydrophobic contacts. This tool predicts the stability of each of the docked structures.

The molecular dynamics (MD) simulations were performed using Gromacs version 2019.1 (Abraham et al., 2015) with the gromos43a1 force field file and plotted using ggplot2. MD simulations were run for each of six different HLA allele structures complexed with their original ligands or with the predicted epitopes (see commands in Supplementary Table S1). Use one of the HLA structures 5XOV as an example. MD simulations were run first using the original 5XOV structure containing the HLA-A*24:02 protein, the beta-2-microglobulin supporting structure, and the original ligand the HIV-1 Nef138-10 peptide. After establishing a baseline, simulation was run on the model, which contains the HLA-A*24:02 protein and the beta-2-microglobulin supporting structure docked with the predicted epitope IYPNASLLI. Simulations were run in the same way over 40 ns by first solvating the molecules in a cube of TIP3P water. Then necessary K + or NA- counterions were added to balance the charge of the system to net zero. Once neutral the model underwent energy minimization using the steepest descent minimization algorithm, to ensure that there were no steric clashes or incorrect geometry. After minimization the simulation solvent and ions were equilibrated, first for temperature by heating the system to 300 K during a 100 ps constant volume simulation with a 2 fs time step. Then the system was equilibrated for pressure at 1 atm during a 100 ps simulation with a 2 fs time step. Both the system’s temperature and pressure were regulated using the Berendsen algorithm. Finally, the simulation production parameters were set, specifying a total simulation runtime of 40 nanoseconds, while the temperature and pressure were held constant at 300 K and 1 atm using the v-rescale temperature and Parrinello − Rahman pressure coupling method.

Once MD simulations were run on each allele epitope and allele original ligand bound structure a variety of methods were used to validate the binding interaction the complexes. RMSD and RMSF information was extracted using the gmx_mpi rms and gmx_mpi rmsf commands and plotted with ggplot2 in R. While further binding interaction verifications were performed on the simulations of the allele epitope bound complexes. PDB structure snapshots of the complexes at timepoints of 1, 20, and 40 ns were obtained from the simulations using the Gromacs command gmx_mpi trjconv -dump 1 -pbc nojump and visualized with UCSC Chimera. The Radius of Gyration (Rg) for each of the simulations was calculated with the gmx_mpi gyrate command and plotted with ggplot2 in R. In addition, the free binding energy of each of the simulations was calculated using the g_mmpbsa tool. G_mmpbsa was run using the g_mmpbsa -pdie 2 -pbsa -decomp command along with the required input files and a production file instructing the program to run calculations of both polar and apolar environments without using the WCA model. The accompanying MmPbSaStat.py python script was then used to calculate the free binding energy from the polar and apolar xvg files generated by g_mmpbsa, which was then plotted with ggplot2 in R.

Results

Sequence retrieval and analysis of CEACAM1, CEACAM5, and CEACAM6 proteins

To design an immunogenic multi-epitope vaccine against cancer-related CEACAM molecules such as CEACAM1, CEACAM5, and CEACAM6, a target region that is highly conserved among these molecules needed to be identified. To this aim, the sequences of CEACAM1, CEACAM5, and CEACAM6 were retrieved from the National Center for Biotechnology Information (NCBI) database. Multiple sequence alignment reveals that the N-terminal domain of these molecules is highly conserved with the similarity >91% (Figure 2(a)). Hence, the domain serves as a perfect region for vaccine design. For the sake of simplicity, this N-terminal domain is referred to as the antigen hereinafter.

Figure 2. Multiple Sequence Alignment of CEACAM molecules and the structures of the conserved N-terminal domain. (a) Domains of CEACAM molecules. The motifs of CEACAM1, 5 and 6 are shown schematically. The sequence of the N-terminal domain, Ig-like V set, of CEACAM1 (P13688), CEACAM5 (P06731) and CEACAM6 (P40199) are aligned, together with the sequence retrieved from the 3 D structure of the CEACAM1 N-terminal domain (PDB: 5DZL Chain A). The predicted MHC class I and II restricted epitopes are highlighted. (b) The protein sequence of the antigen (Chain A in 5DZL). The yellow highlighted regions are the MHC I restricted epitopes obtained from the epitope prediction servers. The red highlighted region represents the MHC II restricted epitopes. The overlapping regions have been colored using yellow text. (c) Three-dimensional structure of the antigen (5DZL, Chain A). The MHC I and II restricted epitopes in the antigen structure are colored with the same coloring scheme as (b).

The ProtParam server (Wilkins et al., 1999) was used for the physiochemical analysis of the CEACAM sequences. For instance, the CEACAM1 sequence contains 526 amino acid residues. Its molecular weight is 57560.38 Da and the theoretical pI for the protein is 5.65. It was found that the theoretical half-life of the molecule in mammalian cells is 30 h and the GRAVY (Grand Average of Hydropathy) value is −0.382, which classifies the protein as mildly hydrophilic. This observation is consistent with the fact that CEACAMs attach to the cell membrane and the N-terminal domain of the molecules protrudes into the extracellular matrix (Beauchemin et al., 1999; Beauchemin & Arabzadeh, 2013).

CTL and HTL epitope prediction

CTLs and HTLs are two subsets of T lymphocytes with CD8+ and CD4+ glycoproteins respectively. CD8+ CTLs interact with the antigen-presenting MHC class I molecules, whereas

CD4+ HTLs interact with the antigen-presenting MHC class II molecules. We predicted antigenic epitopes that can interact with CTLs and HTLs using multiple immunoinformatics tools. Tables 1 and 2 present the consensus epitopes identified by different tools, which shall increase confidence in using these epitopes for designing the vaccine. Note that the MHC alleles were selected based on the ranking of the overall scores provided by the IEDB server and subsequently confirmed by other immunoinformatics tools (Table 1 and 2). Several alleles have PDB structures available (Figure S1, supplementary material), which can be used for molecular docking and molecular dynamics simulation studies.

As a result, 7 CTL epitopes were repeatedly predicted by five different methods (Table 1) and 7 HTL epitopes were detected by three tools (Table 2). Several CTL and HTL epitopes are overlapping (Figure 2(b,c)), indicating that these regions are highly detectable by different HLA alleles. Further analysis of the population coverage of HLA alleles showed that HLA-A*24:02 and HLA-A*01:01 are found in 40.40% and 52.80% of the global population, respectively (Table S2, supplementary material), indicating that a vaccine candidate based on these epitopes may be effective for a large human population. Note that the corresponding population data for HLA-DRB (MHC class II) alleles are not available.

Construction of a multi-epitope vaccine and physiochemical properties assessment

To construct a multi-epitope vaccine, the predicted CTL and HTL epitopes were combined with linkers that play a principal role in the functional and structural features of a protein vaccine (Beauchemin & Arabzadeh, 2013). A tandem fusion of these epitopes without proper linkers may result in the generation of a dysfunctional protein with unknown characteristics. In previous studies, AAY has been used as a linker between CTL epitopes for enhancement of epitope presentation whereas GPGPG has been used to link the HTL epitopes (Hajighahramani et al., 2017). These two linkers were used in designing of the vaccine (Figure 3(a)). One of the fragments, TQNDTGFYTLQVIK, is predicted as both a CTL and a HTL epitope (Tables 1 and 2), suggesting that this fragment is recognized by CD8+ and CD4+ cells.

Figure 3. 1-D, 2-D and 3-D structure of a vaccine candidate. (a) The protein sequence of the vaccine candidate. The MHC I restricted epitopes are highlighted in yellow, while MHC II restricted epitopes are highlighted in red. The yellow text in red regions represents the overlapping epitopes. The MHC I and II epitopes are joined by linker sequences GPGPG and AAY. (b) The second structure of the vaccine candidate. The secondary structure of the vaccine was predicted by the SOPMA server (Geourjon & Deleage, 1995). The letters ‘h’, ‘c’ and ‘e’ stand for α-helix, random coil and extended β-strand respectively. (c–f) The 3-D structure of the vaccine candidate predicted by a homology modeling method with red color for MHC II epitopes (c), yellow for MHC II epitopes (d), and green for interferon gamma epitopes (e) predicted by IFNepitope server (Dhanda et al., 2013).

The presence of IFN-γ produced by CD4+ T cells at the site of infection is important to manage neutrophil recruitment and CXC chemokine production (McLoughlin et al., 2008). To confirm that the vaccine candidate can induce IFN-γ, the IFNepitope server was used (Dhanda et al., 2013) and it was found that the vaccine candidate contains multiple IFN-γ inducer epitopes (Table S3, supplementary material).

B-cell epitope mapping, antigenicity, and allergenicity prediction of designed vaccine

B-lymphocytes are the key player in humoral immunity by antibody production. They are also one of the main types of antigen-presenting cells. To identify B-cell epitopes in the vaccine candidate, the BepiPred Linear Epitope Prediction 2.0 web tool (Jespersen et al., 2017) from IEDB was employed. As a result, 7 peptide regions of variable lengths were predicted to be targeted by B cells. Of the 7 predicted regions, 4 were too short (<10) and therefore removed from further consideration. The remaining 3 regions NVAEGKEVLLL, GPGTQNDTGFYTL, and PNASLLIQNVTQNDTGFYTL were selected not only because they are longer (>10 residues), but also because these regions have a smoothed B-cell epitope likelihood score above the threshold of 0.5 (Jespersen et al., 2017). The binding regions were predicted based on the number and concentration of beta-turns in the sequence calculated according to the Chou and Fasman algorithm (Chou & Fasman, 1979). Successful identification of B-cell epitopes in the peptide vaccine indicates that it may have the ability to enhance humoral immunity as well as cell-mediated immunity.

To examine whether the designed vaccine is immunogenic in nature, its antigenicity was determined by using the ANTIGENpro server (Magnan et al., 2010). It was found that the vaccine has an antigenicity probability of 0.63. Note that the predicted probability of antigenicity score is between 0 and 1, where a higher value indicates a greater likelihood that the input sequence is antigenic. The antigenicity score obtained for this vaccine highlights its antigenic nature and this value exceeds the desired antigenicity value threshold of a 0.6 (Jain et al., 2019).

The designed vaccine was further examined to determine if it is a potential allergen. The AllerTOP online server (Castro et al., 2018) was used to determine its allergenicity. It was found that the vaccine sequence is likely to be a non-allergen as it is most closely related to a nonallergic sequence in its database (UniProt ID: Q8WVR3). This result indicates that the designed vaccine is nonallergic in nature and probably safe for human use.

To characterize the immunogenicity and immune response profile of the designed vaccine, in silico immune simulations were conducted using the C-ImmSim server (Rapin et al., 2010). The levels of antibodies are not elevated in the primary response. The secondary and tertiary responses are characterized by marked increases in levels of IgM, IgG + IgM and IgG₁ + IgG₂ and B-cell populations (Figure S2a, b, supplementary material). This profile indicates the development of immune memory and subsequent clearance of the antigen. A similar pattern is also seen in T_H (helper) and T_c (cytotoxic) cell populations with corresponding memory development (Figure S2C, D, supplementary material). These results suggest that the designed vaccine likely induce immune reactions as evidenced by a marked increase in the generation of secondary responses.

Vaccine structure prediction, refinement, and assessment

The secondary structure of the vaccine candidate was predicted by the SOPMA server (Geourjon & Deleage, 1995), in which most of the structure is covered by β sheets (‘e’) and coils (‘c’) (Figure 3(b)). To model the 3-dimensional structure of the vaccine candidate, predictions were performed using both homology modeling and ab initio modeling methods. The antigen structure 5DZL was used as the template because the vaccine sequence has 57.3% identity with the antigen sequence (Figure S3, supplementary material). Consistent with the predicted secondary structure (Figure 3(b)), the 3-D model is characterized by multiple β sheets and coils, in which the CTL epitopes (Figure 3(c)), HTL epitopes (Figure 3(d)) and IFN-γ epitopes (Figure 3(e)) are highlighted.

The predicted 3-D models were refined by GalaxyRefine (Ko et al., 2012; Shin et al., 2014). Table S4 (supplementary material) presents the top five refined models based on the original homology model. Ramachandran plots were used to illustrate the effect of the refinement. Before the refinement, the homology modeling model showed 93.9% of the residues in the favored regions (Figure 4(a)). After the refinement, the numbers were increased to 100% (Figure 4(b), Table S4, supplementary material), indicating that the quality of the structures was greatly improved after the refinement. Similar improvement has been seen for original ab initio model (Table S5, supplementary material).

Figure 4. Ramachandran plots before and after refinement. Ramachandran plots of the vaccine candidate structure predicted by homology modeling method (a) and after refinement by GalaxyRefine refinement (b).

A detailed examination of other parameters such as MolProbity, clash and poor rotamer (Table S4, supplementary material) revealed that the Model 3 of the homology modeling structures has the best quality with the lowest MolProbity score, clash score, no poor rotamer and 100% of residues occurring in the favored regions in the Ramachandran plot (Table S4, supplementary material). Similarly, the Model 2 of the ab initio modeling structures has the best quality with the lowest MolProbity score, clash score, no poor rotamer and 95.3% of residues occurring in the favored regions in the Ramachandran plot (Table S4, supplementary material). The superposition of these two structures showed a high similarity with the RMSD value of 4.45 Å (Figure S4, supplementary material), which further enhances confidence in the 3-D structure of the vaccine. As such, the Model 3 structure was selected for molecular docking studies.

Molecular docking of epitopes with HLA structures

To understand if the epitopes properly interact with MHC molecules, the interactions of these molecules with the original ligands contained in their PDB structures was first examined. For the six HLA complex structures available in PDB, all ligands are bound to the antigen-binding groove of MHC molecules (Figure S1, supplementary material). These results suggest that the epitopes should appear in the same location for proper interactions with MHC molecules.

The docking of the epitopes to the corresponding HLA allele structures was performed on the ClusPro Server and 39 docked epitope-HLA models were generated. Among them, only the models with the lowest energy score were selected. Detailed analysis of these models showed that all epitopes are bound to the antigen-binding groove (Figure 5(a)–(f)). The lowest energy scores of epitope-HLA models are comparable (supplementary material, Table S6), with the minimal score being −684.1 kCal/mol, indicating that the epitopes interact favorably with the HLA structures.

Figure 5. Docking of epitopes to HLA structures. The epitopes (red) were docked to six HLA 3 D structures available in PDB using Boston University’s ClusPro server. The six complex structures found in PDB include four MHC Class I molecules (a–d, see Table 1) and two MHC Class II molecules (e and f, see Table 2). The epitopes were modeled by PEP-Fold 3.0 that uses ab initio modeling.

As an illustration, a detailed examination of HLA-epitope interactions was performed on HLA-A*24:02 (5XOV) and the corresponding epitope ‘IYPNASLLI’ using the Ligplot Tool (Wallace et al., 1995). It was found that the epitope has numerous hydrogen bonds and hydrophobic interactions with the HLA molecule. Specifically, the epitope region has 6 H-bond interactions where all H-bond length varies from 2.5 − 3.5 Å and multiple hydrophobic interactions between various amino acid residues from the allele (Figure 6). Notably, this particular epitope was predicted to interact with HLA-A*24:02 by several immunoinformatics tools (Table 1). In other words, the docking data confirm the consensus predictions of the immunoinformatics methods (Table 1), illustrating the spatial feasibility of interactions between the epitope regions of the designed vaccine and designated MHC molecules.

Figure 6. Interactions between HLA-A*24:02 allele structure and its epitope IYPNASLLI. The residues in blue color represent the epitope from the vaccine, while the green colored residues represent part of the HLA receptor. The bonds represented with green dashed lines show the hydrogen bond and the relative distance is represented in Angstroms. The comb-like residues are the hydrophobic patches found on the receptor as well as the residue atoms of the vaccine.

Toll-like receptors (TLRs) are key components of the innate immune system, recognizing a variety of microbial products (Medzhitov, 2001). MoreoverTLRs play an important role in tumor progression (Shcheblyakov et al., 2010). It has been reported that the ТLR4 agonists lipopolysaccharide (LPS) from Gram-negative bacteria possess high antitumor activity, when administered intra-tumorally (Okamoto et al., 2006). Moreover, it was found that LPS binds to the central domain of TLR4 protein (Ain et al., 2020). To understand if TLR4 can recognize the six epitopes in the vaccine, the epitopes were docked to the central domain of TLR4 (the B chain of 4G8A) using ClusPro. Top-ranked docking structures have epitopes located on the central domain of TLR4 (Figure S5, Table S7; supplementary material), indicating that TLR4 interacts with the epitopes similar to LPS. These results suggest a potential role of innate immunity against the cancer vaccine.

Molecular dynamics simulation for vaccine-HLA complex

To further study the stability of the epitope-bound HLA structures, molecular dynamics (MD) simulations were performed in GROMACS (Abraham et al., 2015) to compare the stability of the epitope-HLA complexes to that of the HLA complexes bound with their original ligands. The stability was measured by RMSD and RMSF values. RMSD computed along a trajectory is the RMSD averaged over atoms as a function of time, while RMSF computed along a trajectory is the RMSF averaged over time as a function of individual atoms. The RMSD of each complex was plotted over time in nanoseconds. The backbone RMSD value of the original ligand-bound complex averaged to 0.2 − 0.4 Å (red lines in Figure 7(a–f)), while the backbone RMSD value of the epitope bound complex averaged to 0.4 − 0.8 Å (blue lines in Figure 7(a-f)). Note that for both the epitope-docked and original ligand-docked simulations, the RMSD values level out over time, which indicates that the epitope-bound complexes are stable in solution.

Figure 7. Molecular Dynamics represented through RMSD Plots. The red-colored lines plot RMSD values of the HLA molecule from six complex structures bound with the original peptide ligand found in PDB, while the blue-colored lines plot the RMSD values of the HLA molecule in the complex bound with the predicted epitopes. The six complex structures found in PDB include four MHC Class I antigens (a–d, see Table 1) and two MHC Class II antigens (e and f, see Table 2).

In addition, the RMSF values of both the HLA molecule complexed with the original ligand (red lines in Figure 8(a–f)) and the HLA molecule complexed with the epitopes (blue lines in Figure 8(a–f)) were plotted. The results showed a lack of significant RMSF fluctuations between the two complexes. This lack of significant variation suggests that the epitope-HLA complexes have similar stability as the original HLA complexes and that the binding of the epitopes seems not destabilizing to the HLA molecules.

Figure 8. Molecular Dynamics represented through RMSF Plots. The red-colored lines plot RMSF values of the HLA molecules from six complex structures bound with the original peptide ligands found in PDB, while the blue-colored lines plot RMSF values of the HLA molecules in the complexes bound with the predicted epitopes. The six complex structures contain four MHC Class I antigens (a–d, see Table 1) and two MHC Class II antigens (e and f, see Table 2).

PDB structure snapshots of the vaccine peptide complex MD simulations at 1, 20, and 40 ns showed the peptide’s continued interaction and localization to the known ligand binding domain (Figure S6, supplementary material). This result indicates that the epitopes remain stable at the docked site. The Radius of Gyration (Rg) of a complex is a measure of the compactness and can be used to indicate complex stability. Rg was calculated for all 6 epitope bound structures, 3bo8, 3rl1, 4hx1, 5xov, 6biy, and 6cpn at each time step in the MD simulation Figure S7(a–f), supplementary material). The average Rg for each of the simulations was 2.320, 2.325, 2.159, 2.176, 2.392, and 2.363 nm, respectively. The timestep plots showed that the Rg values fluctuate closely to the average Rg values indicating that the overall shape of the protein complex is stable after binding with the epitopes (Yadav et al., 2018). The g_mmpbsa tool was used to calculate the free binding energy of the six MD simulations at each simulation timestep (Figure S8(a–f), supplementary material). The average free binding energies for each simulation were also calculated as −16.658 ± 54.984, −449.296 ± 49.972, −184.390 ± 38.207, −28.661 ± 33.628, 307.072 ± 61.670, and −160.319 ± 77.177 kJ/mol, respectively. The negative free binding energy values observed indicate that five of the six predicted epitopes are strongly bound with their HLA receptor alleles making them promising candidates in a cancer vaccine (Ahmad et al., 2017).

Discussion

The goal of this project was to design a vaccine construct targeting the highly conserved N-terminal domain of the cancer related CEACAMs using a suite of immunoinformatics and molecular modeling tools. This vaccine, if successful, potentially has a universal antitumor effect on all cancers that overexpress CEACAM proteins.

To this goal, six epitopes have been identified from the domain using immunoinformatics tools, which are predicted to interact with MHC class I alleles in CTLs and MHC class II alleles in HTLs. The predicted CTL and HTL epitopes have the lengths of 9-mer (Table 1) and 15-mer (Table 2) respectively. These specific lengths of epitopes were selected based on prior studies on the typical distribution of peptides presented by MHC I and II molecules. Based on the previous study (Bettencourt et al., 2020), MHC-I peptides have a narrow distribution with a predominant peak at 9 mer. The selection of 9-mer epitopes (Table 1) has exactly this length. By contrast, MHC-II peptides show a wide distribution from 12 mer to 18 mer, which is consistent with 11-19 mer from another study (Barra et al., 2018). The selection of 15-mer epitopes (Table 2) is the center of the length distribution.

The selection of the alleles was based on the ranking of the overall scores from the IEDB server. The allele with the highest score is HLA-A*24:02 (Table 1). Interestingly, HLA-A*24:02-restricted CTLs recognize antigens of leukemia (Tawara et al., 2017), lung cancer (Yamada et al., 2003), and stomach cancer (Murahashi et al., 2016), indicating that HLA-A*24:02 is one of the HLA alleles that are critical for cancer antigen recognition. Since CEACAM molecules, especially CEACAM6, are overexpressed in the aforementioned cancers (Hammarstrom, 1999), it is highly likely that HLA-A*24:02 also recognize epitopes from CEACAMs, as predicted in this study (Table 1).

These epitopes were then linked by AAY and GPGPG motifs for proper separation and presentation of the epitopes to the host immune system. AAY and GPGPG linkers are purposely selected and used in the vaccine construct. AAY is used to link together CTL (MHC I) epitopes to enhance epitope presentation. That is, the vaccine is cleaved by the proteasomal and lysosomal degradation systems after the AAY motif in cytoplasm. Then the generated C-terminal of vaccine binds to the transporter associated with antigen processing (TAP) protein complex, which delivers the epitopes to the endoplasmic reticulum (ER), where they bind to nascent MHC I molecules (Bergmann et al., 1996). Thus, binding of epitopes to the TAP transporter, with the help of the AAY linker, is vital for presenting them to MHC I molecules. On the other hand, the GPGPG linker was used because the vaccine contains HTL (MHC II) epitopes and the linker can stimulate HTL response and conserves conformational dependent immunogenicity of the epitopes (Livingston et al., 2002).

The interactions between epitopes and HLA structures were shown by molecular docking and molecular dynamics simulation. It was found that the epitopes can bind to the antigen-binding groove of the MHC molecules and this binding is stable over time in solution. These computational results suggest that the designed vaccine candidate has great potential to become an effective vaccine. Although the candidate has a predicted antigenicity value (0.63) that is higher than the threshold (0.6), follow-up experiments are needed to test if the vaccine is immunogenic in humans. One way to enhance the immunogenicity of the vaccine candidate is to use suitable adjuvants (see below).

One of the most important drawbacks of the subunit vaccine is their relatively low immunogenicity compared to whole cell inactivated or live attenuated vaccines (Khalaj-Hedayati et al., 2020). To solve this problem, adjuvants are often used in the design of a subunit vaccine. Adjuvants have been widely deployed to further increase the effectiveness of vaccines. They increase vaccine effectiveness in several ways from increasing immune system response to aiding in vaccine transport and increasing the duration of its exposure (Temizoz et al., 2016). Vaccine delivery-based adjuvants operate on the principle of creating a deposit of the antigenic compound at the vaccination site that slowly released over time prolonging immune response (Guy, 2007). The water in oil emulsion adjuvant, Montanide ISA-51, has shown promise activating T-cell response in cancer vaccines clinical trials that correlate positively with increased patient survival (Chianese-Bullock et al., 2005;; Neninger Vinageras et al., 2008) . Any clinical trial attempting to further test or derive a cancer vaccine from this work is recommended to consider including an appropriate immunity-enhancing adjuvant such as Montanide ISA-51.

Several CEACAMs such as CEACAM5 and CEACAM6 are highly expressed in primary and metastatic cancers in breast, pancreas, colon, and lung (Blumenthal et al., 2007), which make them potential targets for cancer vaccines. However, CEACAMs are self-antigens that make it difficult to break immune tolerance when part of them is used in a vaccine. Several approaches have been proposed and tested to break the immune tolerance of cancer self-antigens (Makkouk & Weiner, 2015). Additional experiments are required for further development and refinement of the vaccine candidate to break immune tolerance prior to clinical trials.

Conclusions

CEACAMs belong to a family of membrane glycoproteins. Several family members such as CEACAM1, 5 and 6 are highly expressed in a variety of cancers, and thereby they serve as targets for cancer vaccines. Here, a peptide-based vaccine candidate is developed targeting the conserved N-terminal domain of these molecules utilizing a suite of immunoinformatic and molecular modeling tools against cancer-associated CEACAMs. The vaccine candidate contains epitopes predicted to bind to MHC class I and II molecules and has high antigenicity and low allergenicity. The epitopes can bind to the antigen-binding groove of MHC molecules and such binding is stable over time in solution. More experiments are needed to test the immunogenicity of the candidate in humans. Suitable adjuvants can be used for further enhancement of the vaccine. Appropriate strategies to break immune tolerance should be considered prior to clinical trials.

Authors’ contributions

AG and AJR designed and performed the experiments. AG and AJR and FC drafted and edited the manuscript. All authors read and approved the final manuscript.

Disclosure statement

The authors declare that they have no conflicts of interest.

Additional information

Funding

We are thankful for the research support from NIH (grant no. R15GM116102) and the College of Science (F.E.A.D funding) of the Rochester Institute of Technology.