Devendra Sir Cover Page copy.jpg BIBECHANA 17 (2020) 50-57 50 BIBECHANA ISSN 2091-0762 (Print), 2382-5340 (Online) Journal homepage: http://nepjol.info/index.php/BIBECHANA Publisher: Department of Physics, Mahendra Morang A.M. Campus, TU, Biratnagar, Nepal Constant velocity pulling and unfolding of thyroid hormone receptor by steered molecular dynamics Tika Ram Lamichhane, Hari Prasad Lamichhane* Central Department of Physics, Tribhuvan University, Kirtipur, Kathmandu, Nepal *Email: tikaramlamichh@gmail.com Article Information: Received: July 06, 2019 Accepted: October 7, 2019 Keywords: dynamicsmolecularSteered receptorhormoneThyroid Triiodothyronine Ligand binding domain; Protein unfolding ABSTRACT Unfolding pathways of T3 liganded thyroid hormone receptor (THRT3) can be studied by using the protocols of steered molecular dynamics (SMD). Theory of constant velocity pulling has been implemented to the structure of THRT3 in a neutral water-ion solution equilibrated up to 20 ns. The globular form of THRT3 is completely unfolded extending N-C termini from 38 Å to 876 Å at a constant speed of 0.1 Å/ps by means of 8.5 ns long SMD simulations. The peak force measured in the intermediate conformations is related to a burst of backbone H- bonds among -helices and -hairpins. With decrease in H-bonds, electrostatic energy increases by losing gradually the secondary structure and separating  and -strands in solution. The force at the end (t > 8.5 ns) increases steeply with the large increase in bond-angle and bond-length potentials when the system becomes completely unfolded. The hydrophobic ligand binding domain (LBD) of THR- with load bearing H-bonds protects T3 from water attack. Even after complete unfolding of THR- LBD, the position of T3 is not deviated more than 2.5 Å and a large number of water molecules remain in the surrounding of this domain area. This is a strong evidence for the mechanochemical stability of a receptor protein’s LBD towards hormone activated gene expressions followed by ligand binding and dissociation. 1. Introduction Steered molecular dynamics (SMD) is used to unfold the proteins and to study their elastic properties visualizing the different unfolding pathways [1]. The biophysical phenomena behind ligand binding, dissociation and conformational changes of thyroid hormone receptors (THR) are important in triiodothyronine (T3) stimulated gene expressions [2]. The physical properties such as echo dephasing, heat capacity, thermal diffusivity and thermal conductivity of liganded and/or unliganded THR-subtypes in folding states are previously studied [3, 4] by using MD simulations. This work is licensed under the Creative Commons CC BY-NC License. https://creativecommons.org/licenses/by-nc/4.0/ DOI: https://doi.org/10.3126/bibechana.v17i0.25870 http://nepjol.info/index.php/BIBECHANA mailto:tikaramlamichh@gmail.com https://creativecommons.org/licenses/by-nc/4.0/ https://doi.org/10.3126/bibechana.v17i0.25870 Tika Ram Lamichhane, Hari Prasad Lamichhane / BIBECHANA 17 (2020) 50-57 51 A point mutation in THR- gene causes resistance to thyroid hormones. The mutational impacts are observed distinctly on the protein-hormone systems by analyzing conformations and interaction energies through molecular dynamics approach [5]. SMD is a technique to know the structure-function relationships of the protein-hormone complex through unbinding of hormone or unfolding of protein under the application of time-dependent external forces. The elastic properties of the biomolecular systems subjected to deformations by SMD are in close agreement with the experimental results obtained from atomic force microscopy (AFM) and optical tweezers [6-8]. The mechanical stability of thyroid hormone like heavy ligand receptors is governed by protein’s secondary structure and pulling geometry [9]. Thyroid hormone dissociation or unfolding of THR strands proceeds a frictional path with constant velocity along x-direction defined by Langevin’s equation [1] 𝜇�̇� = − 𝑑𝑈 𝑑𝑥 + 𝐹(𝑥, 𝑡) + 𝜎𝑓(𝑡) (1) where  is time dependent frictional coefficient having dimension of [MT -1 ], F(x, t) is deforming force, 𝑈(𝑥) is potential governing ligand dissociation or protein unfolding pathways and 𝑓(𝑡) is the stochastic or fluctuating force term having coupling coefficient 𝜎. In SMD simulation, the SMD atom is attached to a dummy atom through a virtual spring. In one dimensional pulling, the dummy atom moves with constant velocity (�⃗� = 𝑑�⃗�/𝑑𝑡) so that the SMD atom experiences the force vector �⃗�(𝑥, 𝑡) = 𝑘(�⃗�𝑡 − ∆�⃗�) depending on the linear distance between these atoms [10]. So, the external potential energy [11] is given by 𝑈(𝑥, 𝑡) = 1 2 𝑘[(�⃗�𝑡 − ∆�⃗�). �⃗⃗�]2 (2) where 𝑘 is spring constant that specifies the stiffness of the applied harmonic restraining force, ∆�⃗�(𝑡) = �⃗�(𝑡) − �⃗�0 is tagged group displacement with �⃗�(𝑡) and �⃗�0 being actual and initial positions of the SMD atom and �⃗⃗� is the direction of pulling. 2. Methodology The initial structure of T3-liganded THR- isoform of nuclear receptor super family was taken from the protein data bank code 3GWS [12]. The THR- ligand binding domain (LBD) complex has the chain length of -helices and -forms with 259 amino acids and 3895 atoms. In the folding state or the globular form, end to end distance of THR- LBD is 38 Å. The THR- LBD actively binds T3 hormones having 35 atoms including 3 iodine atoms. The simulation packages such as protein structure file (psf) generation and solvation with water (TIP3P model) and ions providing cellular environment were prepared and structural and graphical analysis were performed by using visual molecular dynamics (VMD-1.9.3) [13]. In accordance with nanoscale molecular dynamics (NAMD-2.12) protocols [14], the topologies and parameters required for MD simulations of the THR- LBD complex were obtained from CHARMM force fields for proteins [15, 16]. The T3-hormone was parameterized with the help of Zoete’s force field generation tool [17]. The T3-liganded THR- LBD (THRT3) was fully solvated into a water droplet of radius 37.5 Å consisting of 17245 water molecules neutralized with 26 Na + and 16 Cl - ions in the concentrations of 0.15 mol/L. The system’s energy was minimized up to 3000 conjugate gradient steps and it was equilibrated up to 20 ns with NAMD protocols. Velocity Verlet algorithm [18] was used for the equilibration simulations with the integrator parameter of 2 fs/step, Langevin thermostat at 310 K and barostat at 1-atm and damping coefficient of 1 ps -1 . For the Lenard-Jones interactions, a 12 Å cut-off with smooth switching function starting at 10 Å was applied with 1-4 scaling 1.0. The final coordinates of the solvated THRT3 were extracted from the equilibrated droplet in order to perform SMD simulations. The SMD simulation was conducted setting the C atom of the last residue-460 as the SMD atom and the C atom of the first residue-202 as the fixed Tika Ram Lamichhane, Hari Prasad Lamichhane / BIBECHANA 17 (2020) 50-57 52 atom. The virtual spring between the SMD atom and the dummy atom was 7 kcal/mol/ Å 2 . The one dimensional pulling was performed at constant velocity of 0.1 Å/ps in the direction along a vector connecting the fixed atom and the SMD atom. The THRT3 complex became completely unfolded at the simulation time of about 8.5 ns and the SMD was conducted up to 10 ns. The intermediate conformational states of the THRT3 structures were visualized and the images were generated by using the VMD software. The related physical parameters such as radius of gyration (RG), root mean square deviation (RMSD), extension or end- to-end distance (x), forces and energies were noted down and the graphical analysis was performed with the plotting program XMGRACE. The SMD simulation was repeated three times to check the accuracy of the obtained results. 3. Results and Discussion In the folding state or completely stable globular form, N-C termini or end to end length of the THRT3 system is 38 Å. The receptor protein is unfolded smoothly with simulation time by breaking H-bonds among -helices or -sheets during the constant velocity (0.1 Å/ps) pulling of the SMD atom. The changing extension and RG of THRT3 over the course of simulation are shown in Figures 1-a & 1-b, respectively. The system becomes completely unfolded after the simulation time of 8.5 ns. At t = 8.5 ns, the end-to-end length is 876 Å and RG of about 258 Å. The THRT3 conformations responsible for the peak force are B, C and D as indicated in force vs extension graph (Figure 2). The structure A is the initial folding state in equilibrated form whereas E is the final unfolding state extended to 1023 Å. The force at the end (t > 8.5 ns) of the simulation increases when THRT3 becomes completely unfolded due to stretching of the protein single strand. The structures (A, B, C, D and E) formed by sequential unzipping of H-bonds are shown in Figure 3 and the related physical parameters such as time, length, RG, RMSD, number of H-bonds (at 3 Å internal distance and 20 o angle), force, kinetic energy (KE), potential energy (PE) and electrostatic energy are reported in Table 1. The resistance of H- bonds ruptures simultaneously causing the structural change and rapid extension of THRT3 LBD. Unfolding of native structures needs the maximum pulling force as they represent the bottom of the steep free energy well [19]. (a) (b) Fig. 1: End-to-end distance or extension and radius of gyration (RG) of THRT3 over the course of SMD simulation. The T3 binding domain formed by -helices and - hairpins has been verified to be the most stable region, i.e. the hormone binds strongly in the LBD of THR- because this region surrounded by a large mass of water remains almost unchanged (Figure 3-C) even up to 533 Å extension and 5 ns SMD simulation. Even after the T3 binding pocket is completely unfolded and the protein RMSD is raised up to 300 Å (Figure 4-a), the position of T3 does not shift more from its actual position as indicated by its RMSD (< 2.5 Å shown in Figure 4- Tika Ram Lamichhane, Hari Prasad Lamichhane / BIBECHANA 17 (2020) 50-57 53 b). Abrupt breaking of hydrophobic contacts between helices 8 and 12 is required for the dissociation of T3 hormone from THR- LBD. The ligand binding and dissociation pathways observed in THR-isoforms by using SMD simulations have been explained in the previous studies [10, 20, 21]. RMSD of T3 gets small step-up jump, i.e. the ligand/hormone becomes slightly unstable in the unfolding states of THRT3 responsible for the peak force generation as shown in the Figure 4-b. Along with the elongating system, H-bonds decrease in number (Figure 5-a) whereas electrostatic energy increases up to the complete unfolding state of THRT3 (Figure 5-b). At the time of complete unfolding (t  8.5 ns), H-bonds reduce to 10 and the ranges of KE, PE and electrostatic energy are 13345, -5900 and -69000 kcal/mol, respectively. The actual data are reported in Table 1. After the complete unfolding state (t > 8.5 ns), H-bonds in a small number and electrostatic energy both remain almost constant as shown in Figure 5. Fig. 2: Force vs extension plot showing initial folding state (A), intermediate states (B, C) having greater force associated with the breaking of H-bonds among -helices or -sheets and completely unfolding states (D, E). The force at the end (t > 8.5 ns) increases when THRT3 becomes completely unfolded. Table 1: Physical parameters at different unfolding states of THRT3 protein responsible for peak force Unfolding states Time (ps) End-to-end distance (Å) Force (kcal/mol-Å) No. of H- bonds RG (Å) RMSD (Å) KE (kcal/mol) PE (kcal/mol) Elect. energy (kcal/mol) A 0 38.00 0.00 60 18.86 0.61 13809.01 -65680.75 -75115.80 B 1673 203.64 15.40 54 34.06 23.75 13832.59 -64576.79 -74093.63 C 4964 532.87 17.17 31 125.17 121.28 13446.59 -61668.09 -71620.80 D 8404 875.90 17.10 10 257.22 254.02 13343.50 -58944.43 -69681.59 E 10000 1023.29 98.21 9 300.66 297.56 13273.46 -50967.98 -69695.32 C D Tika Ram Lamichhane, Hari Prasad Lamichhane / BIBECHANA 17 (2020) 50-57 54 (A) t = 0 ps, x = 38.00 Å (C) t = 4964 ps, x = 532.87 Å (D) t = 8404 ps, x = 875.90 Å (E) t =10000ps, x = 1023.29 Å (B) t = 1673 ps, x = 203.64 Å Fig. 3: Different unfolding states responsible for peak force provided with the simulation time (t) and end-to-end distance (x) during the constant velocity (0.1 Å/ps) stretching of THRT3. Fig. 4: Root mean square deviation (RMSD) of (a) receptor protein and (b) T3-hormone during unfolding of THRT3. Tika Ram Lamichhane, Hari Prasad Lamichhane / BIBECHANA 17 (2020) 50-57 55 Fig. 5: Variation of (a) H-bonds and (b) electrostatic energy over the course of 10 ns SMD simulation at constant velocity (0.1 Å/ps) pulling of THRT3. Fig. 6: Energy profile diagram during unfolding of THRT3 by SMD at constant velocity (0.1 Å/ps ) pulling. (a) (b) Fig. 7: Variation of (a) kinetic and (b) internal potential energies over the course of 10 ns SMD simulation at constant velocity (0.1 Å/ps) pulling of energy. Tika Ram Lamichhane, Hari Prasad Lamichhane / BIBECHANA 17 (2020) 50-57 56 The Linard Jones potential is getting raised linearly with small slope and the energy terms: dihedral and improper remain almost constant throughout the SMD simulation. Bond-angle and bond-length energies are almost constant before 8.5 ns and they rise up steeply after 8.5 ns when THRT3 is completely unfolded as shown in Figure 6. Figures 7-a & 7-b show that KE decreases slowly and linearly throughout the simulation, and conversely, PE increases slightly up to the state of complete unfolding of THRT3. When the system is completely unfolded (t > 8.5 ns), PE increases steeply. The cause behind the steeply changed PE is the fast increasing force generated after unwinding the molecular system. The load bearing strands are shielded by water and it interacts with bond-breaking events between such strands [9]. Even if the strand is fully dissociated, water makes H-bonding to the exposed sites as in Figure 3. The hydrophobic LBD of THR- with load bearing H-bonds is protected from water attack so that T3 is not dragged away from its pocket position even if the system is completely unfolded. In order to employ the reversible work comparable to AFM experiments, very low speed pulling is to be implemented while unwinding the molecular system which is practically difficult due to computational limits. However, higher pulling velocities do not influence the reliability of the SMD results [22]. Continuous breaking of H-bonds up to the complete unfolding of THRT3 (Figure 5-a) and force vs extension graph (Figure 2) are comparable to AFM results even at this pulling speed of 0.1 Å/ps. The SMD technique of crack propagation allows the identification of intermediate conformations of THRT3 responsible for the gene transcriptional activities. 4. Conclusion Triiodothyronine nuclear receptor (THRT3) is completely unfolded resulting end-to-end length of 876 Å by 8.5 ns long SMD simulations performed at constant velocity of 0.1 Å/ps. The peak force associated with the intermediate conformations is due to breaking of H-bonds among -helices and - hairpins. Though RMSD of T3 does not exceed 2.5 Å throughout the simulation, T3 is deviated more from its mean position in the conformational states of peak force generation. Hydrophobic LBD, i.e. T3 binding pocket of the receptor that gets surrounded by water is unfolded only in the last of the simulation. It is an evidence for the stability of THR- LBD towards ligand (T3) binding and dissociation. Even after complete stretching of the protein strand LBD, the position of T3 remains almost constant surrounded by water molecules. Along with decrease in H-bonds, electrostatic energy increases gently during unfolding. There are linear changes with slightly decreased kinetic energy and slightly increased van der Waals energy at this constant velocity pulling. However, force or net potential of the system increases rapidly at the end (t > 8.5 ns) of the simulation. Dihedral and improper energies remain almost unchanged throughout the SMD simulations, but bond-angle and bond-length energies increase steeply after complete unfolding (t > 8.5 ns) of the system. Thus, one dimensional mechanical pulling at constant velocity is important technique to better understand the unfolding pathways of THRT3 like nuclear receptors. Acknowledgements We would like to acknowledge the computing resources provided by Prof. Dr. Raju khanal at his Plasma Lab, Central Department of Physics, Tribhuvan University, Kathmandu, Nepal. The partial financial support for this research has been provided by Nepal Academy of Science and Technology as the PhD fellowship to the first author. References [1] S. Izrailev, S. Stepaniants, B. Isralewitz, D. Kosztin, H. Lu, F. Molnar, W. Wriggers, K. Schulten, Steered molecular dynamics, Computational molecular dynamics: challenges, methods, ideas, (1999) 39-65, Springer, Berlin, Heidelberg. [2] T. M. Ortiga-Carvalho, A. R. Sidhaye, F. E. Wondisford, Thyroid hormone receptors and resistance to thyroid hormone disorders, Nature Rev. Endocrinol 10 (2014) 582. doi.org/10.1038/nrendo.2014.143 [3] T. R. Lamichhane, S. Paudel, B. K. Yadav, H. P. Lamichhane, Echo dephasing and heat capacity https://dx.doi.org/10.1038%2Fnrendo.2014.143 Tika Ram Lamichhane, Hari Prasad Lamichhane / BIBECHANA 17 (2020) 50-57 57 from constrained and unconstrained dynamics of triiodothyronine nuclear receptor protein, J. Biol. Phys. 45(2019) 107-135. doi.org/10.1007/s10867-018-9518-3 [4] T. R. Lamichhane, H. P. Lamichhane, Heat conduction by thyroid hormone receptors, AIMS Biophys. 5(2018) 245-256. doi.org/10.3934/biophy.2018.4.245 [5] T. R. Lamichhane, S. Paudel, B. K. Yadav, H. P. Lamichhane, Molecular dynamics approach to the I431V mutational impact on thyroid hormone recep tor -beta, BIBECHANA, 16 (2019) 79-91. doi.org/10.3126/bibechana.v16i0.21109 [6] J. J. Booth, D. V. Shalashilin, Fully atomistic simulations of protein unfolding in low speed atomic force microscope and force clamp experiments with the help of boxed molecular dynamics, J. Phys. Chem. B 120 (2016) 700-708. doi.org/10.1021/acs.jpcb.5b11519 [7] Z. Mártonfalvi, P. Bianco, K. Naftz, G. G. Ferenczy, M. Kellermayer, Force generation by titin folding, Protein Sci 26 (2017) 1380-1390. https://doi.org/10.1002/pro.3117 [8] J. Schönfelder, D. De Sancho, R. Perez-Jimenez, The power of force: insights into the protein folding process using single-molecule fo rce spec troscopy, J. Mol. Biol. 428 (2016) 4245- 4257. doi.org/10.1016/j.jmb.2016.09.006 [9] D. L. Guzmán, A. Randall, P. Baldi, P., Z. Guan, Computational and single-molecule force studies of a macro domain protein reveal a key molecular determinant for mechanical stability. Proc. Natl. Acad. Sci. U. S. A. 107(2010) 1989-1994. doi.org/10.1073/pnas.0905796107 [10] L. Martínez, I. Polikarpov, M. S. Skaf, Only subtle protein conformational adaptations are required for ligand binding to thyroid hormone receptors: simulations using a novel multipoint steered molecular dynamics approach, J. Phys. Chem. B 112(2008) 10741-10751. doi.org/10.1021/jp803403c [11]J. R. Gullingsrud, R. Braun, K. Schulten, Reconstructing potentials of mean force through time series analysis of steered molecular dynamics simulations, J. Comp. Phys. 151(1999) 190-211. doi.org/10.1006/jcph.1999.6218 [12] A.S. Nascimento, S. M. G. Dias, F. M. Nunes, R. Aparício, A. L. Ambrosio, Bleicher, L., A. C. M. Figueira, M. A. M. Santos, M. de Oliveira Neto, H. Fischer, M. Togashi, Structural rearrangements in the thyroid hormone receptor hinge domain and their putative role in the receptor function, J. Mol. Biol. 360 (2006) 586-598. doi.org/10.1016/j.jmb.2006.05.008 [13] W. Humphrey, A. Dalke, K. Schulten, VMD— Visual Molecular Dynamics, J. Mol. Graphics 14 (1996) 33–38. doi.org/10.1016/0263-7855(96)00018-5 [14]W. Humphrey, A. Dalke, K. Schulten, VMD— Visual Molecular Dynamics, J. Mol. Graphics 14 (1996) 33–38. doi.org/10.1016/0263-7855(96)00018-5 [15] A. D. MacKerell, M. Feig, C. L. Brooks, Extending the treatment of backbone energetics in protein force fields: Limitations of gas‐phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations, J. Comp. Chem. 25(2004) 1400-1415. doi.org/10.1002/jcc.20065 [16] A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, All-atom empirical potential for molecular modeling and dynamics studies of proteins, J. Phys. Chem. B 102 (1998) 3586–3616. doi.org/10.1021/jp973084f [17]V. Zoete, M. A. Cuendet, A. Grosdidier, O. Michielin, SwissParam: a fast force field generation tool for small organic molecules, J. Comp. Chem. 32 (2011) 2359-2368. doi.org/ 10.1002/jcc.21816 [18] L. Verlet, Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard- Jones molecules, Phy. Rev. 159 (1967) 98. doi.org/10.1103/PhysRev.159.98 [19] J. J. Booth, D. V. Shalashilin, Fully atomistic simulations of protein unfolding in low speed atomic force microscope and force clamp experiments with the help of boxed molecular dynamics, J. Phys. Chem. B 120(2016) 700-708. doi.org/10.1021/acs.jpcb.5b11519 [20] D. Kosztin, S. Izrailev, K. Schulten, Unbinding of retinoic acid from its receptor studied by steered molecular dynamics, Biophys. J. 76(1999) 188-197. doi.org/10.1016/S0006-3495(99)77188-2 [21] L. Martínez, P. Webb, I. Polikarpov, M. S. Skaf, Molecular dynamics simulations of ligand dissociation from thyroid hormone receptors: evidence of the likeliest escape pathway and its implications for the design of novel ligands, J. Med. Chem. 49(2006) 23-26. doi.org/10.1021/jm050805n [22] J. L. Zhang, Q. C. Zheng, H. X. Zhang, Unbinding of glucose from human pulmonary surfactant protein D studied by steered molecular dynamics simulations, Chem. Phys. Lett. 484(2010) 338-343. doi.org/10.1016/j.cplett.2009.12.022 https://doi.org/10.1007/s10867-018-9518-3 https://doi.org/10.3934/biophy.2018.4.245 https://doi.org/10.3126/bibechana.v16i0.21109 https://doi.org/10.1021/acs.jpcb.5b11519 https://doi.org/10.1002/pro.3117 https://doi.org/10.1016/j.jmb.2016.09.006 https://doi.org/10.1073/pnas.0905796107 https://doi.org/10.1021/jp803403c https://doi.org/10.1006/jcph.1999.6218 https://doi.org/10.1021/acs.jpcb.5b11519 https://doi.org/10.1016/S0006-3495(99)77188-2 https://doi.org/10.1021/jm050805n