Molecular Determinants and Bottlenecks in the Dissociation Dynamics

Subscriber access provided by READING UNIV

Article

Molecular Determinants and Bottlenecks in the Dissociation Dynamics of Biotin-Streptavidin Pratyush Tiwary J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b09510 • Publication Date (Web): 08 Nov 2017 Downloaded from http://pubs.acs.org on November 9, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry B 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 27

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

Molecular Determinants and Bottlenecks in the Dissociation Dynamics of Biotin-Streptavidin Pratyush Tiwary∗ Department of Chemistry and Biochemistry and Institute for Physical Science and Technology, University of Maryland, College Park 20742, USA. E-mail: [email protected]

Abstract Biotin-streptavidin is a very popular system used to gain insight into protein-ligand interactions. In its tetrameric form, it is well-known for its exceptionally high kinetic stability, being one of the strongest known non-covalent interactions in nature, and is heavily used across the biotechnological industry. In this work we gain understanding of the molecular determinants and bottlenecks in the dissociation of the dimeric biotinstreptavidin system in wild type and with a point mutation. Using recently proposed enhanced sampling methods with full atomistic resolution, we reproduce the experimentally reported effect of the mutation on the dissociation rate. We also answer a longstanding question regarding cause/effect in the coupled events of bond stretching and bond hydration during dissociation and establish that in this system, it is the bond stretching and not hydration which forms the bottleneck in the early parts of the dissociation process. We believe these calculations represent a step forward in the use of atomistic simulations to study pharmacokinetics. An improved understanding of biotinstreptavidin dissociation dynamics should also have direct benefits in biotechnological

1

ACS Paragon Plus Environment


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

and nanobiotechnological applications.

Introduction The interaction of streptavidin with biotin in the homotetramer form is one of the strongest known non-covalent interactions in nature. 1–4 This system forms a central component of countless biological, biotechnological and nano-biotechnological applications, and has long been used as a model system to understand the molecular basis of high affinity binding. 5,6 Designing specific mutations in streptavidin allows tailoring a wide range of affinities and rate constants 7 as needed for diverse practical applications. Further, given that many of the structural motifs found in biotin-streptavidin are utilized in numerous other protein-ligand interactions, it is very desirable to map out the structural and dynamical motifs underlying the slow dissociation kinetics in this system. Experiments have reported both the binding affinity (KD ) and the association/dissociation rate constants (kon , kof f respectively, with KD = kof f /kon ) for this system with and without mutations in the protein. Of these numbers, the dissociation rate constant kof f (also often represented through its inverse, the residence time) has recently received a lot of attention as possibly being a much better predictor of drug efficacy compared to KD . 8–10 While experiments can provide a direct measurement of equilibrium dissociation constant kof f , it is not so easy to derive from experiments direct insight into the molecular determinants of the residence time. These determinants could come from a diverse range – protein-ligand contacts, solvation, protein/ligand flexibilities and many others. Armed with an atomistic resolution understanding of the role played by these determinants in the dissociation dynamics, one could then propose structural modifications that would assist in rational design of drugs with desired residence times. With the advent of petaflop computing and reliable forcefields, all-atom molecular dynamics (MD) simulations seem ideally poised to gain atomic resolution insight into the process of ligand dissociation. However, the reported residence

2


Page 2 of 27

Page 3 of 27


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 ACS Paragon Plus Environment


50

1

χ

χ

trial

trial

χoptimized

χ 40

30

optimized

0.8

2

1

0.6

4

3

e−λi

F(χ) in kJ/mol

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 4 of 27

20

0.4

10

0.2 (b)

(a) 0 0

0.2

0.4

0.6

0.8

0 0

1

2

χ

4

6 8 process index i

10

12

Figure 2: (a) Approximate free energies in kJ/mol as function of the trial and optimized RC (red and black respectively) for WT. To help with comparison, in this plot both RCs have been scaled by their respective maximum recorded values (reported in SR). 4 states are visible in this free energy profile, where state 1 is bound and state 4 is unbound. The other two are metastable intermediates. The very high free energy in the end of state 4 is because the runs were stopped when the trajectories reached state 4 and committed to it - hence this state is not fully sampled. (b) Eigenvalue spectrum corresponding to trial and optimized RC (red circles and black asterisks respectively) for WT. Consistent with the 3 barriers seen in (a) for the optimized RC, we mark the spectral gap as the difference between the eigenvalues for 4th and 5th eigenmodes (since the first eigenmode with λ0 = 0 is just the stationary state). It is clear how the SGOOP-identified RC has sharper and higher barriers in (a), and a much larger spectral gap in (b). Analogous plots for free energies and eigenvalue spectrum for the mutant N23A have been provided in the SR.

4


Page 5 of 27

Spectral Gap 1.6

1

1.4 0.8 1.2 0.6 5

1

c

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


0.8

0.4

0.6 0.2 0.4 0 0

0.2

0.4

c2

0.6

0.8

1

Figure 3: Here we show the relative contributions of D128-biotin bond stretching (c2 ) and D128 hydration (c5 ) to the spectral gap of the dissociation dynamics when viewed as function of different trial RC χ = c1 Ψ1 + c2 Ψ2 + c5 Ψ5 setting c1 = 1. Since Ψ3 and Ψ4 were found not to be a part of the optimized RC, their corresponding weights c3 and c4 are both kept 0 for this illustration. The full dictionary {Ψi } of collective variables is defined in Table 1. i=5 P Similar spectral gap profiles for different χ = ci Ψi for WT and for N23A are provided i=1

in the SR. The spectral gap is in arbitrary units and has been scaled so that the RC with χ = {1, 0, 0, 0, 0}, which was the trial RC, has a spectral gap of 1 unit. The optimized RC with highest spectral gap is χ = Ψ1 + 0.75Ψ2 . is about 6 orders of magnitude lower than for the full tetrameric system. 21 With recently proposed enhanced sampling methods 11,12 we study the biotin-streptavidin system (PDB ID: 3ry2) in its wild-type (WT) and with an asparagine-alanine mutation (N23A). 22 We expect and find that considering the dimer instead of tetramer brings down the residence time from more than a day to about a second, which is a timescale regime where the enhanced sampling methods of this work have been shown to be reliable. 14–16 Note that there is no fundamental reason why these methods should not work for much longer timescales, and this will be the subject of future investigations. For the tetrameric biotin-streptavidin system the WT and N23A have been reported to have kof f values of 4.4 ± 0.3 × 10−5 s−1 and 1030 ± 220 × 10−5 s−1 respectively. 22 Thus, a single point mutation diminishes the residence time by a staggering nearly 250 times. Our simulations reproduce this behavior in the dimeric biotin-streptavidin and give crucial insight into its molecular origins. Our calculated residence times are also in the rough order of magnitude range expected from the reported difference in binding affinity 5



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

for the tetrameric versus dimeric systems. 21 Furthermore, we settle a long-standing question about the cause versus effect nature of two events, namely water movement into the binding pocket and lengthening of a specific biotin-streptavidin hydrogen bonding interaction (D128-biotin bond). Previous simulations have demonstrated that these two events occur in a coupled manner before dissociation can proceed. 2,3 However the authors could not comment on the relative importance of the two events. In this work, using the recently proposed method “Spectral gap optimization of order parameters (SGOOP)" proposed by Tiwary and Berne, 11 we optimize the reaction coordinate (RC) for biotin dissociation which quantifies the relative importance of these two and other molecular determinants in dissociation. We establish that it is the hydrogen bond stretching that forms the bottleneck for the dissociation to proceed. Using the optimized RC for the WT and N23A forms, we then perform multiple independent “infrequent metadynamics" simulations. 12,13 These give direct estimates of the various dissociation timescales and pathways. All kinetic observables are validated through statistical self-consistency tests. Table 1: One-dimensional RC was constructed through SGOOP as a linear combination of the following 5 collective variables. See relevant residues marked in Fig. 4. Label Ψ1 Ψ2 Ψ3 Ψ4 Ψ5

Collective variable (CV) Distance between centre of mass of S45:OG and N49:N from centre of mass of Biotin:C11 and Biotin:N2 Distance between D128:CG and Biotin:N1 Hydration state of S45:OG Hydration state of N49:N Hydration state of D128:CG

METHODS System specifications The initial structure of biotin bound to streptavidin in its dimeric form was obtained from protein data bank (PDB:3RY2). Amber ff99SB*-ILDN all-atom force-field and TIP4P water 6


Page 6 of 27

Page 7 of 27

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


model were used. 23,24 All the non-protein molecules were removed, and protein preparation wizard within Schrodinger’s software package Maestro was used to post process the structures of protein for simulation. The residue mutation N23A was also carried out using Maestro. Protonation states of the amino acid residues were assigned assuming the systems are at pH 7.0. The ligand charges were parametrized using the general amber forcefield and AM1-BCC charges using the software package Chimera. 25,26 The protein and ligand structures were concatenated and then solvated with TIP4P water. Counterions in the form of sodium chloride (NaCl) salt were added at a concentration of 150 mM to maintain physiological conditions. The full systems, which comprised of 82,793 atoms for each of WT and N23A, including protein, solvent and charge-compensating counter-ions, were prepared and equilibrated. The full structure’s energy was minimized, and then further equilibration MD simulations was run first in the NVT (at temperature 300 K) and then in the NPT ensembles (at pressure 1 bar and temperature 300 K) for 1 ns each, keeping the position of ligand restrained in its binding pose. Nose Hoover Thermostat 27 and Parinello-Rahman barostat 28 were employed. Periodic boundary conditions were applied in all three directions. The equations of motion were evolved with a time step of 2 fs. All our subsequent simulations were performed using constant temperature and pressure ensemble (with isotropic pressure) at temperature T = 300 K using the software GROMACS 5.0.7 patched with the plugin PLUMED 2.1. 29,30

Infrequent Metadynamics for kinetics In this work, we use the infrequent metadynamics approach 12,31 to obtain dissociation trajectories and kinetic rate constants. This involves periodic but infrequent biasing of a lowdimensional RC in order to increase the probability of escape from metastable states where the system would ordinarily be trapped for extended periods of time. 32 If the RC can demarcate all relevant stable states of interest, and if the time interval between biasing events is infrequent compared to the transition path time spent in the ephemeral transition state (TS)

7



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 8 of 27

regions, then one increases the likelihood of not adding bias in the TS regions and thereby keeping pristine the dynamics during barrier crossing. This in turn preserves the sequence of transitions between stable states that the unbiased trajectory would have taken. 12 The acceleration α of transition rates through biasing which directly yields the true unbiased rates, is then through generalized transition state theory: 12,33–35 α = heβV (s,t) it

(1)

Here s is the RC being biased, β = 1/kB T is the inverse temperature, V (s, t) is the bias experienced at time t and the subscript t indicates averaging under the time-dependent potential. One can verify a posteriori if the requirements of Ref. 12,31 were satisfied by checking if the cumulative distribution function for the transition times is Poisson. 31 Note that all rate constants are calculated through the use of the accelerated time, 12,34 without having to converge the free energy. Here we use SGOOP 11,36 for the construction of a desirable RC (next section). Total 5 collective variables (CVs) were used in this work (Table 1) and an optimized 1-dimensional RC χ was constructed as a linear combination of these 5 CVS using SGOOP (next section). As can be seen in Table 1, out of these 5 CVs 3 were distance variables d between the ligand and different parts of the protein. The other 3 were hydration CVs proportional to the number of water molecules in the vicinity of specific protein residues. Since the CVs need to be smoothly differentiable in order to be biased, the hydration variable w was implemented through the function

w=

XX

f (ri,j )

(2)

i∈A j∈B

Here the set A denotes the respective protein residues from Table 1. The set B comprises all the oxygen atoms of various water molecules. ri,j is the distance in nm between atoms i

8


Page 9 of 27

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


and j coming from these two sets respectively. The function

f (r) ≡

(1 − (r/0.6)6 ) (1 − (r/0.6)10 )

(3)

is a switching function that makes w smoothly differentiable. As described in the following section and in SI, the optimized RC was found through SGOOP to be χ = d1 + 0.75d3 . The bias was added once every 5 picoseconds as a Gaussian function of χ: V (χ, t) = h0 e

V (χt ) B ∆T

−k

e−(χ−χt )

2 /2σ 2

. Here χt is the value of the RC at a certain

time t, h0 is the initial Gaussian amplitude which we took as 2 kJ/mol. ∆T , which we took to be ∆T = 14 T is the so-called tempering parameter which smoothly modulates the amplitude of the Gaussian each time a certain point in the RC space is revisited. 32,37 The Gaussian width σ was taken as 0.02 units. The free energy as a function of χ, defined below, can be estimated from the bias as (within an additive constant 38 ): F (χ) = −kB T logP (χ) = −

T + ∆T V (χ) ∆T

(4)

where P (χ) is the stationary probability density as a function of χ. We define the unbound state as states with χ > 3, beyond which we find that the ligand diffuses in the bulk. All 14 independent simulations for both WT and N23A were ran until they gave full dissociation. This took between 2 nanosecond of unscaled MD time for the fastest exit event to as much as 11 nanoseconds for the slowest event. Through each of the simulations we found no deterioration of the protein structure as measured by the RMSD value of the heavy atoms, which stayed within 3 of the native pose.

SGOOP for reaction coordinate optimization Spectral gap optimization of order parameters (SGOOP) 11,36 is a method to optimize a RC as a combination (linear or non-linear) of a set of many given candidate collective variables 9



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

Ψ = (Ψ1 , Ψ2 , ..., Ψd ). SGOOP considers the optimal RC to have the following characteristics: (1) it can demarcate between various metastable states that the full system possesses, and (2) it can maximize separation of timescales between visible slow and hidden fast processes. This timescale separation is calculated as the spectral gap between the slow and fast eigenvalues of the transition probability matrix on a grid along a trial RC. 11,36 The transition probability matrix is calculated in SGOOP using an approximate kinetic model that can be derived for example through the principle of Maximum Caliber. 36,39,40 Let {λ} denote the set of eigenvalues of this matrix, with λ0 = 0 < λ1 ≤ λ2 ≤ .... We refer to e−λn − e−λn+1 as the spectral gap, where n is the number of barriers apparent from the free energy estimate that are higher than a user-defined threshold (typically ≥ kB T ). Note that this calculation is done without any pre-knowledge of the number of metastable states the system might have, and as such makes it easy to use SGOOP for systems with an arbitrary number of metastable states as in this work. The key input to SGOOP is an estimate of the stationary probability density P (χ) (or equivalently the free energy) of the system, accumulated through a biased simulation performed along a sub-optimal trial RC χ given by some linear or non-linear function χ = f0 (Ψ), where Ψ denotes the larger set of candidate CVs. This is done by performing 5 independent frequent metadynamics runs (biasing interval = 2 ps) 32,38 and using Eq. 4 to average over the estimate from these runs. The bias deposition rate is higher for these preliminary runs since the objective in this part is to sample the free energy landscape and not necessarily preserve kinetic pathways. Metadynamics 38 allows to easily reweight for the stationary probability density for any RC in post-processing without having to repeat simulations. This stationary density estimate is complemented with short unbiased MD runs (10 ns long) in the bound pose, which provide the dynamical observables for SGOOP, namely the average number of transitions per unit time along any trial RC (see SR). These serve as a proxy for the explicit knowledge of the diffusivity tensor. Given these pieces of information we use the principle of Maximum Caliber 36,39,40 to set up an unbiased master equation for the

10


Page 10 of 27

Page 11 of 27

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


dynamics of various trial CVs f (Ψ). Through a post-processing optimization procedure we then find the optimal RC as the f (Ψ) which gives the maximal spectral gap of the associated transfer matrix, and hence the optimized RC. We refer to Ref. 36 and to SI for further details of the theory behind SGOOP and its application in this work.

Results We start with a trial RC (Ψ1 in Table 1) for both WT and N23A forms, with which we perform preliminary metadynamics dissociation simulations to explore the free energy. Given the rough nature of the RC, these runs are only for preliminary exploration of the configuration space. We also perform short unbiased MD runs in the bound pose to calculate dynamical observables for the Maximum Caliber framework 39,40 needed by SGOOP (Methods). These metadynamics and unbiased MD runs are input into SGOOP which processes them to optimize the RC as a linear combination of 5 order parameters (Table 1, Fig. 1(b)). We chose these particular 5 order parameters to compose the dictionary from which the RC is optimized because (a) they capture the stabilizing interactions in the bound pose, and (b) we were specifically interested in the mechanistic hypothesis of solvation versus bond stretch. A different dictionary could very well have been used in SGOOP. Note that SGOOP as implemented in this work through the Maximum Caliber framework, takes into explicit account the diffusion anisotropy between different collective variables. 36 Roughly speaking, this anisotropy captures the relative rates at which different collective degrees of freedom relax in response to change in the environment (see Ref. 36 for more details). The importance of taking such anisotropy into account while studying ligand unbinding, and the profound effects it can have on reaction pathways has been highlighted in several papers including but not limited to Ref. 41–44 A second set of “infrequent metadynamics" runs is performed using this optimized RC, with much slower bias deposition frequency in order to avoid perturbing the transition states

11



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: State-to-state transitions rates and various residence times for WT biotinstreptavidin dimer.

12


Page 12 of 27

Page 13 of 27




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 14 of 27

and meet the requirements of Ref. 12,31 (Methods). 14 such independent runs each were performed for the WT and N23A forms. All sets of metadynamics runs (i.e. for input into SGOOP as well as infrequent metadynamics) were started in the bound x-ray pose and stopped when the ligand was fully unbound and freely diffusing. Both sets of dissociation times pass the Kolmogorov-Smirnoff test of Ref. 31 as well as other statistical analyses reported in this section and in Supporting Information (SI). Our calculated kof f values are 2.1 ± 1.5 s−1 and 218 ± 82 s−1 respectively for the WT and N23A - thus we find the mutant to be about 100 times faster at dissociation than the wild type protein. The error bars here and elsewhere in this work denote standard errors corresponding to 95% confidence intervals. Note that these are in rough agreement with the kof f values for the full tetrameric system reported in the introduction taking into account that the dimeric system has 5-6 order of magnitude smaller binding affinity KD , and that if one assumes kon to be diffusion limited, much of the variation in KD = kof f /kon should come from kof f . We would however like to point out that we make such an assumption with great reservation, and that often the measured association constant can be significantly slower than diffusion-controlled values. We make this assumption here only for the purpose of ascertaining if the dissociation timescales of the dimeric system are in the rough ballpark relative to the tetrameric system which has been studied experimentally.

Dissociation dynamics for WT and N23A can be mapped onto a one-dimensional reaction coordinate Several groups over the decades have asked the question: does dissociation of biotin-streptavidin proceeds through well-defined ligand exit pathways that can be described by a reactioncoordinate (RC) model? 2,3 For this system, the experimental evidence is in favor of this statement, given well-defined exponential statistics for the ligand residence times. 3 Here, we explicitly construct such a RC using the method SGOOP, which allows optimizing a lowdimensional RC from a dictionary of possible collective variables or order parameters. 11 The 14


Page 15 of 27

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


key idea in SGOOP is that the best RC will show a higher gap between the timescales associated with slow and fast processes, referred to in the literature as spectral gap (Methods). The dictionary of 5 such collective variables used for the systems here is detailed in Table 1, which comprises various protein-ligand contacts and hydration states of specific protein residues. For the dimeric biotin-streptavidin considered here, we find that both the WT and N23A dissociation dynamics can be reliably mapped into a 1-d RC picture (Fig. 2) with very similar optimized RCs for WT and N23A. By reliability here we allude to the notion introduced in Ref., 11 namely that by optimizing the spectral gap of the RC we have reduced the likelihood of the existence of orthogonal slow degrees of freedom relevant to the dissociation process. Fig. 2(a) shows the approximate 1-d free energy profile as a function of the trial RC used by SGOOP (see Methods) and as a function of the optimized RC. Fig. 2(b) shows that the spectral gap is indeed larger for the optimized RC. Corresponding results for N23A are given in SI. Note that these free energies were generated using a sub-optimal RC and can not be expected to be fully converged. We have shown in the past that SGOOP works remarkably well with such sub-optimal free energies. 18 Another question which has received much attention concerns the nature of the coupling between protein hydration and the lengthening of a specific biotin hydrogen bonding interaction. Stayton et al and others have demonstrated through simulations that an early mechanistic event in biotin-streptavidin dissociation involves these two events (viz. D128 hydration and D128-biotin bond stretching) occuring in a concerted manner. 2,3 Here, we find through SGOOP that the optimized RC involves contribution from the D128-biotin bond stretching but not from the hydration state of D128 (Fig. 3), which is unequivocal evidence that, at least for this system, it is the D128-biotin bond stretching which contributes to the bottleneck. Of course this picture will likely change for other systems - what is more important here is the demonstrated ability to answer such questions with quantitative confidence. In SI, we give similar spectral gap plots for the full exploration of the 5-dimensional order parameter space for WT and N23A. Note that our findings do not imply that the D128

15



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

hydration is not important to the unbinding process. It is however slaved to the D128-biotin bond stretching. Thus, driving the D128 hydration is not likely to lead to D128-biotin bond stretching, while driving the D128-biotin bond stretching is likely to cause D128 hydration. This is what we find in our simulations as reported here using the optimized RC.

Dissociation dynamics for WT and N23A can be mapped into 4 states and sequential movement between them When viewed as a function of the optimized reaction coordinate, we find that the free energy landscape Fig. 2 for both the WT and N23A maps into 4 states, namely bound, unbound and 2 intermediates. Through infrequent metadynamics runs with the optimized RC (Methods) we find these to have diverse lifetimes across several order of magnitudes. The representative structures illustrated in Fig. 4 were picked from various unbinding trajectories as per the protocol described in SI. Undoubtedly there will be many structures in each state that differ along orthogonal degrees of freedom not captured by our optimized RC. Here we do not report them due to the ansatz in this work as in Ref. 11 that the only configurations relevant to the dissociation pathway are the ones that can be distinguished using the optimized RC. Here we describe in details the states for the WT (Fig. 4 and Fig. S1 for WT and N23A respectively) - analogous description for N23A is given in the SI. State 1 is the bound state with a lifetime of 0.45±0.32 sec. Some of the distinguishing and stabilizing interactions here include direct hydrogen bonds between biotin and the residues D128, N23 and S45. To leave state 1, the interactions of biotin with D128 and with N23 distort from being direct to water-mediated, and are finally broken. State 2 is a key intermediate state, with a lifetime of 11±8 msec. The characterizing interaction here is a water mediated bond between biotin and S45, and interactions between biotin and N49. Once the S45 interaction is broken, the system reaches state 3 which is the last intermediate state and has a lifetime of 0.1±0.02 µsec. Here the biotin interacts mainly with N49 and adjoining residues on the periphery of the binding pocket. This state is composed of several shallow energy basins with fast 16


Page 16 of 27

Page 17 of 27

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


interconversion. The biotin is almost fully solvent exposed now and finally it reaches state 4, which is the freely diffusing unbound state.

N23A dissociation has the same rate-determining step as WT, only 100 times faster We extracted the residence times in various states and rate constants for movements between them by pooling the various infrequent metadynamics runs (Methods). This gives a matrix of transition rates and a discrete-time master equation for movement between the various states. The eigenvalues and eigenvectors of this matrix carry precise and useful information about the slow steps in dissociation and their time scales. We calculate this transition matrix and its eigenvalues for different values of the minimum time the trajectory must spend in a certain state before being labelled as being committed to that state. Such a metric, introduced previously by Tiwary, Mondal and Berne, 14 is inspired by the notion of friction-induced spurious recrossing of the barrier between two states, and attempts to discount such recrossings from genuine transition events. This is akin to the lag time in markov models, 45 and the dynamics can be trusted only for those eigenvalues that converge beyond a certain minimum commitment time. As can be seen in Figs. 5(a-b) the eigenvalues for both the WT and N23A are well-converged with respect to this metric. Figs. 5(c-d) show the corresponding eigenvectors. Of special importance is the slowest eigenvalue and eigenvector - which corresponds to the dynamical bottleneck in dissociation of biotin-streptavidin. Interestingly, both for WT and N23A the exit from state 1 forms the bottleneck, as can be seen from the node in the eigenvector i.e. movement happens between states with maximally positive and negative eigenvectors. For N23A this step is around 100 times faster than for WT, which is explained by the missing N23-biotin interaction in the mutant. Once state 1 is exited, the rest of the dissociation process is similar for both forms. Note that the eigenvectors for λ2 and for λ3 in Fig. 5 are flipped in sign with respect to each other when comparing the WT and the N23A forms, but keeping approximately the same 17



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

magnitudes. This is because given an eigenvector v, its negative −v is also an equally valid solution with the same eigenvalue and we can not distinguish between these two. In Fig. 5, it is thus the location of the node and the magnitude of the eigenvector in respective states that carries relevant physical information.

Discussion Biotin-streptavidin is a classic system used to gain insight into protein-ligand interactions. It is heavily used across the biotechnological and nanobiotechnological industry, and is wellknown for its extremely long biotin residence times. In this work we studied the dimeric biotin-streptavidin system, which is slightly more tractable than the full tetrameric system, in its wild type and N23A mutated forms. Using recently proposed enhanced sampling methods 11,12 we demonstrated that atomistic simulations could reproduce the experimentally reported variation in residence time. Our results are significant because as of yet there are only a few examples of simulations that can capture mechanistic details of ligand release all of which use state-of-the-art simulation methods. 13–15,46–48 Our simulations go far beyond just being able to calculate residence times - we also identify the various molecular determinants and bottlenecks in the dissociation process and quantify their roles. One of these is the longstanding question regarding cause/effect in the coupled events of D128-biotin and bond stretching and water entry into the binding pocket to hydrate D128. We establish that at least in the current dimeric system (which itself has been the subject of previous investigations 21 ), it is the bond stretching that forms the bottleneck in the early parts of the dissociation. Thus, we can speculate on the basis of our findings that the effect on residence time would be more profound if one was to restrict/ease the movement of D128, compared to easing the process of hydration of that residue. Our reported residence times for the dimeric system are in rough agreement with what one expects after taking into account the difference in binding affinities between the dimeric and the tetrameric systems.

18


Page 18 of 27

Page 19 of 27

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 this work we were able to map the dynamics on a one-dimensional RC as a combination of given collective variables. This is not entirely surprising given that the experimentally reported residence times for the tetrameric biotin-streptavidin system follow a mono-exponential distribution. However, in more complicated systems, there is no guarantee that a one-dimensional RC exists or can be found. We are however hopeful that the methods and statistical tests used in this work will signal in a self-consistent manner if either the dictionary of trial collective variables needs to be expanded, or if the one-dimensional RC framework is simply not sufficient. We have in fact found this to be the case for dissociation in still more complex systems such as some of those in Ref. 49 Encouraged by the success in this system, we are now in the process of looking at the full tetrameric system with a larger number of mutations. One crucial difference between these two systems is that we anticipate the need to incorporate movement between the two dimeric units as an additional member of the dictionary of collective variables from which the RC is optimized. Even though Grubmüller et al comment in their original paper using a biotin-streptavidin monomer instead of tetramer 5 that they do not expect this simplification to change the results, we wonder if we will see surprises when we simulate the full tetramer. Specifically, given the considerable protein fluctuations in specific surface loop regions (see Fig. 1 for example), having a second dimer unit which clamps these fluctuations could have profound affects on the dissociation pathways. Furthermore, other systems involving biotin are now known to display considerable allosteric affects as well - often involving surface loops, 50,51 and we wonder if we might find similar affects at work in the full tetramer.

Conclusions In this work we used recently proposed enhanced sampling methods to gain understanding into the molecular determinants and bottlenecks in the dissociation of the dimeric biotinstreptavidin system in wild type and with a point mutation. We were able to reproduce

19



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 experimentally reported effect of the mutation on the dissociation rate. We also gained insight into a longstanding question regarding cause/effect in the coupled events of bond stretching and bond hydration during dissociation. Our results demonstrate that in this system and for the choice of molecular dynamics force-fields made here, it is the bond stretching and not hydration which forms the bottleneck in the early parts of the dissociation process. We conclude by sharing our optimism that recent and continuing developments in atomistic simulations and enhanced sampling are now taking us into an exciting age of reliable and routine calculations of the dynamics of protein-ligand interactions.

Supporting Information Available The supporting information file contains further details of SGOOP, representative spectral gap landscapes for the wild type and the mutant, information on the reliability of timescales, and other technical details.

ACKNOWLEDGMENTS The author is grateful to Dorothy Beckett for suggesting this system, for several helpful discussions, and for reading the manuscript very carefully. The author is also grateful to Bruce Berne and Michele Parrinello for countless inspiring discussions. The authors acknowledge the University of Maryland supercomputing resources (http://hpcc.umd.edu) and resources from Extreme Science and Engineering Discovery Environment (XSEDE) [TGMCA08X002] made available for conducting the research reported in this paper. The author also acknowledges financial support from University of Maryland start-up grant. The author thanks Anthony Clark for help with system preparation, and Purushottam Dixit for explaining to him the principle of Maximum Caliber. Author contributions: P.T. designed the research, performed calculations, analyzed 20


Page 20 of 27

Page 21 of 27

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


data and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the SI. Additional data related to this paper may be requested from the author.

References (1) Weber, P. C.; Ohlendorf, D.; Wendoloski, J.; Salemme, F. Structural origins of highaffinity biotin binding to streptavidin. Science 1989, 243, 85–88. (2) Hyre, D. E.; Amon, L. M.; Penzotti, J. E.; Le Trong, I.; Stenkamp, R. E.; Lybrand, T. P.; Stayton, P. S. Early mechanistic events in biotin dissociation from streptavidin. Nat. Struc. Mol. Bio. 2002, 9, 582–585. (3) Freitag, S.; Chu, V.; Penzotti, J. E.; Klumb, L. A.; To, R.; Hyre, D.; Le Trong, I.; Lybrand, T. P.; Stenkamp, R. E.; Stayton, P. S. A structural snapshot of an intermediate on the streptavidin–biotin dissociation pathway. Proc. Natl. Acad. Sci. 1999, 96, 8384–8389. (4) Caswell, K.; Wilson, J. N.; Bunz, U. H.; Murphy, C. J. Preferential end-to-end assembly of gold nanorods by biotin- streptavidin connectors. J. Amer. Chem. Soc. 2003, 125, 13914–13915. (5) Grubmüller, H.; Heymann, B.; Tavan, P. Ligand binding: molecular mechanics calculation of the streptavidin-biotin rupture force. Science 1996, 271, 997–999. (6) Young, T.; Abel, R.; Kim, B.; Berne, B. J.; Friesner, R. A. Motifs for molecular recognition exploiting hydrophobic enclosure in protein–ligand binding. Proc. Natl. Acad. Sci. 2007, 104, 808–813.

21



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

(7) Chilkoti, A.; Stayton, P. S. Molecular origins of the slow streptavidin-biotin dissociation kinetics. J. Amer. Chem. Soc. 1995, 117, 10622–10628. (8) Copeland, R. A. The drug-target residence time model: a 10-year retrospective. Nat. Rev. Drug. Discov. 2016, 15, 87–95. (9) Tautermann, C. S. Impact, determination and prediction of drug–receptor residence times for GPCRs. Curr. Opin. Pharm. 2016, 30, 22–26. (10) Dickson, A.; Tiwary, P.; Vashisth, H. Kinetics of ligand binding through advanced computational approaches: a review. Curr. Top. Med. Chem. 2017, 17, 2626–2641. (11) Tiwary, P.; Berne, B. J. Spectral gap optimization of order parameters for sampling complex molecular systems. Proc. Natl. Acad. Sci. 2016, 113, 2839–2844. (12) Tiwary, P.; Parrinello, M. From metadynamics to dynamics. Phys. Rev. Lett. 2013, 111, 230602–230606. (13) 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, E386–E391. (14) Tiwary, P.; Mondal, J.; Berne, B. J. How and when does an anticancer drug leave its binding site? Science Advances 2017, 3, e1700014–e1700020. (15) Casasnovas, R.; Limongelli, V.; Tiwary, P.; Carloni, P.; Parrinello, M. Unbinding kinetics of a p38 MAP kinase type II inhibitor from metadynamics simulations. J. Amer. Chem. Soc. 2017, 139, 4780–4788. (16) Wang, Y.; Papaleo, E.; Lindorff-Larsen, K. Mapping transiently formed and sparsely populated conformations on a complex energy landscape. eLife 2016, 5, e17505–e17516. (17) Stuckey, J. I.; Dickson, B. M.; Cheng, N.; Liu, Y.; Norris, J. L.; Cholensky, S. H.; Tempel, W.; Qin, S.; Huber, K. G.; Sagum, C.; et al, A cellular chemical probe targeting 22


Page 22 of 27

Page 23 of 27

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 chromodomains of Polycomb repressive complex 1. Nature Chem. Bio. 2016, 12, 180–187. (18) Tiwary, P.; Berne, B. J. How wet should be the reaction coordinate for ligand unbinding? J. Chem. Phys. 2016, 145, 054113–054118. (19) Fleming, K. L.; Tiwary, P.; Pfaendtner, J. New approach for investigating reaction dynamics and rates with ab initio calculations. J. Phys. Chem. A 2016, 120, 299–305. (20) Bochicchio, D.; Salvalaglio, M.; Pavan, G. M. Into the dynamics of a supramolecular polymer at submolecular resolution. Nature Comm. 2017, 8, 147–157. (21) Aslan, F. M.; Yu, Y.; Vajda, S.; Mohr, S. C.; Cantor, C. R. Engineering a novel, stable dimeric streptavidin with lower isoelectric point. J. Biotech. 2007, 128, 213–225. (22) Klumb, L. A.; Chu, V.; Stayton, P. S. Energetic roles of hydrogen bonds at the ureido oxygen binding pocket in the streptavidin- biotin complex. Biochemistry 1998, 37, 7657–7663. (23) 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. (24) 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. (25) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comp. Chem. 2004, 25, 1157–1174. (26) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF ChimeraâĂŤa visualization system for exploratory research and analysis. J. Comp. Chem. 2004, 25, 1605–1612. 23



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

(27) Evans, D. J.; Holian, B. L. The nose–hoover thermostat. J. Chem. Phys. 1985, 83, 4069–4074. (28) Parrinello, M.; Rahman, A. Crystal structure and pair potentials: A moleculardynamics study. Phys. Rev. Lett. 1980, 45, 1196–1199. (29) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theor. Comp. 2008, 4, 435–447. (30) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comp. Phys. Comm. 2014, 185, 604–613. (31) Salvalaglio, M.; Tiwary, P.; Parrinello, M. Assessing the reliability of the dynamics reconstructed from metadynamics. J. Chem. Theor. Comp. 2014, 10, 1420–1425. (32) Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing important fluctuations: rare events and metadynamics from a conceptual viewpoint. Ann. Rev. Phys. Chem. 2016, 67, 159–184. (33) Berne, B. J.; Borkovec, M.; Straub, J. E. Classical and modern methods in reaction rate theory. J. Phys. Chem. 1988, 92, 3711–3725. (34) Voter, A. F. Hyperdynamics: Accelerated Molecular Dynamics of Infrequent Events. Phys. Rev. Lett. 1997, 78, 3908–3911. (35) Grubmüller, H. Predicting slow structural transitions in macromolecular systems: Conformational flooding. Phys. Rev. E 1995, 52, 2893–2899. (36) Tiwary, P.; Berne, B. J. Predicting reaction coordinates in energy landscapes with diffusion anisotropy. J. Chem. Phys. 2017, 147, 152701–152707. (37) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603–020606. 24


Page 24 of 27

Page 25 of 27

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


(38) Tiwary, P.; Parrinello, M. A time-independent free energy estimator for metadynamics. J. Phys. Chem. B 2014, 119, 736–742. (39) Pressé, S.; Ghosh, K.; Lee, J.; Dill, K. A. Principles of maximum entropy and maximum caliber in statistical physics. Rev. Mod. Phys. 2013, 85, 1115–1132. (40) Dixit, P. D.; Jain, A.; Stock, G.; Dill, K. A. Inferring transition rates of networks from populations in continuous-time Markov processes. J. Chem. Theor. Comp. 2015, 11, 5464–5472. (41) Agmon, N.; Hopfield, J. CO binding to heme proteins: a model for barrier height distributions and slow conformational changes. J. Chem. Phys. 1983, 79, 2042–2053. (42) Poon, G. G.; Peters, B. Accelerated nucleation due to trace additives: a fluctuating coverage model. J. Phys. Chem. B 2015, 120, 1679–1684. (43) Moritz, C.; TrÃűster, A.; Dellago, C. Interplay of fast and slow dynamics in rare transition pathways: the disk-to-slab transition in the 2d Ising model. J. Chem. Phys. 2017, 147, 152714–152721. (44) Berezhkovskii, A.; Szabo, A. One-dimensional reaction coordinates for diffusive activated rate processes in many dimensions. J. Chem. Phys. 2005, 122, 014503–014506. (45) Prinz, J.-H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J. D.; SchÃĳtte, C.; NoÃľ, F. Markov models of molecular kinetics: generation and validation. J. Chem. Phys. 2011, 134, 174105–174111. (46) Plattner, N.; Noe, F. Protein conformational plasticity and complex ligand-binding kinetics explored by atomistic simulations and Markov models. Nat. Comm. 2015, 6 . (47) Morando, M. A.; Saladino, G.; DâĂŹAmelio, N.; Pucheta-Martinez, E.; Lovera, S.; Lelli, M.; López-Méndez, B.; Marenchino, M.; Campos-Olivas, R.; Gervasio, F. L. Con-

25



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

formational selection and induced fit mechanisms in the binding of an anticancer drug to the c-Src kinase. Sci. Rep. 2016, 6, 24439–24447. (48) Dickson, A.; Lotz, S. D. Multiple ligand unbinding pathways and ligand-induced destabilization revealed by WExplore. Biophys. Jour. 2017, 112, 620–629. (49) Kwarcinski, F. E.; Brandvold, K. R.; Phadke, S.; Beleh, O. M.; Johnson, T. K.; Meagher, J. L.; Seeliger, M. A.; Stuckey, J. A.; Soellner, M. B. Conformation-selective analogues of Dasatinib reveal insight into kinase inhibitor binding and selectivity. ACS Chem. Bio. 2016, 11, 1296–1304. (50) Weaver, L. H.; Kwon, K.; Beckett, D.; Matthews, B. W. Corepressor-induced organization and assembly of the biotin repressor: a model for allosteric activation of a transcriptional regulator. Proc. Natl. Acad. Sci. 2001, 98, 6045–6050. (51) Streaker, E. D.; Beckett, D. Ligand-linked structural changes in the Escherichia coli biotin repressor: the significance of surface loops for binding and allostery. J. Mol. Bio. 1999, 292, 619–632.

26


Page 26 of 27

Page 27 of 27



Molecular Determinants and Bottlenecks in the Dissociation Dynamics

Recommend Documents