Construction of IL15-NGR chimeric protein

In order to construct a chimeric IL-15-NGR protein, amino acid sequence of IL15 and NGR was retrived from the from NCBI database. Both the sequence were fused through the rigid linker that results in a chimeric protein with a single chain of 133 Amino acids (Fig. 1). A rigid linker was used to ensure the cleavage of fusion protein in two functional peptides and to prevent the disulphide bridge formation between peptides.

Fig. 1

Schematic representation of IL-15 fusion with NGR

Secondary structure prediction

The amino acid sequence dictates the shape, structure, and function of a protein. GOR IV analysis predicted the secondary structure of the submitted chimeric protein, specifying the distribution of alpha helices and beta sheets. The predicted structure unveiled that the construct comprises 50.38% (67 residues) organized into alpha helices, 12.03% (16 residues) forming extended strands, and 37.59% (50 residues) adopting random coil conformations (Fig. 2). The higher percentage of alpha helix in chimeric construct ensures structural stability of IL15-NGR (Sakurai et al. 2005).

Fig. 2

Secondary structure prediction of IL-15/NGR via GOR IV

Evaluation of tertiary structure of chimeric IL-15-NGR protein

Tertiary structure of designed chimeric protein was predicted by homology modeling using I-TASSER and trRosetta online servers. The top 5 models predicted by both servers were evaluated for selecting the most accurate configuration of designed chimeric protein on the basis of highest C-score and TM-score. The finalized model among I-TASSER was having C-score − 1.50 and TM-score 0.53 ± 0.15 while the finalized model from trRosetta was having TM-score of 0.915. The selected most stable models were evaluated further for making comparison among the models predicted by two softwares I-TASSER and trRosetta.

Suitable 3D structure selection, validation, and quality assessment

To determine the stereochemistry of designed chimeric protein, the final IL15-NGR 3D structure obtained from I-TASSER and trRosetta was uploaded on SAVES online server. The most reliable structure among the two selected models from different softwares (I-TASSER and trRosetta) was selected on the basis of Verify 3D, Ramachandran plot and ERRAT value (Table 2).

Table 2 Comparison of the quality of 3D models predicted by different software

The results indicate that the most stable and reliable structure was predicted by trRosetta (as shown in Fig. 3) that shows 97.6% residues in most favored region and ERRAT value of 96.72 as shown in Fig. 4B. Through GalaxyRefine tool, trRosetta predicted IL15-NGR 3D structure was refined and validated by Ramachandran plot, ERRAT analysis. Ramachandran plot analysis was done by PROCHECK which shows that the IL15-NGR model contains 99 residues in Ramachandran most favored region (Fig. 4A). For analyzing the non-bonded interactions statistics among various atom typesIL15- NGR 3D structure was uploaded on ERRAT program. Higher ERRAT value indicates high quality of 3D structure. ERRAT value of final model was predicted to be 96.72 (Fig. 4C).

Fig. 3

PYMOL illustration of IL-15/NGR fusion protein

Fig. 4

Structure validation and quality index of IL15-NGR fusion protein

Figure 5A illustrates the model quality, with black dots representing the IL-15/NGR fusion protein’s data points. These dots, situated among dense clusters of X-ray and NMR data, signify a negative Z-score of -6.2 which is consistent with high-quality structures, affirming the fusion protein’s reliable folding. The broad distribution of data points justify structural variability, yet the fusion protein’s proximity to these clusters reinforces confidence in its computationally generated model. In Fig. 5B, the energy profile graph from ProSA-web for the IL-15/NGR fusion protein indicates predominantly negative values across most sequence positions, suggesting structural stability. Fluctuations, particularly with a smaller window size, reveal localized variations in stability, smoothed out in larger window sizes, implying overall structural integrity.

Fig. 5

Quality assessment graphs and solubility prediction of IL-15/NGR fusion protein

Prediction of the physiochemical properties of chimeric protein

Using the Protparam web service, the physical and chemical characteristics of the protein were estimated after submitting amino acid sequence. Table 3 summarized all the Parameters of chimeric protein. The chimeric protein has calculated theoretical PI value 4.89 and estimated molecular weight 14503.49. The ratio of negatively charged residues (Asp + Glu) in chimeric protein is 19 and positively charged residues (Arg + lys) in chimeric protein is 11, this imparts that this chimeric protein (IL-15/NGR) has Acidic properties. For protein-protein interaction studies, the extinction coefficient value of chimeric protein plays a very important role. The estimated half-life of chimeric protein in vitro among mammalian reticulocytes is 1.4 h and in E.coli the estimated half-life is more than 10 h. The value of instability index estimated is 50.38 and aliphatic index value is 96.09 which indicates the protein stability over wide range of temperature. The value of gravity index indicates the reactivity of chimeric protein. The Protparam server generated − 0.047 GRAVY value of chimeric protein. This fusion protein’s extremely low GRAVY index suggests that IL-15 NGR may lead to a more effective interaction with water.

Table 3 Physiochemical properties of IL15-NGR fusion peptideSolubility prediction in chimeric protein

The FASTA amino acid sequence was processed using the protein-sol programme in order to evaluate solube expression of fusion protein in E.coli. The prosol web server determines the solubility of proteins by contrasting the average solubility of the population dataset in E. coli with the anticipated solubility of the query sequence (Hebditch et al. 2017). The threshold value for protein solubility set by protein-sol server is 0.45. As value of predicted chimeric protein is 0.844 (Fig. 5C) which is greater than threshold value, this indicates that designed protein is highly soluble. The PI value of the designed protein is 4.880.

Evaluation of toxicology, allergenicity, and antigenicity

According to predictions made by Toxinpred, Algpred Vaxigen servers, Algpred, Vaxigen servers, the engineered chimeric protein IL-15-NGR is nontoxic. The protective antigen score was 0.5 overall. The finding imply that chimeric protein is non-immunogenic and can be used as a promising therapeutic target for cancer.

Docking analysis

In order to predict interaction between designed chimeric protein and its IL-15Rα/IL-15Rβ/γc receptor, molecular docking was performed by using the Cluspro online server. ClusPro operates on the Fast Fourier Transform Correlation technique, it assess the docked complexes in three stages. First, the rigid docking is carried out by analysing the conformations in billions with the aid of PIPER.Second, it predicts the biggest cluster and displays the most probable docked complex by creating clusters of one thousand lowest energy structures based on RMSD (root mean square deviation). The third stage involves energy reduction to stabilise the chosen docked model (Kozakov et al. 2017). Among 10 models generated by ClusPro server (Table 4), the selection of the best structure for protein-protein interaction studies was guided by factors such as model scores, electrostatic energies, and van der Waals attraction as described previously (Muhammad Rehman et al. 2023). Cluster 0 stands out with 189 members, the highest among all clusters, indicating a robust and consistent set of docking poses that suggest a reliable interaction model. The weighted score of -617.2 and the lowest energy score of − 801.0 for Cluster 0 highlight a highly favorable and stable binding affinity, essential for effective therapeutic targeting. The lowest energy score indicates the most stable conformation of the protein complex, signifying that Cluster 0 (Fig. 6B) can achieve a highly stable interaction with the receptor, which is crucial for the therapeutic efficacy of the fusion protein. In comparison, Cluster 1, with 58 members and a weighted score of − 676.0, and Cluster 2, with 53 members, also show strong interactions with a weighted score of -815 and a lowest energy of − 833. A more negative docking score indicates a higher likelihood of binding, while a higher cluster number enhances confidence in the binding model. Despite the favorable scores of Cluster 2, the significantly lower member count reduces its reliability. Clusters 3 to 10 have progressively fewer members and varying energy scores, which, while sometimes favorable, do not match the reliability of Cluster 0 due to their lower frequency of occurrence. The selection of Cluster 0 is further justified by the balance between a high member count and excellent energy scores, highlighting its potential as the most stable and reproducible interaction between the fusion protein and the IL-15Rα/IL-15Rβ/γc receptor.

Fig. 6

Docked complex and interacting residues involved hydrogen bonds (blue lines), salt bridges (red line) and disulphide bonds (yellow lines), in binding of IL15-NGR fusion peptide with its cognate receptor

Table 4 Cluster score of docked complexes

The docking of the IL-15-NGR fusion protein with the IL15Rα/IL15Rβ/γc receptor using Hdock yielded notable results. The top ten models displayed docking scores ranging from − 245.18 to − 217.88 (Table 5), with the highest confidence score being 0.8703. These negative docking scores suggest favorable binding interactions between the fusion protein and the receptor. The ligand RMSD values ranged from 26.01 Å to 53.14 Å, indicating variability in the predicted binding poses and suggesting flexibility in the interaction sites. The ProQ quality assessment of the input models showed an LGscore of 2.288 and a MaxSub score of 0.205 for the receptor, and an LGscore of 2.944 and a MaxSub score of 0.341 for the ligand. These scores fall into the “good” to “very good” range, indicating that the structural models used for docking were of high quality. The selection of ClusPro-generated docked model was based on binding energy, the number of interactions, and the stability of the docked complex. These factors made the ClusPro generated complex superior, prompting its selection for further interaction analysis and molecular dynamics simulation studies.

Table 5 Top 10 docked models summary on Hdock serverProtein-protein interaction studies

Estimating the theoretical binding energies among amino acid residues at protein-protein interfaces presents a challenging task. With advancements in bioinformatics computations, a variety of software and online servers are now available, employing empirical and analytical approaches. However, it’s important to note that results from each module may vary, influenced by the submitted 3D structure of the protein under investigation (Lubkowski et al. 2018). We utilized the PDBsum program to gain insights into the docked complexes and analyze the forces of attraction at the molecular level. The PDBsum server facilitated the analysis of interactions within the docked complex. The PDB file of the docked complex was submitted to the PDBsum online server for analysis. The length of the interface area and the number of residues therein determine the strength of binding; a higher percentage of amino acids with a longer interface area indicates greater binding strength and a stable conformation (Khan et al. 2021). In the docked complex, a total of 21 hydrogen bonds and 6 salt bridges were identified (Table 6). In the Fig. 6A and B, C, chain D represents the newly designed fusion peptide, chain A represents the IL-2/IL-15Rβ receptor subunit, chain C represents the IL-15Rα receptor subunit, and chain B represents the IL-2/1L-15γc receptor subunit. Between chain D and chain B, hydrogen bonds of chain D are contributed by Glu64, His60, Asp61, Ser83, Thr81, Asp56, Ser51, Ser54, Glu53, Leu52, Gln48, Leu45, Glu89, and Glu87. Hydrogen bonds of chain B are contributed by Arg167, Asn198, Asn20, Tyr18, Asn203, Lys201, Arg193, Arg171, Arg173, and Trp189. Additionally, hydrogen bonds of chain D are contributed by Asp22, while the hydrogen bond of chain C is contributed by Lys51. Furthermore, between chain D and chain B, salt bridges of chain D are contributed by Glu64, Glu89, Glu87, and Lys41, while salt bridges of chain B are contributed by Arg167, Arg171, Arg173, and Glu188. These hydrogen bonds and salt bridges play central roles in stabilizing the interaction between the fusion peptide and the receptor subunits, providing insights into the molecular mechanisms underlying their binding affinity. The binding energy calculation of the IL-15-NGR fusion protein with the IL15Rα/IL15Rβ/γc receptor, as predicted using the HawkDock online server, is − 2.98 kcal/mol. This negative binding energy indicates a favorable interaction between the fusion protein and the receptor, suggesting that the IL-15-NGR fusion protein can effectively bind to the IL15Rα/IL15Rβ/γc receptor. The magnitude of the binding energy, while modest, supports the stability of the complex formed. This binding energy is crucial for the therapeutic potential of the fusion protein, as it underlines its ability to target and engage the receptor with sufficient affinity. However, while computational predictions provide valuable insights, these results should be further validated through experimental binding assays to confirm the interaction strength and to refine the binding efficiency. Finally the results of binding affinity prediction on PRODIGY server suggest that the IL-15/NGR fusion protein binds to its receptor with high affinity at 37 °C. The strongly negative ΔG value of − 12.6 kcal mol−1 and the low Kd (constant of dissociation) of 1.2e-09 M, both indicate a strong and stable interaction, suggesting that the fusion protein is likely to effectively engage its receptor under physiological conditions.

Table 6 Interface statistics of chain A, B, C, DMD simulation analysis

Molecular Dynamics (MD) simulations are performed to refine docking results by allowing for conformational flexibility and dynamic adjustments. It helps to assess the stability of the complex over time, investigate key binding interactions, and observe conformational changes that affect binding affinity. The graph in Fig. 7A depicts the radius of gyration (Rg) of a docked complex over a 100-nanosecond (ns) simulation period, indicating changes in the compactness of the molecular structure. Initially, from 0 to 20 ns, the Rg increases from approximately 3.1 nm to around 3.3 nm, suggesting an initial expansion of the complex. Between 20 and 60 ns, the Rg continues to rise and fluctuates around 3.4 nm, indicating further expansion and structural adjustments within the complex. In the later phase, from 60 to 100 ns, the Rg stabilizes and fluctuates between 3.2 nm and 3.4 nm, reflecting a balance between expansion and contraction as the complex reaches a more equilibrated state. These observations suggest that the complex undergoes significant initial conformational changes and expansion, followed by a period of relative stability with minor fluctuations in its compactness, indicative of the complex’s dynamic adaptation and stabilization during the simulation. The Fig. 7B illustrates the Root Mean Square Deviation (RMSD) of molecular system over a simulation period of 100 nanoseconds (ns), depicting the structural changes relative to the initial configuration. Initially, from 0 to 20 ns, the RMSD increases sharply from close to 0 to about 0.8 nm, indicating significant structural adjustments as the system moves away from its starting conformation. Between 20 and 60 ns, the RMSD continues to rise but at a slower rate, suggesting ongoing but less dramatic changes in the structure. In the final phase, from 60 to 100 ns, the RMSD stabilizes, fluctuating around 1.0 to 1.2 nm, which signifies that the system has reached a more equilibrated state with only minor fluctuations. This trend reflects the typical behavior in MD simulation where system undergoes an initial period of rapid structural changes before settling into a stable conformation. The graph in Fig. 7C illustrates the area of the interaction interface (in nm²) between two molecules over time (ns). Initially, during the first 20 ns, the interface area fluctuates between approximately 350 and 370 nm², indicating that the complex is undergoing initial adjustments and conformational changes to find a more stable binding orientation. From 20 to 60 ns, the interface area shows a slight downward trend, fluctuating around 350 to 360 nm², suggesting that the complex is becoming more stable as the interacting surfaces of the molecules adjust to form a tighter interface. In the final phase of the simulation, from 60 to 100 ns, the interface area continues to decrease slightly, with fluctuations between 330 and 350 nm², indicating that the complex is reaching a more stable conformation. The overall decrease in the interface area over time suggests that the molecules are coming closer together, enhancing their interactions, while the fluctuations highlight the dynamic nature of the molecular interactions. The graph in Fig. 7D illustrates the Root Mean Square Fluctuation (RMSF) of a docked complex over a 600-nanosecond (ns) simulation period, providing insights into the flexibility of different regions within the complex. Initially, from 0 to 100 ns, the RMSF values show significant fluctuations, reaching up to approximately 0.8 nm, indicating that certain regions of the complex exhibit considerable movement. In the middle phase, between 100 and 400 ns, the RMSF values display periodic peaks and troughs, with a pronounced spike around 200 ns where the RMSF peaks at about 1.4 nm, suggesting that specific regions of the complex experience increased flexibility, possibly due to conformational adjustments or interactions with the receptor. In the later phase, from 400 to 600 ns, the RMSF values generally remain below 0.6 nm with occasional fluctuations, except for a notable increase towards the end, peaking again around 1.4 nm, indicating that while most regions of the complex remain relatively stable, certain areas experience heightened flexibility. These observations suggest that the docked complex undergoes dynamic structural changes, with specific regions exhibiting significant flexibility, reflecting its interaction and adaptation processes during the simulation. The MD simulation results indicate that the docked complex undergoes significant conformational adjustments initially, achieving greater stability over time. Overall, these observations reflect the complex’s dynamic adaptation and eventual stabilization during the simulation.

Fig. 7

MD simulation graphs of RMSD, Rg, RMSF and SASA generated by GROMACS for IL15-NGR fusion peptide