Structure prediction of transferrin receptor protein 1 (TfR1) by homology modelling, docking, and molecular dynamics simulation studies.

Transferrin receptor protein 1 (TfR1) is an important molecule in anti-cancer therapy. Targeted delivery of such therapeutic compounds improves their cellular uptake and circulation time, thereby enhancing therapeutic efficacy. Drug designing is therefore used to engineer molecules with structures that facilitate specific interactions. However, this process requires a thorough knowledge of all the interactions, including the three-dimensional (3D) and quaternary structures (QS) of the interacting molecules. Since structural information is available for only a part of the full TfR1 sequence, in the present study, we predicted the whole structure of TfR1 using homology modelling, docking, and molecular dynamics simulations. Homology modelling is used to generate 3D structures of TfR1 using MODELLER, I-TASSER, and RaptorX programs. Verify3D and Rampage server evaluated the quality of the resultant models. According to this evaluation, the model built by the RaptorX server and validated by Verify3D (compatibility: 83.82%) had the highest number of residues (95.5%) within the favoured regions of the Ramachandran plot, making it the most reliable 3D protein structure for TfR1 compared with others. The QS of TfR1 was built using HADDOCK and SymmDock docking software, and the results were evaluated by the ligand root mean square deviation (l-RMSD) value computed using the ProFit software. This showed that both HADDOCK and SymmDock gave acceptable results. However, the HADDOCK result was more stable and closest to the native complex structure with disulfide bonds. Therefore, the HADDOCK complex was further refined using both SymmRef and GalaxyRefineComplex until the medium l-RMSD rank was reached. This QS was successfully verified using nanoscale molecular dynamics (NAMD) energy minimization. This model could pave the way for further functional, structural, and therapeutic studies on TfR1.

Introduction

Transferrin receptor protein 1 (TfR1), is the primary receptor responsible for regulating cellular uptake of iron from transferrin (Ponka and Lok, 1999). TfR1 is expressed in all nucleated cells of the body but at different levels (Qian et al., 2002). Furthermore, it usually highly expressed in proliferating cells and cancers of the pancreas, colon, lungs, breasts, bladder, and lymphocytes. This reflects the need for iron as a cofactor in the enzymes involved in DNA synthesis resulting in cell proliferation and rapid cell division. Thus, this makes TfR1 an attractive target for developing anticancer strategies that can be used in combination with other therapeutic drugs (Daniels et al., 2006; Jeong et al., 2016; Peer et al., 2007; Qian et al., 2002; Richardson et al., 2009; Tortorella and Karagiannis, 2014). Additionally, TfR1 is the most widely studied receptor for targeted drug delivery to the brain because the transferrin-transferrin receptor complex is stable and does not undergo endosomal degradation (Patel and Patel, 2017). TfR1 is a homodimeric type II transmembrane protein consisting of two identical glycosylated subunits with an approximate mass of 95 kDa, each linked by two disulfide bonds to form a dimer. Each of its polypeptide subunits contain 760 amino acids, that are made up of a short N-terminal cytoplasmic domain (residues 1 to 67), a hydrophobic transmembrane domain (residues 68 to 88), and a large, globular extracellular C-terminal domain (residues 89 to 760) that contains the binding site for transferrin. Crystallographic studies of the ectodomain region have revealed that it possesses a butterfly-like shape with three subdomains (protease-like domain, apical domain, and helical domain) (Lawrence et al., 1999). Although previous studies have demonstrated the TfR1 structure of residues 121 to 760, the remaining structure has not been studied owing to the difficulties associated with the overexpression and subsequent crystallisation of transmembrane (TM) proteins. This also makes the X-ray diffraction or nuclear magnetic resonance (NMR) spectroscopy studies of the remaining residues difficult (Käll, 2010; Lawrence et al., 1999; Qian et al., 2002).

The unmodeled segment of TfR1 includes a short N-terminal cytoplasmic domain (residues 1–67) containing the internalization motif 20YTRF23 that was a part of the first natural variant (p.Tyr20His) caused by a missense mutation in TFRC, which codes for TfR1. The mutation impaired receptor endocytosis, resulting in a genetically undiagnosed form of combined immunodeficiency. This finding suggests a new role of TFR1 in immunity (Rothenberger et al., 1987; Collawn et al., 1993; Jabara et al., 2015). Furthermore, the cytoplasmic domain includes post-translational modifications such as acylation at Cys62 and phosphorylation at Ser24, which contribute to the regulation of endocytosis (Davis et al., 1986; Jing and Trowbridge, 1990). Additionally, it interacts with various proteins, including transferrin receptor trafficking protein, c-aminobutyric acid type A receptor-associated protein, adapter complex-2, ADP-ribosylation factor GTPase-activating protein with a coiled-coil, programmed cell death 6 interacting protein, ankyrin repeat and Pleckstrin homology domains, and heat shock protein HSC70 (Dai et al., 2004; Geminard et al., 2001; Green et al., 2002; Nesterov et al., 1999; Tosoni et al., 2005). The second part of the unmodeled segment of TfR1 is the transmembrane domain (residues 68–88), which anchors TfR1 to the membrane and guides signal polypeptide insertion into the membrane (Zerial et al., 1987, 1986). The final part is the stalk (residues 89–120) of about 30Ao that immediately follows the transmembrane pass and connects the ectodomain with the cytoplasmic domain, creating a sufficiently large gap that allows Tf molecules to bind. The stalk contains two intermolecular disulfide bonds—one formed by Cys89 and one formed by Cys98—covalently linking the monomers as well as the 0-linked glycan at Thr104 and thus contributing to receptor stability and activity. Additionally, a trypsin-mediated cleavage of the stalk at Arg100 can produce the soluble form of TfR1 (Dukovski et al., 2009; Hayes et al., 1992; Jingl and Trowbridgel, 1987; Rutledge et al., 1994; Williams and Ennss, 1991). Therefore, this study aimed to predict the whole structure of TfR1 for better understanding its structural, functional, and therapeutic behaviour.

Materials and methods

Target sequence

The TfR1 protein sequence used in this study, that contains 760 residues, was obtained from the Universal Protein Knowledgebase (UniProtKB, http://www.uniprot.org/) database in the Fasta format.

Physiochemical properties of TfR1

Protparam (https://web.expasy.org/protparam/) was employed for the calculation of different properties such as theoretical isoelectric point (pI), molecular weight, the total number of positively and negatively charged residues, and instability and aliphatic indices (Gasteiger et al., 2006).

Homology modelling

To generate 3D structure of TfR1 we used MODELLER 9.22 (Webb and Sali, 2016) (https://salilab.org/modeller/download_installation.html), I-TASSER (https://zhanglab.ccmb.med.umich.edu/I-TASSER/) (Yang and Zhang, 2015), and RaptorX software (http://raptorx.uchicago.edu/StructurePrediction/predict/) (Källberg et al., 2012; Peng and Xu, 2011). Next, the predicted models were validated using various validation techniques such as Verify3D (http://services.mbi.ucla.edu/Verify_3D/) (Bowie et al., 1991), which determines the compatibility of an atomic model (3D) with its amino acid sequence (1D) and Rampage server (http://mordred.bioc.cam.ac.ukrapper/rampage.php) to analyse the Ramachandran Plot (Lovell et al., 2003). Additionally, ModRefiner server (https://zhanglab.ccmb.med.umich.edu/ModRefiner/) (Xu and Zhang, 2011) was used to perform energy minimization and structure refinement of the modelled protein to remove unrealistically close steric clashes and significant deviations from ideal geometry resulting from the modelling process.

Molecular docking simulations

TfR1 is an asymmetrical homodimeric protein that contains a pair of identical subunits, but in the protein structure prediction step of the current study, only one subunit was predicted. To confirm its structure, an asymmetric docking simulation to build the TfR1 quaternary structure was performed using SymmDock (http://bioinfo3d.cs.tau.ac.il/SymmDock/) (Schneidman-Duhovny et al., 2005), and HADDOCK multi-body interface that requires a guru-level access (http://haddock.science.uu.nl/services/HADDOCK2.2/haddockserver-multi.html) (de Vries et al., 2010; Karaca et al., 2010). The interface information that we got from the TfR1 ectodomain crystal structure (2nsu.pdb) using PyMOL InterfaceResidues script (Schrodinger LLC, 2015), was converted to TaBLe information (TBL) file via the GenTBL server (http://milou.science.uu.nl/services/GenTBL/). The generated TBL file together with the TfR1 structure, were supplied to the multibody server as an input for the docking through the unambiguous distance restraints option, and the center of mass restraints was enabled. The number of structures was increased to 10000 for rigid body docking, 400 for semi-flexible refinement, and 400 for the explicit solvent refinement. Also, the C2 symmetry segment pair 1 was defined as the following symmetry segment: A and B (starts at residue 1 and ends at residue 760). Finally, all the other parameters were left at their default settings, and the docking was submitted (Karaca et al., 2010). The quality of the docking was computed according to the ligand Root Mean Square Deviation (l-RMSD), by using the ProFit software (http://www.bioinf.org.uk/profit/) (Bioinf.org.uk, 2018; McLachlan, 1982). In the C2 dimer of TfR1, its two subunits were arbitrarily labelled as A and B, and the l-RMSD was calculated by fitting the backbone atoms of the A subunit (fit zones 124–760) for both the reference structure (2nsu.pdb) as well as the mobile structure, i.e. the resultant models. Next, the l-RMSD was calculated with the backbone atoms of the B subunit (RMS zones 124–760) for both the reference structure and mobile structure. The results of the l-RMSD obtained were ranked according to the critical assessment of predicted interactions (CAPRI) criteria, as follows: incorrect if l-RMSD>10.0, acceptable if l-RMSD = 5.0–10.0, medium if l-RMSD = 1.0–5.0, or high quality if l-RMSD<1.0 (Lensink et al., 2007). The acceptable docking results from SymmDock and HADDOCK were selected to add the disulfide bonds. The formation of the disulfide bonds in the 3D model of TfR1 provides extra stability to the overall protein structure (Abkevich and Shakhnovich, 2000; Darby and Creighton, 1997; Ventura, 2008). TfR1 has about six disulfide bonds, two (Cys353-Cys363 and Cys556–Cys558) between the ectodomains of each monomer and two (Cys89–Cys98) between the two stalks of the dimeric receptor (Jingl and Trowbridgel, 1987; Lawrence et al., 1999). The best way to include the disulfide bonds in the TfR1 model is to use the model-disulfide.py script of MODELLER 9.22 software (Webb and Sali, 2016). Finally, the most stable complex was selected for refinement by SymmRef V1.2 and GalaxyRefineComplex that used to remove steric clashes and allow backbone and sidechain flexibility (Heo et al., 2016; Mashiach-Farkash et al., 2011).

Molecular dynamics (MD) simulations

The best structural model of TfR1 complex obtained from the docking refinement procedure was minimized and subjected to MD simulation using NAMD membrane protein tutorial (Aksimentiev et al., 2016). The production run was performed in NPT ensemble (constant number of particles, pressure, and temperature) with no restrictions for 80 ns. MD simulation of all TfR1 complexes in the 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane was performed in a solvent environment to reach an equilibrate state. Energy minimization for TfR1 was based on RMSD variation of Cα and the energy plot. The temperature of the simulation system was set to 310 K. The MD simulation was performed by nanoscale molecular dynamics (NAMD) software (http://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=NAMD) version 2.13 (Phillips et al., 2005), while the visual molecular dynamics (VMD) software version 1.9.3 (http://www.ks.uiuc.edu/Development/Download/download.cgi?PackageName=VMD) was used to display, analyze, and animate trajectories visually (Humphrey et al., 1996).

Results and discussion

In the present study, TfR1 (UniProtKB - P02786) were selected and characterized according to physicochemical properties by using Protparam (Table 1) (Gasteiger et al., 2006).

Table 1

Physical and chemical parameters of TfR1.

Physical and chemical parameters	TfR1
Number of amino acids	760 aa
Theoretical pI	6.18
Molecular weight	84871.38
Instability index	21.79
Aliphatic index	82.12
Total number of positively charged residues	(Arg + Lys): 83
Total number of negatively charged residues	(Asp + Glu): 89
Grand average of hydropathicity (GRAVY)	-0.243

The protein is classified as stable since the instability index computed to be 21.79 (Table 1).

The final predicted models of TfR1 that were generated through homology modelling were visualized using chimeraX (Figure 1) (Goddard et al., 2018). For the template structure for MODELLER (Basic and Advanced Modeling), we selected PDB files with the greatest identity to our query sequence in the BLAST analysis. We selected four PDB files, 3S9LA (PDB 3s9l, chain A) (99.69%), 1DE4C (99.84%), 6D03A (99.84%) and 1CX8A (99.53%). I-TASSER uses the TM-align structural alignment program to match the first model to all structures in the PDB library and the highest TM-score indicates the closest structural similarity. We obtained 3kasA with the highest TM-score (0.832). RaptorX uses the score to show the alignment score between 0 (the worst) and the domain length. RaptorX predicted 3 domains (domain1 = 643 residues, domain2 = 52 residues, and domain3 = 65 residues). The first one was created using the protein structure 3s9lA, with a score of about 701, which was slightly beyond the sequence length due to estimation error. The second one was created using the protein structure (4apsA, 2witA, 6an7C, 6cseM, 4cadC), with a mean score of about 30.2, and the final domain was created using the protein structure (4ckbA, 2xfaA, 4bedB, 5hstA, 5o60O), with a mean score of about 15.4. If high-quality templates are not available, RaptorX models the alignment using an inhouse free modeling program to generate final models. Unaligned portions of the template are also folded by free modeling (Källberg et al., 2012). Validation techniques such as Verify3D demonstrated that the models generated with I-TASSER and RaptorX had compatibility of 86.05% and 83.82%, respectively, indicating good quality models (Table 2). Additionally, the model built by RaptorX server (Figure 1) had the highest number of residues (95.5%) within the favoured regions of the Ramachandran plot (described in Table 3) (Bowie et al., 1991; Lovell et al., 2003). Therefore, this structure was found to be the most reliable among all the predicted models and used for the next step.

Figure 1

Homology modelling of TfR1 structure using basic MODELLER (a), advanced MODELLER (b), I-TASSER (c) and RaptorX (d). The red coloured region is the predicted structure. The images of the structures were visualised in ChimeraX.

Table 2

Verify3D results of the models.

Model	(Verify 3D) % residues with averaged 3D-1D score >0.2
MODELLER (Basic Modelling)	75.26% Fail
MODELLER (Advanced Modelling)	79.47% Fail
I-TASSER	86.05% Pass
RaptorX	83.82% Pass

Table 3

Summary of results from Ramachandran plot of the models.

Model	Region
	Favoured	Allowed	Outlier
Basic MODELLER	91.4%	5.7%	2.9%
Advanced MODELLER	93.8%	3.8%	2.4%
I-TASSER	82.7%	10.4%	6.9%
RaptorX	95.5%	3.7%	0.8%

Table 4 shows the rank and l-RMSD of the obtained top 12 complexes using HADDOCK and Symmdock docking software. According to the CAPRI ranks, the complexes (2,3,4,6,7,8,9,10,12) were incorrect since the l-RMSD was >10.0, while the complexes (1,5,11) were acceptable since the l-RMSD was between 5.0 and 10. To increase the stability of the quaternary structure, disulfide bonds were added using the model-disulfide.py script of MODELLER software for both acceptable docking results from SymmDock and HADDOCK. The HADDOCK result has a more stable structure with disulfide bonds; this model was subsequently refined using SymmRef and GalaxyRefineComplex, and then reranked by ProFit to obtain the final complex with a medium rank with an l-RMSD equal to 4.226 (between 1.0 and 5.0). This complex was selected to be the best docked complex for subsequent analysis (Figure 2) (Supplementary Video I) and was successfully refined using NAMD energy minimization (Karaca et al., 2010; Lensink et al., 2017; Phillips et al., 2005; Schneidman-Duhovny et al., 2005).

Table 4

Summary of results of l-RMSD from Profit for the models.

Protein-protein complex	Software	l-RMSD
complex -1	SymmDock	5.413
complex -2	SymmDock	42.781
complex -3	SymmDock	71.620
complex -4	SymmDock	37.164
complex -5	SymmDock	7.745
complex -6	SymmDock	34.934
complex -7	SymmDock	43.379
complex -8	SymmDock	77.290
complex -9	SymmDock	56.557
complex -10	SymmDock	46.238
complex -11	HADDOCK	6.731
complex -12	HADDOCK	46.328

Figure 2

The best structure of the TfR1 docked complex after refinement and formation of disulfide bonds. It is shown in molecular surface representation that generated in chimeraX. The two subunits are colored in red and green.

Supplementary content related to this article has been published online at https://doi.org/10.1016/j.heliyon.2020.e03221.

The following are the supplementary data related to this article:

TfR1

Microplastics ingestion process of goldfish when food was present.1

Molecular dynamics simulations

The TfR1 complex was energy minimized at -300000 kcal/mol. Moreover, Xmgrace software (Humphrey et al., 1996) displayed that the RMSD curves of the system were balanced after 40 ns, maintaining fluctuations between 6Å and 7Å. This suggested the TfR1 complex had reached the equilibrium state and achieved structural stability (Figure 3).

Figure 3

Represents the plotting of the backbone of root mean squared deviation (RMSD) versus time in nanosecond (ns) generated in Xmgrace software.

Conclusion

In this study, we used homology modelling, molecular docking simulations, and molecular dynamics simulations to determine the whole structure of TfR1. The results obtained indicate that the predicted model of TfR1 complex displays high quality and lacks any serious steric clashes. Therefore, this model could pave the way for further structural, functional, and therapeutic studies, specifically concerned with in silico studies. Additionally, this study can help in future mutational studies to pinpoint the roles of various amino acid residues and to predict the effects of mutations on TfR1 structure, and its biological functions.

Declarations

Author contribution statement

Maha Ateeq AL-Refaei: Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Rania Marwan Makki: Analyzed and interpreted the data; Wrote the paper.

Hani Mohammed Ali: Conceived and designed the experiments; Analyzed and interpreted the data.

Funding statement

This work was supported by King Abdulaziz City for Science and Technology (KACST) (Proposal number: 1-18-01-009-0155).

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.