In silico designing of a new cysteine analogue of hirudin variant 3 for site specific PEGylation.

Hirudin is an anticoagulant agent of the salivary glands of the medicinal leech. Recombinant hirudin (r-Hir) displays certain drawbacks including bleeding and immunogenicity. To solve these problems, cysteine-specific PEGylation has been proposed as a successful technique. However, proper selection of the appropriate cysteine residue for substitution is a critical step. This study has, for the first time, used a computational approach aimed at identifying a single potential PEGylation site for replacement by cysteine residue in the hirudin variant 3 (HV3). Homology modeling (HM) was performed using MODELLER. All non-cysteine residues of the HV3 were replaced with the cysteine. The best model was selected based on the results of discrete optimized protein energy score, PROCHECK software, and Verify3D. The receptor binding was investigated using protein-protein docking by ClusPro web tool which was then visualized using LigPlot+ software and PyMOL. Finally, multiple sequence alignment (MSA) using ClustalW software and disulfide bond prediction were performed. According to the results of HM and docking, Q33C, which was located on the surface of the protein, was the best site for PEGylation. Furthermore, MSA showed that Q33 was not a conserved residue and LigPlot+ software showed that it is not involved in the hirudin-thrombin binding pocket. Moreover, prediction softwares established that it is not involved in disulfide bond formation. In this study, for the first time, the utility of the in silico approach for creating a cysteine analogue of HV3 was introduced. Our study demonstrated that the substitution of Q33 by cysteine probably has no effect on the biological activity of the HV3. However, experimental analyses are required to confirm the results.

INTRODUCTION

Hirudin is an anticoagulant agent of the salivary glands of the medicinal leech; Hirudomedicinalis. This protein is a single polypeptide that contains 65–66 amino acids and has a molecular weight of approximately 7KDa. It also has three intra-molecular disulfide bands between the residues 6–14, 16–28, and 22–39(12). Because of the pharmaceutical usefulness of hirudin and its limited availability in a natural form, the production of recombinant hirudin (r-Hir) is of great interest(34). Nevertheless, recombinant hirudin displays a short plasma half-life of approximately 60–100 min; therefore, several injections are needed, which makes it more expensive for patients(5). Furthermore, it has certain drawbacks, including bleeding and immunogenicity(6).

To resolve these issues, PEGylation has been proposed as a successful technique to improve the half-life and immunogenicity(78).

While PEGylation can increase the size, drug solubility, and bioavailability of r-Hir, it also decreases the dose frequency and immunogenicity(9).

Among the different PEGylation methods, lysine PEGylation is most commonly used; however, due to the high prevalence of this amino acid in the molecule and the low chemical selectivity of the amine group, positional isomers occur(57).

A potential result of non-specific PEGylation is a decrease in biologic activity(10). Thus, site specific PEGylation is an effective method to obtain a mono-PEGylated form. Because the number of free cysteine residues on the protein's surface is lower than that of other conjugable groups (i.e. lysine residues), free cysteine PEGylation is considered to be an efficient tool for the enhancement of the pharmacokinetic properties of therapeutic proteins(11).

The three-dimensional structure of protein must not be altered by the insertion of cysteine residue(12) ; however, increasing the potential for disulfide bond formation and protein dimerization are drawbacks of the addition of free cysteines via genetic engineering(13). Therefore, the proper selection of the appropriate sites for substitution with cysteine residue is a critical step in the process(1012). The in silico design of new variants with the proper site for PEGylation is currently considered a useful method, and it is therefore frequently used in studies of folding, kinetics, and different interactions of proteins.

The results of computational studies have revealed insights into the mechanism of protein dynamics. These studies have also helped identify a variety of structural and dynamical characteristics of protein interactions(141516).

For the first time, a computational approach was used in this study that aimed to identify a single potential PEGylation site for replacement by cysteine residue in the hirudin variant 3 (HV3; hirPA).

MATERIALS AND METHODS

Homology modeling and design mutant library

To the best of the authors’ knowledge, a three-dimensional structure for HV3 has not been identified; therefore, HM was performed using MODELLER version 9.11 (http://www.salilab.org) (17). Sequence-structure alignment against the protein data bank (PDB) database was carried out using the HHpred toolkit(18) to determine the most appropriate template. Ten thousand models were generated using MODELLER version 9.11(17), and a crystal structure of hirudin variant 1 (HV1) (PDB code: 1HRT, chain I) from the protein data bank (http://www.rcsb.org) was used as a template. To identify the best location for cysteine-specific PEGylation, all non-cysteine residues of the HV3 (60 amino acids) were replaced with cysteine. Among the 10,000 generated molecules of the cysteine analogs of HV3, the five best models of each mutant were selected based on their discrete optimized protein energy (DOPE) score(19). The quality of the models was then checked by a torsion angles analysis of the protein backbone with PROCHECK software(20) and Verify3D(21).

Receptor binding

To investigate the receptor (thrombin) binding of new mutants, protein-protein docking was applied by using the ClusPro web tool(22), which functions based on the surfaces’ complementarities. Subsequently, the resulting models were sorted by clustering based on their docking score, as described by Comeau, et al.(23). A two-dimensional interacting residues diagram was visualized using LigPlot+ software(24) and the 3D structure was drawn by means of PyMOL(25).

Multiple sequence alignment

The FASTA format of hirudin variants from the universal protein resource (http://www.uniprot.org)(26) was used in multiple sequence alignment (MSA) using ClustalW software with the blocks substitution matrix (BLOSUM62) scoring matrix and default parameter settings. For the alignment visualization and the secondary structure elements’ predictions, JalView and Jpred softwares were applied, respectively(2728).

Disulfide bond prediction

To predict the bonding state of cysteines in hirudin, online software that used the different algorithms, including Dipro(29), CYSPRED(30), Disulfind(31) and Metadetector(32), were employed.

RESULTS

For this study, the computational methods were employed to determine the best site for cysteine specific PEGylation in HV3. HM was applied to predict the 3D structure of HV3 and its binding to the receptor.

The best template for HM was obtained by sequence-structure alignment using the HHpred toolkit. The HV1 (PDB code: 1HRT), which had an 85% identity and 1.518 similarity, was selected as the best template for HM Table 1. Among 10,000 generated models, the best five models were selected based on their DOPE score (Table 2). The best model (Fig. 1A) was then subjected to further analyses, including stereochemical conformation checking (Ramachandran plot) and Verify 3D, to identify the association between the 3D structure and the sequence. ProsaII was then applied for structural quality validation. It was found that the covalent bonds and the related angles were acceptable (Fig. 1B and 1C).

Table 1

HHpred toolkit analysis

Table 2

The best models of HV3 based on the DOPE score and the related analyses

Fig. 1

(A) 3D structure of the best model of HV3. The structure contains a 3–10 helix in the C-terminal, an anti-parallel beta sheet, and an unstructured region in the C-terminal. No alpha helix is observed. (B) Stereochemical analysis of the model HV3-normal is drawn from PROCHECK. The plot confirms the stereochemical quality with 90.2% of residues fall in the most favored region and 7.8% in allowed region and 2% in generously allowed region. No residue was located in the disallowed region. (C) Verify 3D analysis of HV3. As the scores were higher than zero, the model correctness was confirmed.

PyMOL software was employed to visualize the selected model. As a result, a folded globular anti-parallel beta sheet was shown at the N-terminal domain, and an unstructured tail with a 3–10 helix was found at the C-terminal region (Fig. 1A).

To design the cysteine variants, single amino acid replacement was performed using HM. Then, the best models were introduced to Cluspro to investigate the binding affinity of HV3 to the receptor. Finally, the best mutant for cysteine PEGylation was selected based on structural quality and binding energy in docking. The HM and docking showed that the Q33C was the best site for PEGylation. Because the Q33 is located on the surface of the protein, it would be an appropriate site for the PEGylation (Fig. 2A). The surface accessibility for Q33 was defined using the visualization with PyMOL. After further analyses using Ramachandran plot and Verify3D, the accuracy of the model was confirmed; the quality and binding pattern of this variant were similar to the normal pattern (Fig. 2B and 2C).

Fig. 2

(A) 3D structure of the HV3-Q33C. Cys33 is located on the surface of the protein so, it would be accessible for the PEGylation. (B) Stereochemical property analysis of the model HV3-Q33C is drawn from PROCHECK. The plot confirms the stereochemical quality of the model, with 92.2% of residues fall in the most favored region and 7.8% in allowed region. No residue is located in the disallowed region. (C) Verify 3D analysis of HV3-Q33C. As the scores were higher than zero, the model correctness was confirmed.

The 2D visualization of hirudin-thrombin using LigPlot+ software showed that the binding pocket had not changed. In other words, Q33C was not involved in ligand-receptor binding.

For further confirmation, MSA was carried out using ClustalW software. Different variants of hirudin were subjected to MSA. As can be found from the Fig. 3, the Q33 residue (code: 09944) was not conserved; therefore, it was not likely to be essential to the protein function and could be changed into another amino acid.

Fig. 3

Multiple sequence alignment using ClustalW software. The Q33 residue is shown not to be a conserve residue (code: 09944).

Tha disulfide bond formation prediction was performed using different online software. The results showed that the hirudin did not use the Q33C for disulfide bond formation. Hence, the native disulfide bonds of hirudin did not change after the cysteine insertion.

DISCUSSION

In silico designing has been used for different proteins(333435). In this regard, we applied this approach for site specific PEGylation, a well-known method used to prolong the circulating half-lives of protein therapeutics and to improve both their stability and efficiency(363738). The primary challenge is determining the best conformational location for PEGylation. Cysteine-specific PEGylation includes introducing a ‘free’ cysteine amino acid into a target protein that has not participated in a disulfide bond formation. The PEGylation reaction is carried out with a cysteine-specific PEG reagent. A major technical problem with this method is the identification of a proper site in the protein where the substituted or inserted cysteine residue can be easily PEGylated and thus have no effect on the biological activity of the protein(39). To perform enzymatic PEGylation by microbial transglutaminase in G-CSF, two different mutations, P132Q and Q134N, were selected using bioinformatic approach. It was shown that an in silico study is a useful approach for identifying specific sites for PEGylation(15). Cysteine specific PEGylation is performed for different proteins, such as interferon α-2(39) and GM-CSF(40). Hirudin, a recently FDA-approved drug, has been PEGylated with different PEGylation reagents, including PEG-NHS ester disulfide(41), SC-mPEG(4243), and methoxyl-polyethylene glycol (mPEG)(44). For the first time, the utility of the in silico approach for creating a cysteine analogue of hirudin variant 3 has been introduced in this study. The results of this study showed that Q33 would be the best residue for substitution with cysteine. As can be observed from the Ramachandrn plot, the mutant model's quality is quite similar to the normal variant. Furthermore, because Q33 is located on the surface of hirudin (Fig. 2A), it would be easily accessible for the cysteine specific PEGylation reagent. In addition, based on the MSA analysis, Q33 is not a conserved residue of the hirudin molecule; therefore, it probably has no significant effect on the biological activity of the protein. Further confirmation of the appropriateness of selecting the Q33 residue was provided by the results of the online software used including DIpro(29), CYSPRED(30), Disulfind(31), and Metadetector(32), which showed the binding state of cysteines in the proteins. These software programs employed different algorithms. In this regard, Q33 was not involved in disulfide bond formation. Moreover, there has been another study in which the Asp33 in the hirudin variant 1(HV1) was changed to lysine, and PEGylation was successfully completed(8). The LigPlot+ also revealed that the Q33 was not involved in the hirudin-thrombin binding pocket. This is in accordance with previously reported findings(45).

CONCLUSION

This was the first study in which the in silico approach was applied for creating a cysteine analogue of HV3. Our study demonstrated that the substitution of Q33 by cysteine probably has no effect on the biological activity of the HV3. However, experimental analyses are required to confirm the results. In conclusion, an in silico study can be a fast and useful approach for predicting a proper site for the cysteine-specific PEGylation of recombinant hirudin.