記住我
ORIGINAL ARTICLE
Year : 2023 | Volume
: 18
| Issue : 2 | Page : 159-176
A comparative study of the arazyme-based fusion proteins with various ligands for more effective targeting cancer therapy: an in-silico analysis
Rezvan Mehrab1, Hamid Sedighian2, Fattah Sotoodehnejadnematalahi1, Raheleh Halabian2, Abbas Ali Imani Fooladi2
1 Department of Biology, Science and Research Branch, Islamic Azad University, Tehran, I.R. Iran
2 Applied Microbiology Research Center, Systems Biology and Poisonings Institute, Baqiyatallah University of Medical Sciences, Tehran, I.R. Iran
Correspondence Address:
Abbas Ali Imani Fooladi
Applied Microbiology Research Center, Systems Biology and Poisonings Institute, Baqiyatallah University of Medical Sciences, Tehran, Tel & Fax: +98-2188068924
I.R. Iran
Source of Support: None, Conflict of Interest: None
CheckDOI: 10.4103/1735-5362.367795
Background and purpose: Recently, the use of immunotoxins for targeted cancer therapy has been proposed, to find new anticancer drugs with high efficacy on tumor cells with minimal side effects on normal cells. we designed and compared several arazyme (AraA)-based fusion proteins with different ligands to choose the best-targeted therapy for interleukin 13 receptor alpha 2 (IL13Rα2)-overexpressed cancer cells. For this purpose, IL13Rα2 was selected as a receptor and IL13 and IL13.E13K were evaluated as native and mutant ligands, respectively. In addition, Pep-1 and A2b11 were chosen as the peptide ligands for targeted cancer therapy.
Experimental approach: Several bioinformatics servers were used for designing constructs and optimization. The structures of the chimeric proteins were predicted and verified by I-TASSER, Q-Mean, ProSA, Ramachandran plot, and Verify3D program. Physicochemical properties, toxicity, and antigenicity were predicted by ProtParam, ToxinPred, and VaxiJen. HawkDock, LigPlot+, and GROMACS software were used for docking and molecular dynamics simulation of the ligand-receptor interaction.
Findings/Results: The in silico results showed AraA-A2b11 has higher values of confidence score and Q-mean score was obtained for high-resolution crystal structures. All chimeric proteins were stable, non-toxic, and non-antigenic. AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 retained its natural structure and based on ligand-receptor docking and molecular dynamic analysis, the binding ability of AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 to IL13Rα2 was sufficiently strong.
Conclusion and implications: Based on the bioinformatics result AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 was a stable fusion protein with two separate domains and high affinity with the IL13Rα2 receptor. Therefore, AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 fusion protein could be a new potent candidate for target cancer therapy.
Keywords: Arazyme, Cancer; In silico; Interleukin 13; Targeting therapy; Peptide ligand.
An applied approach in the clinical microbiology field is the study of the effect of microbial metabolites, especially microbial toxins and enzymes for the treatment of various types of cancer [1],[2]. Moxetumomab pasudotox is the generic name of an anticancer drug, which is a fusion protein consisting of the catalytic fragment of Pseudomonas exotoxin-A (PE38) and the fragment variable (Fv) portion of the anti-CD22 antibody, approved by the United States Food and Drug Administration (FDA) for the hairy cell leukemia treatment in 2018 [3].
Another drug for which the FDA has made recommendations specifically for additional clinical/statistical data and analysis, in addition to chemistry, manufacturing, and control issues related to a recent pre-approval inspection and product quality, is oportuzumab montox (Vicineum™), contains the PE and anti-EpCAM single chain antibodies to treat bladder cancer [4],[5]. MT-3724 is an engineered form of Shiga-like toxin A subunit in phase II clinical trials for B-cell lymphoma [6]. E7777 is a recombinant cytotoxic protein composed of the human interleukin (IL) 2 and diphtheria toxin, and E7777 is a modified-processes Ontak (Denileukin diftitox). This modified-processes toxin to a reduced level of aggregated and/or misfolded protein impurities and the percentage of active proteins increased. In phase II clinical trials in Japanese patients, sufficient efficacy with a manageable safety profile of E7777 was confirmed for the treatment of cutaneous and peripheral T-cell lymphoma [7].
Proteases are microbial metabolites with many biological functions. Arazyme (AraA), an extracellular metalloproteinase secreted by Serratiaproteamaculans, has strong proteolytic activity on various protein substrates such as albumin, elastin, creatinine, and collagen in a wide pH and temperature range [8],[9]. Several studies have reported the strong anticancer activity of AraA, especially its antimetastatic effect on various cancers such as colon, ovarian, and melanoma cancer cell lines [8],[10],[11]. Ghadaksaz et al. have investigated a fusion protein consisting of transforming growth factor alpha third loop (TGFαL3)-targeted AraA for breast cancer therapy [12]. The main mechanism of AraA for cancer treatment is its protease-dependent action. In addition, this enzyme has non-proteolytic antitumor effects which are mostly produced by the induction of cytotoxic T lymphocytes [10]. The hepatoprotective anti-inflammatory activity [13] and the inhibitory effect of this enzyme on allergic inflammation have also been demonstrated [14]. Moreover, AraA induces antioxidant signaling and inhibits the release of T helper cell type 2 cytokine [15].
For many years, ligand-targeted therapeutics have been studied to improve the efficacy of cancer treatments, but in recent years these strategies have changed significantly [16]. In general, the density of the target receptor on the surface of tumor cells must be higher than normal cells or overexpressed [17]. IL13 is secreted by T-helper 2 cells, mast cells, and CD4 cells such as macrophages [18]. This cytokine plays an important role in immune responses to cancer, inflammation, and infection and binds to both IL13 receptor alpha 1 (IL13Rα1) and IL13Rα2 receptors [19]. The IL13Rα2 has a high affinity for IL13 and is overexpressed in various cancer, compared to low or no expression in natural tissues. High-affinity binding of IL13Rα2 for IL13 induces signals via an activator protein 1-dependent pathway independent of STAT-6, which results in increased activity of TGF-β, but it does not induce signal transmission through the JAK/STAT6 pathway like IL13Rα1 [20]. IL13.E13K is a modified IL13 with a specific high affinity for IL13Rα2 containing a modified glutamic acid in the 13th amino acid instead of lysine [21].
Recently, linear peptide ligands have been considered due to their small size, low molecular weight, high flexibility, high specificity, low antigenicity, and non-immunogenicity [22],[23]. One of the peptide ligands which binds specifically to the IL-13Rα2 is Pep-1, a linear peptide with CGEMGWVRC synthetic sequence [22]. The A2b11 peptide (WALRVKAG), which was a T7 random peptide phage display, was shown to bind to IL13Ra2 with the highest specificity [24].
Based on the findings of the in silico study, to obtain an effective fusion protein, comparison of the designed proteins is easier and rapid analysis before starting the laboratory experiments [25],[26]. Moreover, the pre-experimental study utilizing bioinformatics tools is cost-effective and fast in contrast to the experimental methods. In a recent study, several novels engineered analogs of serratiopeptidase were computationally designed and molecular dynamically simulated. Based on the results of the in silico study, two engineered analogous of serratiopeptidase (T344 and T380) were chosen as appropriate candidates for the experimental study [26].
Bioinformatics tools (such as servers, software, etc.) are designed and developed to extract information such as sequence or retrieve data from genomic sequence databases as well as predict biomolecule structures including DNA, RNA, proteins, etc. Also, this field of science is widely used in different studies of biological research such as the design and discovery of new drugs and vaccines plus structure and function prediction [27]. In addition, these tools provide access to scientific databases in various fields of the life sciences, including genomics, proteomics, systems biology, phylogeny, population genetics, etc.
Since it is necessary to find anticancer agents with high efficacy on tumor cells and minimal side effects on normal cells, here, we have tried to design and compare several AraA-based fusion proteins with different ligands to choose more effective targeted therapy for some cancers with overexpression of IL13Rα2 using bioinformatics tools. Consequently, IL13Rα2 was selected as the receptor and IL13 as well as IL13.E13K were considered native and mutant ligands, respectively, for targeted cancer therapy. In addition to IL13 and IL13.E13K, we also studied Pep-1 and A2b11 as peptide ligands for targeted drug delivery to IL13Ra2-overexpressed cancer cells. Our next study will confirm these in silico results using in vitro and in vivo studies.
Materials and methodsDesigning the chimeric constructs
The protein sequence of AraA (Serratia proteamaculans) and IL13 protein were obtained from Uniprot (with accession numbers Q2VU40 and P35225, respectively) and GenBank database (with accession number AAX21094.1 and X69079.1, respectively). In this study, the main component of chimeric constructs was the AraA protein as a functional site that fused with various ligands including IL13, IL13.E13K, Pep-1 and A2b11 for targeted drug delivery to IL13Rα2-overexpressed cancer cells. Several linkers were used for fusion protein design and tested by the I-TASSER server to obtain the best linker for separating two functional domains of the fusion proteins that do not interfere with each other or are only minimally disturbed, compared to their native 3D structure. Finally, the sequences of the chimeric proteins were codon optimized according to E. coli codon usage using the Java Codon Adaptation tool (JCat, http://www.jcat.de/) to increase the protein expression. In addition, GC content and codon adaptation index (CAI) were calculated by this server.
The physicochemical property and secondary-structure prediction
In order to compute the physiochemical characterization of the designed constructs, including molecular weight, theoretical isoelectric point, the total number of positive and negative residues, half-life, instability index, extinction coefficient, aliphatic index, and grand average of hydropathicity, the Expasy ProtParam server (https://web.expasy.org/protparam/) was used. To predict the percentage of secondary structures in the fusion proteins, the GOR IV online database was employed [28].
3D-Structure prediction of fusion proteins
The 3D model of the chimeric proteins was predicted using the I-TASSER (https://zhanglab.ccmb.med.umich.edu/I-TASSER/), which generates 3D models along with a confidence score (C-Score) for the quality of the predicted structure. I-TASSER predicts the 3D structure of the fusion protein and function based on a hierarchical approach from sequence to structural and functional model [28]. To investigate the conformational state of the domains in fusion proteins, 3D structures aligned with protein data bank (PDB) file 1SAT (based on the template, Q2VU40) for AraA and 1IK0 solution structure of human IL13 using the template modeling (TM)-Align [29] and CHIMERA software V1.10.2. Finally, the 3D structure of the superior fusion protein was also predicted by Robetta (http://new.robetta.org/) [30] and RaptorX (http://raptorx.uchicago.edu/) [31] online servers and the predicted models were compared by different servers to check the similarities. RaptorX is a server for predicting the 3D structure of the fusion protein without using templates, including solvent accessibility, secondary structure and disordered regions [30]. The Robetta server generates 3D structural models for fusion protein via either comparative modeling or de novo structure prediction techniques [30].
Refinement and validation of fusion protein structures
In order to improve the 3D protein structure predicted to be closer to the native structure, all fusion proteins were refined using the 3D refine server (http://sysbio.rnet.missouri.edu/3Drefine/index.html) [32].
The refined 3D structure of all fusion proteins validated by both QMEAN (https://swissmodel.expasy.org/qmean/) [33] and ProSA-web (https://prosa.services.came.sbg.ac.at/prosa.php) [28]. The quality measured as Z-score was computed by ProSA and indicated that the predicted structures were in the score range of the native proteins with similar size.
For more validation, the predicted structures were checked using UCLA-DOE LAB - SAVES v6.0 server for several parameters including ERRAT, VERIFY, PROVE and PROCHECK (https://saves.mbi.ucla.edu/). By ERRAT, nonbonded interactions between different types of atoms were analyzed and showed the value of the error function in comparison with highly refined structures. The overall quality factor for good structure with high resolution and lower resolutions is valued at around ≥ 95% and 91%, respectively. The compatibility of the 3D model with its amino acid sequence was determined by assigning a structural class based on its environment and location (polarity, non-polarity, alpha, beta, loop, etc.) and comparing the results to good structures by the VERIFY3D program. In verified structures, at least 80% of the residues have a 3D/1D score ≥ 0.2.
The PROVE program calculated volumes of atoms in the macromolecules and the Z-score deviation for the predicted 3D structures from highly resolved and refined PDB structures were calculated. If the total buried outlier protein atoms of the predicted structure, was > 5%, the 3D model failed. The PROCHECK program was applied to evaluate the stereochemical quality of the 3D structure. In this program, residue-by-residue geometry and overall structural geometry were analyzed and then compared to the stereochemical parameters of the high-resolution and well-refined PDB structures.
Prediction of cleavage sites
Proteasome cleavage sites in the fusion protein have been predicted by Netchop 3.1. This server is the best neural network for cleavage site predictions of the human proteasome and approximately 75% of the cleavage sites have been correctly predicted with false positives of about 15% [17]. The C-term 3.0 network was used for predicting the cytotoxic T lymphocyte epitopes. This network is a database containing 1260 class I major histocompatibility complexes ligands.
B-Cell epitopes prediction
Linear B cell epitopes of the chimeric proteins were predicted by the ABCpred server (https://webs.iiitd.edu.in/raghava/abcpred/ABC_submission.html) with a threshold of 0.8 and 16 mers epitope length [34]. Continuous B cell epitopes of fusion proteins were predicted using the BcePred server (https://webs.iiitd.edu.in/raghava/bcepred/bcepred_submission.html) with the 1.9 threshold based on the antigenic propensity [35].
Toxicity and antigenicity prediction
The toxicity of the chimeric proteins has been predicted by ToxinPred. ToxinPred is an in silico database server for checking the presence of toxic peptides in proteins. The foremost dataset used in this method consists of 1805 toxic peptides (≤ 35 residues) [36]. To identify the antigenicity of the chimeric proteins the VaxiJen V2.0 server was employed. VaxiJen is an in silico approach to the prediction of subunit vaccines and protective antigens [28].
Ligand-receptor docking
In this study, AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 fusion protein was selected and further investigated for its natural ligand. Molecular docking simulation is performed using HawkDock web server. HawkDock is a web server to predict and analyze the protein-protein complex based on computational docking and molecular mechanics-generalized born surface area (MM/GBSA). Our proteins were docked with the HawkDock web server [37] and ranked according to binding free energy scores. The hydrogen bond and hydrophobic interactions between the IL13 and IL13Rα2 of the docked complexes were drawn using the LigPlot+ software.
Among the molecular species docked by the HawkDock server, the best molecular species were extracted and subsequently, were selected as the primary structure to enter the molecular dynamics (MD) simulation step.
MD simulation of complexes
The MD simulation of the best docking pose of top-ranked complexes was accomplished using the GROMACS simulation package (2020) [38]. The all-atom optimized potentials for liquid simulations force field were used for simulating our system [39]. The protein-protein complexes were solvated in the MD simulation box with transferable intermolecular potential 3P (TIP3P) solvent molecules and then neutralized by adding 0.15 mol/L Na+/Cl- ions. Energy minimization of all systems was performed by the steepest descent algorithm. The systems were then equilibrated using NVT and NPT with 400 ps steps, respectively, utilizing a V-rescale Berendsen thermostat and a Parrinello-Rahman [38]. The heating of the systems was gradually increased from 0 to 310 K, and the pressure of the systems was set to 1 atm for the NVT and NPT ensembles, respectively. The particle-mesh Ewald (PME) and the LINCS algorithms were applied to assess all electrostatic connections and to restrain all bond lengths in the protein, respectively. Moreover, periodic boundary condition was utilized during the simulation [37].
The MM/GBSA method was also, applied to calculate the binding free energy using the following equations:
In these equations, Greceptor, Gligand, and Gcomplex refer to free energies of the receptor, ligand, and complex, respectively. A sum of changes in molecular mechanical (MM) gas-phase binding energy (ΔEmm), the solvation-free energy (ΔGsolv), and entropic (-TΔS) contributions were applied to compute the binding free energy (ΔGbind). Van der Waals (ΔEvdW) and gas-phase electrostatic (ΔEele) energies create the molecular mechanical term (ΔEmm). Nonpolar (ΔGnonpol) and polar (ΔGpol) energies create ΔGsolv. The ΔGnonpol was calculated using ΔGnonpol = γSASA + β [40].
ResultsDesigned and optimized codon chimeric constructs
The major component of the chimeric constructs was the Arazyme protein (AraA, Serratia proteamaculans) located at the functional domain of the fusion proteins and considered to be an anticancer agent. Several ligands were investigated for targeting cancer therapy to obtain the best ligand with high affinity to IL13ralphα2, including IL13 as a native ligand, IL13.E13K as a modified ligand, Pep-1 and A2b11 as peptide ligands. The Schematic model of designed chimeric proteins is shown in [Figure 1].
Figure 1: Schematic model of designed chimeric constructs. (A) AraA + linker + native IL13; (B) AraA + linker + IL13.E13; (C) AraA + linker + Pep1; (D) AraA + linker + A2b11.Several hydrophobic linkers (34 linkers listed in [Table 1]) used in the last studies, were tested by the I-TASSER server to obtain two separate functional domains of the fusion proteins do not interfere with each other or are only minimally disturbed, compared to their native 3D structure of the 33 linkers, only linker-31 (A(EAAAK)4ALEA(EAAAK)4A)2 was selected as the appropriate hydrophobic linker for AraA-IL13 and AraA-IL13.E13K chimeric constructs with two distinct functional domains. Ultimately, the restriction enzyme BamH1 site at the 5' end and HindIII site at the 3' end of the construct were regarded to facilitate the cloning in the pET28a+ vector for further cloning purposes [41]. Improved GC content and CAI to near 50% and 1 (GC% of E. coli is about 50%), respectively. Therefore, the stability of mRNA in the chimeric gene is expected to increase.
Table 1: List of linkers investigated to obtain fusion proteins with separate domains.The physicochemical properties and secondary structure prediction
The results of the physicochemical properties of the chimeric proteins from the Expasy ProtParam server are presented in [Table 2]. The instability index of all designed fusion proteins is less than 40, indicating stability. AraA-(G4S)1-Pep1 and AraA-(G4S)1-A2b11 fusion proteins have better interaction with water due to their lower grand average of hydropathicity. AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL 13 and AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13.E13K with high aliphatic index (73.12) have thermal stability for a wide range of temperatures.
The results of the secondary structure of fusion proteins using the GOR IV server are shown in [Table 3]. The most common secondary structural elements are the random coils and then alpha helix. The beta sheets are considered to be the third common secondary structure.
Table 3: The secondary structure of fusion proteins predicted by GOR IV.The tertiary structure of designed fusion proteins
The 3D models of the designed fusion proteins were predicted using the I-TASSER (https://zhanglab.ccmb.med.umich.edu/I-TASSER/) online server, which generates 3D models along with a C-score showing the quality of the predicted structure. I-TASSER-predicted Model.1 of five probable tertiary structures is shown in [Figure 2]. Afterward, the 3D model of the superior fusion protein was predicted using Robetta and RaptorX to confirm the predicted models of I-TASSER [Figure 3].
Figure 3: 3D structures of the selected fusion protein were predicted by (A) Robetta and (B) RaptorX; (C) Robetta model aligned with I-TASSER model; (D) model RaptorX aligned with I-TASSER model; and (E) the results obtained from TM-align server for aligned AraA-IL13 fusion protein with protein data bank files of AraA and IL 13 using Chimera software 1.10.2. AraA, Arazyme; IL, interleukin.3D Models of chimeric proteins predicted by I-TASSER showed the chimeric proteins with two separate main domains linked together with the linker. The C-score of model 1 was the maximum when compared to other models. The C-score is generally reported in the range of-5 to 2 and a higher C-score indicates higher confidence in the model. In addition, the TM-score and root-mean-square deviation (RMSD) for model 1 were reported in [Figure 2].
In the next step, model 1 of the predicted 3D structure of the selected fusion protein by I-TASSER was aligned with the PDB file of AraA and IL13 using the TM-Align server. TM-Score for AraA and IL13 aligned with AraA-IL13 were 0.887 and 0.866, respectively (0.0 < TM-score < 0.30, random structural similarity and 0.5 < TM-score < 1.00, in about the same fold) [28] which indicates that the protein fusion components are similar to their native structure [Figure 3]E.
Refinement and validation of fusion protein structures
To improve the predicted 3D protein structure to be similar to the native structure, all fusion proteins were refined using the 3D refine server. Then, the refined 3D structure of all fusion proteins was confirmed with the output of both QMEAN (https://swissmodel.expasy.org/qmean/) and ProSA (https://prosa.services.came.sbg.ac.at/prosa.php). The higher value of the C-score signifies a higher confidence score and vice versa. AraA-(G4S)1-Pep1 and AraA-A2b11 (without linker) have higher C-score and Q-mean scores. The 3D structures of all designed fusion proteins were analyzed by ProSA-web to detect errors in the 3D model of fusion proteins. Fusion proteins with (A(EAAAK)4ALEA(EAAAK)4A)2 linker are out of the native protein range [Table 4].
Table 4: Structure refinement and scores of Q-mean and ProSA for the designed fusion proteins.For more validation, the predicted structures were checked using UCLA-DOE LAB - SAVES v6.0 server with several parameters including ERRAT, VERIFY, PROVE, and PROCHECK (https://saves.mbi.ucla.edu/). Unfortunately, all overall quality factors for all designed fusion proteins were less than 90%, while overall quality factors for good structure with high resolution and lower resolutions were calculated to be about ≥ 95% and 91%, respectively.
In structures verified by the Verify3D program, at least 80% of the residues have a 3D/1D score ≥ 0.2. Designed fusion proteins with (A(EAAAK)4ALEA(EAAAK)4A)2 linker have less than 80% of the residues with the averaged 3D-1D score ≥ 0.2, except AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 that has 82.77% of the residues with scored ≥ 0.2 in the 3D/1D profile. In the PROVE program, buried outlier protein atoms total for the 3D model of AraA-(EAAAK)1-A2b11, AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 fusion proteins were calculated to be > 5%, therefore these structures have failed.
In the Ramachandran plot presented by the PROCHECK program, a good-quality model would be expected to have over 90% amino acid residues in the most favored regions. Among the designed structures, AraA-A2B1 (without linker) fusion protein with the highest score has 82.9% residues in the most favored regions with the lowest errors/warnings [Table 5].
Table 5: The assessment of a 3D refined model of designed fusion proteinsCleavage sites predicted in designed fusion proteins
Proteasome cleavage sites inside the fusion protein were predicted by Netchop 3.1 server. The Net Chop neural network-based method was the most effective presently-available system for the predictions of the cleavage site. The new version of Net Chop has predicted approximately 75% of cleavage sites correctly with false positives near 15%. The number of the cleavage sites on the fusion proteins computed with NetChop 3.1 server using version C-term and threshold 0.800000 has been shown in [Table 6].
Table 6: The predicted cleavage sites of the human proteasome for designed fusion proteins by Netchop 3.1.Cleavage sites of the human proteasome by AraA in full sequence are shown in [Table 6] and AraA sequence for other designed fusion proteins is not shown.
B-Cell epitopes
The linear B cell epitopes predicted for AraA and AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 have been shown in [Table 7]. The predicted B-cell epitopes for IL13.E13K (modified IL13) were similar to the native IL13. The continuous B-cell epitopes for other residues in the designed protein were not identified by Bcepreds, except, the TDFIVGGGGSWALRVK sequence in start position 500 of AraA-(G4S)1-A2b11 fusion protein with the score 0.87 predicted as B-cell epitopes. The continuous B cell epitopes predicted for AraA, the main component of fusion proteins, and IL13 are shown in [Table 8].
Table 7: The linear B-cell epitopes predicted in designed fusion proteins using the ABCpred prediction server.
Table 8: The continuous B-cell epitopes predicted in chimeric proteins using Bcepred server.The predicted continuous B-cell epitope for IL13.E13K (modified IL13) was similar to the native IL13. The continuous B-cell epitopes for other residues in the designed protein were not identified by Bcepreds, except for AraA-Pep1 (without linker) fusion protein that predicted TDFIVCGE sequence as B-cell epitopes.
Toxicity and antigenicity of chimeric proteins
For toxicity prediction, the designed chimeric proteins were scanned with ToxinPred server. This server predicts whether overlapping peptide/analog is toxic or not. Based on the protein scanning using the ToxinPred server, all designed chimeric proteins are considered non-toxic for peptide fragment length 10. To identify the antigenicity, VaxiJen V2.0 was employed. All designed fusion proteins predicted probably are not antigenic for the tumor model with threshold 1 using the VaxiJen server.
Ligand-receptor docking
In this study, AraA-(A(EAAAK)4ALEA (EAAAK)4A)2-IL13 fusion protein was selected and further investigated due to its natural structure of components. Protein-protein docking was done by the HawkDock web server [37] and ranked according to binding free energy scores [Table 9]. The schematic 2D view of the IL13-IL13Rα2 and chimeric IL13-IL13Rα2 docked complexes were generated by LigPlot+ software. The hydrogen bond and hydrophobic interactions between complexes are presented in [Figure 4]. To validate the docking results, 100 ns MD simulations on the top-ranked binding pose are performed.
Figure 4: (A) Interaction between IL13 and IL13Rα2; (B) interaction between chimeric protein and IL13Rα2 receptor; (C) the 2D view of the docked complexes were analyzed using LigPlot+ software. The hydrogen bond (green dashed lines) and hydrophobic interactions (spoked arc) between the IL13-IL13Rα2; (D) chimeric IL13-IL13Rα2. IL, Interleukin.
Table 9: Protein-protein docking results of chimeric protein with IL13Rα2.Molecular dynamics simulation of complexes
The RMSD values for the IL13-Rα2 and chimeric-IL13-Rα2 were calculated as a time-dependent parameter for clarification of the displacement of Cα atoms in the two molecular types during the MD simulation [Figure 5]A. According to the RMSD curve, The IL13-Rα2 structure has a small fluctuation compared with chimeric-IL13-Rα2. The high value of the RMSD's chimeric-IL13-Rα2 can be arising from the big size of the chimeric-IL13. Also, the RMSD values of the IL13 in the two complexes were estimated separately [Figure 5]B. According to [Figure 5]B, IL13 has a similar pattern of RMSD values during simulations gradually reaching a stable state with a smaller fluctuation range after 5000 ps. These findings mean that the IL13 has good stability in the two complexes [Figure 5].
Figure 5: (A) The Cα-RMSD graph for the IL13Rα2 and chimeric-IL13Rα2 in the complex state; (B) the Cα-RMSD graph for the IL13 and chimeric-IL13 at the complex state; (C) the Cα-RMSF graph for the IL13 in the complexes state. IL13Rα2, Interleukin 13 receptor alpha 2.In general, based on the results obtained from the analysis of structure fluctuations, the protein-protein complexes have good structural stability during simulation times. Also, to estimate the strength of bonds for the protein-protein complexes and to clarify the residues which play the important role in protein-protein interactions, free energy calculations and binding energy decompositions per residue were done based on the 100 frames retrieved from the last 10 ns of the MD simulation trajectories. MM/GBSA approach was used to calculate components of the binding energy and to assess the energetic contribution of residue to the binding free energy by energy decomposition scheme [Table 10].
Table 10: Contribution of each energy component associated with binding free energy (kcal mol-1).According to [Table 10], the non-polar or hydrophobic energies (ΔEnonpolar) have a large contribution to the binding energy. It can be said that these hydrophobic energies are the major force in ligand-receptor binding and electrostatic interactions (ΔEpolar) in comparison with hydrophobic interactions, which appear to be unfavorable for ligand binding to the receptor.
Direct intermolecular interactions (ΔEele) in the IL13Κα2 complex are favorable, but the contribution of these interactions is neutralized by the large component (ΔGGB) which plays a major role in the electrostatic interactions. These interactions (ΔEele) play an inverse role in the chimeric protein complex with IL13Rα2. However, polar interactions generally play a negative role in the interactions between the two species complexes. To obtain details of the intermolecular interactions, the contribution of each residue in the binding energy between different species of proteins was calculated (data are not shown). Considering the contribution of each component in the free binding energy, the major and desired energies involved in the ligand binding are Van der Waals (vdW) and electrostatic (ele) interactions, while the polar solvation energy (ΔGgb) in binding the ligand to the protein is undesirable. vdW interactions and non-polar solvent energies are responsible for covering hydrophobic groups of proteins.
Therefore, it can be stated that hydrophobic residues play an important role in forming the binding envelope between the studied protein complexes and assure the stability of the complex system during the simulation. These residues also prevent the breakdown of the protein species complexes [Figure 6].
Figure 6: The 3D view of the final protein-protein complexes. (A) IL13Rα2-IL13; (B) IL13Rα2-IL13 in chimeric state. Residues referring to protein-protein interface sites (cyan and magenta sticks) are labeled. The IL13Rα2 is represented as the red cartoon and IL13 as the yellow cartoon. IL13Rα2, Interleukin 13 receptor alpha 2. Discussion Common treatments for cancer are chemotherapeutic agents, radiation therapy and etc. that kill both tumor cells and normal cells. Therefore, it may lead to some very unbearable side effects [42]. In different types of cancers, several receptors and enzymes are overexpressed. Therefore, in the last decade, one of the new strategies for target cancer therapy is the identification of overexpressed receptors [42]. Techniques and applications of bioinformatics provide computational methods for designing new drugs with specific ligands for overexpressed receptors. Previous studies have shown that IL13Κα2 is overexpressed on several cancer cells [43]. In addition, in recent studies, AraA has been introduced as a promising new antitumor chemotherapy agent [8],[15]. Ghadaksaz et al. used AraA as an effective domain in AraA-linker-TGFαL3 chimeric proteins [12]. In the current study, we have also selected AraA as the catalytic domain in the designed fusion proteins. These fusion proteins were designed to achieve the highest efficiency of the AraA for targeted cancer therapy and minimal side effect on normal cells [10]. Afterward, various ligands with different linkers were investigated to obtain the best ligand with high affinity to IL13Rα2, including IL13 as native ligand, IL13.E13K as modified ligand, Pep-1 and A2b11 as peptide ligands for targeting IL-13Rα2-overexpressed cancer cell.
At the beginning of the in silico analysis, physicochemical parameters were computed by ProtParam server. The instability index for all designed fusion proteins was less than 40, therefore predicted as stable. The high aliphatic index in the range of 67.10 to 73.12 indicates the thermostability of all designed fusion proteins with a high number of hydrophobic amino acids. In order to increase the mRNA stability and high expression of the fusion proteins, nucleic acid codons were optimized based on the E. coli host.
After predicting the refined 3D structure of all designed fusion proteins, the structure was validated by both QMEAN and ProSA. For more validation predicted structures were checked by UCLA-DOE LAB - SAVES v6.0 server using several parameters including ERRAT, VERIFY, PROVE, and PROCHECK. Among designed fusion proteins, AraA-(G4S)1- Pep1 and AraA-A2b11 (without linker) have higher values of C-score model 1 of 3D structure predicted by I-TASSER. Moreover, these fusion proteins have high Q-mean scores with high-resolution crystal structures. By more validation of predicted structures, it was determined that AraA-A2b11 (without linker) had highly favored regions with 82.9% in the Ramachandran plot and had verified structure by Verify3D program with 92.38% of the residues and averaged 3D-1D score ≥ 0.2. Based on the protein scanning, all designed chimeric proteins are considered non-toxic and non-antigenic.
The best linker was considered to create the maximum distance between the two fused parts so that the folding of the two structures is close to the original structure. Based on I-TASSER results, among 33 linkers, only (A(EAAAK)4ALEA(EAAAK)4A)2 linker was selected as an appropriate rigid linker in AraA-IL13 and AraA-IL13.E13K chimeric constructs. It is due to the huge structure of the AraA that requires maximum distance from the other parts of the molecule for the best folding and for the activity with enhanced stability. Furthermore, due to their small size and high flexibility, peptide ligands may lead to folding in the structure of AraA by fusing to this large enzyme. It may lead to reduced accessibility and prevent their binding to the receptor. Therefore, we selected the AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 fusion protein for further studies, because the native ligand (IL13) and two separate domains including (A(EAAAK)4ALEA(EAAAK)4A)2 rigid hydrophobic linker improve the biological activity of the protein, increase the stability of the structure and lead to a high expression due to the correct protein folding [33],[44]. Tertiary structures of AraA-(A(EAAAK)4 ALEA(EAAAK)4A)2-IL13 chimeric protein have two separate functional domains including IL13 as a specific ligand and AraA as an anticancer agent. The 3D structure of this chimeric protein was predicted by the I-TASSER, Robetta, and RaptorX servers with (A(EAAAK)4ALEA(EAAAK)4A)2 linker that were approximately similar and each domain was analogous to the native component, indicating the preservation of the functional structure. Ghadaksaz et al. used the A(EAAAK)4ALEA(EAAAK)4A linker for their construct [33], but we used two series of this linker so that two separate domains do not interfere with each other or are only minimally disturbed, compared to their native 3D structure.
The analysis of physicochemical parameters showed that this fusion protein had an acidic nature with an isoelectric point of 5.01. Our fusion protein had 708 amino acids and its molecular weight was estimated to be 74.76 KDa. Furthermore, this chimeric protein is probably non-antigenic and non-toxic.
The HawkDock server was used as a molecular docking method to model the interactions between proteins at the atomic level, which allows us to clarify the behavior of interface residues in the binding pocket of target proteins [37]. Protein-protein complex with a high-affinity score was selected as the best pose mode. The MD simulations were used to evaluate the stability of the ligand (IL13) in the binding state. The stability of the ligand binding state in the binding pocket indicates the accuracy and stability of the ligand binding state. To determine the stability of the ligand binding state in the binding pocket, parameters such as the Cα-RMSD and Cα-root-mean-square fluctuation (RMSF) values of the protein-protein complexes were calculated. Simulation can be performed to determine the ligand bindings to the receptor over time. RMSD analysis can reveal Cα atoms replacement in ligand binding states. By examining this parameter, it is possible to observe the changes in the complex state gradually during MD simulation times, and therefore, the structure corresponding to this change from the simulation film may be extracted.
The current study showed that after an increase in the RMSD graph, the structural fluctuations reach a plateau state which was confirmed by the evaluation of complexes up to 60 ns showing that the pattern of changes at 30 and 60 ns is similar. Thus, based on a comparison between the general and minor fluctuations of the ligand (IL13) structure in the monomeric and chimeric states incorporated into the IL13Rα2 junction envelope, it can be clearly seen that the ligand in the junction envelope exhibits a steady state and completely adapts to the binding site and does not dissociate. Also, the partial structural fluctuations of the IL13 in native and recombinant molecules showed that the pattern of each molecular species is almost similar, indicating the stability of the related molecular species during the simulation. Therefore, based on the results of the present in silico study, we predicted that AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 fusion protein tightly binds to IL13Rα2 through IL13 domain and can exhibit potent anticancer function through the AraA domain. In the next study phase, in vivo, the anticancer activity of this fusion protein will be investigated on different cell lines.
ConclusionThe main purpose of this study is to introduce a new potent fusion protein, which can be a candidate for cancer therapy. In order to find a new anti-cancer drug, we investigated the structure and function of the several fusion proteins that contain AraA for the elimination of cancer cells and a ligand that can be attached to the IL13Rα2 receptor by various bioinformatics tools. The results of this study demonstrated that AraA-IL13 fusion protein with native ligand was a stable chimeric protein with two separate domains and high affinity to IL13Rα2 receptor that is overexpressed in various human tumor cells. Therefore, the bioinformatics results of the present study demonstrated although some considerate fusions protein also show some desirable characteristics as well but AraA-(A(EAAAK)4ALEA(EAAAK)4A)2-IL13 fusion protein could be a new potent candidate for in vitro and in vivo studies in order to target cancer therapy.
Acknowledgments
The authors would like to thank Dr. Abdolamir Ghadaksaz for his technical help.
Conflict of interest statement
The authors declared no conflict of interest in this study.
Authors’ contribution
A. A. Imanifooladi and R. Halabian supervised the project; A. A Imanifooladi and H. Sedighian conceived the methodology and validation; R. Mehrab performed the experimental parts of the study; R. Mehrab, H. Sedighian, and R. Halabian analyzed the data; R. Mehrab wrote the original draft of the manuscript and provided grammatical revisions to the manuscript by F. Sotoodehnejadnematalah. The final version of the manuscript was approved by all authors.
References
留言 (0)