Chimeric proteins (CP) included ricin_A, a linker, and an approved ErbB2-specific penetrating peptide (Figure 2A). First, the amino acid composition of ricin_A was retrieved from PDB ID 2AAI in FASTA format (Rutenber et al., 1991). Then, either a flexible (GGGGS) or rigid (EAAAK) linker, and a peptide (KCCYSL (Ringhieri et al., 2017), MARAKE (Houimel et al., 2001), WYSWLL (Karasseva et al., 2002), MARSGL (Sugo et al., 2013), MSRTMS (Houimel et al., 2001), and WYAWML (Karasseva et al., 2002)) were added to the end of ricin_A sequence manually in a text editor (Figure 2B). The peptide sequences were obtained from the Immunet BDB database (http://immunet.cn/bdb).

3D structures of all templates were predicted by the MODELLER software 15.9 (Eswar et al., 2006) based on the most similar available conformations. The best final models were selected after their qualities were checked by the ERRAT (Colovos and Yeates, 1993) and VERIFY3D servers. Moreover, Ramachandran diagrams were obtained using the PROCHECK server (Laskowski et al., 1993) (Supplementary Figure S3).

Molecular dynamics (MD) simulations were employed to investigate whether the ricin_A moiety of CPs retains the conformation it had in ricin since any significant change in ricin_A may affect the ribosome-inactivating activity of CPs. MD simulations were performed using the GROMACS package, version 2020 (Van Der Spoel et al., 2005), and Charmm36 forcefield for making topology files. All systems were solvated in a cubic box with a minimum distance of 1.0 nm from the edges of the box and filled with SPC water molecules (Mark and Nilsson, 2001) as a reliable water model for aqueous solutions of biomolecules (Zielkiewicz, 2005).

Energies of all systems were minimized for 50,000 steps, followed by a thermal equilibrium step using the Berendsen thermostat at 310 K. Finally, systems were subjected to 100-ns production simulations. The MD trajectories were saved every 10 ps. The pressure was equilibrated for 1 ns to achieve the pressure of 1 bar using the Berendsen barostat. LINCS (Hess et al., 1997) and PME mesh (Darden et al., 1993) algorithms were used to constrain the bond parameters and the calculation of long-range electrostatic interactions, respectively. During the simulations, the Fourier grid spacing and Coulomb radius were set at 0.16 and 1.2 nm, respectively. The cutoff radius for van der Waals interactions was 1.2 nm.

It is proved that a protein’s activity depends on its dynamic behavior (Daniel et al., 2003; Torgeson et al., 2022). Therefore, CPs whose ricin_A moieties showed the greatest similarity in behavior during the MD simulations to that of chain A were considered to have a minimal change in their ribosome-inactivating activity and were subjected to molecular docking simulations against ErbB2. Docking was carried out by the HADDOCK server (Dominguez et al., 2003) which is a fully automated server designed for protein-protein docking simulations. It requires PDB files of the proteins as input. The docking score of this server is a linear sum of energy terms such as van der Waals, electrostatic, desolvation, and restraint violation energies the higher quantities of which indicate greater energy constraints of complex formation.

To validate the docking studies, trastuzumab was docked against ErbB2. The X-ray crystallized PDB ID 6OGE (Hao et al., 2019) was obtained from the RCSB data bank. It included trastuzumab, Fab-Pertuzumab, and ErbB2. After making sure that HADDOCK produced the same interactions and orientations of trastuzumab, CPs were docked and the trastuzumab-binding region of ErbB2 was introduced (Hao et al., 2019) as the active residues to the server while the remained parameters remaining as default. Then, the complexes with the lowest binding scores were selected. Finally, we used PRODIGY server (Xue et al., 2016) to decipher the potential binding affinity of each CP for the receptor.

The toxic potential of CPs was examined by the ToxDL server which is devised to predict toxic domains in protein structures by deep learning (http://www.csbio.sjtu.edu.cn/bioinf/ToxDL/) (Pan et al., 2021).

Six peptides with the capability of binding to the trastuzumab-binding pocket were connected to ricin_A by either a rigid or a flexible linker (Figure 2) which resulted in 12 initial templates. The rigid linker (EAAAK) is an α-helix maker previously used to design fusion proteins against various cancers such as breast (Thiesen et al., 1987; Epler et al., 2012) and cervical (Park et al., 2022) cancer. It contains a salt bridge between its glutamic acid and lysine which keeps a fixed distance between the peptide and ricin_A, leading to the maintenance of independent functions of these compartments by creating a stable helix structure (Zielkiewicz, 2005).

The flexible linker included smaller GGGGS residues (Klein et al., 2014). Flexible linkers are usually applied when the connected domains require a certain degree of movement or interaction. This linker has been designed for recombinant fusion proteins and is used to treat some cancers (Valiyari et al., 2020).

Two main considerations of designing proteins are their foldability and functionality which in most cases are related to each other. Here, although the length of peptides and linkers are shorter compared with ricin_A, it is not unlikely they can disturb the foldability of CPs and hence negatively impact their ribosome inactivating function. To elucidate any potential effect of adding such sequences, we used homology modelling and molecular dynamics simulations.

To elucidate the 3D conformations, all of twelve initial templates were subjected to homology modeling. Models with the best dope scores (the lowest scores) were selected as final 3D structures. Moreover, the quality of models was determined by three programs ERRAT (Colovos and Yeates, 1993), VERIFY3D (Eisenberg et al., 1997), and PROCHECK (Laskowski et al., 2006) on the SAVES server (Tables 1, 2; Supplementary Figures S1, S2). ERRAT evaluates non-bonded atomic interactions, and its higher scores indicate higher quality. The generally accepted score of >50 indicates a high-quality model (Table 1). VERIFY3D determines the compatibility of an atomic model (3D) with its amino acid sequence. It assigns a structural class to each residue (alpha, beta, loop, polar, non-polar, etc.) and compares the results to other structures. VERIFY3D scores of higher than 80% means the acceptable quality of the model (Table 1). Ramachandran plots of the models (Supplementary Figure S3) were depicted, and their statistics are shown in Table 2. It can be seen that at least 83.7% of residues were in the most favorable regions and maximally 1.3% in disallowed regions, suggesting an acceptable quality of the models.

The quality of the homology model results indicates that the addition of linkers and peptide sequences had no negative effect on foldability of ricin_A.

To make sure that the integrity, structure, and subsequent inhibitory activity of ricin_A moiety is not affected by the addition of linkers and peptides, we performed MD simulations. CPs whose ricinA moiety’s behavior was the most similar to that of free ricinA were then subjected to molecular docking.

The RMSD of ricin’s backbone atoms compared to the initial structure as a reference is shown in Figures 3A, B, and their average values are shown in Figure 4A. In the CP12 structure, ricin_A has the least stable structure as seen from the severe fluctuations in its plot. The ricin_A moiety of CP12 also showed the greatest divergence from ricin_A in ricin compared to other CPs, with an average RMSD of 0.8 (Figure 4A). Moreover, except for CPs 2, 5, 7, 8, and 10, which had the closest plots and average values (0.35, 0.33, 0.38, 0.36, and 0.36 nm, respectively) to ricin, other CPs diverged from ricin dramatically, suggesting a greater structural change in ricin_A moiety when it is joint with peptides and linkers in these CPs.

To further explore ricin’s significant structural changes in each CP, further analyzes were carried out. It is shown in Figures 3C, D that ricinA‘s radius of gyration (Rg) in CPs 1, 3, 4, and 12 (with the average values of 1.99 nm for CPs 1, 3, and 4 and 2.03 nm for CP12) diverged the most from that of ricin (with the average Rg of 1.89 nm) (Figure 4B) suggesting that the addition of their relevant recombinant moiety may induce significant structural change in ricin_A structure. This matter increases the possibility of losing the ricin_A inhibitory effect in these CPs. On the other hand, the remaining CPs had favorable Rg values (Figure 4B).

SASA (solvent-accessible surface area) analysis is known as an indicator of the surface area of the protein. Increased values of SASA suggest the expansion of the protein structure. SASA plots of all systems are shown in Figures 3E, F and their average quantities are depicted in Figure 4C. Compared with free ricin (SASA = 136 nm), it can be seen that CPs 2, 5, 6, 7, 8, and 10 have the closest SASA values to free ricin. This suggests that the ricin moiety of these CPs has a similar surface to free ricin. The Rg and SASA plots show consistency so that CPs with a significant increase in ricin’s Rg exhibited a more divergent SASA than the free ricin.

Furthermore, RMSF analysis was performed to understand better the mobility of the ricin part in each CP. Figures 3G, H show that residues 175–260 in the ricin part of CPs 4, 7, and 12 have the most significant movements during the simulations. This suggests that adding a linker and peptide to ricin_A increased the flexibility of this region compared to single ricin. Ricin_A has three main regions: domain I consist of β -sheets, while domain II is an α -helical structure. Domain III plays a major role in dimer formation by binding to chain B. The active site of ricin contains highly conserved residues (Tyr80, Tyr123, Glu177, Trp211, and Arg180 (Eswar et al., 2006); see Figure 3I). As seen from the RMSF plots in Figure 3, the most fluctuating region of CPs 4, 7, and 12 is located within ricin’s active site, which can diminish the catalytic activity of these complexes. Therefore, this matter can be another reason for removing these CPs from our test cases.

Furthermore, principal component analysis (PCA) of Cα atoms of ricin was employed to understand how the protein backbone in CPs behaves during the simulations (Figure 5). The results show that CPs 1, 3, 4, 11, and 12 have greater overall motions than ricin_A since they cover a broader part of their conformational space. This implies that adding linkers and peptides increases collective motions in the ricin compartment, which in turn may affect ricin’s functionality. The remaining CPs that showed similar behaviors to ricin in the previous MD analyses had smaller conformational space coverage, indicating their limited collective motions. This suggests that adding the linker and peptides in these CPs makes ricin_A rather rigid.

According to the analyzes performed in the MD section (RMSD, RMSF, Rg, SASA), CP2, 5, 7, 8, and 10 had the most stable structure and were the most similar to free ricin_A. Thus, we chose them for the docking step. It is worthy of note that none of the designed CPs had toxic potential according to the ToxDL server results (Supplementary Table S1). Keeping a balance between having toxicity for tumors and safety for normal cells might be a challenge in protein design. However, the proposed CPs make the common solutions like designing a carrier unnecessary because while they are safe, they may still have their ribosome inactivating potential. Another challenge that needs to be addressed is the specificity and internalization of CPs. Instead of using a carrier, we hypothesized that trastuzumab binding cavity of ErbB-2 can be exploited as a natural carrier for designed CPs. Moreover, this affords the opportunity of providing specificity for CPs. We examined this idea with molecular docking.

MD results indicated that ricin_A of CPs 2, 5, 7, 8, and 10 may most likely retain its inhibitory activity. Using molecular docking, we aimed to examine whether the selected CPs were capable of binding to trastuzumab-binding pocket on ErbB-2 since this region internalizes trastuzumab upon its binding (Austin et al., 2004) in a Caveolae/Lipid-Raft Mediated mechanism (Liang et al., 2021). The last frame of each MD trajectory entered the docking step. We used HADDOCK and PRODIGY servers to find best potential binding modes and ΔG for each CP, respectively. The trastuzumab-binding pocket of ErbB-2 includes its domain IV residues Pro579, Glu580, Asp582, Gln583, Lys591, Asp592, Pro593, Pro594, Phe595, Asp607, Leu608, Tyr610, Lys615, Gln624, Cys626, and Pro627 (Figure 6). The results indicated that except CP8, other CPs showed higher docking scores compared to the well-known anti-cancer agent, trastuzumab (with the binding score of −90.), suggesting that they may have lower energy constraints for binding to the receptor (Table 3). Regarding the binding affinity (ΔG), CPs 2 and 10 had the maximum quantities suggesting their highest potential for binding (Table 3). This can be supported by the highest numbers of hydrogen bonds these CPs established (15 and 10 hydrogen bonds, respectively; Figures 7C, D). Moreover, it can be seen that CPs 2 and 10 had higher binding affinity compared with trastuzumab (ΔG = −10.3 (kcal.mol−¹)).