Unbinding Pathway by Enhanced Sampling

Subscriber access provided by Columbia University Libraries

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Exploring Ligand Binding/Unbinding Pathway by Enhanced Sampling Selectively of Ligand in Protein-Ligand Complex Qiang Shao, and Weiliang Zhu J. Phys. Chem. B, Just Accepted Manuscript • Publication Date (Web): 03 Sep 2019 Downloaded from pubs.acs.org on September 3, 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 40 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

Exploring Ligand Binding/Unbinding Pathway by Enhanced Sampling Selectively of Ligand in ProteinLigand Complex Qiang Shao1,2,3*, Weiliang Zhu1,2,4 1Drug

Discovery and Design Center, CAS Key Laboratory of Receptor Research, Shanghai Institute

of Materia Medica, Chinese Academy of Sciences, 555 Zuchongzhi Road, Shanghai, 201203, China 2University 3Beijing

of Chinese Academy of Sciences, Beijing, 100049, China

National Laboratory for Molecular Sciences, 1st North Street, Zhongguancun, Beijing,

100080, China 4Open

Studio for Druggability Research of Marine Natural Products, Pilot National Laboratory for

Marine Science and Technology, 1 Wenhai Road, Aoshanwei, Jimo, Qingdao, 266237, China *To

whom correspondence should be addressed. Qiang Shao, Tel: +86 21 50806600-1304, 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 40

Abstract Understanding the protein-ligand binding is of fundamental biological interest and essential for structure-based drug design. The difficulty in capturing the dynamic process, however, poses a great challenge for current experimental and theoretical simulation techniques. A selective integrated-tempering-sampling molecular dynamics (SITSMD) method offering an option for enhanced sampling selectively of the ligand in protein-ligand complex was utilized to quantitatively illuminate the binding of benzamidine to the wild-type trypsin protease and its two mutants (S214E and S214K). The SITSMD simulations could produce consistent results as the extensive conventional MD simulation and gave additional insights into the binding pathway for test protein-ligand complex system using significantly saved computational resource and time, indicating the potential of such a method in investigating protein-ligand binding. Additionally, the simulations identified the different roles of trypsin-benzamidine vdW and electrostatic interaction in the binding: the former interaction works as the driving force for dragging the benzamidine close to the native binding pocket, and the latter interaction mainly contributes in stabilizing the benzamidine inside the pocket. The S214E mutation introduces more favorable electrostatic interactions and as a result both vdW and electrostatic interactions drive the benzamidine binding, lowering the binding and unbinding free energy barrier. In contrast, the S214K mutation prohibits the binding of the benzamidine to the native ligand-binding pocket by introducing disliked charge-charge interactions. In summary, these findings suggest that the change in specific residues could modify the protein druggability, including the binding kinetics and thermodynamics.

2


Page 3 of 40 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


Introduction The binding of ligands, including various endogenic substrates and exogenetic drugs, underpins all biomolecular processes. For decades, binding affinity that quantitatively describes the protein-ligand interaction strength for a ligand to be docked in a target pocket has been often used as an important indicator for drug efficacy, with a high binding affinity implying a better drug efficacy.1-5 Lately, it is realized that compared to binding affinity, the binding kinetics displays a better correlation with drug efficacy and thus should be more predictive for in vivo drug action.6-8 From the landscape point of view, the ligand binding might be simply treated as a two-state process between the unbound and bound states that form the stable basins of attraction.9 While the relative free energy difference between the two states governs the binding affinity, the kinetics of the binding is supposedly determined by the height of the relevant free energy barrier. That is to say, the binding affinity can be solely governed by the endpoint states whereas the binding kinetics is influenced by the detailed pathway connecting the endpoints. Therefore, understanding the complete process of ligand binding transition and underlying mechanism rather than simply predicting the binding affinity data is of fundamental biological interest and essential for structure-based drug design.9-12 The difficulty in measuring the dynamic binding process is, however, noticeable and poses a great challenge for both experimental and theoretical simulation techniques. The development of specialized computer hardware, e.g., Anton supercomputer, allows for the performance of biomolecular simulations with the lengths up to hundreds of microseconds or even milliseconds.13-15 Extensive MD simulations thus have been used to explore the complete binding pathways for multiple systems including the binding of the cancer drug dasatinib or the kinase inhibitor PP1 to Src kinase,16 three antagonists and one agonist to the β2-adrenergic receptor (β2AR).10 Despite

3



Page 4 of 40

these successes, direct observation of binding/unbinding events using brute-force MD simulations on general-purpose computers is still often out of the question. As a matter of fact, while only a few binding/unbinding events can be captured by conventional MD (coupled with the analysis technique of Markov state model, network analysis, or machine learning algorithm et al.),17-21 most binding events have been reported using various enhanced sampling algorithms.9,22 For instance, metadynamics and the variants have been widely used to evaluate the kinetic quantities such as the dissociation rate (koff) and association rate (kon).9,22 The random acceleration molecular dynamics (RAMD) and steered molecular dynamics (SMD),23 memory enhanced random acceleration (MERA) MD and immune algorithm (IA),24 and WExplore technique25 have been developed to explore the unbinding pathways of various protein-ligand systems. In comparison to the sampling challenges in unbinding that focuses on disengaging protein-ligand strong interactions, the sampling challenges in binding correspond to the diffusion and translocation of the ligand through a multitude of possible states before searching out the exact binding pocket for docking. Multiple enhanced sampling methods have been reported for binding simulations, e.g., the modeling employing limited data (MELD) accelerated MD,26,27 Gaussian accelerated molecular dynamics (GaMD),28 et al. Selective integrated-tempering-sampling molecular dynamics (SITSMD) is an enhanced sampling method offering an option for enhanced sampling selectively of predefined object(s) in a simulation system (e.g., often used for protein conformational sampling in a protein-solvent system).29-32 In a ligand binding simulation, SITSMD can be hopefully used to make an enhanced sampling on ligand only to enforce a random walking for the ligand molecule in aqueous solution and on the protein surface to facilitate its searching for the protein binding pocket(s). The remaining part of the system, e.g., the (structurally rigid) protein target and solvent environment, can be treated with a non-enhanced

4


Page 5 of 40 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


conventional MD form to avoid unnecessary sampling enhancement. The subsequent statistical analysis on the trace of the ligand allows for the quantitative exploration of the ligand binding/unbinding pathways on the free energy surface. Such a method should be essentially suitable for offering a computationally efficient service for protein-ligand binding investigation, particularly for structurally rigid protein targets. With this aim, SITSMD simulations were tested in the present study to explore the binding/unbinding pathway for a trypsin-benzamidine complex, a system often studied by various simulation methods.20,21,25,33-36 Additionally, to further evaluate the reliability of SITSMD method in ligand binding pathway exploration and understand the mechanism of trypsin-benzamidine binding, not only wild-type trypsin (TrypsinWT) but also two residue mutants (TrypsinS214E and TrypsinS214K) were simulated to demonstrate the influence of the residue composition of protein target on ligand binding/unbinding. Ser214 is called the fourth member of the catalytic triad on the basis of its conservation throughout the serine protease family, which contributes hydrogen bonding to Asp102, a residue regulating the catalytic power of trypsin (the replacement of Asp102 with Asn or Ala reduces kcat by approximately 104 fold).37 The free energy landscape analysis on the sampled conformational space revealed the detailed binding/unbinding pathways of all test complex systems. The structure characterizations of the binding pathway were compared to the counterparts reported in the previous simulation studies to evaluate the reliability of SITSMD simulation results, indicating the potential of such a method in exploring ligand binding/unbinding pathway and investigating binding/unbinding mechanism.

5



Page 6 of 40

Figure 1. (A) Crystallographic structure of wide-type trypsin-benzamidine complex system (PDB code: 3PTB). The native ligand binding pocket is shown by purple mesh. (B) Structure of benzamidine molecule. (C) Magnified view of the position of Ser214 relative to the native binding pocket (left), and the stabilization of the benzamidine in the pocket by the hydrogen bonding with the side-chain carboxylic group of Asp189 and the backbone carbonyl group of Gly219 (right). Computational Details Trypsin-benzamidine Simulation System Preparation. The atomic coordinates of the complexes under study were retrieved from protein data bank (PDB code: 3PTB (TrypsinWT), 1ANB (TrypsinS214E), 1ANC (TrypsinS214K)).37,38 The Ca2+ cation and all the crystal water in PDB structures were maintained whereas the benzamidine was relocated from the binding pocket to a position remote to protein (>20 Å, as shown in Figure S1). AMBER 14 suite of program39 was employed for simulation in which the underlying force fields are FF99SBildn force field40 for protein and TIP3P model41 for water. The benzamidine was parameterized using generalized Amber force field (GAFF).42 For each complex system, the experimental structure was solvated in a cubic box (detailed size shown in Table S1 in the Supporting Information). Multiple Na+ anions were added to neutralize the charges on protein and ligand and a total of ~11000 water molecules were subsequently added to fill the empty position of the cubic box.

6


Page 7 of 40 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


SITSMD Simulation Details. All simulations have been carried out using the CPU implementation of the parallelized pmemd program in AMBER14. The detailed simulation speed depends on the simulation system size and the number of CPUs used. The modification of the pmemd program in SITSMD doesn’t significantly alter the computational speed of the molecular simulation. Each constructed trypsin-benzamidine complex system was initially minimized for 2500 steps with the protein being fixed using harmonic restraints (force constant = 500.0 kcal mol-1 Å-2). Then the simulation system was heated to 300 K and equilibrated for 20 ns with harmonic restrains applied on protein backbone atoms (force constant = 10.0 kcal mol-1 Å-2). The equilibrated systems with varied coordinates and velocities were used for multiple independent long-time production runs yielding simulation data for statistical analysis. Production runs were performed at constant temperature of 300 K and constant pressure of 1 atm (NPT ensemble), with the potential energy modified via SITSMD method to encourage thorough sampling over a large potential energy range. The SHAKE algorithm43 was used to fix all covalent bonds involving hydrogen atoms and periodic boundary conditions were used to avoid edge effects. The Particle Mesh Ewald method44 was applied to treat long-range electrostatic interactions and the cutoff distance for long-range terms (electrostatic and van der Waals energies) was set as 10.0 Å. The simulation data in each product run was collected every 2.0 ps. The details of SITSMD method have been described in our previous literatures.29-32,45-47 Generally, P-SITSMD method yields the effective potential energy (Eeff) so as to modify the energy distribution in a simulation system through the introduction of a sum-over-temperature non-Boltzmann distribution factor:

E eff  

1 β E ln  n k e k β0 k

(1)

where E is the original potential energy of the system,  0  1 / k BT0 , (kB is the Boltzmann constant and 7



Page 8 of 40

T0 the temperature of the system), 𝛽𝑘’s are a series of temperatures that cover both low and high temperatures around T0, 𝑛𝑘’s are parameters with their values being determined in a preliminary iteration process30 (see the converged 𝑛𝑘 values determined for trypsinWT-benzamidine system as an example in Figure S2). As a result of the proper adjustment of the values of 𝑛𝑘’s, the potential energy is supposed to be sampled evenly in a large desired range corresponding to the integration over a large temperature range  k  given in Eq. 1, leading to a sufficient sampling on potential energy surface. For a protein-ligand complex, to selectively enhance the sampling of the ligand so as to remove the unnecessary enhancement for protein and solvent, the total system potential energy can be separated as follows: E  E L  E P  E Sol  E L - P  E L - Sol  E P  Sol

(2)

where EL: the energy of the ligand, EP: the intra-protein energy, ESol: the energy of solvent (including water molecules and all necessary ions involved),

EL-P: ligand-protein interaction energy, EL-Sol:

ligand-solvent interaction energy, and EP-Sol: protein-solvent interaction energy. The effective potential can be rewritten as the following form:

E eff  E P  E Sol  E P  Sol 

1  β (E  E E ) ln  n k e k L L - P L - Sol β0 k

(3)

The ligand relevant components of the potential energy can be sampled in enhanced manner intentionally while leaving other components sampled with conventional forms. In the present SITSMD simulation, a total of 60 discrete temperatures  k  at evenly spaced intervals in the range of 280K-550 K were utilized in Eq. 2. As a result, while the kinetic energy and temperature of the whole simulation system, and the non-enhanced potential energy part (EP+ESol+EP-Sol) are all maintained in similar ranges as in the conventional MD (Figure S3A-C), the target energy part (EL+EPSol+EL-Sol)

can be sampled more randomly in a much more broadened range (Figure S3D). 8


Page 9 of 40 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 converged calculation of SITSMD simulation yields a biased distribution in the configuration space. The unbiased distribution function at simulation temperature can be recovered as: ρ(r)  ρ eff (r)e

 β 0 (E(r)  E eff (r))

(4)

where ρ eff (r) is measured in SITSMD simulation. Technically, in the simulation trajectories, each snapshot (conformation) is collected with a corresponding potential energy Eeff recorded. While the configuration space of the simulation system is sampled under a biased force which produces a biased distribution function in the configuration space, the visiting probability of each conformation is not unity as that in non-enhanced conventional MD simulation. Therefore, in the statistical analysis after the simulation data is collected, the “real” visiting probability of each snapshot needs to be recovered using the reweighting factor ( W  e

-β 0 (E(r)  E eff (r))

 e

/ n k e

 β 0 (E L  E L - P  E L - W )

 β k (E L  E L - P  E L - W )

). That is

k

to say, for each snapshot, there is a specific visiting probability which is given by the reweighting factor. The reweighted ensemble average of any thermodynamic quantity A can be calculated as A 

AW W

ΔG(X)  

. And the free energy landscape along any reaction coordinate(s) X is calculated by:

1 P(X) ln β Pmax

(5)

where P(X) is the recovered probability distribution along X, which can be recovered from peff (r, x), the probability of observing x with Eeff (r, x): P(x) 



p(r, x)dr 

p

eff

(r, x)  W

W

.

It is noteworthy that for an enhanced sampling simulation with modified energy, the energy noise affect in the implemented reweighting protocol could be a potential flaw for the accurate evaluation of thermodynamic quantities. As shown in Figure S4A-B, the biased energy is exponentially distributed in a relatively narrow range, making the small biased energies dominating the contribution. Therefore, 9



Page 10 of 40

the efficiency lost due to the reweighting should be not severe in the SITSMD simulations. Protein-ligand Binding Free Energy Calculation by MM/PBSA Simulations. The SITSMD was compared with the standard approach of molecular mechanics/Poisson Boltzmann surface area (MM/PBSA) 48 in evaulating the tryspin-benzamidine binding free energy. The simulation systems of MM/PBSA measurement were prepared from the crystal structures of 3PTB, 1ANB, and 1ANC, respectively. The minimization, heating, and constrained equilibrium processes were performed similarly as in SITSMD preprocess, and subsequently the equilibrium simulations were run for ~2 ns at 300 K with conventional MD. For each system, the MM/GBSA energy was calculated using MMPBSA.py.MPI program in AMBER14 on 200 snapshots collected from the equilibrium simulation. The van deer Waals (vdW) and electrostatic interaction energies between trypsin and benzamidine were extracted for data analysis, respectively. Results Reversible Binding and Unbinding Events Observed in Individual Simulation Trajectories. 4~5 independent SITSMD trajectories with varied initial coordinates and velocities were performed for each protein-ligand complex system, giving the accumulated simulation length of ~1 μs/system. The root-mean-square deviation (RMSD) of ligand relative to the crystal structure (ligand RMSD) is utilized as a suitable order parameter to describe the translocation of benzamidine relative to trypsin. As shown in Figure 2, in each single trajectory, under the SISTMD sampling technique, the ligand RMSD can be sampled in a large range. The remaining trajectories can be seen in Figures S5-S7. Additionally, the two-dimensional free energy landscapes as the function of ligand RMSD and the coordinates of benzamidine relative to the pocket in the x/y/z directions, respectively, further indicate that the benzamidine can translocate in a large scope in the aqueous solution and on the protein surface

10


Page 11 of 40 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 S8). Intriguingly, the benzamidine can not only spontaneously reach the experimentally identified (native) binding pocket but also unbind from the pocket repeatedly for wild-type trypsin (Figure 2 and Table 1), yielding enough binding/unbinding events for the quantitative evaluation of the binding/unbinding pathway. The repeated binding and unbinding events can be also seen in the case of S214E mutant but not S214K, reflecting the residue mutation influence on the trypsinbenzamidine binding.

Protein

Ligand

Ntraja

ttrajb (μs)

Nbindingc

Nunbindingc

TrypsinWT

Benzamidine

5

1.2

10

7

TrypsinS214E Benzamidine

5

1.6

9

8

TrypsinS214K Benzamidine

4

1.4

n/a

n/a

Table 1. Parameter summary for SITSMD simulations. aNtraj: the number of trajectories; bttraj: the total simulation time; cNbinding/Nunbinding: the binding/unbinding events observed from the trajectories.

Figure 2. Representative SITSMD simulation trajectories shown by the time-series of ligand RMSD for (A) wild-type trypsin, (B) S214E, and (C) S214K mutants. Cyan colored arrows indicate the 11



Page 12 of 40

occurrence of the binding of benzamidine to the native binding pocket. Estimation of Binding/Unbinding Thermodynamics. After the reweighting approach (Eq. 4), the potential of mean force (PMF) for the binding transition is depicted along the order parameter of ligand RMSD for all three complex systems. One can see from Figure 3 that the ligand binding is a multiplestate transition: As the native binding pocket observed in the crystal structure (e.g., at the ligand RMSD of ~3.5 Å for wild-type trypsin) has the lowest free energy for benzamidine binding, multiple local energy minima also exist corresponding to the trapping of the benzamidine in non-native metastable binding poses. It is noteworthy that previous unbiased MD simulations on trypsin and G-proteincoupled receptor (GPCR) systems also suggested that the ligands enter the ligand-binding pockets by passing through specific metastable intermediate states.10,20 Therefore, the native binding pocket of the wide-type trypsin is thermodynamically favorable for benzamidine binding. The S214E mutated trypsin reduces the binding/unbinding free energy barrier and thus further promotes the binding and unbinding of benzamidine to trypsin native binding pocket (Figure 3). In contrast, the S214K mutation prohibits the benzamidine to bind to the native binding pocket but meanwhile allows it to be trapped at non-native binding poses.

12


Page 13 of 40 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. Potential of mean force (PMF) for the binding/unbinding of benzamidine along the coordinate of ligand RMSD for wild-type trypsin (TrypsinWT), S214E mutant (TrypsinS214E), and S214K mutant (TrypsinS214K), respectively. It is noteworthy that the binding free energy of trypsinWT-benzamidine was experimentally identified to be -6.20 kcal/mol.49 Doudou et al. proposed a practical approach to calculate the standard binding free energy from one-dimensional PMF profile:50 ΔG  ΔG PMF  ΔG V , where ΔG PMF is the free energy change between the bound and unbound states and ΔG V  k B Tln(Vu /V0 ) corresponds to the free energy of taking the ligand from a volume V0 = 1661 Å3 (1M concentration) to the unbound region sampled by the ligand. Using such a method, the trypsinWT-benzamidine binding free energy calculated on the basis of Figure 3 is roughly -5.06 kcal/mol, differing by 1.14 kcal/mol from the experimental data. The same method was also used in the extensive conventional MD simulations of trypsinWT-benzamidine binding,20,33 giving similar simulated binding free energy values as here. In contrast, the MM/PBSA calculation with either FF99SBildn or another often used force filed, FF03,51 produced apparently deviated binding free energy from the experimental data (Table S2). Therefore, the SITSMD could give more reasonable binding free energy than MM/PBSA does. Characterization of Binding/Unbinding Pathway and the Influence from S214 Mutations. Ligand poses corresponding to the local minima in the PMF profiles were structurally identified using the hierarchical clustering analysis by the MMTSB toolset.52 Note that each simulation consists of a single ligand and protein, and ligands are superimposed for visualization. One can see for the binding of benzamidine to the wild-type trypsin (Figure 4A), a wide well with relatively high free energy at ligand RMSD of 35~45 Å is presented in the PMF profile, corresponding to a “diffusion” (D) state of benzamidine in which the ligand molecule diffuses in the bulk with the popular positions remote from

13



Page 14 of 40

the protein surface (Figure S9A). As the benzamidine is proximal to trypsin, it could make the first structural contact with the protein at an “a” pose (Figure 4) where the ligand is stabilized by hydrophobic interactions with Tyr59, Tyr94, and Val90 residues (Figure 5). Such a pose has been also observed in previous MD simulation and Markov mode analysis by Buch et al. (see the defined S1 state in Figure 3B of Ref. 20). In addition, an “c” binding pose observed here is also similar to the S2 state in Ref. 20 which consists of stronger π-π stacking interactions between the benzamidine with Tyr39 and Tyr151 (see Figure 5 here and Figure 3B in Ref. 20). Accordingly, the free energy of the “c” state is significantly lower than that of “a” state. The translocation of the benzamidine between the two states (a→c) passes through two additional poses of “b1” and “b2”. In comparison to the “a” state, the “b1” state maintains the hydrophobic interactions with Tyr59 and Tyr94 but loses the interaction with Val90 and meanwhile, recruits the hydrogen bonding from the backbone carbonyl group of Ser96. The “b2” state, on the other hand, is more biased to the “c” state with the formation of a π-π stacking only between the benzamidine and Tyr39 (Figure 5). Before the benzamidine arrives the final native binding pocket (N state), it needs to undergo another high-free-energy-level “d” intermediate which has favorable hydrophobic interactions from Leu99 and Trp215 and unfavorable charge-charge interaction from the positively charged side-chain of His57. Such a distinct state has not been observed in previous trypsin-benzamidine binding studies. Finally, in the N state, the benzamidine molecule is stabilized by the hydrogen binding from Asp189 and Gly129 (as identified in the crystal structure) and an additional hydrophobic interaction from Val243. Therefore, combining Figures 4 and 5, a lowest free energy binding pathway of benzamidine towards the wild-type trypsin is quantitatively explored as a→b1→b2→c→d→N and the unbinding is reversed, which is accompanied by the continuous formation and breaking of protein-ligand hydrophobic interactions in the early binding stage and the

14


Page 15 of 40 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


formation and breaking of electrostatic interactions later.

Figure 4. (A) Metastable states along the binding/unbinding pathway of the benzamidine towards the wild-type trypsin and (B) the visual inspection of their structure characterization. The benzamidine is represented by licorice mode and colored by free energy level in a red-white-blue (RWB) color scale.

Figure 5. Magnified views displaying the trypsin-benzamidine interactions in all metastable states as shown in the PMF profile for the benzamidine binding the wild-type trypsin (Figure 4). The benzamidine and the protein residues involving in the interactions are represented by licorice mode. The S214E mutation alters not only the binding thermodynamics of benzamidine but also the detailed binding pathway. As shown in Figure 6, diffusing from the bulk (Figure S9B), the benzamidine reaches the first structural contact with the S214E mutated trypsin at either the “b′1” or “b′2” pose. In the former state, the ligand is stabilized by the hydrogen binds from the side-chain 15



Page 16 of 40

hydroxyl group of Ser171 and the backbone carbonyl group of Asp223 and the hydrophobic interaction from Leu185. The latter state is similar to the “b2” state in the wild-type tryspin that is characterized by the π-π stacking between the benzamidine and Tyr39 residue (see the superposition of the two states in Figure 7A). As the benzamidine further approaches to the N state, it needs to undergo an “c′” and subsequently “d′” state. The “c′” state is introduced by the S214E mutation and thus not found in the case of wild-type trypsin. In such a state, the benzamidine is mainly attracted by the favorable chargecharge interaction between its positively charged amidine moiety and the negatively charged carboxylic side-chains of Glu214 and Asp102 as well as the hydrophobic interactions from neighboring Leu99 and Val227 residues. The “d′” state shares the same disliked charge-charge interaction from the His57 as in the “d” state of wild-type trypsin whereas the presence of Glu214 modify the surrounding environment of the benzamidine and introduces more electrostatic interactions and thus reduces the free energy level of such a state compared to the wild-type trypsin (see Figures 5 and 6 and the pose superposition in Figure 7B). Finally, the N state is slightly biased to the experimentally identified position and as a result the hydrogen bond from Gly129 to the benzamidine is replaced by the one from the backbone of Trp215.

16


Page 17 of 40 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 6. Metastable states along the binding/unbinding pathway of the benzamidine towards the S214E mutated trypsin and (B) the visual inspection of their structure characterization. (C) Magnified views displaying the trypsin-benzamidine interactions in all metastable states as shown in the PMF profile for the benzamidine binding the S214E mutated trypsin.

Figure 7. The superposition of the binding poses observed among the wild-type trypsin and its S214E 17



Page 18 of 40

and S214K mutants. (A) The “b2” pose of wild-type trypsin (black) and “b′2” pose of S214E mutant (red). (B) The “d” pose of wild-type trypsin (black) and “d′” pose of S214E mutant (red). (C) The “b′1” pose of S214E mutant (red) and “b′′1” pose of S214K mutant (green). (D) The “c′” pose of S214E mutant (red) and “c′′” pose of S214K mutant (green). Although the S214K mutated trypsin does not attract the benzamidine into the native binding pocket, its non-native binding poses are quite similar to those in S214E (Figures 7C and D). For instance, while the “c′′” pose in S214K is slightly deviated from the “c′” pose in S214E, the “b′1” pose in S214E shares the similar hydrophobic interaction with Leu185 and hydrogen bond with the backbone of Asp223 as at the “b′′1” pose or the similar hydrophobic interaction with Leu185 and hydrogen bond with the backbone of Ser171 as at the “b′′2” pose in S214K (Figure 8). Thus, the S214K mutation adjacent to the native binding pocket might exert more or less similar influence on the remote binding poses as in S214E mutation but meanwhile, introduces disliked charge-charge interaction (S214Kbenzamidine), rejecting the benzamidine to reach the native binding pocket.

Figure 8. Metastable states along the binding/unbinding pathway of the benzamidine towards the S214K mutated trypsin and (B) the visual inspection of their structure characterization. (C) Magnified 18


Page 19 of 40 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


views displaying the trypsin-benzamidine interactions in all metastable states as shown in the PMF profile for the benzamidine binding the S214K mutated trypsin. Protein Conformational Plasticity vs Ligand Binding. Proteins often undergo conformational changes during the ligand binding process.12 Previous extensive MD simulations have observed different binding pocket structure for wild-type trypsin.20,33 The conserved three loops have been used for evaluating the structure of the native binding pocket: residues 188-196, 214-220, and 226-230.33 To see the protein conformational plasticity, we first visualized the structure changes of the entire wildtype trypsin and its local binding pocket loops, respectively, in single trajectory. As shown in Figure S10, the RMSD of the entire protein structure fluctuates within the range of 0 to 2 Å, with the value increasing corresponding to the RMSD increasing of the binding pocket loops. Next, we measured the two-dimensional free energy landscape as the function of ligand RMSD and the pocket loop RMSD. One can see from Figure 9 that the pocket loops keep in the crystal structure with small RMSDs in the binding process but could deviate from the experimental structure (RMSD increasing to ~2.0 Å) as the ligand reaches the native binding pocket. Such a deviated pocket structure is induced by the opening of the loop of residues 214-220. Intriguingly, the same opened pocket structure has been also observed in the previous extensive MD simulation by Plattner and Noe (see the structure in Figure 9 here and Supplementary Figure 7 in Ref. 33). Therefore, the present SITSMD simulation could reflect the conformational plasticity of trypsin protease.

19



Page 20 of 40

Figure 9. Two-dimensional free energy landscapes as the function of ligand RMSD and the protein pocket loop RMSD (residues 188-196, 214-220, and 226-230), along with the representative pocket loop structures in the distinct states. The contours are spaced at intervals of 0.5 kBT. Protein-ligand Interaction Energy along the Binding Pathway. During the trypsin-benzamidine binding process, the benzamidine undergo the interaction from the trypsin and solvent. We analyzed in details the interaction energy from trypsin and solvent to benzamidine (protein/solvent-ligand) and decomposed into van der Waals (vdW) and electrostatic terms. As shown in the two-dimensional free energy landscapes as the function of ligand RMSD and protein/solvent-ligand interaction energy (Figure 10A), the total interaction energy decreases from the very beginning of the binding and drops significantly near the native binding pocket of the wild-type trypsin (see the minimum free energy path (MFEP) achieved by scanning the lowest free energy points of the free energy surface using a sophisticated minimum energy path surface analysis (MEPSA) tool53). Specifically, the decrease of the total interaction energy in the early binding stage is mainly induced by the decrease of vdW interaction energy (Figure 10B) whereas the significant dropping of the total interaction energy 20


Page 21 of 40 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 the late stage is mainly a result of the decease of electrostatic interaction energy (Figure 10C). Besides, we also analyzed the interaction from trypsin to benzamidine (protein-ligand). One can see from Figure 10D that the wild-type trypsin-benzamidine interaction energy is positive during most of the binding process, which is contributed from the favorable vdW interaction energy (Figure 10E) and unfavorable electrostatic interaction energy (Figure 10F). Therefore, the trypsin-benzamidine vdW interaction works as the main driving force for dragging the benzamidine close to the native binding pocket of the wild-type trypsin, and finally the trypsin-benzamidine both vdW and electrostatic interactions stabilize the benzamidine inside the pocket. These two kinds of interactions play different roles in the entire binding pathway. Accordingly, one can see from Figure 5 that the hydrophobic interactions (vdW energy) from the trypsin residues towards benzamidine are throughout the binding pathway whereas the electrostatic interactions like hydrogen binding participate into the protein-ligand interaction mainly in the late stage of the ligand binding process for wild-type trypsin.

Figure 10. (A-C) Two-dimensional free energy landscapes as the function of ligand RMSD and the protein/solvent-ligand interaction energy as well as the vdW or electrostatic component for wild-type trypsin system. (D-F) Two-dimensional free energy landscapes for protein-ligand interaction energy. 21



Page 22 of 40

The contours are spaced at intervals of 0.5 kBT. The purple triangles indicate the minimum free energy paths (MFEPs) calculated by MEPSA tool.53 In the case of S214E mutant, both vdW and electrostatic components of the protein/solvent-ligand interaction energy are decreased along the binding pathway, and such a cooperative decease of two kinds of energies results in the decrease of the total protein-ligand interaction energy (Figure 11A-C). Additionally, for the protein-ligand interaction energy, the detailed value of electrostatic energy is significantly decreased as comparison to the case of wild-type trypsin, and both vdW and electrostatic energies decrease from the beginning of the binding process (Figure 11D-F). As a result, the binding of benzamidine towards S214E mutated trypsin is under favorable vdW and electrostatic energies from trypsin. As shown in Figure 6C, not only the hydrophobic interactions but also electrostatic interactions (including charge-charge interaction and hydrogen binding) are involved throughout the entire binding pathway, working as the driving force for the binding of benzamidine to the S214E mutated trypsin.

Figure 11. (A-C) Two-dimensional free energy landscapes as the function of ligand RMSD and the protein/solvent-ligand interaction energy as well as the vdW or electrostatic component for trypsin S214E mutant system. (D-F) Two-dimensional free energy landscapes for protein-ligand interaction 22


Page 23 of 40 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


energy. The contours are spaced at intervals of 0.5 kBT. Discussion and Conclusions The dynamic binding process of a drug to a target pocket contains the relevant kinetics and thermodynamics information that is correlated to the drug efficacy and should be predictive for in vivo drug action. To visualize the detailed binding transition at high-resolution is still experimentally difficult. On the other hand, while MD simulation is an attractive alternative approach by producing a simultaneous resolution in time, space, and energy, its application in biomolecular systems often remains beyond the assessable computational time-scale and requires enhanced sampling techniques to accelerate the simulation. To this end, a selective integrated-tempering-sampling molecular dynamics (SITSMD) method allowing an option for enhanced sampling selectively of the ligand in protein-ligand complex could be an optimal choice to explore the binding/unbinding transition pathway with economized computational resource and time. In this work, we performed SITSMD simulation to investigate the binding of benzamidine with wild-type trypsin protease and its S214E and S214K mutants. These complex systems were chosen because of two reasons: 1) Trypsin-benzamidine has been often used as a simulation model for the theoretical studies of the ligand binding and/or unbinding, allowing for the evaluation of the present SITSMD simulation results, 2) The studies of the influence of S214 mutation, a residue affecting the catalytic power of trypsin, on the binding/unbinding pathway of benzamidine might shed light on the detailed binding/unbinding mechanism and evaluate the sensitivity of SITSMD in sampling ligand binding/unbinding. The detailed free energy landscape analysis indicates that the trypsin-benzamidine binding is not a two-state transition and metastable intermediate states are involved, consistent with the observations 23



Page 24 of 40

in previous non-enhanced sampling conventional MD simulations.10,20 Specifically, the observed intermediate states could be very similar to those in the previous generalized ensemble simulation and Markov state model with the same force field on the benzamidine but different force field on protein (FF99SBildn vs. FF99SB).20 For instance, the “a” and “c” binding poses here (Figure 5) are the same as the defined S1 and S2 poses in Ref. 20, respectively, for the wild-type trypsin. Meanwhile, the consumed simulation time is extremely reduced (1.2 μs vs. 49.5 μs). Additionally, additional states are observed here to connect the path between “a” to “c” state and the subsequent “c” to the native binding pocket (N) state, yielding a complete pathway of a→b1→b2→c→d→N for the binding of benzamidine to wild-type trypsin and the unbinding is in reverse. In summary, using significantly saved computational resource and time, the SITSMD simulation can produce consistent results as the extensive MD simulation and give additional insights into the binding pathway for test protein-ligand complex system. The S214E mutation still allows the benzamidine to reach the native binding pocket but meanwhile alters the binding/unbinding pathway by passing through different binding poses and accordingly lowers the binding and unbinding free energy barrier. In contrast, although the S214K mutant could share similar binding poses as S214E does, it prohibits the binding of the benzamidine to the native ligand-binding pocket. It is then observed that the trypsin-benzamidine vdW and electrostatic interactions play different roles in the binding pathway: while the vdW interaction works as the driving force for dragging the benzamidine close to the native binding pocket of the wild-type trypsin, the electrostatic interaction stabilizes the benzamidine inside the native binding pocket. The S214E mutation introduces more electrostatic interactions and as a result both vdW and electrostatic interactions are involved throughout the entire binding pathway to drive the benzamidine cooperatively.

24


Page 25 of 40 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


These findings suggest that the change in specific residues could modify the protein druggability, including the binding kinetics and thermodynamics, and demonstrate that the SITSMD is capable of capturing the residue mutation induced change in binding pathway.

Supporting Information Table indicating the simulation parameter summary for the trypsinbenzamidine complex systems under study; front and side views of the initial positions of trypsin and benzamidine in the simulations; the variance of the values of individual terms of nk’s along the preliminary iteration process for a SISTMD simulation; representative trajectories indicating the comparison of the kinetic energy, the temperature, and the non-enhanced and enhanced sampling terms of the potential energy in SITSMD and conventional MD simulations; the variance of the biased potential energy (Eeff-E) and the corresponding distribution in SITSMD; the remaining trajectories for the binding of benzamidine to wild-type trypsin and the two mutants;

two-dimensional free energy

landscapes as the function of ligand RMSD and the coordinates of benzamidine relative to the pocket in the x/y/z directions; table showing the components of the trypsin-benzamidine binding free energy calculated with the protein force fields of FF99SBildn or FF03 by MM/PBSA; popular positions of benzamidine relative to trypsin in the “diffusion” (D) state for wild-type and S214E mutated trypsin; the variance of the RMSD of the entire protein and the RMSD of protein pocket loops in a representative trajectory for wild-type trypsin. This information is available free of charge via the Internet at http://pubs.acs.org.

Notes There are no conflicts to declare.

25



Page 26 of 40

Acknowledgements This work was supported by the National Key Research and Development Program (Grant No. 2016YFA0502301), and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase) under Grant No. U1501501. The simulations were run at TianHe 1 supercomputer in Tianjin and TianHe II supercomputer in Guangzhou, China.

References (1) Wang, L.; Wu, Y. J.; Deng, Y. Q.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; et al. Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by Way of a Modern Free-Energy Calculation Protocol and Force Field. J Am Chem Soc 2015, 137, 2695-2703. (2) Borea, P. A.; Dalpiaz, A.; Varani, K.; Gilli, P.; Gilli, G. Can Thermodynamic Measurements of Receptor Binding Yield Information on Drug Affinity and Efficacy? Biochem Pharmacol 2000, 60, 1549-1556. (3) Baron, R.; McCammon, J. A. Molecular Recognition and Ligand Association. Annu Rev Phys Chem 2013, 64, 151-175. (4) Wang, J.; Shao, Q.; Cossins, B. P.; Shi, J.; Chen, K.; Zhu, W. Thermodynamics Calculation of Protein-Ligand Interactions by QM/MM Polarizable Charge Parameters. J Biomol Struct Dyn 2016, 34, 163-176. (5) Yu, Y.; Wang, J.; Shao, Q.; Shi, J.; Zhu, W. Effects of Drug-Resistant Mutations on the Dynamic Properties of HIV-1 Protease and Inhibition by Amprenavir and Darunavir. Sci Rep 2015, 5, 10517. (6) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drug-Target Residence Time and Its Implications for Lead Optimization. Nat Rev Drug Discov 2006, 5, 730-739. (7) Schuetz, D. A.; de Witte, W. E. A.; Wong, Y. C.; Knasmueller, B.; Richter, L.; Kokh, D. B.; Sadiq, S. K.; Bosma, R.; Nederpelt, I.; Heitman, L. H.; et al. Kinetics for Drug Discovery: An IndustryDriven Effort to Target Drug Residence Time. Drug Discov Today 2017, 22, 896-911. (8) Zhang, R. M.; Monsma, F. Binding Kinetics and Mechanism of Action: Toward the Discovery and Development of Better and Best in Class Drugs. Expert Opin Drug Dis 2010, 5, 1023-1029. (9) Dickson, A.; Tiwary, P.; Vashisth, H. Kinetics of Ligand Binding through Advanced Computational Approaches: A Review. Curr Top Med Chem 2017, 17, 2626-2641. (10)Dror, R. O.; Pan, A. C.; Arlow, D. H.; Borhani, D. W.; Maragakis, P.; Shan, Y. B.; Xu, H. F.; Shaw, D. E. Pathway and Mechanism of Drug Binding to G-Protein-Coupled Receptors. Proc Natl Acad Sci USA 2011, 108, 13118-13123. (11)Pang, X. D.; Zhou, H. X. Rate Constants and Mechanisms of Protein-Ligand Binding. Annu Rev Biophys 2017, 46, 105-130. (12)Skjaerven, L.; Reuter, N.; Martinez, A. Dynamics, Flexibility and Ligand-Induced 26


Page 27 of 40 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


Conformational Changes in Biological Macromolecules: A Computational Approach. Future Med Chem 2011, 3, 2079-2100. (13)Shaw, D. E.; Dror, R. O.; Salmon, J. K.; Grossman, J. P.; Mackenzie, K. M.; Bank, J. A.; Young, C.; Deneroff, M. M.; Batson, B.; Bowers, K. J.; et al. Millisecond-Scale Molecular Dynamics Simulations on Anton. Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis (SC09) ACM, New York, 2009. (14)Sultan, M. M.; Kiss, G.; Pande, V. S. Towards Simple Kinetic Models of Functional Dynamics for a Kinase Subfamily. Nat Chem 2018, 10, 903-909. (15)Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H. F.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annu Rev Biophys 2012, 41, 429-452. (16)Shan, Y. B.; Kim, E. T.; Eastwood, M. P.; Dror, R. O.; Seeliger, M. A.; Shaw, D. E. How Does a Drug Molecule Find Its Target Binding Site? J Am Chem Soc 2011, 133, 9181-9183. (17)Huang, Y. M. M.; Raymundo, M. A. V.; Chen, W.; Chang, C. E. A. Mechanism of the Association Pathways for a Pair of Fast and Slow Binding Ligands of HIV-1 Protease. Biochemistry 2017, 56, 1311-1323. (18)Decherchi, S.; Berteotti, A.; Bottegoni, G.; Rocchia, W.; Cavalli, A. The Ligand Binding Mechanism to Purine Nucleoside Phosphorylase Elucidated via Molecular Dynamics and Machine Learning. Nat Commun 2015, 6, 6155. (19)Huang, D. Z.; Caflisch, A. The Free Energy Landscape of Small Molecule Unbinding. Plos Comput Biol 2011, 7, e1002002. (20)Buch, I.; Giorgino, T.; De Fabritiis, G. Complete Reconstruction of an Enzyme-Inhibitor Binding Process by Molecular Dynamics Simulations. Proc Natl Acad Sci USA 2011, 108, 1018410189. (21)Doerr, S.; De Fabritiis, G. On-the-Fly Learning and Sampling of Ligand Binding by HighThroughput Molecular Simulations. J Chem Theory Comput 2014, 10, 2064-2069. (22)Bruce, N. J.; Ganotra, G. K.; Kokh, D. B.; Sadiq, S. K.; Wadel, R. C. New Approaches for Computing Ligand-Receptor Binding Kinetics. Curr Opin Struct Biol 2018, 49, 1-10. (23)Niu, Y. Z.; Li, S. Y.; Pan, D. B.; Liu, H. X.; Yao, X. J. 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. (24)Rydzewski, J.; Nowak, W. Memetic Algorithms for Ligand Expulsion from Protein Cavities. J Chem Phys 2015, 143, 124101. (25)Dickson, A.; Lotz, S. D. Multiple Ligand Unbinding Pathways and Ligand-Induced Destabilization Revealed by Wexplore. Biophys J 2017, 112, 620-629. (26)Morrone, J. A.; Perez, A.; MacCallum, J.; Dill, K. A. Computed Binding of Peptides to Proteins with Meld-Accelerated Molecular Dynamics. J Chem Theory Comput 2017, 13, 870-876. (27)Morrone, J. A.; Perez, A.; Dene, Q.; Ha, S. N.; Holloway, K.; Sawyer, T. K.; Sherborne, B. S.; Brown, F. K.; Dill, K. A. Molecular Simulations Identify Binding Poses and Approximate Affinities of Stapled α-Helical Peptides to MDM2 and MDMX. J Chem Theory Comput 2017, 13, 863-869. (28)Miao, Y. L.; Huang, Y. M. M.; Walker, R. C.; McCammon, J. A.; Chang, C. E. A. Ligand Binding Pathways and Conformational Transitions of the HIV Protease. Biochemistry 2018, 57, 15331541. (29)Yang, L.; Gao, Y. Q. A Selective Integrated Tempering Method J Chem Phys 2009, 131, 27



Page 28 of 40

214109. (30)Yang, L.; Liu, C.; Shao, Q.; Zhang, J.; Gao, Y. Q. From Thermodynamics to Kinetics: Enhanced Sampling of Rare Events. Acc Chem Res 2015, 48, 947-955. (31)Yang, L. J.; Shao, Q.; Gao, Y. Q. Enhanced Sampling Method in Molecular Simulations. Prog Chem 2012, 24, 1199-1213. (32)Shao, Q.; Xu, Z.; Wang, J.; Shi, J.; Zhu, W. Energetics and Structural Characterization of the "DFG-Flip'' Conformational Transition of B-Raf Kinase: A SITS Molecular Dynamics Study. Phys Chem Chem Phys 2017, 19, 1257-1267. (33)Plattner, N.; Noe, F. Protein Conformational Plasticity and Complex Ligand-Binding Kinetics Explored by Atomistic Simulations and Markov Models. Nat Commun 2015, 6, 7653. (34)Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of Protein-Ligand Unbinding: Predicting Pathways, Rates, and Rate-Limiting Steps. Proc Natl Acad Sci USA 2015, 112, E386-E391. (35)Teo, I.; Mayne, C. G.; Schulten, K.; Lelievre, T. Adaptive Multilevel Splitting Method for Molecular Dynamics Calculation of Benzamidine-Trypsin Dissociation Time. J Chem Theory Comput 2016, 12, 2983-2989. (36)Votapka, L. W.; Jagger, B. R.; Heyneman, A. L.; Amaro, R. E. Seekr: Simulation Enabled Estimation of Kinetic Rates, a Computational Tool to Estimate Molecular Kinetics and Its Application to Trypsin-Benzamidine Binding. J Phys Chem B 2017, 121, 3597-3606. (37)Mcgrath, M. E.; Vasquez, J. R.; Craik, C. S.; Yang, A. S.; Honig, B.; Fletterick, R. J. Perturbing the Polar Environment of Asp102 in Trypsin: Consequences of Replacing Conserved Ser214. Biochemistry 1992, 31, 3059-3064. (38)Marquart, M.; Walter, J.; Deisenhofer, J.; Bode, W.; Huber, R. The Geometry of the Reactive Site and of the Peptide Groups in Trypsin, Trypsinogen and Its Complexes with Inhibitors. Acta Crystallogr B 1983, 39, 480-490. (39)Case, D.A.; Babin, V.; Berryman, J. T.; Betz, R.M.; Cai, Q.; Cerutti, D.S.; Cheatham, T.E. III; Darden, T.A.; Duke, R.E.; Gohlke, H.; et al. AMBER 14; University of California: San Francisco, CA, 2014. (40)Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber FF99SB Protein Force Field. Proteins 2010, 78, 1950-1958. (41)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, 926-935. (42)Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J Comput Chem 2004, 25, 1157-1174. (43)Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical-Integration of Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J Comput Phys 1977, 23, 327-341. (44)Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald - an N.Log(N) Method for Ewald Sums in Large Systems. J Chem Phys 1993, 98, 10089-10092. (45)Shao, Q. Enhanced Conformational Sampling Technique Provides an Energy Landscape View of Large-Scale Protein Conformational Transitions. Phys Chem Chem Phys 2016, 18, 2917029182. (46)Shao, Q.; Yang, L.; Zhu, W. Selective Enhanced Sampling in Dihedral Energy Facilitates 28


Page 29 of 40 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


Overcoming the Dihedral Energy Increase in Protein Folding and Accelerates the Searching for Protein Native Structure. Phys Chem Chem Phys 2019, 21, 10423-10435. (47)Shao, Q.; Zhu, W. Ligand Binding Effects on the Activation of the EGFR Extracellular Domain. Phys Chem Chem Phys 2019, 21, 8141-8151. (48)Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S. H.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; et al. Calculating Structures and Free Energies of Complex Molecules: Combining Molecular Mechanics and Continuum Models. Acc Chem Res 2000, 33, 889-897. (49)Mares-Guia, M.; Shaw, E. Studies on the Active Center of Trypsin. J Biol Chem 1965, 240, 1579-1585. (50)Doudou, S.; Burton, N. A.; Henchman, R. H. Standard Free Energy of Binding from a OneDimensional Potential of Mean Force. J Chem Theory Comput 2009, 5, 909-918. (51)Xu, L.; Sun, H. Y.; Li, Y. Y.; Wang, J. M.; Hou, T. J. Assessing the Performance of MM/PBSA and MM/GBSA Methods. 3. The Impact of Force Fields and Ligand Charge Models. J Phys Chem B 2013, 117, 8408-8421. (52)Feig, M.; Karanicolas, J.; Brooks, C. L. MMTSB Tool Set: Enhanced Sampling and Multiscale Modeling Methods for Applications in Structural Biology. J Mol Graph Model 2004, 22, 377-395. (53)Marcos-Alcalde, I.; Setoain, J.; Mendieta-Moreno, J. I.; Mendieta, J.; Gomez-Puertas, P. MEPSA: Minimum Energy Pathway Analysis for Energy Landscapes. Bioinformatics 2015, 31, 38533855.

TOC Graphic

29



Figure 1. (A) Crystallographic structure of wide-type trypsin-benzamidine complex system (PDB code: 3PTB). The native ligand binding pocket is shown by purple mesh. (B) Structure of benzamidine molecule. (C) Magnified view of the position of Ser214 relative to the native binding pocket (left), and the stabilization of the benzamidine in the pocket by the hydrogen bonding with the side-chain carboxylic group of Asp189 and the backbone carbonyl group of Gly219 (right). 104x60mm (300 x 300 DPI)


Page 30 of 40

Page 31 of 40 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 2. Representative SITSMD simulation trajectories shown by the time-series of ligand RMSD for (A) wild-type trypsin, (B) S214E, and (C) S214K mutants. Cyan colored arrows indicate the occurrence of the binding of benzamidine to the native binding pocket. 124x95mm (300 x 300 DPI)



Figure 3. Potential of mean force (PMF) for the binding/unbinding of benzamidine along the coordinate of ligand RMSD for wild-type trypsin (TrypsinWT), S214E mutant (TrypsinS214E), and S214K mutant (TrypsinS214K), respectively. 124x87mm (300 x 300 DPI)


Page 32 of 40

Page 33 of 40 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 4. (A) Metastable states along the binding/unbinding pathway of the benzamidine towards the wildtype trypsin and (B) the visual inspection of their structure characterization. The benzamidine is represented by licorice mode and colored by free energy level in a red-white-blue (RWB) color scale. 124x67mm (300 x 300 DPI)



Figure 5. Magnified views displaying the trypsin-benzamidine interactions in all metastable states as shown in the PMF profile for the benzamidine binding the wild-type trypsin (Figure 4). The benzamidine and the protein residues involving in the interactions are represented by licorice mode. 124x88mm (300 x 300 DPI)


Page 34 of 40

Page 35 of 40 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 6. Metastable states along the binding/unbinding pathway of the benzamidine towards the S214E mutated trypsin and (B) the visual inspection of their structure characterization. (C) Magnified views displaying the trypsin-benzamidine interactions in all metastable states as shown in the PMF profile for the benzamidine binding the S214E mutated trypsin. 124x97mm (300 x 300 DPI)



Figure 7. The superposition of the binding poses observed among the wild-type trypsin and its S214E and S214K mutants. (A) The “b2” pose of wild-type trypsin (black) and “b′2” pose of S214E mutant (red). (B) The “d” pose of wild-type trypsin (black) and “d′” pose of S214E mutant (red). (C) The “b′1” pose of S214E mutant (red) and “b′′1” pose of S214K mutant (green). (D) The “c′” pose of S214E mutant (red) and “c′′” pose of S214K mutant (green). 104x123mm (300 x 300 DPI)


Page 36 of 40

Page 37 of 40 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 8. Metastable states along the binding/unbinding pathway of the benzamidine towards the S214K mutated trypsin and (B) the visual inspection of their structure characterization. (C) Magnified views displaying the trypsin-benzamidine interactions in all metastable states as shown in the PMF profile for the benzamidine binding the S214K mutated trypsin. 124x105mm (300 x 300 DPI)



Figure 9. Two-dimensional free energy landscapes as the function of ligand RMSD and the protein pocket loop RMSD (residues 188-196, 214-220, and 226-230), along with the representative pocket loop structures in the distinct states. The contours are spaced at intervals of 0.5 kBT. 124x108mm (300 x 300 DPI)


Page 38 of 40

Page 39 of 40 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 10. (A-C) Two-dimensional free energy landscapes as the function of ligand RMSD and the protein/solvent-ligand interaction energy as well as the vdW or electrostatic component for wild-type trypsin system. (D-F) Two-dimensional free energy landscapes for protein-ligand interaction energy. The contours are spaced at intervals of 0.5 kBT. The purple triangles indicate the minimum free energy paths (MFEPs) calculated by MEPSA tool.49 145x97mm (300 x 300 DPI)



Figure 11. (A-C) Two-dimensional free energy landscapes as the function of ligand RMSD and the protein/solvent-ligand interaction energy as well as the vdW or electrostatic component for trypsin S214E mutant system. (D-F) Two-dimensional free energy landscapes for protein-ligand interaction energy. The contours are spaced at intervals of 0.5 kBT. 145x96mm (300 x 300 DPI)


Page 40 of 40

Unbinding Pathway by Enhanced Sampling

Recommend Documents