Effects of Deformation Rate on the Unbinding Pathway of the MMP8-Aggrecan_IGD Complex in Cartilage

Deng Li and Shuwei Chang,

1 Department of Civil Engineering,National Taiwan University,Taipei 10617,Taiwan.

Abstract:Mechanical force plays a critical role in the remodeling and degradation of cartilage tissues.The cartilage tissue generates,absorbs,and transmits mechanical force,enabling specific biological processes in our body.A moderate intensity mechanical force is necessary for cartilage tissue remodeling and the adaptation of biomechanical properties,but a high intensity mechanical force can lead to pathological degradation of cartilage tissue.However,the molecular mechanism of cartilage degradation is still unclear.We use full atomistic simulations with SMD simulations to investigate whether the magnitude of mechanical force affects the unbinding pathway of the MMP8-Aggrecan_IGD complex.We find that when the pulling velocity is slow,the mechanical force required to unbind the Aggrecan_IGD from MMP8 is higher,and a three-step unbinding pathway is observed.On the other hand,when the pulling velocity is fast,the mechanical force required to unbind the Aggrecan_IGD from MMP8 is lower,and a two-step unbinding pathway is observed.Our results help us to understand how the magnitude of the mechanical force affects the unbinding pathway of the enzyme-ligand complex in cartilage tissue at the molecular level.

Keywords:Mechanical force,steered molecular dynamics (SMD),unbinding pathway,potential of mean force (PMF).

1 Introduction

Cartilage tissue is an important structural component in our body,providing biomechanical properties for joints [Fosang,Rogerson,East et al.(2008);Pratta,Yao,Decicco et al.(2003)].Mechanical force plays a critical role in the remodeling and degradation of cartilage tissues [Grodzinsky,Levenston,Jin et al.(2000)].The cartilage tissue generates,absorbs and transmits mechanical force,enabling specific biological process in our body [Julkunen,Harjula,Iivarinen et al.(2009);Borrelli,Torzilli,Grigiene et al.(1997)].These biological processes are strongly influenced by biomechanical signals mediated by mechanical forces [Pufe,Lemke,Kurz et al.(2004)].A moderate intensity mechanical force is necessary for cartilage tissue remodeling and the adaptation of biomechanical properties,but a high intensity mechanical force can lead to pathological degradation of cartilage tissue and cause related diseases such as developmental dysplasia of the hip (DDH) and osteoarthritis (OA) [Ewers,Jayaraman,Banglmaier et al.(2003);O’Kane,Hutchinson,Atley et al.(2006);Bassey,Rothwell,Littlewood et al.(1998);Ewers,Dvoracek-Driksna,Orth et al.(2001);Muller and Seddon (1953)].

Aggrecan is one of the dominant extracellular matrix (ECM) constituents in cartilage tissues providing relatively low friction for smooth articular movement [Venn and Maroudas (1977)].Aggrecan degradation in cartilage tissue occurs predominantly through proteolysis and has been attributed to the action of ADAMTS (a disintegrin and metalloproteinase with thrombospondin motifs) and MMPs (matrix metalloproteinases) [Durigova,Nagase,Mort et al.(2011);Malemud (2017)].Many diseases are associated with pathological degradation of aggrecan caused by high intensity mechanical force [Buckwalter (2012);Felson (2013);Huang and Wu (2008)].During the degradation of aggrecan,the enzymes bind to specific sites in the interglobular domain (IGD) of the aggrecan core protein and form a protein-ligand complex for further hydrolyzing [Huang and Wu (2008);Durigova,Nagase,Mort et al.(2011)].However,it remains unclear on how the mechanical forces may affect the molecular structures of the enzyme and the ECM in cartilage.

Full atomistic simulations with steered molecular dynamics (SMD) simulations have been successfully used to investigate the unbinding pathway of many protein-ligand complexes (Tab.1) [Dellago and Hummer (2013);Ferreira,Franca and Leite (2017);Patel,Berteotti,Ronsisvalle et al.(2014);Shen,Guan,Xu et al.(2012);Tekpinar (2018);Vashisth and Abrams (2008);Zhang,Santos,Zhou et al.(2016)].In the SMD simulations,an external force or a constant velocity is applied to the system,and the responses of the system are recorded and analyzed.SMD simulations allow us to apply mechanical forces on the ligands to investigate the unbinding pathway and the conformation changes of the proteinligand complex from its equilibrium state to accelerate the transitions between different local energy minima [Aprodu,Redaelli and Soncini (2008)].

The objective of this study is to understand the effects of mechanical force on the unbinding pathway of the MMP8-Aggrecan_IGD complex.We combine the peptide of a representative cleavage site in the aggrecan IGD with MMP8 to form MMP8-Aggrecan_IGD complex [Humphrey,Dalke and Schulten (1996)].Full atomistic simulations with SMD simulations are performed to investigate whether the magnitude of mechanical force affects the unbinding pathway of the MMP8-Aggrecan_IGD complex.Simulations with different SMD parameters,including different force constants and pulling velocities,are performed to explore the effects of different mechanical force and potential of mean force (PMF) is calculated to identify the difference of energy profile.

2 Computational details

The real sequences of aggrecan core protein in the interglobular domain of human are used to generate the aggrecan core protein peptide.The sequences of the aggrecan core protein in the interglobular domain (Aggrecan_IGD) are obtained from the UniProt protein database (http://www.uniprot.org/uniprot/P16112) [Renaux and Consortium (2018)].The chosen sequences are PDMELPLPRNITEGE～ARGSVILTVKPIFEV.The scissile bond is depicted with a～sign.The homology model of Aggrecan_IGD is constructed using Modeller 9.19 [Fiser and Sali (2003)].NAMD 2.12 with CHARMM27 force field is used for all the simulations [Phillips,Braun,Wang et al.(2005);MacKerell,Bashford,Bellott et al.(1998);Jo,Cheng,Islam et al.(2014)].The homology model is energy minimized with conjugate gradient scheme,and the initial structure of the MMP8 catalytic domain is adapted from the protein data bank of PDB ID:2OY4 [Bertini,Calderone,Fragai et al.(2006)].The MMP8 and Aggrecan_IGD are combined to form MMP8-Aggrecan_IGD complex by using Visual Molecular Dynamics program [Humphrey,Dalke and Schulten (1996)].The MMP8-Aggrecan_IGD complex is solvated by using the TIP3P water molecules.Ions are added to the model to neutralize the simulation system (Fig.1).After the model construction,the complex is first energy minimized with conjugate gradient algorithm,followed by NPT equilibration for 80 ns at a constant temperature of 310 K and with 1 bar pressure.Simulation trajectories from 50 ns to 70 ns are used to analyze.

Figure 1:MMP8-Aggrecan_IGD complex solvated in water box which rendered in molecular surface renderer called MSMS [Sanner,Olson and Spehner (1996)].The complex used for the SMD simulation consists of MMP8 which rendered in NewCartoon and Aggrecan_IGD ligand which rendered in CPK.Chloride and sodium ions are shown in light blue and yellow.The constant velocity is applied at the Cα atom of the first residue of the Aggrecan_IGD ligand

To study the effects of different mechanical forces,we perform SMD simulations with different pulling velocities.In the constant velocity SMD,pulling velocity has a significant influence on the magnitude of the unbinding force.In the previous studies,there have been various velocities adapted in different SMD studies (Tab.1).In order to mimic the atomic force microscopy experiment,the pulling velocity used in the SMD simulation should be as slow as possible.On the other hand,the force constant should be large enough to be able to measure the local dissociation potential;however,it should not be too large as the fluctuation of the force-displacement curve increases as the force constant increases.The choice of pulling velocities and force constant are based on a balance between computer resource and accuracy.We chose two different pulling velocities:0.01 ?/ps and 0.05 ?/ps,both with a force constant of 7 kcal/mol/?2,to investigate the effects of the magnitude of the mechanical forces.The MMP8 is fixed during the SMD simulations.To obtain reliable PMF profiles of the Aggrecan_IGD dissociation from MMP8,a number of additional simulations are conducted by repeating the SMD trajectories.Ten SMD trajectories are generated for each PMF calculation.Based on the protocol of Christoph Dellago et al.[Dellago and Hummer (2013)],we use second order cumulant expansion of Jarzynski’s equality to calculate the PMF [Jarzynski (1997);Crooks (2000);West,Olmsted and Paci (2006);Dellago and Hummer (2013)].

Table 1:The choice of force constant and velocity in different SMD study

Table 1:Continued

3 Result and discussion

3.1 Molecular structure of the MMP8-Aggrecan_IGD complex

To analyze the structural stability of the MMP8-Aggrecan_IGD complex,the root-meansquare deviations (RMSD) are calculated.The RMSD from the initial equilibrium state of the MMP8,Aggrecan_IGD,and MMP8-Aggrecan_IGD complex are shown as a function of simulation time in Fig.2.The RMSD profiles show that the structural conformation of the complex reaches a stable conformation at 25 ns.The average RMSD of the MMP8-Aggrecan_IGD complex is 2.61 ?,with a standard deviation of 0.26 ?.The average RMSD of the MMP8 is 2.44 ?,with a standard deviation of 0.25 ?.The differences of RMSD between the MMP8-Aggrecan_IGD complex and the MMP8 is only 0.17 ? suggesting that the binding of Aggrecan_IGD does not alter the molecular structure of MMP8 significantly.

Figure 2:Root-mean-square deviations (RMSD) with respect to the initial equilibrium state for MMP8-Aggrecan_IGD complex

We analyze the hydrogen bonds of MMP8-Aggrecan_IGD complex with a distance cutoff of 3.5 ? and an angle cut-off of 30°.Aggrecan_IGD bound MMP8 mainly through five hydrogen bonds:Tyr219-O…HN-Ile21,Glu180-O…HN-Arg17,Gly158-O…HN-Arg17,Leu160-NH…O-Glu15and Ala161-O…HN-Glu15.The number of hydrogen bonds between the MMP8 and the Aggrecan_IGD is plotted as a function of simulation time in Fig.3.It is shown that the number of hydrogen bonds has a large fluctuation in the first 25 ns simulations.After 25 ns,the number of hydrogen bonds approaches a stable value of 4.7815.The results of the RMSD and the number of hydrogen bonds confirms that our simulation reaches a stable structure for SMD simulations.

Figure 3:The number of hydrogen bonds (blue line) between the MMP8 and the Aggrecan_IGD as a function of the simulation time

3.2 Effects of pulling velocity on the unbinding pathway

SMD simulations are performed to investigate the unbinding pathway to provide molecular insights into the interactions between the MMP8 and the Aggrecan_IGD.Force-reaction coordinate curves and the number of hydrogen bonds are analyzed to investigate the unbinding pathway of MMP8-Aggrecan_IGD complex at different pulling velocity [Dellago and Hummer (2013);Ferreira,Franca and Leite (2017);Patel,Berteotti,Ronsisvalle et al.(2014);Shen,Guan,Xu et al.(2012);Tekpinar (2018);Vashisth and Abrams (2008);Zhang,Santos,Zhou et al.(2016)].

The force-reaction coordinate curve of the MMP8-Aggrecan_IGD complex at a pulling velocity of 0.01 ?/ps is shown in Fig.4.There are three peaks on the force-reaction coordinate curve during the SMD simulation.The first peak occurs at a reaction coordinate of 24 ?.The second peak is at the reaction coordinate of 35 ?.The maximum peak appears at 46 ?.At about 55 ?,the Aggrecan_IGD is completely pulled out from the MMP8,and thus the pulling forces go to zero.

Fig.5 shows the unbinding pathway of the MMP8-Aggrecan_IGD complex at the pulling velocity of 0.01 ?/ps.Fig.5(A) shows the structure of MMP8-Aggrecan_IGD complex before applying mechanical force.When pulling Aggrecan_IGD out of the MMP8 catalytic domain,the Aggrecan_IGD first lost two hydrogen bonds:Ala161-O…HN-Glu15and Leu160-NH…O-Glu15,as shown in Fig.5(B) at the reaction coordinate of 24 ?,which is corresponding to the first peak in Fig.4.The hydrogen bond Gly158-O…HN-Arg17is broken (Fig.5(C)) when the force reaches the second peak in Fig.3.When the reaction coordinate reaches 41 ?,Aggrecan_IGD formed a new hydrogen bond with MMP8 (Gly158-NH…OAla16) (Fig.5(D)).At the reaction coordinate of 44 ?,the newly formed hydrogen bond is broken (Fig.5(E)).Then the hydrogen bond (Glu180-O…HN-Arg17) breaks right after the maximum force.Finally,the last relatively weaker hydrogen bond (occupancy less than 50%) is immediately broken (Tyr219-O…HN-Ile21) (Fig.5(F)).These results reveal a threestep unbinding pathway for the MMP8-Aggrecan_IGD complex.

The force-reaction coordinate curve of the MMP8-Aggrecan_IGD complex at a pulling velocity of 0.05 ?/ps is shown in Fig.6.There are totally two peak forces during the SMD simulation.The first peak,which is also the maximum peak occurs at a reaction coordinate of 30 ?.The second peak appears at a reaction coordinate of 46 ?.At about 55 ?,the pulling force goes to zero.

Fig.7 reports the unbinding pathway of the MMP8-Aggrecan_IGD complex at the pulling velocity of 0.05 ?/ps.The initial structure of MMP8-Aggrecan_IGD complex is showing in Fig.7(A).Aggrecan_IGD loses two hydrogen bonds (Ala161-O…HN-Glu15and Leu160-NH…O-Glu15) (seen Fig.7(B)) at the reaction coordinate of 30 ? where the maximum force appears.When the reaction coordinate reaches 46 ?,the other two hydrogen bonds (Glu180-O…HN-Arg17and Gly158-O…HN-Arg17) are broke (Fig.7(C)) corresponding to the second peak in Fig.6.The last hydrogen bond is also broken immediately (Fig.7(D)).These results reveal a two-step unbinding pathway for the MMP8-Aggrecan_IGD complex at a higher pulling velocity.

In summary,the unbinding pathway of the MMP8-Aggrecan_IGD complex is altered by the deformation rate.When the pulling velocity is lower,the mechanical force required to unbind the Aggrecan_IGD from MMP8 is higher,and a three-step unbinding pathway is observed.On the other hand,when the pulling velocity is higher,the mechanical force required to unbind the Aggrecan_IGD from MMP8 is lower,and a two-step unbinding pathway is observed.The specific transition state (Fig.5(D)) in the unbinding pathway is observed at a lower pulling velocity.

Figure 5:Unbinding pathway of the MMP8-Aggrecan_IGD complex at the pulling velocity of 0.01 ? /ps.MMP8 is rendered in white NewCartoon,and Aggrecan_IGD is rendered in licorice.Dashed lines mean the hydrogen bonds between MMP8 and Aggrecan_IGD [Humphrey,Dalke and Schulten (1996)]

Figure 6:Force-reaction coordinate curve during the SMD simulations at the pulling velocity of 0.05 ?/ps.The smoothed force-reaction coordinate curve at the pulling velocity of 0.01 ?/ps is represented in red color

Figure 7:Unbinding pathway of MMP8-Aggrecan_IGD complex at the pulling velocity of 0.05 ?/ps.MMP8 is rendered in white NewCartoon,and Aggrecan_IGD is rendered in licorice.Dashed lines mean hydrogen bonds between MMP8 and Aggrecan_IGD [Humphrey,Dalke and Schulten (1996)]

3.3 Effects of different force constant on the potential of mean force (PMF)

As the discussion above,the unbinding pathway of MMP8-Aggrecan_IGD complex is controlled by the deformation rate.The transition state (Fig.5(D)) is only observed at lower pulling velocity.To further study the influence of the force constants as Aggrecan_IGD unbinding MMP8 at lower pulling velocity (0.01 ?/ps),PMF is calculated to estimate the free energy during the unbinding pathway with different force constants (K=0.597 kcal/mol/?2,K=7 kcal/mol/?2).The PMF in Fig.8 shows the energy change along the unbinding pathway of the MMP8-Aggrecan_IGD complex.When K=0.597 kcal/mol/?2,PMF increases slightly at the beginning mainly due to the hydrogen bonds of Ala161-O…HN-Glu15and Leu160-NH…O-Glu15until 35 ?.Then the PMF increases sharply until 55 ? due to the hydrogen bonds of Glu180-O…HN-Arg17and Gly158-O…HN-Arg17.Finally,the PMF is converged after 95 ? displacement.As for K=7 kcal/mol/?2,PMF increases slightly at the beginning also because of the hydrogen bonds of Ala161-O…HN-Glu15and Leu160-NH…O-Glu15until 31 ? displacement.Then the PMF increased sharply until 48 ? displacement also due to the hydrogen bonds of Glu180-O…HN-Arg17and Gly158-O…HNArg17.Finally,the PMF is converged after 96 ? displacement.Both PMF profiles have similar free energy change during the unbinding pathway of the MMP8-Aggrecan_IGD complex,suggesting that force constants do not affect the PMF significantly at a relatively lower pulling velocity.

Figure 8:Potential of Mean Force of the SMD simulations at different force constant

4 Conclusion

By using the full atomistic simulation approach,we investigate the effects of mechanical force on the unbinding pathway of the MMP8-Aggrecan_IGD complex from the atomistic scale.We find that the pulling velocity alters the unbinding pathway of MMP8-Aggrecan_IGD complex.When the pulling velocity is slow,the mechanical force required to unbind the Aggrecan_IGD from MMP8 is higher,and a three-step unbinding pathway is observed.On the other hand,when the pulling velocity is fast,the mechanical force required to unbind the Aggrecan_IGD from MMP8 is lower,and a two-step unbinding pathway is observed.An additional transition state with a new hydrogen bond forming between the MMP8 and the Aggrecan_IGD along the unbinding pathway is observed only at the lower pulling velocity.We also investigate the effects of the force constant on the PMF estimation along the unbinding pathway.We found that the PMF profiles obtained with a different force constant have similar free energy change along the unbinding pathway of MMP8-Aggrecan_IGD complex.This indicates that force constant has a relatively low influence for the PMF calculation.This study helps us to understand how the pulling velocity and force constant of SMD affects the unbinding pathway of the MMP8-Aggrecan_IGD complex at the molecular level.

Acknowledgments:The authors appreciate the financial support from the Ministry of Science and Technology (107-2218-E-002-064and 108-2636-E-002-016),Taiwan (R.O.C).We are grateful to the National Center for High-performance Computing for computer time and facilities.