p53 Complex Investigated by Parallel

Subscriber access provided by Iowa State University | Library

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Dissociation Process of MDM2/p53 Complex Investigated by Parallel Cascade Selection Molecular Dynamics and Markov State Model Duy Phuoc Tran, and Akio Kitao J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b10309 • Publication Date (Web): 15 Jan 2019 Downloaded from http://pubs.acs.org on January 15, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Dissociation Process of MDM2/p53 Complex Investigated by Parallel Cascade Selection Molecular Dynamics and Markov State Model Duy Phuoc Tran, Akio Kitao* School of Life Sciences and Technology, Tokyo Institute of Technology

*Corresponding author: Akio Kitao School of Life Science and Technology, Tokyo Institute of Technology, 2-12-1, Ookayama, Meguro-ku, Tokyo, 152-8550, Japan Tel: +81-3-5734-3373 Fax: +81-3-5734-3372 e-mail: [email protected]

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 30

Abstract Recently, we efficiently generated dissociation pathways of a protein-ligand complex without applying force-bias with parallel cascade selection molecular dynamics (PaCS-MD) and showed that PaCS-MD in combination with the Markov state model (MSM) yielded a binding free energy comparable to experimental values. In this work, we applied the same procedure to the complex of MDM2 protein and the transactivation domain of p53 protein (TAD-p53), the latter of which is known to be very flexible in the unbound state. Using 30 independent MD simulations in PaCS-MD, we successfully generated 25 dissociation pathways of the complex, which showed complete or partial unfolding of the helical region of TAD-p53 during the dissociation process within the average simulation time of 154.8 ± 46.4 ns. The standard binding free energy obtained by the combination with 1D-, 3D- or Cα-MSM was in good agreement with those determined experimentally. Using 3D-MSM based on the center of mass position of TAD-p53 relative to MDM2, the dissociation rate constant was calculated, which was comparable to those measured experimentally. Cα-MSM based on all Cα coordinates of TAD-p53 reproduced experimentally-measured standard binding free energy, and dissociation and association rate constants. We conclude that the combination of PaCS-MD and MSM offers an efficient computational procedure to calculate binding free energies and kinetic rates.


Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


I. Introduction The dissociation of protein-ligand complexes is an essential process in biological systems; however, observation of these processes is challenging. Molecular simulations are useful methods to investigate molecular mechanisms of protein dissociation in silico. For example, Paul et al. performed six independent Hamiltonian replica exchange molecular dynamics (MD) simulations of the complex of MDM2 and PMI peptide for 1 s and simulated dissociation events. By combining the results of conventional MD simulations of the bound state for 500 s in total, they constructed the Markov State Model (MSM) and reproduced residence time of the peptide in agreement with experimental values.1 Rydzeski et al. studied the dissociation of hyperzine A from Torpedo California acetylcholinesterase using 4 s unbiased allatom MD simulations and biased metadynamics with memetic sampling for determination of the reaction pathway2. They reproduced values of experimentally determined kinetic rates and explained the major dissociation pathways from a similar free energy cost based on the residence time. These examples demonstrated the power of molecular simulations to calculate kinetic rates, but also show the difficulty in sampling dissociation events by only using unbiased conventional MD. In addition to conventional MD, more advanced sampling techniques, e.g., steered molecular dynamics3, τ-random acceleration MD4, metadynamics and the extended versions5–7, string method8,9, mile-stoning10, adaptive biased force11,12, accelerated MD13,14, weighted ensemble15–17, WExplore18,19, and other methods not listed here have been used to investigate dissociation and association processes. Recently, we applied parallel cascade selection molecular dynamics (PaCS-MD)20,21 to simulate the dissociation process of the complex of hen egg white lysozyme (LYZ) and tri-N-acetyl-D-glucosamine (triNAG)22. PaCS-MD is an unbiased enhanced conformational sampling method, which comprises cycles of multiple unbiased MD simulations combined with a selection of MD snapshots as the initial structures for the next cycle. Since the probability of rare event occurrence toward dissociation was drastically enhanced by the selection and re-randomization of the initial velocities, PaCS-MD easily induced the dissociation of the LYZ/triNAG complex within the order of 100 – 101 ns and a total computational cost of 102 ns MD simulation time, despite the fact that no dissociation was observed during 1 𝜇𝑠 of conventional MD. Because one trial of PaCS-MD generates a dissociation pathway as a combination of multiple short MD trajectories that mutually overlap in conformational space, the generated trajectories can be utilized to construct a MSM. We showed that PaCS-MD/MSM yielded a binding free energy value comparable to the experimental values. 3 ACS Paragon Plus Environment


Page 4 of 30

PaCS-MD has been proven to be a powerful tool to estimate the binding free energy of a protein-ligand complex in combination with MSM22. Consequently, we were motivated to extend the method to estimate the binding free energy and dissociation rate constant of a more flexible protein-peptide complex. In this work, we investigate murine double minute clone 2 protein (MDM2) bound to the transactivation domain of p53 protein (TAD-p53)23. P53 is well-known for its tumor suppressor function, in which p53 responds to stress by blocking cell cycle progression or instigating programmed cell death, while overexpression of the MDM2 gene blocks the p53-mediated transactivation of a reporter gene bearing a p53-responsive element24. TAD-p53 has been proven to be a very flexible peptide that can adopt multiple conformations17. The formation of the MDM2/TAD-p53 complex performs the following three roles: 1) MDM2 ubiquitinates p53 upon binding; 2) MDM2 interacts with p53 to block its binding with target DNA; and 3) MDM2 promotes the export of p53 out of the cell nucleus, preventing it from interacting with target DNA25. We selected this complex for our study because of this biological importance as well as the very flexible nature of free TAD-p53. TAD-p53 is well-structured in the complex form but is expected to be less structured during the dissociation process. In addition, the standard binding free energy, and dissociation and association rate constants of this complex have been measured experimentally, which enables us to compare with calculated results. In this paper, we use PaCS-MD to generate dissociation pathways of the MDM2/TAD-p53 complex without applying force-bias within a short MD simulation time. By analyzing the generated trajectories by the following three methods: one-dimensional MSM based on the inter-center of mass (inter-COM) distance between the two molecules (1D-MSM); three-dimensional MSM based on the COM position of TAD-p53 relative to MDM2 (3D-MSM); and by Cα coordinates of TAD-p53 relative to MDM2 (CαMSM), the standard binding free energy is calculated, which is in good agreement with those determined experimentally. Using 3D-MSM and Cα-MSM, the dissociation rate constant was calculated, which is comparable to those measured experimentally. Cα-MSM reproduced the associate rate constant as well. We show that the combination of PaCS-MD and MSM (PaCS-MD/MSM) can reproduce binding free energy and kinetic rates which were measure experimentally.

2. Methods


Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.1. Parallel Cascade Selection Molecular Dynamics PaCS-MD generates conformational transition pathways using cycles of multiple independent MD simulations20. A simulated system for each MD is referred to as a replica hereafter. Typically, the number of replicas (nrep) is 10 – 100. We selected nrep = 30, which we confirmed to be the optimum number to construct MSM with a single trial of PaCS-MD in the case of the MDM2/TAD-p53 complex. After each cycle of PaCS-MD, the top nrep structures were selected from the trajectories of the previous cycle as the initial structures for the next cycle from the ranking of the generated snapshots based on longer interCOM distances between MDM2 and p53, d. Each MD trajectory (0.1 ns) contains 200 snapshots with 0.5 ps interval. Therefore, 30 initial structures for each cycle were selected from 6,000 snapshots. The new MD cycle was restarted after re-randomization of the initial velocities with the Maxwell-Boltzmann distribution. The previous study showed that the selection and re-randomization significantly accelerated the dissociation of the LYZ/triNAG complex22. The cycle of MD and selection were repeated until TADp53 completely dissociated from MDM2 up to d = 7.0 nm. We conducted 25 trials of PaCS-MD to generate multiple dissociation pathways.

2.2. Markov State Model and Kinetic Properties TAD-p53 dissociation pathways generated by PaCS-MD were analyzed by MSM. As shown in the previous paper22, one trial of PaCS-MD generates a set of unbiased trajectories that highly overlap in conformational space along the dissociation pathway because new MDs are started from snapshots of the previous cycle. Therefore, a set of the generated trajectories can be analyzed by MSM. In the case of the LYZ/triNAG dissociation22, we constructed one-dimensional MSM using the inter-COM distance d (1DMSM), which we also performed for the MSM2/TAD-p53 dissociation in this study. If a sufficient number of PaCS-MD trials are conducted, trajectories from different trials mutually overlap, which enables the construction of MSM from merged PaCS-MD trajectories of the multiple trials. In this case, we employed three-dimensional COM positions of TAD-p53 relative to MDM2 to construct MSM, which we call 3DMSM. As the third method to construct MSM, we considered all C coordinates of TAD-p53 relative to MDM2.



Page 6 of 30

To construct the MSM, it is necessary to cluster the MD snapshots into microstates via clustering methods. The transition probability matrix T is then estimated as described in the literature26. The stationary distribution of  can be easily obtained by solving the equation, 𝝅 = 𝝅𝑻.

(1)

Since each microstate represents a certain inter-COM distance 𝑑𝑖 in 1D-MSM (i: index of microstate), the potential of mean force (PMF) can be obtained as a natural logarithm of the stationary distribution as, 1

𝑊(𝑑𝑖) = ― 𝛽ln 𝜋(𝑑𝑖),

(2)

where β is thermodynamic beta. In the case of 3D-MSM, the 3D stationary probability 𝜋(𝒓𝑗) of state j is mapped into 1D-PMF by the following equation: 1

𝑊(𝑑𝑖) = ― 𝛽𝑙𝑛∑𝑑

𝑖

― 𝛿𝑑/2 ≤ 𝑑(𝒓𝑗) < 𝑑𝑖 + 𝛿𝑑/2

𝜋(𝒓𝑗)

(3)

where 𝒓𝑗 indicates the 3D position of microstate j and 𝛿𝑑 is a bin size (0.04 nm in this work). In the plots shown in the Results section, we set min(𝑊(𝑑𝑖)) = 0. The PMF difference between the bound and unbound states, ∆𝑊 = 𝑊(𝑑𝑢𝑛𝑏𝑜𝑢𝑛𝑑) ―𝑊(𝑑𝑏𝑜𝑢𝑛𝑑),

(4)

can be used to calculate the binding free energy of the simulated system. The calculation of Δ W is conducted based on the expectation that the PMF converges to a certain value if the distance between two molecules is far enough and inter-molecular interactions are negligible, independent of dissociation pathways. In Results section, we will show that ΔW converges very well independent of pathways. This type of approaches were used to further calculate binding free energy27–30. To obtain the standard binding free energy Δ Gº from the binding free energy calculated from PMF ∆𝐺𝑃𝑀𝐹, we should consider the correction term, so called “standard state correction31”, originated from the difference between the sampled unbound volume Vu and standard volume 𝑉𝑜( = 1661 Å) as27, 1

𝑉𝑢

∆𝐺 ∘ = ∆𝐺𝑃𝑀𝐹 ― 𝛽ln 𝑉𝑜.

(5)

∆𝐺𝑃𝑀𝐹 is given by, 1

𝑄𝑏

𝑉𝑏

1

∆𝐺𝑃𝑀𝐹 = ― 𝛽ln 𝑄𝑢= ― 𝛽ln 𝑉 𝑒 ―𝛽∆𝑊

(6)

𝑢


Page 7 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where 𝑄𝑏 = ∫𝑏𝑒 ―𝛽𝑊(𝑟)𝑑𝑟 is integrated over the bound state. Eq. (6) is obtained also considering 𝑄𝑢 = ∫𝑢𝑒 ―𝛽𝑊(𝑟)𝑑𝑟 = 𝑉𝑢𝑒 ―𝛽∆𝑊 integrated over the unbound state because W(r) converges to ∆𝑊 in the unbound state. From Eqs (5) and (6), we obtain, ∆𝐺 ∘ = ―∆𝑊 + ∆𝐺𝑣,

(6)

where the second term is the so-called volume correction of the free energy difference, 1

𝑉𝑏

∆𝐺𝑣 = ― 𝛽ln 𝑉𝑜,

(7)

Comparing to Eq. 7 of the literature27, Eq. (5) in this work does not require a correction term for the restrained simulation because we used unbiased simulations. Procacci and Chelli recently assessed the methods to evaluate the standard state correction31, and discussed the differences in the definition of the bound volume. In the method used in the literature28, 𝑉𝑏 is obtained by 𝑉𝑏(𝑑) = ∫𝑑(𝒓) ≤ 𝑑𝑒 ―𝛽𝑊(𝒓)𝑑𝒓.

(8)

This method was successfully used to estimate the standard binding free energies28,30,40,41,32–39. In this treatment, 𝑉𝑏 is defined as the weighted volume over all the range in which the interactions between the two molecules are non-negligible. For large enough d, Vb(d) was shown to be constant28. We also examined another method in which the integration was made in the range W(r) = 0 as, 𝑉𝑏 = ∫𝑊(𝒓) = 0𝑒 ―𝛽𝑊(𝒓)𝑑𝒓.

(9)

In this case, 𝑉𝑏 only considers the physical volume in the tightly bound state. By using Eq. (8), 𝑉𝑏 should be estimated to be relatively large, which gives the upper limit of the volume estimation in this paper. In contrast, 𝑉𝑏 obtained by Eq. (9) shows the lower limit. The results from both treatments will be shown in Results section. In order to make comparison of the dimensional effect in MSM, we also examined the third model to construct the MSM based on the Cα coordinates of TAD-p53 relative to MDM2, namely Cα-MSM. To generate the time series trajectories of the Cα coordinates relative to MDM2, we first performed a least squared fitting of all snapshots to the MDM2 backbone atoms in the initial snapshot of the PaCS-MD, which yielded 51 dimensional space (17 residues ×3). We took advantage of the time lagged independent component analysis (TICA)42–44 to reduce the dimension to 45, which maintains 90% of the original information. We also examined the results without the reduction of the dimension in one case. To overcome the high dimensional integral in the volume correction, the standard binding free energy is estimated by the equation 5 in the literature45 as follows. 7 ACS Paragon Plus Environment


∆𝐺 ∘ = ― 𝑘𝐵𝑇ln

( ) ― 𝑘 𝑇ln ( ) 𝑝𝑏𝑜𝑢𝑛𝑑 𝑒𝑞 𝑝𝑏𝑢𝑙𝑘 𝑒𝑞

𝑉𝑢

𝐵

𝑉0

Page 8 of 30

(10)

where the first and second terms correspond to ‒ΔW and ‒ΔGv, respectively. Vu is the unbound volume, which is defined as a volume of the actual sampling volume of the unbound state which will be defined in Section 3.2.

The association and dissociation rate constants, kon and koff, are the inverse of the Mean First Passage Time (MFPT) of the association and dissociation as, 1

(11)

𝑘𝑜𝑛 = 𝑀𝐹𝑃𝑇𝑜𝑛.𝐶𝑐𝑜𝑚𝑝 1

(12)

𝑘𝑜𝑓𝑓 = 𝑀𝐹𝑃𝑇𝑜𝑓𝑓

in which 𝐶𝑐𝑜𝑚𝑝 is the ligand concentration in the MD simulation. MFPTs can easily be computed using the transition matrix T28. We used MSM Builder package46 and pyEMMA package47 for constructing COM-based MSMs (1D and 3D) and C- MSM, respectively.

2.3. Model set-up and simulation procedure We adopted the crystal structure of the MDM2/TAD-p53 complex (PDB ID: 1YCQ23). Although residues 1329 of TAD-p53 were included in the sample, the structure of residues 17-27 were only determined by X-ray crystallography (hereafter ‘crystal region’) and those of Pro13, Leu14, Ser15, Gln16, Glu28, and Asn29 were not solved. We generated the missing residues and hydrogen atoms, and capped the C- and N-termini with acetyl and N-methyl groups, respectively, as shown in Fig. 1. The complex was solvated into a box of 18.8 × 9.9 × 8.9 nm3 with TIP3P water molecules48. The initial inter-COM vector from MDM2 to TAD-p53 was oriented parallel to the x axis. Sodium and chloride ions were added to realize an ion concentration of 0.15 M and achieve charge neutrality. AMBER ff99SB-ILDN force field49 was used for the protein complex. All simulations were performed using GROMACS 5.1.250. The simulations to realize the equilibrated state were carried out as follows: 1) The solvated models were energy-minimized by the steepest descent method followed by the conjugate gradient method with positional restraints applied on heavy atoms of the crystal region (force constant: 1000 kJ/mol.nm2). 2) The system with the same restraints was heated from 0 to 400 K within 1 ns and equilibrated at 400 K for 500 ps, then cooled down to 300 K in the next 500 ps, and equilibrated at 300 K for 1 ns. This step was conducted to ensure that the modeled residues were in the lower energy configuration. In this step, the NVT ensemble simulation was employed. 3) MD 8 ACS Paragon Plus Environment

Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulation with the NPT ensemble was conducted for the next 100 ps at 300 K with a relaxation time of 0.1 ps for the heat bath coupling, and at 1.0 atm with a relaxation time of 2.0 ps for the isotropic pressure coupling. 4) The force constant of positional restraints was reduced by 100 kJ/mol.nm2 every 100 ps until it vanished (total 0.9 ns). 5) The system was equilibrated for 10 ns using the NPT ensemble with weak positional restraints on the backbone atoms of MDM2 in order to prevent the translation and rotation of MDM2 during the simulation. To suppress the translation and rotation of MDM2, the positional restraint (force constant: 100 kJ/mol.nm2) was imposed also during PaCS-MD on the backbone atoms of MDM2 excluding the binding interface region (binding interface is defined via CAPRI criteria51). The equation of motion was integrated by the leap-frog scheme52 with bond constraints by the LINCS method (steps 2-4)53 and by velocity Verlet methods54 without bond constraints (step 5 and production runs). The isothermal condition was established by the velocity rescaling (steps 2-3)55 and Nosé-Hoover methods (steps 4,5 and production runs)56,57. The isobaric condition was achieved using the Berendsen (step 3)58, ParrinelloRahman (step 4)59, and MTTK barostats (step 5 and production runs)60. After the equilibration, we performed 1 µs MD simulation, which was used to examine the stability of the complex. We also selected 25 different structures every 10 ns from the beginning of the 1 µs MD trajectory. From each of these structures, 1 ns MD was conducted (hereafter cycle 0), and 30 snapshots of the first PaCS-MD cycle were selected. Starting from these configurations, we simulated dissociation of the MDM2/TAD-p53 complex by PaCS-MD, similar to the case of the protein/ligand complex22. We used 30 replicas and 0.1 ns MD simulations for each cycle. Snapshots were recorded every 0.5 ps and used for the analysis. 25 PaCS-MD trials were conducted and the total simulation time of PaCS-MD was 3.895 μs including cycle 0.

3. Results 3.1. Interactions between MDM2/TAD-p53 and stability of the complex The binding interface between MDM2 and TAD-p53 consists of one helix and two turns of MDM2, and one helix of TAD-p53 (Fig. 1a and b). A few positively charged residues (mainly lysines) are situated on the binding interface of MDM2, while TAD-p53 has two negatively charged glutamic acids (GLU17 and GLU28) (Fig. 1c). There are five hydrogen bonds (SER15–LYS66, GLU17–LYS90, PHE19–GLN68, TRP23–ILE50, PRO27–TYR96. The orders are: MDM2 residue–TAD-p53 residue.) and 17 hydrophobic interactions in the binding interface as shown (Fig. 1d). These key interactions were maintained as in the crystal structure during the equilibration simulation. PHE19 and TRP23, which are thought to be key amino acids23, stayed inside the hydrophobic cleft of MDM2. These results indicate that weak but specific binding of the complex is maintained, consistent with the NMR observations by Chi et al.61 Computational mutagenesis by Moreira et al. suggested that MET20 and TYR22 are new hotspots, in addition to PHE19, TRP23 and LEU2662. Moreover, they also indicated 9 ACS Paragon Plus Environment


Page 10 of 30

that backbone hydrogen bonds of p53 with MDM2 are important. The helical region of TAD-p53 spans from PHE19 to LEU25 as found in the work by Joseph et al.63 To examine the stability of the complex during the typical time scale of MD, we carried out 1 µs MD after equilibration. As shown in Fig. S1, the inter-COM distance was 1.4 ± 0.6 nm for the entire region of the complex and 1.3 ± 0.4 nm for the residues of the crystal region. The values after ‘±’ indicate the standard deviation hereafter unless noted. The root-mean-square deviation (RMSD) of TAD-p53 from the structure at time t = 0 shows the stability of the crystal region, while the RMSD for the whole region indicates large fluctuations. The latter consistently indicates flexibility of the unsolved regions of TAD-p53 in the crystal structure.

3.2. Dissociation pathways generated by PaCS-MD PaCS-MD easily generated dissociation pathways within short MD runs, which took 51.6 ± 15.5 cycles (average ± standard deviation) of 0.1 ns MDs and a total of 154.8 ± 46.4 ns MD until reaching the complete dissociation defined by the inter-COM distance d = 7.0 nm. Similar to the case of the LYZ/triNAG dissociation22, the dissociation of MDM2/TAD-p53 can be split into 3 steps, corresponding to 3 macrostates, i.e., bound, partiallybound, and unbound states as in Fig. 2. In the bound state (d ≤ 1.90 nm), most of the key interactions between MDM2 and TAD-p53 were maintained, while a part of TAD-p53 dissociated from MDM2 in the partially bound state (3.45 ≥ d > 1.90 nm, gray shaded area in Fig. 2). In the unbound state (d > 3.45 nm), TAD-p53 completely dissociated and d increased linearly. To visualize the TAD-p53 dissociation pathways, we generated the reactive trajectories from each trial, defined as the concatenated trajectory that connects the initial bound and the final unbound state along the dissociation pathway20,22 as shown in Fig. 3a. By looking into the details of the reactive trajectories and observing the solvation of the binding interface, we found that dissociation and solvation of the interface residues occurred first from either PHE19 or LEU26 of TAD-p53 (12 among 25 trials for each case). Interestingly, PHE19, TRP23 and LEU26 dissociated altogether in one case. These results indicated that dissociation gradually evolves, starting from either side of the TAD-p53 helix with almost equal probability, which likely minimizes the energy change along the dissociation pathway. Solvation of the binding interface took place from the MDM2 helix (THR45-GLN61) where water molecules mediated and broke the hydrophobic interactions between TAD-p53 and MDM2, leading to complete dissociation of TAD-p53. Distribution of the TAD-p53 COM positions from all the trajectories from one trial (shaded light purple region spheres in Fig. 3a) and all the trials (Fig. 3b) were then investigated. As seen in Fig. 3a, a reactive trajectory is surrounded by the other trajectories in the trial, forming a continuous conformational space connecting the bound and unbound states, similar to the previous case22. To quantify the sampled conformational space, we examined the 10 ACS Paragon Plus Environment

Page 11 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


effective diameter of the conformational surface, σ(d), as a cross section of trajectories as a function of d. σ(d) is defined as σ(d) =(4S(d)/π)1/2 where S(d) is a cross section of the sampled TAD-p53 COM positions at d. After the dissociation, the value of σ(d) for 1D-MSM averaged over the trials gradually increased in the partially bound state and reached a converged value of 0.6 – 0.8 nm at around d = 3.5 nm (Fig. 3c). In the case of 3D-MSM constructed from merged trajectories of all PaCS-MD trials, the sampled conformational space expanded into a cone-like shape (Fig. 3b), which resulted in a linear increase in σ(d) (Fig. 3c). Since TAD-p53 is expected to be intrinsically disordered, we monitored its fraction of helical content along dissociation pathways (Fig. 3d). On average, the fraction was around 0.3 in the bound state and decreased to 0.2 upon dissociation, which is consistent with the disordered nature of TAD-p53. Interestingly, individual fractions largely depended on the trial. In two cases (thin red and black lines), the helical structure was completely lost, whereas in one case (thin yellow line), helical structure was completely lost in the partially-bound state but partly refolded in the unbound state. Since PaCS-MD does not apply any force-bias during the simulation, we expect that more natural dissociation pathways were generated.

3.3. Free energy analysis by 1D-MSM To calculate PMF in 1D-MSM, we employed the K-Means method64 with the K-Mean++ algorithm65 to select the initial guess of the cluster centers. As mentioned in the Method section, the MD trajectories were recorded every 0.5 ps. We chose 50 clusters, which is a reasonable number for determining smooth PMF with an average intercluster distance of 0.1 nm as shown in our previous work22. We carefully monitored the convergence of the cluster centers throughout multiple trials of K-Means and confirmed that the cluster centers were widely distributed. To construct better MSM, we performed an estimation of the implied timescale versus different lag times (τ) in the reversible condition and determined the optimal lag time as shown in Fig. S2a. Although evolution of the implied timescales versus τ does not directly guarantee quality of the constructed MSM, it does suggest the optimal lag time66. We found that the flat distribution of the implied time spans from 10 to 40 ps in all of the ten slowest states and that disconnected states only occur when τ > 40 ps. Although the plateau regime of the slowest motion still continues from τ = 45 to 50 ps (the largest implied timescale), MSM contains more than 2 strongly connected components, of which the largest population spans from 85% to 97%, depending on the trial and chosen lag time. Doerr and De Fabritiis suggested that the lag time should be chosen as the minimum of the maximum lag time before disconnecting the states and half of the trajectory length.45 Therefore, we judged τ = 40 ps to be the optimum value for 1D-MSM. To further analyze the MSM constructed with different conditions, we first examined the reversible MSMs of fully connected states, which are denoted as rMSM and default MSM methods unless mentioned, at 3 different lag 11 ACS Paragon Plus Environment


Page 12 of 30

times of 20, 30 and 40 ps. Since all the processes should be reversible, the use of the detailed balance in constructing MSM is reasonable. However, it is also interesting to examine the effect of the detailed balance condition. For this purpose, we also built the irreversible MSM (irMSM) without using the detailed balance condition. Unfortunately, the irMSM can be constructed only up to τ = 20 ps, and longer lag times yielded negative eigenvectors of the transition matrices. This suggests that more statistics are needed for irMSM since the longer τ leads to less statistics in the calculation of the transition matrix. Figure 3e shows the PMFs that were obtained as the average over 25 trials. The error bars in the figure indicate standard error of the mean. PMFs show linear increments in the bound state (d ≤ 1.90 nm) and convergence in the unbound state (d > 3.45 nm) in all four cases, which is essential for accurate calculation of ΔW. To calculate ΔW by Eq. (4), W(dbound) was obtained as the average in the range 4.00 ≥ d > 3.50 nm and W(dunbound) = 0. The PMF value of rMSM in this region increases as the lag time becomes longer. Compared to rMSM with the optimal τ = 40 ps, ΔW in irMSM with τ = 20 ps was slightly larger. Standard binding free energy ΔGº was obtained by Eq. (5). Note that the volume integral in Eq. (8) is calculated using W(r) and σ(d) as shown in Fig. S3. The volume correction term ΔGv completely converged at the end of the bound state d = 1.90 nm, showing good agreement with the previous work by Buch et al.28, also indicating that our choice of the bound state was reasonable. It should be noted that d is ~ 1.3 nm in the crystal complex structure. We also examined another volume correction method by Eq. (9) in some cases. Calculated standard binding free energies after the volume correction are shown in Table 1. The correction terms calculated by Eq. (8) better reproduced the standard binding free energy. The calculated binding free energy with the optimum lag time (τ = 40 ps) shows the best agreement with the experimental values obtained by isothermal titration calorimetry (ITC)67, surface plasmon resonance (SPR)68, and dynamic force spectroscopy (DFS)69. Δ Gº obtained by irMSM was also comparable to the experimental ΔGº. As expected from the previous work on LYZ/triNAG, 1D-MSM is sufficient to calculate the binding free energy. Interestingly, the binding free energy does not depend greatly on the residues considered (Table 1), suggesting that MDM2 only binds to the core region of TAD-p53.

3.4. Free energy and kinetic rate analysis by 3D- and Cα-MSM Although 1D-MSM has the power to calculate the standard binding free energy accurately, each trial only considers a straight dissociation pathway from the bound to unbound states because PaCS-MD enhances sampling only in the radial direction by selecting snapshots with longer inter-COM distances. To consider movements toward other directions, we merged all the trajectories obtained by all the PaCS-MD trials. As shown in Fig. 3b, the trajectories of the 3D COM positions overlap greatly in conformational space, which enables us to construct the 3D-MSM. For this purpose, the COM positions of TAD-p53 were shifted such that the backbone atoms of MDM2 12 ACS Paragon Plus Environment

Page 13 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


were in the best-fit position. To determine microstates, the spatial clustering method70 with two cluster distances, dclus = 0.10 and 0.12 nm, was used, which resulted in 10098 and 6606 clusters, respectively (Table 1). The sampling frequency in 3D-MSM was every 1.0 ps to reduce memory usage during the MSM calculation, and the snapshots with inter-COM distances d ≤ 4.0 nm were used. By carefully checking the relationship between lag time and implied timescale (Fig. S2b and c), τ = 40 ps was identified as the best choice, which was the same as in the case of 1D-MSM. The PMFs obtained by 3D-MSM with τ = 20, 30, and 40 ps are shown in Fig. 4. The standard error of the mean of the PMF was calculated by the bootstrap method, by randomly choosing 90% of the snapshots for 1,000 times in the case of dclus = 0.12 nm. We did not complete the bootstrap method with dclus = 0.10 nm because it was too time consuming. With τ = 40 ps, the values of ‒∆W in 3D-MSM were higher than that in 1D-MSM, however, after the volume correction, the calculated ∆Go values were mutually in good agreement, and comparable to the experimental values (Table 1). We also calculated the dissociation and association rate constants, koff and kon, by Eqs (8) and (9) (Table 1). Figure 5 shows the lag time dependence of koff and kon. It also shows the convergence tendency of both koff and kon with the increment of lag time. As far as dclus ≤ 0.15 nm, the calculated koff exhibits agreement with the range of the experimental values (shown in gray shaded area in Figure 5) whereas kon was overestimated by two orders of magnitude. This is probably due to the difference between dissociation and association processes as discussed below. In this work, we generated many dissociation pathways, which formed a cone-like shape as shown in Fig. 3b. The agreement between the calculated and experimental koff values suggests that major dissociation events should happen within this cone-like volume and contributions of the pathways passing through other regions to koff are expected to be negligible. This is a reasonable assumption because dissociation always starts from the same binding pose, and once the ligand dissociates far enough, its contribution should be small. Our method is more efficient in estimating the kinetics rates; the current calculations only used 3.895 μs in total MD time (the total number of PaCSMD cycles × 0.1 ns × nrep), whose computational cost was 2/3 of that in the previous study by Metadynamics7, and 1/15 of that in extended WE sampling17. kon calculated by 3D-MSM was clearly overestimated. Although this 3D-MSM cannot directly estimate kon from the dissociation PaCS-MD trajectories, kon (kon,eff) can be effectively calculated from ∆Go and koff from the relation, ∆Go = −β-1 ln (Co kon,eff / koff) where Co = 1 M is the standard concentration. As shown in Table 1, kon,eff are comparable to the experimental kon values. Since direct accurate calculation of kon should require additional computational investment, indirect calculation of kon from ∆Go and koff can be used as a method to drastically reduce the simulation time.



Page 14 of 30

To yield deeper insight into the overestimation of kon, we also conducted Cα-MSM in which the number of dimensions were drastically increased by using the Cα coordinates of TAD-p53 after performing least-squares fitting of MDM2 and the TICA transformation (details in calculation section). The lag time evolution of implied timescale shows the convergence from 20 ps (Fig S4). We found that the bound state (d > 3.5 nm) is in the region of PMF lower than -10.0 kcal/mol while the unbound state (d < 1.9 nm) corresponds to the region of PMF > -4.5 kcal/mol (Fig S5). We further estimated the binding free energy and kinetic rates as shown in Table 1. Interestingly, the results with τ = 50 ps show good agreements not only in ∆Go and koff. but also in kon with those from the ITC and stopped-flow experiments67. We also examined the results without the reduction of the TICA dimension in one case, but the results were comparable with those with the reduction.

4. Discussion and Conclusions In this work, we used PaCS-MD to generate multiple dissociation pathways for the MDM2/TAD-p53 complex. Despite its simplicity, PaCS-MD provided a thorough sampling of the pathways without applying force-bias to the system and easily achieved complete dissociation in an average simulation time of 154.8 ± 46.4 ns. It should be noted that PaCS-MD does not use any biasing force to guide the trajectories because conventional MD simulation is conducted for each replica. Overall accelerated increase of inter-COM distance observed during PaCS-MD is only caused by the selection of the naturally-occurred fluctuations with longer inter-COM distances and by their use in the next cycle as the initial structures. Therefore, PaCS-MD does not deform distribution of any parameter in each MD. For example, in the previous study22, we confirmed that the diffusion constant of a dissociated ligand was equivalent to that of the free diffusion constant. In combination with MSM, PaCS-MD was shown to be very efficient in calculating the binding free energy, not only for the protein-ligand complex22 but also for a more flexible protein-peptide complex. The binding free energy ∆Go obtained by both 1D- and 3D-MSM was in good agreement with those determined experimentally. Using 3D-MSM based on the COM position of TAD-p53 relative to MDM2, the dissociation rate constant koff was calculated, which is comparable to those measured experimentally. The optimum lag time (τ = 40 ps) was chosen based on the lag time−implied time relationship and the connection of the states. At the same time, this choice provided the best agreement of the calculated ∆Go and koff as shown in Fig. 3 and Table 1. Interestingly, PaCS-MD generated dissociation pathways that contained complete or partial unfolding of the helical region of the MDM2-bound TAD-p53 as expected in the literature17 and references therein. This suggests the ability of PaCS-MD to generate more natural dissociation pathways since it consists of unbiased MD simulations. In this work, we showed that, with a very simple reaction coordinate (inter-COM distance), ∆Go can be estimated in agreement with experimental values within acceptable range of error. Given that sampling is enough 14 ACS Paragon Plus Environment

Page 15 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in statistics, the agreement with experimental data might be related to the fact that the free energy difference is independent of the choice of pathway and only depends on the bound state and the unbound state, which can be any state far enough from the bound state. However, this is not the case with kinetic rates where thorough sampling of all possible states are needed. In fact, we made efforts to estimate the kinetic rates with 1D-MSM without complete success, which showed instability of the calculation. Some trials gave values very close to the experimental koff (for example, 0.83.2 s-1), while others were three-orders higher than the experimental value (1019102.4 s-1). In contrast to 1D-MSM, the results of 3D-MSM were in good agreement with the experimental values. To further examine the stability of 3D-MSM, we performed a bootstrapping calculation, in which 17 PaCS-MD trials (half of the dataset) from 25 trials were randomly selected in the calculation of koff with dclus = 0.12 nm and τ = 40 ps. The obtained koff was 45.8 ± 1.2 s-1, which was one-order greater than that with a small standard error, indicating the importance of thorough sampling along the directions perpendicular to the dissociation directions. While 3D-MSM overestimated kon, Cα-MSM estimated it in agreement with the experimental values. This implies that the calculation of kon requires higher resolution in MSM in contrast to koff where 3D-MSM is sufficient. TAD-p53 was shown to make conformational changes upon the binding into the cleft as suggested in the literatures17,71, which might causes the overestimation of the association event because the same COM positions can contain multiple conformations of TAD-p53 which probably should be distinguished as distinct microstates. During the process of the dissociation, TAD-p53 cannot make large conformational change in the bound and partially-bound states, which may be related to the results that 3D-MSM was sufficient to calculate koff. These results indicate that the resolution of MSM affects the quantities which can be estimated.

Acknowledgement This research was supported by MEXT/JSPS KAKENHI (No. JP25104002 and JP15H04357) to A.K., and by MEXT as "Priority Issue on Post-K Computer” (Building Innovative Drug Discovery Infrastructure through Functional Control of Biomolecular Systems) to A.K. The computations were partly performed using the supercomputers at the RCCS, The National Institute of Natural Science, and ISSP, The University of Tokyo. This research also used computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID: hp150270, hp160207, hp170254, and hp180201). D.T. acknowledges financial support from Hitachi Scholarship Program No. S175 of the Hitachi Global Foundation.

Supporting Information



Page 16 of 30

Stability of MDM2/TAD-p53 during 1 μs conventional MD, implied timescales versus lag time for 1D- and 3DMSM, volume correction term of the binding free energy, implied timescales versus lag time for Cα-MSM, and free energy landscape of the MDM2/TAD-p53 complex.

References (1)

Paul, F.; Wehmeyer, C.; Abualrous, E. T.; Wu, H.; Crabtree, M. D.; Schöneberg, J.; Clarke, J.; Freund, C.; Weikl, T. R.; Noé, F. Protein-Peptide Association Kinetics beyond the Seconds Timescale from Atomistic Simulations. Nat. Commun. 2017, 8 (1), 1095.

(2)

Rydzewski, J.; Jakubowski, R.; Nowak, W.; Grubmüller, H. Kinetics of Huperzine A Dissociation from Acetylcholinesterase via Multiple Unbinding Pathways. J. Chem. Theory Comput. 2018, 14 (6), 2843– 2851.

(3)

Izrailev, S.; Stepaniants, S.; Isralewitz, B.; Kosztin, D.; Lu, H.; Molnar, F.; Wriggers, W.; Schulten, K. Steered Molecular Dynamics. Comput. Mol. Dyn. Challenges, Methods, Ideas SE - 2 1999, 4, 39–65.

(4)

Kokh, D. B.; Amaral, M.; Bomke, J.; Grädler, U.; Musil, D.; Buchstaller, H. P.; Dreyer, M. K.; Frech, M.; Lowinski, M.; Vallee, F.; et al. Estimation of Drug-Target Residence Times by τ-Random Acceleration Molecular Dynamics Simulations. J. Chem. Theory Comput. 2018, 14 (7), 3859–3869.

(5)

Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90 (23), 238302.

(6)

Casasnovas, R.; Limongelli, V.; Tiwary, P.; Carloni, P.; Parrinello, M. Unbinding Kinetics of a P38 MAP Kinase Type II Inhibitor from Metadynamics Simulations. J. Am. Chem. Soc. 2017, 139 (13), 4780–4788.

(7)

Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of Protein–Ligand Unbinding: Predicting Pathways, Rates, and Rate-Limiting Steps. Proc. Natl. Acad. Sci. 2015, 112 (5), E386–E391.

(8)

E, W.; Ren, W.; Vanden-Eijnden, E. String Method for the Study of Rare Events. Phys. Rev. B 2002, 66 (5), 052301.

(9)

Maragliano, L.; Roux, B.; Vanden-Eijnden, E. Comparison between Mean Forces and Swarms-ofTrajectories String Methods. J. Chem. Theory Comput. 2014, 10 (2), 524–533.

(10)

Májek, P.; Elber, R. Milestoning without a Reaction Coordinate. J. Chem. Theory Comput. 2010, 6 (6), 1805–1817.

(11)

Fu, H.; Shao, X.; Chipot, C.; Cai, W. Extended Adaptive Biasing Force Algorithm. An On-the-Fly Implementation for Accurate Free-Energy Calculations. J. Chem. Theory Comput. 2016, 12 (8), 3506–


Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3513. (12)

Zhao, T.; Fu, H.; Lelièvre, T.; Shao, X.; Chipot, C.; Cai, W. The Extended Generalized Adaptive Biasing Force Algorithm for Multidimensional Free-Energy Calculations. J. Chem. Theory Comput. 2017, 13 (4), 1566–1576.

(13)

Hamelberg, D.; Mongan, J.; McCammon, J. A. Accelerated Molecular Dynamics: A Promising and Efficient Simulation Method for Biomolecules. J. Chem. Phys. 2004, 120 (24), 11919–11929.

(14)

Tang, Z.; Chang, C. A. Systematic Dissociation Pathway Searches Guided by Principal Component Modes. J. Chem. Theory Comput. 2017, 13 (5), 2230–2244.

(15)

Huber, G. A.; Kim, S. Weighted-Ensemble Brownian Dynamics Simulations for Protein Association Reactions. Biophys. J. 1996, 70 (1), 97–110.

(16)

Bhatt, D.; Zhang, B. W.; Zuckerman, D. M. Steady-State Simulations Using Weighted Ensemble Path Sampling. J. Chem. Phys. 2010, 133 (1), 014110.

(17)

Zwier, M. C.; Pratt, A. J.; Adelman, J. L.; Kaus, J. W.; Zuckerman, D. M.; Chong, L. T. Efficient Atomistic Simulation of Pathways and Calculation of Rate Constants for a Protein-Peptide Binding Process: Application to the MDM2 Protein and an Intrinsically Disordered P53 Peptide. J. Phys. Chem. Lett. 2016, 7 (17), 3440–3445.

(18)

Dickson, A.; Lotz, S. D. Multiple Ligand Unbinding Pathways and Ligand-Induced Destabilization Revealed by WExplore. Biophys. J. 2017, 112 (4), 620–629.

(19)

Lotz, S. D.; Dickson, A. Unbiased Molecular Dynamics of 11 Min Timescale Drug Unbinding Reveals Transition State Stabilizing Interactions. J. Am. Chem. Soc. 2018, 140 (2), 618–628.

(20)

Harada, R.; Kitao, A. Parallel Cascade Selection Molecular Dynamics (PaCS-MD) to Generate Conformational Transition Pathway. J. Chem. Phys. 2013, 139 (3), 035103.

(21)

Harada, R.; Kitao, A. Nontargeted Parallel Cascade Selection Molecular Dynamics for Enhancing the Conformational Sampling of Proteins. J. Chem. Theory Comput. 2015, 11 (11), 5493–5502.

(22)

Tran, D. P.; Takemura, K.; Kuwata, K.; Kitao, A. Protein–Ligand Dissociation Simulated by Parallel Cascade Selection Molecular Dynamics. J. Chem. Theory Comput. 2018, 14 (1), 404–417.

(23)

Kussie, P. H.; Gorina, S.; Marechal, V.; Elenbaas, B.; Moreau, J.; Levine, A. J.; Pavletich, N. P. Structure of the MDM2 Oncoprotein Bound to the P53 Tumor Suppressor Transactivation Domain. Science (80-. ). 1996, 274 (5289), 948–953.

(24)

Momand, J.; Wu, H.-H.; Dasgupta, G. MDM2 — Master Regulator of the P53 Tumor Suppressor Protein.



Page 18 of 30

Gene 2000, 242 (1), 15–29. (25)

Zhao, Y.; Aguilar, A.; Bernard, D.; Wang, S. Small-Molecule Inhibitors of the MDM2-P53 ProteinProtein Interaction (MDM2 Inhibitors) in Clinical Trials for Cancer Treatment. J Med Chem 2015, 58 (3), 1038–1052.

(26)

Bowman, G. R.; Pande, V.; Noé, F. An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation; 2014; Vol. 797.

(27)

Doudou, S.; Burton, N. A.; Henchman, R. H. Standard Free Energy of Binding from a One-Dimensional Potential of Mean Force. J. Chem. Theory Comput. 2009, 5 (4), 909–918.

(28)

Buch, I.; Giorgino, T.; De Fabritiis, G. Complete Reconstruction of an Enzyme-Inhibitor Binding Process by Molecular Dynamics Simulations. Proc. Natl. Acad. Sci. U. S. A. 2011, 108 (25), 10184–10189.

(29)

Gumbart, J. C.; Roux, B.; Chipot, C. Efficient Determination of Protein–Protein Standard Binding Free Energies from First Principles. J. Chem. Theory Comput. 2013, 9 (8), 3789–3798.

(30)

Limongelli, V.; Bonomi, M.; Parrinello, M. Funnel Metadynamics as Accurate Binding Free-Energy Method. Proc. Natl. Acad. Sci. 2013, 110 (16), 6358–6363.

(31)

Procacci, P.; Chelli, R. Statistical Mechanics of Ligand-Receptor Noncovalent Association, Revisited: Binding Site and Standard State Volumes in Modern Alchemical Theories. J. Chem. Theory Comput. 2017, 13 (5), 1924–1933.

(32)

Buch, I.; Harvey, M. J.; Giorgino, T.; Anderson, D. P.; De Fabritiis, G. High-Throughput All-Atom Molecular Dynamics Simulations Using Distributed Computing. J. Chem. Inf. Model. 2010, 50 (3), 397– 403.

(33)

Gumbart, J. C.; Roux, B.; Chipot, C. Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy? J. Chem. Theory Comput. 2013, 9 (1), 794–802.

(34)

De Vivo, M.; Masetti, M.; Bottegoni, G.; Cavalli, A. Role of Molecular Dynamics and Related Methods in Drug Discovery. J. Med. Chem. 2016, 59 (9), 4035–4061.

(35)

Sun, H.; Li, Y.; Li, D.; Hou, T. Insight into Crizotinib Resistance Mechanisms Caused by Three Mutations in ALK Tyrosine Kinase Using Free Energy Calculation Approaches. J. Chem. Inf. Model. 2013, 53 (9), 2376–2389.

(36)

Biedermannová, L.; Prokop, Z.; Gora, A.; Chovancová, E.; Kovács, M.; Damborsky, J.; Wade, R. C. A Single Mutation in a Tunnel to the Active Site Changes the Mechanism and Kinetics of Product Release in Haloalkane Dehalogenase LinB. J. Biol. Chem. 2012, 287 (34), 29062–29074.


Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(37)

Chen, L. Y. Hybrid Steered Molecular Dynamics Approach to Computing Absolute Binding Free Energy of Ligand–Protein Complexes: A Brute Force Approach That Is Fast and Accurate. J. Chem. Theory Comput. 2015, 11 (4), 1928–1938.

(38)

Doudou, S.; Sharma, R.; Henchman, R. H.; Sheppard, D. W.; Burton, N. A. Inhibitors of PIM-1 Kinase: A Computational Analysis of the Binding Free Energies of a Range of Imidazo [1,2-b] Pyridazines. J. Chem. Inf. Model. 2010, 50 (3), 368–379.

(39)

Heinzelmann, G.; Bastug, T.; Kuyucak, S. Mechanism and Energetics of Ligand Release in the Aspartate Transporter GltPh. J. Phys. Chem. B 2013, 117 (18), 5486–5496.

(40)

Li, Y.; Li, X.; Dong, Z. Exploration of Gated Ligand Binding Recognizes an Allosteric Site for Blocking FABP4-Protein Interaction. Phys. Chem. Chem. Phys 2015, 17, 32257.

(41)

Niu, Y.; Li, S.; Pan, D.; Liu, H.; Yao, X. Computational Study on the Unbinding Pathways of B-RAF Inhibitors and Its Implication for the Difference of Residence Time: Insight from Random Acceleration and Steered Molecular Dynamics Simulations. Phys. Chem. Chem. Phys 2016, 18, 5622–5629.

(42)

Molgedey, L.; Schuster, H. G. Separation of a Mixture of Independent Signals Using Time Delayed Correlations. Phys. Rev. Lett. 1994, 72 (23), 3634–3637.

(43)

Pérez-Hernández, G.; Paul, F.; Giorgino, T.; De Fabritiis, G.; Noé, F. Identification of Slow Molecular Order Parameters for Markov Model Construction. J. Chem. Phys. 2013, 139 (1), 7653.

(44)

Schwantes, C. R.; Pande, V. S. Improvements in Markov State Model Construction Reveal Many NonNative Interactions in the Folding of NTL9. J. Chem. Theory Comput. 2013, 9 (4), 2000–2009.

(45)

Doerr, S.; de Fabritiis, G. On-the-Fly Learning and Sampling of Ligand Binding by High-Throughput Molecular Simulations. J. Chem. Theory Comput. 2014, 10 (5), 2064–2069.

(46)

Beauchamp, K. A.; Bowman, G. R.; Lane, T. J.; Maibaum, L.; Haque, I. S.; Pande, V. S. MSMBuilder2: Modeling Conformational Dynamics on the Picosecond to Millisecond Scale. J. Chem. Theory Comput. 2011, 7 (10), 3412–3419.

(47)

Scherer, M. K.; Trendelkamp-Schroer, B.; Paul, F.; Pérez-Hernández, G.; Hoffmann, M.; Plattner, N.; Wehmeyer, C.; Prinz, J.-H.; Noé, F. PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models. J. Chem. Theory Comput. 2015, 11 (11), 5525–5542.

(48)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926–935.

(49)

Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E.



Page 20 of 30

Improved Side-Chain Torsion Potentials for the Amber Ff99SB Protein Force Field. Proteins Struct. Funct. Bioinforma. 2010, 78 (8), 1950–1958. (50)

Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29 (7), 845–854.

(51)

Lensink, M. F.; Méndez, R.; Wodak, S. J. Docking and Scoring Protein Complexes: CAPRI 3rd Edition. Proteins Struct. Funct. Bioinforma. 2007, 69 (4), 704–718.

(52)

Van Gunsteren, W. F.; Berendsen, H. J. C. A Leap-Frog Algorithm for Stochastic Dynamics. Mol. Simul. 1988, 1 (3), 173–185.

(53)

Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18 (12), 1463–1472.

(54)

Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. J. Chem. Phys. 1982, 76 (1), 637–649.

(55)

Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126 (1), 014101.

(56)

Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81 (1), 511–519.

(57)

Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31 (3), 1695–1697.

(58)

Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, a; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 (1984), 3684–3690.

(59)

Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52 (12), 7182–7190.

(60)

Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87 (5), 1117–1157.

(61)

Chi, S. W.; Lee, S. H.; Kim, D. H.; Ahn, M. J.; Kim, J. S.; Woo, J. Y.; Torizawa, T.; Kainosho, M.; Han, K. H. Structural Details on Mdm2-P53 Interaction. J. Biol. Chem. 2005, 280 (46), 38795–38802.

(62)

Moreira, I. S.; Fernandes, P. A.; Ramos, M. J. Protein–Protein Recognition: A Computational Mutagenesis Study of the MDM2–P53 Complex. Theor. Chem. Acc. 2008, 120, 533–542.


Page 21 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(63)

Joseph, T. L.; Lane, D.; Verma, C. Stapled Peptides in the P53 Pathway: Computer Simulations Reveal Novel Interactions of the Staples with the Target Protein. Cell Cycle 2010, 9 (22), 4560–4568.

(64)

MacQueen, J. Some Methods for Classification and Analysis of Multivariate Observations. In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics; Fifth Berkeley Symposium on Mathematical Statistics and Probability; University of California Press: Berkeley, Calif., 1967; pp 281–297.

(65)

Arthur, D.; Vassilvitskii, S. K-Means++: The Advantages of Careful Seeding. In Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms; 2007; pp 1027–1035.

(66)

Pande, V. S.; Beauchamp, K.; Bowman, G. R. Everything You Wanted to Know about Markov State Models but Were Afraid to Ask. Methods 2010, 52 (1), 99–105.

(67)

Schon, O.; Friedler, A.; Bycroft, M.; Freund, S. M. V.; Fersht, A. R. Molecular Mechanism of the Interaction between MDM2 and P53. J. Mol. Biol. 2002, 323 (3), 491–501.

(68)

Domenici, F.; Frasconi, M.; Mazzei, F.; D’Orazi, G.; Bizzarri, A. R.; Cannistraro, S. Azurin Modulates the Association of Mdm2 with P53: SPR Evidence from Interaction of the Full-Length Proteins. J. Mol. Recognit. 2011, 24 (4), 707–714.

(69)

Bizzarri, A. R.; Cannistraro, S. Free Energy Evaluation of the P53-Mdm2 Complex from Unbinding Work Measured by Dynamic Force Spectroscopy. Phys. Chem. Chem. Phys. 2011, 13 (7), 2738–2743.

(70)

Senne, M.; Trendelkamp-Schroer, B.; Mey, A. S. J. S.; Schütte, C.; Noé, F. EMMA: A Software Package for Markov Model Building and Analysis. J. Chem. Theory Comput. 2012, 8 (7), 2223–2238.

(71)

Chen, H. F.; Luo, R. Binding Induced Folding in P53-MDM2 Complex. J. Am. Chem. Soc. 2007, 129 (10), 2930–2937.

(72)

Laskowski, R. A.; Swindells, M. B. LigPlot+: Multiple Ligand–Protein Interaction Diagrams for Drug Discovery. J. Chem. Inf. Model. 2011, 51 (10), 2778–2786.

(73)

Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14 (1), 33–38.



Page 22 of 30

Table 1. Standard binding free energy (∆Go) and dissociation and association rate constants (koff and kon) of MDM2/TAD-p53 estimated by using PaCS-MD/MSM. The rows shown in bold indicate the results obtained by the optimum lag time. The values after ‘±’ indicate the standard error of the mean. More detailed information on the quantities is described in the main text.

Method

p53 residues

Reversible MSM

τ

# clusters

(ps)

-∆W

∆Gv

∆Go

koff

kon

kon,eff

(kcal/mol)

(kcal/mol)

(kcal/mol)

(s-1)

(107M-1s-1)

(107M-1s-1)

13-29

Y

20

50

-9.0±0.1

3.1±0.2

-5.9±0.2

13-29

Y

30

50

-10.3±0.1

3.1±0.2

-7.2±0.2

13-29

Y

40

50

-11.4±0.2

3.1±0.2

-8.3±0.3

13-29

N

20

50

-12.1±0.1

3.1±0.2

-9.0±0.2

13-29

Y

20

10098

-9.1±0.7

2.2

-6.9±0.7

56.47

904.86

0.6

13-29

Y

30

10098

-9.8±0.7

2.2

-7.6±0.7

9.98

622.6

0.34

13-29

Y

40

10098

-10.2±0.7

2.2

-7.9±0.7

3.27

477.93

0.19

13-29 a

Y

20

6606

-9.0±0.7

694.43

Y

30

6606

-10.0±0.7

7.49

512.99

13-29 a

Y

40

6606

-10.7±0.7

-6.8±0.7 (-5.8±0.7) -7.8±0.7 (-6.2±0.7) -8.5±0.7 (-7.0±0.7)

68.85

13-29 a

2.2 (3.8) 2.2 (3.8) 2.2 (3.7)

1.57

409.18

13-29

Y

30

3000

-4.63±0.4

-2.93

-7.56±0.4

334.958.8

3.90.3

13-29

Y

40

3000

-5.94±0.3

-2.93

-8.9±0.3

28.059.5

2.440.3

13-29 b

Y

50

3000

-5.5±0.7 (-5.6±0.5)

-2.93 (-2.93)

-8.4±0.7 (-8.5±0.7)

3.121.3 (1.790.7)

1.430.2 (3.370.3)

15-29

-9.1±0.1

2.06±0.09

0.92±0.004

17-26

-9.9±0.1

1.43±0.24

2.2±0.02

SPR68

Full length

-8.8±0.2

0.8±0.3

0.21±0.02

DFS69

Full length

-8.4±0.5

1D-MSM

3D-MSM

C-MSM

ITC and stoppedflow67

0.62 0.36 0.24

The correction term calculated by Eq. (8). The values in the parentheses were obtained by the correction in Eq. (9). b The values in the parentheses were obtained without dimensional reduction in TICA.

a


Page 23 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure Legends Figure 1. Visualization of the equilibrated MDM2/TAD-p53 complexes solvated in the simulation box: a. Overall arrangement of the simulation box; b. a close-up view and three important amino acids inside the binding cleft; c. View from the direction of the x-axis to show the electrostatic potential on the MDM2 binding cleft; d. LigPlot+ diagram72 to show the interactions between MDM2 and TAD-p53. The images in a-c were created by VMD73.

Figure 2. Inter-COM distance between MDM2 and TAD-p53 (d) as a function of the number of cycles in 25 PaCS-MD trials. The largest inter-COM distance in each cycle was plotted.

Figure 3. Features of TAD-p53 dissociation from MDM2. a. TAD-p53 dissociation pathways in 25 PaCS-MD trials represented by the COM positions of the reactive trajectories (small spheres in different colors). All the COM positions from one PaCS-MD trial are shown by shaded light purple around the corresponding reactive trajectory shown in dark blue. b. All the TAD-p53 COM positions in all the trials. Color differences denote different trials. c. Effective diameter of sampled COM area (σ(d)) plotted as a function of inter-COM distance d. Thick line indicates σ(d) averaged over the trials (error bars: standard deviation) and broken line shows the merged trajectories of all the trials. d. Fraction of helical content of TAD-p53. e. Potential of mean force (W(d)) as a function of d. In d, the values for average (thickest black line) and each trial (thinner black, red, yellow and gray lines) are shown. The images in a,b were created by VMD73.

Figure 4. Potential of mean force (W(d)) obtained by 3D-MSM with τ = 20 (red), 30 (green), and 40 ps (black). a. dclus = 0.10 and b. 0.12 nm. The error bars indicate standard error of the mean.



Page 24 of 30

Figure 5. Lag time dependence of a. association (kon) and b. dissociation rate constants (koff) in 3D-MSM. The shaded gray areas indicate the range of the experimental values shown in Table 1.


Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figures Figure 1.



Figure 2.


Page 26 of 30

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3.



Figure 4.


Page 28 of 30

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5.



TOC Graphic


Page 30 of 30

p53 Complex Investigated by Parallel

Recommend Documents