Method for screening candidate CTL epitopes based on combined free energy calculation

Abstract

The invention discloses a method for screening candidate CTL epitopes based on free energy calculation. The method comprises the following steps: constructing a three-dimensional structure of a target MHC I molecule as an initial compound model of pMHC I through homologous modeling; respectively mutating residues of at least one preset anchoring site of polypeptide in the initial compound model into a plurality of natural amino acids, performing molecular dynamics simulation on compounds in the compound model set, and quantifying an energy contribution value of each amino acid type on the preset anchoring site; based on the energy contribution value, constructing a quantitative binding motif of the target MHC I molecule, the quantitative binding motif comprising amino acid types preferred by different anchoring sites and corresponding binding free energy weights; and screening out candidate polypeptides conforming to motif characteristics as CTL epitopes by using the quantitative binding motifs. According to the method, very few high-potential candidate epitopes can be quickly locked from massive antigen sequences, and the cost and time of subsequent experimental verification are greatly reduced.

Description

Technical Field

[0001] This invention relates to the fields of immunology and computational biology, and particularly to a method for screening candidate CTL epitopes based on binding free energy calculation. Background Technology

[0002] Cytotoxic T lymphocytes (CTLs) play a central role in adaptive immunity, precisely eliminating infected or cancerous cells by recognizing antigenic peptides (CTL epitopes) presented on the cell surface by Major Histocompatibility Complex Class I (MHC I) molecules. Therefore, accurate identification of CTL epitopes has become a crucial prerequisite for the development of novel vaccines, tumor immunotherapies, and immunodiagnostic reagents. However, given the vast genomes of pathogens and tumors, relying solely on experimental techniques for epitope screening presents significant challenges due to the enormous workload, lengthy process, and high costs. This necessitates the development of efficient and reliable computer-aided screening methods, a pressing need in this field. This challenge is particularly prominent in veterinary science and animal disease control. For example, developing effective T-cell vaccines against important animal pathogens such as porcine reproductive and respiratory syndrome virus (PRRSV), avian influenza virus, or canine and feline coronaviruses is key to controlling disease transmission. However, research on the polymorphism of animal MHC I molecules (such as DLA in dogs, FLA in cats, and SLA in pigs) and corresponding epitope data are extremely scarce compared to human HLA molecules, which seriously restricts the development of animal-specific cell-based immune vaccines.

[0003] Currently, the development of this technology field presents a situation where multiple paths coexist. While experimental sequencing-based binding peptidomics technology can provide the most direct and reliable evidence, its throughput is limited, its cost is high, and it is insensitive to low-abundance peptides, making its application even more difficult in animal studies where resources are limited and model systems are imperfect. Subsequently, bioinformatics prediction methods based on sequence characteristics have emerged. These tools analyze known binding peptide sequence patterns using machine learning algorithms, achieving extremely high prediction throughput. However, their prediction accuracy is severely limited by the scale and quality of training data. For animal MHC type I, where data is already scarce, the reliability of prediction results is often significantly reduced, or even impossible to predict effectively. More importantly, these models are often black boxes, making it difficult to accurately quantify binding affinity, resulting in a high false positive rate in screening results. Ultimately, extensive experimental verification is still required, failing to fundamentally reduce research and development costs and timelines.

[0004] To overcome the precision limitations of sequence-based methods, structure-based computational simulations have emerged. Among these, molecular dynamics simulations combined with free energy calculation techniques such as MM-PBSA can simulate the dynamic behavior of protein-peptide complexes at the atomic level and accurately calculate their binding free energies, providing unparalleled precision and physical insight. This method, in principle, does not rely on known binding data and has theoretical application potential for animal MHC I molecules lacking training data. However, this high precision comes at the cost of extremely high computational costs. Particularly critical is that systematically exploring the binding preferences of an MHC I molecule (whether human or animal) requires scanning all possible combinations of the peptide sequence, leading to an unsolvable combinatorial explosion problem. For a typical nonapeptide, the possible sequence combinations reach hundreds of billions, making this high-precision method computationally infeasible. This limits its application to post-hoc validation of a very small number of candidate molecules, preventing its use for pre-hoc systematic motif identification and large-scale screening. This constitutes an insurmountable obstacle in data-starved animal MHC research. Summary of the Invention

[0005] This invention provides a method for screening candidate CTL epitopes based on binding free energy calculation, which solves the problem of insufficient accuracy due to data scarcity in the prior art. It realizes the quantitative analysis of the interaction between peptides and MHC I molecules through molecular dynamics simulation and MM-PBSA free energy calculation, thereby constructing their binding motifs and achieving high-throughput and high-precision screening of candidate cytotoxic T lymphocyte epitopes.

[0006] This invention includes the following steps: (a) The three-dimensional structure of the target MHC I molecule was constructed by homology modeling as the initial complex model of pMHC I; (b) For at least one preset anchoring site residue of the peptide in the initial complex model, mutate it into a variety of natural amino acids, and set the non-target site of the peptide to a non-side chain amino acid - alanine, thereby generating a complex model set of single-point mutated pMHC I. (c) Perform molecular dynamics simulations on the complexes in the complex model set, calculate the binding free energy of each mutant peptide to MHC I molecule using the MM-PBSA method, and simultaneously decompose the binding free energy to quantify the energy contribution value of each amino acid type at the preset anchoring site. (d) Based on the energy contribution value, construct the quantitative binding motif of the target MHC I molecule, wherein the quantitative binding motif includes the amino acid type preferred by different anchor sites and their corresponding binding free energy weights; (e) Using the quantified binding motif, candidate protein sequences are scanned, and candidate peptides that meet the motif characteristics are selected as CTL epitopes.

[0007] Furthermore, in step (a), when obtaining the three-dimensional structure of the target MHC I molecule, if there is an experimentally resolved crystal structure, the crystal structure is directly used as the initial model. If no experimentally resolved crystal structure exists, the three-dimensional structure of the target MHC I molecule is constructed by artificially mutagenesis using the crystal structure of a highly homologous MHC I molecule as a template and homology modeling software.

[0008] Further, in step (b), the preset anchoring point includes one or more of the 2nd, 3rd, 5th, and 9th positions of the nonapeptide; the non-targeting site is set to alanine.

[0009] Furthermore, in step (b), mutation operations are performed independently for each preset anchor site. Each mutation changes only the amino acid type of one anchor site, while other anchor sites and non-target sites remain alanine, keeping the polypeptide backbone conformation unchanged. Only the side chain is replaced with an alanine side chain, without involving combined mutations of two or more sites.

[0010] Further, in step (c), the duration of the molecular dynamics simulation is not less than 50 nanoseconds; the energy decomposition is used to obtain the energy contribution of van der Waals interactions, electrostatic interactions and solvation effects for each anchored residue.

[0011] Further, the quantitative binding motif described in step (d) is represented as an energy weight matrix, where the rows of the matrix represent the anchoring sites of the polypeptide, the columns represent 20 natural amino acids, and the matrix element values ​​are the average energy contribution value of the corresponding amino acid at the corresponding site or the score derived therefrom.

[0012] Further, the quantification of binding motifs in step (d) is represented by an energy weight matrix. The matrix represents the anchor sites of the peptide, with 20 natural amino acids listed. The matrix element values ​​are the relative binding preference index (RBPI) of the corresponding amino acid at the corresponding site. The formula for the RBPI is: E(i,j)=μ[ΔG contrib (i,j)] - μ[ΔG contrib (i,j)] Where μ[ ] represents the arithmetic mean of the energy contribution values ​​from at least 3 repeated simulations, i is the site, and j is the amino acid.

[0013] Further, the amino acid sequences of the candidate CTL epitope peptides obtained in step (e) are shown in SEQ ID NO.1-4.

[0014] One or more technical solutions provided in the embodiments of the present invention have at least the following technical effects or advantages: The energy weight matrix constructed in this invention is directly derived from physical force fields and free energy calculations. Its results have a clear physicochemical interpretation and can accurately quantify the energy contribution of amino acid residues. Its prediction accuracy is far higher than that of black-box bioinformatics models based on statistical laws. It is particularly suitable for MHC I alleles with extremely scarce training data, providing a powerful and universal tool for the rapid development of animal disease vaccines.

[0015] This invention compresses the enormous computational load to a controllable range, and combined with a rapid virtual screening step, it can quickly identify a very small number of high-potential candidate epitopes from a massive number of antigen sequences, greatly reducing the cost and time of subsequent experimental verification. Attached Figure Description

[0016] Figure 1 The flowchart shows the method for screening candidate CTL epitopes based on combined free energy calculations. Figure 2 The conformational stability results of peptides corresponding to the 20 natural amino acid residue mutation models for the nonapeptide anchoring sites (P2, P3, and P9) of the feline FLA-H*0401 molecule are shown in Figures A, B, and C, respectively, representing the RMSD trajectory changes of the bound peptides in the 20 natural amino acid residue mutation models corresponding to P2, P3, and P9 during a 50 ns molecular dynamics simulation. Figure 3 The MM-PBSA binding free energy decomposition and sequencing results of 20 natural amino acid residues at the nonapeptide anchoring sites (P2, P3, P9) of the cat FLA-H*0401 molecule are shown in Figures A, B, and C, respectively. Figure 4 Predicted binding motifs for the nonapeptide of the feline FLA-H*0401 molecule; Figure 5 The results of expression and inclusion body purification of feline FLA-H*0401 heavy chain and β2m light chain recombinant protein are shown in Figures A and B, respectively, which are the SDS-PAGE electrophoresis results of FLA-H*0401 heavy chain and β2m light chain recombinant protein expression and inclusion body purification. Figure 6 To verify the binding of VII9 peptide to FLA-H*0401 molecule in vitro refolding experiments; A shows the results of gel chromatography and the corresponding SDS-PAGE electrophoresis; B shows the results of anion exchange chromatography and the corresponding SDS-PAGE electrophoresis. Figure 7To verify the binding of YII9 peptide to FLA-H*0401 molecule in vitro refolding experiments; A shows the results of gel chromatography and the corresponding SDS-PAGE electrophoresis; B shows the results of anion exchange chromatography and the corresponding SDS-PAGE electrophoresis. Figure 8 To verify the binding of SLY9 peptide to FLA-H*0401 molecule in vitro refolding experiments; A shows the results of gel chromatography and the corresponding SDS-PAGE electrophoresis; B shows the results of anion exchange chromatography and the corresponding SDS-PAGE electrophoresis. Figure 9 To verify the binding of VLI9 peptide to FLA-H*0401 molecule in vitro refolding experiments; A shows the results of gel chromatography and the corresponding SDS-PAGE electrophoresis; B shows the results of anion exchange chromatography and the corresponding SDS-PAGE electrophoresis. Detailed Implementation

[0017] This invention provides a method for screening candidate CTL epitopes based on binding free energy calculation. To better understand the above technical solutions, the following will provide a detailed explanation of the technical solutions in conjunction with the accompanying drawings and specific implementation methods.

[0018] The preferred embodiments of the present invention will now be described in detail with reference to specific examples. It should be understood that the following examples are given for illustrative purposes only and are not intended to limit the scope of the invention. Those skilled in the art can make various modifications and substitutions to the present invention without departing from its spirit and essence.

[0019] like Figure 1 As shown, (1) Initial structure preparation steps: The three-dimensional structure of the target MHC I molecule is obtained. If an experimentally resolved pMHC complex crystal structure exists for this molecule, it is directly used as the initial structure. Otherwise, a known MHC I molecule crystal structure with high homology is selected as a template, and a three-dimensional structural model of the target molecule is constructed using homology modeling software (such as Modeller) through artificial mutation. This step ensures the universality of the method and is not limited by whether the structure of a specific MHC I molecule has been experimentally resolved.

[0020] (2) Steps for constructing a systematic single-point mutation model library: Based on the initial structure prepared earlier, the peptide sequence is used as a template. First, all non-critical anchor sites of the peptide are uniformly mutated to alanine (Ala), a simple structure with inert side chains, thereby fixing the skeletal structure of the peptide and eliminating interference from non-critical sites. Subsequently, for each pre-defined critical anchor site (typically P2, P3, P5, P9, etc. for a typical nonapeptide), it is independently and systematically mutated sequentially to 20 natural amino acids. During this process, all other sites (including other anchor sites) are strictly kept to be alanine. Using this strategy, for a peptide containing n anchor sites, only 20 × n pMHC complex models need to be constructed, instead of 20 ^n For example, for four anchor points, only 80 models need to be built, thus reducing the computational load by several orders of magnitude and making systematic high-precision calculations possible.

[0021] (3) Molecular dynamics simulation and MM-PBSA energy calculation and decomposition steps: For each single-point mutation pMHC I complex model constructed in step (2), molecular dynamics simulations were performed in an explicit solvent environment. To ensure that the system reaches sufficient equilibrium and obtains convergent statistical properties, the simulation time for each system was no less than 50 nanoseconds. To evaluate the convergence of the results, at least three repeatable simulations were performed for each independent complex system; each repetition used the same simulation parameters but different initial atomic velocity random seeds.

[0022] Based on the equilibrium molecular dynamics trajectory, the total binding free energy (ΔG_bind) of each mutant peptide with the MHC I molecule was calculated using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method. Furthermore, the total binding free energy was decomposed at the residue level to precisely quantify the van der Waals energy, electrostatic interaction energy, polar solvation energy, and nonpolar solvation energy contributed by the amino acid residues at the mutant site in their interaction with the MHC I molecule's binding pocket. This energy decomposition process was performed separately for each repeated simulation trajectory, and the results were recorded.

[0023] (4) Steps for constructing the quantified sequence "energy weight matrix": Based on all the energy decomposition data obtained in step (3) after error evaluation, a quantized binding motif is constructed, which is presented in the form of an energy weight matrix.

[0024] Matrix construction method: The specific anchoring sites of the peptides are used as rows of the matrix. The 20 natural amino acids are used as columns. Each element in the matrix (E...) i,jThe ) represents a quantified binding preference value, called the Relative Binding Preference Index (RBPI). This index is calculated using the following formula:

[0025] Where i represents the anchoring residue site; j represents 19 natural amino acid residues other than alanine; μ[ ] represents the arithmetic mean of the energy contribution values ​​obtained from multiple repeated simulations. This calculation method effectively eliminates background differences between different systems, allowing the preferences of different amino acids to be compared under a unified standard. To enhance the scientific rigor and reliability of the results, each RBPI value (E i,j An uncertainty metric (such as standard deviation) should be associated. This standard deviation is calculated using an error propagation formula from the standard deviations of repeated simulations of each amino acid mutant and the reference alanine mutant. The constructed energy weight matrix is ​​the quantified binding motif defined in this invention. This matrix provides direct physical insight into the binding preference of MHC molecules: a negative RBPI value indicates that the amino acid is more favorable for binding at the corresponding anchor site compared to alanine; the more negative the value, the stronger the preference. A positive RBPI value indicates that the amino acid is unfavorable for binding.

[0026] (5) Steps for screening candidate CTL epitopes: The quantized binding motif (energy weight matrix) constructed based on the aforementioned steps includes the following steps: Generate a candidate peptide library: Bioinformatics analysis is performed on the whole genome sequence of the pathogen to be screened, translating all possible open reading frames to obtain all its antigenic protein sequences. Subsequently, using a sliding window method, based on the peptide length (e.g., 9 amino acids) corresponding to the binding motif of the target MHC class I molecule, each antigenic protein sequence is systematically traversed to generate a set of overlapping peptide candidate libraries with specified lengths covering all sequences (e.g., for a 300-amino acid protein, nonapeptide sequences from 1-9, 2-10, 3-11... up to 292-300 will be generated).

[0027] Define the optimal anchoring mode: From the "energy weight matrix" constructed in step (4), extract the top 5 amino acids with the highest relative binding preference index for each key anchoring point; Pattern matching filtering: Each nonapeptide in the candidate nonapeptide library is matched with the "preferred amino acid set" defined above for each anchor site. Only nonapeptides whose amino acids at each key anchor site belong to the preferred amino acid set of the corresponding site are retained, forming candidate CTL epitopes.

[0028] (6) In vitro experimental verification steps To verify the reliability of the computational screening, the candidate peptides screened in step (5) were chemically synthesized. Subsequently, using in vitro refolding assembly technology, the synthesized peptides, the heavy chain inclusion body protein of the target MHC I molecule, and the β2-microglobulin (β2m) inclusion body protein were refolded together to form the correct pMHC I ternary complex. The complex was separated by gel filtration chromatography and ion exchange chromatography, and the target peak was detected by SDS-PAGE protein electrophoresis, thus directly and quickly confirming whether the candidate peptide could be effectively bound by the MHC I molecule, completing the process from "… in silico "arrive" in vitro "Closed-loop verification of "

[0029] Example 1: Identification of the peptide-binding motif of the FLA-H*0401 allele in domestic cats and screening of candidate CTL epitopes. (1) Model building and simulation calculation First, we obtained the amino acid sequence of the extracellular region of FLA-H*0401 (GenBank ID: NP_001041626). Using the resolved high-resolution crystal structure (PDB ID: 5XMM) of the typical nonapeptide-binding complex FLA-E*01801 as a template, we performed site-directed mutagenesis using the Modeller tool to make its sequence identical to that of FLA-H*0401. Simultaneously, P2, P3, and P9 sites were defined as key anchoring sites for FLA-H*0401 binding to the nonapeptide.

[0030] For each anchor site (P2 as an example), we mutated it to 20 natural amino acids. While performing single-site mutations, the other three anchor sites (P3 and P9) were fixed to alanine (Ala), thus constructing 20 different pMHC complex models. This operation was repeated for P3 and P9 sites, resulting in a total of 60 single-point mutation models.

[0031] Molecular dynamics simulations were performed on each of the 60 model systems described above. The Gromacs AMBERGS force field and the SPC / E aqueous solvent model were selected to obtain the protein topology file. A periodic cubic box was defined under the aforementioned force field, and an aqueous solvent model was added using the Solvate module in Gromacs software, while the Genion program was used to add Na+ and Cl- equilibrium systems to achieve an electrically neutral state. Before the formal molecular dynamics (MD) simulations, the systems were optimized using the steepest descent (SD) method for 5000 steps to eliminate unreasonable van der Waals position conflicts. Then, positional constraints were added to the ligands in the systems, and isothermal-isochoric ensembles (NVT) and isothermal-isobaric ensembles (NPT) were used to achieve equilibration at 300 K using the Parrinello-Rahman method for 1 ns and 5 ns, respectively. Finally, unrestricted molecular dynamics simulations for more than 50 ns were performed under pre-equilibration conditions. Different initial atomic velocity random seeds were used for each simulation to ensure statistical independence of the sampling. After the simulations were completed, the molecular dynamic trajectories of all 60 model systems were analyzed to evaluate the stability and convergence of the simulations. The root mean square deviation (RMSD) of each peptide in the system over the entire simulation timescale was calculated to determine whether the systems had reached equilibrium. Results are shown below. Figure 2 As shown in Figures A to 2C, the RMSD values ​​of all systems tend to stabilize relatively in the later stages of the simulation, indicating that the system has reached full equilibrium.

[0032] (2) Energy calculation and decomposition Based on the equilibrium trajectory of the last 50 ns of each simulated trajectory, the binding free energy was calculated using the MM-PBSA method, and residue energy decomposition was further performed to accurately extract the energy contribution value (ΔG_contrib) of amino acid residues at each mutation site and sort them. The results are shown in […]. Figure 3 A to 3C.

[0033] (3) Calculation of relative binding preference index and prediction of peptide binding motifs Taking the P2 site as an example, the calculation process of RBPI is demonstrated. For the 20 amino acids at the P2 site, we first calculate the mean of their energy contribution. For example, for P2-valine (Val): (a) From the three repeated simulations, obtain three ΔG_contrib(P2, Val) values, calculate their arithmetic mean, and denot it as μ[ΔG_contrib(P2, Val)]; (b) Similarly, for the reference system P2-alanine (Ala), its mean μ[ΔG_contrib(P2, Ala)] is calculated; (c) Subsequently, the relative associative preference index of P2-Val is calculated according to the formula: E(P2, Val) = μ[ΔG_contrib(P2, Val)] - μ[ΔG_contrib(P2, Ala)] Simultaneously calculate the standard deviation (σ) of this RBPI value. For example, σ[ΔG_contrib(P2, Val)] and σ[ΔG_contrib(P2, Ala)] are the standard deviations of its three repetitions, respectively, calculated using the error propagation formula σ_RBPI = sqrt(σ_Val). 2 + σ_Ala 2 The uncertainty measure of RBPI(P2, Val) was obtained. Based on the above method, the relative binding preference index (RBPI) of FLA-H*0401 molecule at anchor sites P2, P3, and P9 was calculated and converted into peptide binding motifs. The results are shown in […]. Figure 4 .

[0034] (5) Screening of FLA-H*0401 restricted candidate CTL epitopes: (a) Generation of candidate nonapeptide libraries The complete genome sequences of feline parvovirus FPV (AMPV2020) (GenBank: MZ712026.1) and feline herpesvirus FAV-1 (PHIL_01) (GenBank: MH070336.1) were obtained from the NCBI database. All encoded proteins in the genome annotations of the two viruses were selected as candidate antigens, and their full-length sequences were obtained.

[0035] Table 1 Candidate feline viral CTL epitopes ; (b) Define the optimal anchoring mode Based on the relative binding preference index (RBPI) calculation results of the FLA-H*0401 molecule at anchor sites P2, P3, and P9, the top 5 amino acids with the highest relative binding preference index were extracted as the "optimal anchoring mode".

[0036] (c) Pattern matching filtering and results A Python script was written to perform rapid pattern matching between all candidate antigens from step (1) and the defined set of "optimal anchoring patterns". The filtering logic was as follows: only nonapeptides that simultaneously satisfy the following conditions are retained: their P2 amino acid belongs to the set {L, M, V, T, Q}, their P3 amino acid belongs to the set {Y, F, W, L, V}, and their P9 amino acid belongs to the set {V, L, I, M, A}. After this efficient pattern matching filter, four candidate CTL epitopes were quickly selected. Table 1 shows the relevant information of the candidate CTL epitopes.

[0037] (6) Expression and inclusion body purification of FLA-H*0401 heavy chain and β2m light chain recombinant protein (a) Induction: Prepare several portions of ampicillin LB culture medium, and add 200 μL of positive bacterial culture (using PET-21a as the prokaryotic expression vector, in...) Nde I as well as Xho I Two positive bacterial cultures (with FLA-H*0401α and β2m chain fragments inserted between the restriction sites) were added to 100 mL of A+LB solution and shaken for 6 h. Then, they were added to 1 L of LB solution and cultured in a shaker for 2 h until the OD600 value was around 0.6. Subsequently, 1 mL of IPTG (1:1000) was added and the culture was induced in a shaker at 37 °C for 5 h. (b) Collection of bacterial cells: The induced bacterial culture was centrifuged at low temperature to obtain bacterial cell precipitate. The precipitate was then resuspended in 25 mL of ultrapure water and sonicated on ice for 12 s, 18 s interval, and 45 min. The lysed bacterial culture was centrifuged at low temperature at 6000 rpm for 15 min. The supernatant was discarded, and bacterial cell impurities on the surface of the inclusion bodies were gently removed with a pipette tip. (c) Washing: Vortex rinse the inclusion bodies with the prepared washing solution, centrifuge at low temperature and discard the supernatant, repeat washing until the inclusion bodies appear relatively pure; (d) Resuspension: Vortex inclusion bodies in an appropriate amount of resuspension solution until they are in a suspended state. Take a portion of the sample for identification by SDS-PAGE electrophoresis, then centrifuge at low temperature and discard the supernatant as much as possible. (e) Dissolution: Weigh the net weight of the inclusion bodies and calculate the amount required to dissolve them. Adjust the final concentration to 30 mg / mL, add the dissolving solution, and stir at low temperature until no obvious solids remain. Centrifuge multiple times to remove bottom impurities, then aliquot, label, and freeze for storage. Results are shown in [see table below]. Figure 5 A and 5B.

[0038] (7) Refolding and concentration of FLA-H*0401 with antigenic peptides: (a) Preparation of refolding system: Prepare an appropriate amount of refolding solution according to the number of candidate CTL epitopes. Take a beaker of appropriate size, add 50 mL of refolding solution, add a stir bar, cover the mouth of the bottle with a plastic film, and fix the syringe vertically on the film. (b) Adding β2m light chains: Add 300 μL of β2m light chain inclusion body solution (refer to the preparation process of patent CN109669043A) to the syringe, and slowly drop it into the beaker. Place the beaker on a magnetic rack and stir at 4°C for 8 h. (c) Adding peptide: Dissolve the synthesized peptide powder in DMSO, vortex mix, add about 2 mg of peptide to the refolding solution, and stir for 5 min; (d) Adding α chain: Add 1 mL of FLA-H*0401 heavy chain inclusion body solution (refer to the preparation process of patent CN109669043A) to the syringe, place it on a magnetic rack and stir at 4°C for more than 24 h to refold; (e) Concentration: Clean and install the 10 kDa concentration cup in advance, transfer the refolded liquid into the concentration cup, and concentrate it at 4°C; (f) Liquid replacement: After the refolded liquid volume is concentrated to about 15 mL, add 60 mL of molecular sieve solution and continue to concentrate to about 10 mL-15 mL. Transfer to a 15 mL centrifuge tube, centrifuge at 4℃ for 7000 rpm for 5 min and take the supernatant. (g) Concentration again: Prepare a 15mL, 10kDa ultrafiltration tube in advance, centrifuge at low temperature to concentrate the liquid to less than 1 mL, transfer to an EP tube, centrifuge again to collect the supernatant, filter with a 0.22μm filter membrane, centrifuge at 12000 rpm to remove air bubbles, and store in a 4℃ refrigerator for later use.

[0039] (8) Gel chromatography and anion exchange chromatography purification of the FLA-H*0401 antigen peptide complex The Superdex 200 Increase 10 / 300 GL gel chromatography column was pre-equilibrated with one column volume of molecular sieve buffer (20 mM Tris-HCl pH 8.0, 50 mM NaCl) until the salt concentration stabilized at approximately 5.7%. The following parameters were then set: maximum column pressure of 3 MPa, flow rate of 0.8 mL / min, and collection parameters of UV = 20 mAU and peak volume of 1 mL. After UV zeroing, the sample was slowly injected into the loop to begin purification. The sample loading process was monitored periodically, and samples were collected based on peak position and size. The sample at the peak tip was then collected for SDS-PAGE analysis of the binding of the FLA-H*0401 molecule to the candidate CTL epitope.

[0040] After refolding the FLA-H*0401 heavy chain, β2m light chain inclusion body protein, and candidate CTL epitope peptide in vitro, the sample was purified by gel chromatography. Under normal circumstances, the sample showed three main protein peaks in sequence: heavy chain polymer, heavy chain-β2m light chain-peptide complex, and β2m light chain. The heavy chain-β2m light chain-peptide complex appeared at around 15 mL. The position of this peak can be used to determine the binding status of the FLA-H*0401 molecule to the candidate CTL epitope.

[0041] The heavy chain-β2m light chain-peptide complex protein was collected by gel chromatography. Further verification of the relative affinity of the candidate CTL epitope for the FLA-H*0401 molecule was performed using anion exchange chromatography. First, the RESOURCE Q anion exchange column was equilibrated with one column volume of molecular sieve buffer until the salt concentration stabilized (approximately 5.7%). Then, the following parameters were set: maximum column pressure of 5 MPa, flow rate of 1 mL / min, and collection parameters of UV = 20 mAU and peak volume of 1.2 mL. After UV zeroing, the sample was slowly injected into the loop for purification. Once the sample binding to the column was confirmed, linear elution was performed under the following conditions: the proportion of buffer B (10 mM Tris-HCl pH 8.0, 1 M NaCl) was gradually increased from the initial value to 100% over 30 minutes. During elution, the sample loading was monitored periodically, and the sample was collected based on the peak position and peak size. Finally, samples from the peak tip were taken and analyzed by SDS-PAGE to evaluate the relative affinity between the FLA-H*0401 molecule and the candidate CTL epitope.

[0042] The collected heavy chain-β2m light chain-peptide complex protein was purified by anion exchange chromatography. Generally, if a target peak with a good peak shape is formed during elution, it can be determined that the FLA-H*0401 molecule has a relatively strong affinity for the candidate CTL epitope.

[0043] Based on the above, the four predicted candidate CTL epitope peptides were validated, and the results are shown in [the table below]. Figures 6 to 9 The results showed that the four predicted candidate CTL epitope peptides bound to FLA-H*0401 to form a stable MHC I-peptide complex, exhibiting strong binding ability.

[0044] This invention discloses a method for efficiently predicting MHC I molecule peptide binding motifs and identifying candidate CTL epitopes based on MM-PBSA binding free energy calculations. The method systematically constructs a pMHC I complex model library using a single-point scanning strategy, and combines molecular dynamics simulations and MM-PBSA calculations to construct a quantified "energy weight matrix" as the binding motif of the MHC I molecule, which is then used to rapidly screen high-affinity candidate CTL epitopes. This invention effectively avoids the combinatorial explosion problem, providing a precise and efficient tool for candidate epitope discovery in vaccine design and tumor immunotherapy.

[0045] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A method for screening candidate CTL epitopes based on combined free energy calculations, characterized in that, Includes the following steps: (a) The three-dimensional structure of the target MHC I molecule was constructed by homology modeling as the initial complex model of pMHC I; (b) Mutate at least one preset anchoring site residue of the peptide in the initial complex model to a variety of natural amino acids, and set the non-target site of the peptide to a non-side chain amino acid - alanine, thereby generating a complex model set of single-point mutated pMHC I. (c) Perform molecular dynamics simulations on the complexes in the complex model set, calculate the binding free energy of each mutant peptide to MHC I molecule using the MM-PBSA method, and simultaneously decompose the binding free energy to quantify the energy contribution value of each amino acid type at the preset anchoring site. (d) Based on the energy contribution value, construct the quantitative binding motif of the target MHC I molecule, wherein the quantitative binding motif includes the amino acid type preferred by different anchor sites and their corresponding binding free energy weights; (e) Using the quantified binding motif, candidate protein sequences are scanned, and candidate peptides that meet the motif characteristics are selected as CTL epitopes.

2. The method for screening candidate CTL epitopes based on binding free energy calculation according to claim 1, characterized in that, In step (a), when obtaining the three-dimensional structure of the target MHC I molecule, if there is an experimentally resolved crystal structure, the crystal structure is directly used as the initial model. If no experimentally resolved crystal structure exists, the three-dimensional structure of the target MHC I molecule is constructed by artificially mutagenesis using the crystal structure of a highly homologous MHC I molecule as a template and homology modeling software.

3. The method for screening candidate CTL epitopes based on binding free energy calculation according to claim 1, characterized in that, In step (b), the preset anchoring point includes one or more of the 2nd, 3rd, 5th, and 9th positions of the nonapeptide; the non-targeting site is set to alanine.

4. The method for screening candidate CTL epitopes based on binding free energy calculation according to claim 1, characterized in that, In step (b), mutation operations are performed independently for each preset anchor site. Each mutation changes only the amino acid type of one anchor site, while other anchor sites and non-target sites remain alanine, keeping the polypeptide backbone conformation unchanged. Only the side chain is replaced with an alanine side chain, without involving combined mutations of two or more sites.

5. The method for screening candidate CTL epitopes based on binding free energy calculation according to claim 1, characterized in that, In step (c), the duration of the molecular dynamics simulation is not less than 50 nanoseconds; the energy decomposition is used to obtain the energy contributions of van der Waals interactions, electrostatic interactions and solvation effects for each anchored residue.

6. The method for screening candidate CTL epitopes based on binding free energy calculation according to claim 1, characterized in that, The quantitative binding motif described in step (d) is represented by an energy weight matrix, where the rows of the matrix represent the anchoring sites of the peptide, the columns represent 20 natural amino acids, and the matrix element values ​​are the average energy contribution value of the corresponding amino acid at the corresponding site or the score derived therefrom.

7. The method for screening candidate CTL epitopes based on binding free energy calculation according to claim 1, characterized in that, The quantification of binding motifs in step (d) is represented by an energy weight matrix. The matrix represents the anchor sites of the peptide, with columns for 20 natural amino acids. The matrix element values ​​are the relative binding preference index (RBPI) of the corresponding amino acid at the corresponding site. The formula for the RBPI is: E(i,j) = μ[ΔG contrib (i,j)] - μ[ΔG contrib (i,j)]; Where μ[ ] represents the arithmetic mean of the energy contribution values ​​from at least 3 repeated simulations, i is the site, and j is the amino acid.

8. The method for screening candidate CTL epitopes based on binding free energy calculation according to claim 1, characterized in that, The amino acid sequences of the candidate CTL epitope peptides obtained in step (e) are shown in SEQ ID NO.1-4.

Images