Time-Dependent Markov State Models for Single ... - ACS Publications

Feb 22, 2017 - Department of Chemical Engineering, Indian Institute of Technology Bombay, Mumbai 400076, India. •S Supporting Information. ABSTRACT:...
1 downloads 10 Views 7MB Size
Subscriber access provided by Fudan University

Letter

Time-Dependent Markov State Models for Single Molecule Force Spectroscopy Susmita Ghosh, Abhijit Chatterjee, and Swati Bhattacharya J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01094 • Publication Date (Web): 22 Feb 2017 Downloaded from http://pubs.acs.org on February 23, 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.

Journal of Chemical Theory and Computation 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 26

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

Journal of Chemical Theory and Computation

Time-Dependent Markov State Models for Single Molecule Force Spectroscopy

Susmita Ghosh1, Abhijit Chatterjee*2 and Swati Bhattacharya*1,2 1

Department of Physics, Indian Institute of Technology Guwahati, Guwahati 781039, India.

2

Department of Chemical Engineering, Indian Institute of Technology Bombay, Mumbai,

400076, India. KEYWORDS Markov state models, Single molecule Force Spectroscopy, Molecular Dynamics Simulations. Corresponding Author *E-mail: [email protected], [email protected] Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

ABSTRACT

This Letter demonstrates that using time-dependent Markov state models (TD-MSMs) one can obtain molecular-scale insights into force-extension curves for a variety of stretching experiments. A master-MSM constructed at a reference extension forms the basis for generating the required TD-MSM, i.e., the TD-MSM that is appropriate for the stretching experiment can be

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation

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 26

constructed from a single master-MSM. In addition, the availability of state-specific force models enable calculation of force-extension behaviour in a variety of ensembles. Changes in the network topology upon stretching is related through a thermodynamic quantity termed the mechanical disposition. Proof-of-principle is provided using a stretched alanine decapeptide under a time-varying pulling force.

Introduction A stretching force applied on a biomolecule can induce structural changes in the molecule that would have been otherwise rarely sampled in the absence of the applied force1. This has motivated a variety of single molecule force spectroscopy (SMFS) experiments2–6 employed in constant force, constant extension, force-ramp and force-jump modes. SMFS experiments are used to probe the structure, kinetic rates and interactions of molecules. Here we explore whether a similar idea can be exploited for rapidly constructing a detailed kinetic network model of the biomolecule using computer simulations, i.e., kinetic information pertaining to rare transitions can be recovered by stretching the biomolecule at conditions where transitions are more frequent. One potential beneficiary of such an approach would be SMFS experiments themselves. The usual analysis of SMFS experiments involves an excessively coarse-grained view7–14 . For instance, a 1-D reaction coordinate between the “folded” and “unfolded” states invoked in the interpretation of SMFS experiments hides the complexity of a rough free-energy landscape with multiple states and kinetic pathways15. Resolving the individual states can avoid complicated force-dependent rate maps arising from incorrectly lumping multiple intermediate states and competing pathways. Unfortunately, the inability to experimentally observe microscopic structural changes in a molecule presents an obstacle towards gaining higher resolution. Kinetic

ACS Paragon Plus Environment

2

Page 3 of 26

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

Journal of Chemical Theory and Computation

models derived bottom-up from molecular scimulations may be able to fill gaps in our understanding of the experimental force-extension curve. While this may provide unprecedented insights, it also raises broader questions whether i) detailed microscopic kinetic model can be sufficiently versatile to encompass a wide range of stretching experiments and ii) it is possible to convert between network models applicable for each type of SMFS experiment, e.g., constant force or extension, without the need to generate a new model for each new experiment. As a first step in this endeavour, we apply continuous-time Markov state models (MSM)16–19 to pulling experiments. A MSM, written as a master equation  

=   ,

(1)

approximates the dynamics in terms of state-to-state transitions while providing control over the model resolution. Here   is the occupation for a pre-specified list of states as a function of the pulling time and  is the rate matrix describing the rates of interconversion between the states. Generally,  is a constant. MSMs, including discrete-time MSMs20,21, have been used widely to study folding/unfolding/dynamics of peptides in the absence of external pulling forces20,21,16,22. However, extension of MSMs to SMFS experiments is not straight-forward due to the nonequilibrium nature of the force experiments. The main challenge is to relate two topologicallydistinct MSMs generated at separate stretching conditions. The key steps involved includes specifying the range of stretching conditions where the MSMs are valid, determining the list of relevant states/pathways over the entire range, and allowing  to be a function of time. In our time-dependent MSM (TD-MSM) formalism, we assume that vibrational timescales ≪ relaxation timescales within a Markov state  ≪ network relaxation timescales and timescales associated with pulling. This enables one to calculate state-specific elastic, kinetic and thermodynamic properties.

ACS Paragon Plus Environment

3

Journal of Chemical Theory and Computation

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 26

Theoretical Basis

The proof-of-principle of the TD-MSM formalism is provided by our study on a deca-alanine molecule in vacuum tethered at both ends by harmonic springs to two anchor points mimicking a dual optical trap setup (see Figure 1a.). Deca-alanine is a classic example of an alpha-helix with minimal interference from side-chain interactions. The spring constant of the device-molecule interaction is 0.86 kcal.mol-1 Å-2. Temperature  is 300 K. First, we consider a position-clamp setup, i.e., the anchor separation d is the control parameter such that the instantaneous force varies as in the absence of a feedback loop. States and kinetic pathways are identified using state-constrained MD (SCMD) calculations. SCMD entails detecting a transition from a state in question, noting the new state, and returning the system to the same state to record additional transitions23. In effect, this enhanced-sampling technique prevents the system from freely diffusing over the energy landscape and confines it to a particular state. A list of states of deca-alanine is determined by comparing the backbone atoms after aligning the molecule using the Kabsch algorithm using a tolerance of 3 Å. Standard practices like assessing the quality of the Markovian approximation23 (see Section E of Supporting Information for more details) were employed besides verifying the presence of a sharp single-peaked distribution for force and end-to-end distance for each Markov state. Stateconstrained MD snapshots were collected every 0.2  due to the rapid interconversion between states. A transition is said to have occurred when the system continues to reside in the new state for at least 1.2  after the transition is detected in the MD trajectory. Kinetic rates were determined using a maximum likelihood analysis. In general, MSMs constructed piecewise using MD have an associated uncertainty measured in terms of the probability that a pathway missing in the MSM would have been selected in the dynamics23,24. We define the completeness of a

ACS Paragon Plus Environment

4

Page 5 of 26

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

Journal of Chemical Theory and Computation

MSM in terms of one minus this probability. An adaptive algorithm, called programmed SCMD, was used to initiate trajectories from rare states until the transitions from the state have been mapped adequately, ensuring that the MSMs are complete with 90% confidence.

Figure 1. a) SMFS setup for deca-alanine. b) MSM-0 showing the dominant states of decaalanine at constant anchor separation  between 16-26 Å.

ACS Paragon Plus Environment

5

Journal of Chemical Theory and Computation

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

Page 6 of 26

a) Equilibrium occupations for different states in Figure 1b from individual MSMs

constructed for anchor separations  between 16-26 Å. State indices are labelled in the plot. Lines are quadratic fits. b) Corresponding kinetic rates (in symbols) for anchor separations  between 16-26 Å. Initial and final state indices are labelled in the plot. Lines correspond to Equation 3.

Regular MSMs of deca-alanine constructed using MD calculations at constant separations between 16 to 26 Å demonstrate a large topological variation with . These topological differences arise due to missing states and pathways as they become relevant to the dynamical evolution at the new stretching condition. A database of 810 states was compiled that can be utilized between =16-26 Å. Only 10 states have a combined occupation 95 % or more for the range of separations studied. For simplicity, we focus only on these states (see Figure 1b). As  evident from the equilibrium occupations   from the MSMs in Figure 2a, the MSM at 16 Å

has only two dominant states, namely, 1 and 288. Stretched configurations dominate at larger . Past studies on unravelling of the helical structure have not reported any metastable configurations. The MSMs reveal that folding/unfolding proceeds at the C-terminal via states 1, 2, 6 and 17 involving local readjustments (see Figure1b). Fluctuations in the middle residues can also frequently cause the helix to undergo deformation (state 1 to 3) that sometimes may engender an elongated configuration (state 3 to 6). However, this pathway is less preferred. Figure 2b shows the diverse behaviour of the kinetic rates as they span several orders of magnitude between 16-26Å. Some rates like 1-3 (9-1) increase (decrease), some like 2-3 have a maximum, while others like 1-113 have a relatively small variation.

ACS Paragon Plus Environment

6

Page 7 of 26

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

Journal of Chemical Theory and Computation

Consider a kinetic pathway  ⇋ ′ between states  and ′. We develop a kinetic theory relating the rates to the free energies of these states. Using transition state theory, the forward and backward rates are given by   =  exp

!"# 

and   =  exp

!"$ 

, respectively,

where the pre-exponential factors  and  are assumed to be constant, Δ&  and Δ&  are the free energy barriers, ' ( = )  and ) is the Boltzmann constant. To be thermodynamically * 

consistent, we require *#  = exp −'Δ&-  such that Δ&-  = &-  − &  is the $

free energy change along the forward pathway and &  denotes free energy of state . Following the Bell-Evans-Polyani principle25, we assume that the free energy barrier Δ& at trap separations  and . are related through Δ&  = Δ& .  + 0 1 where 1 = Δ&-  − Δ&- . .

(2)

1, termed as the mechanical disposition at separation , measures the thermodynamic preference for states  and ′ upon application of a force relative to a reference condition. Note that each pair of states possesses its own mechanical disposition, which has to be evaluated separately. The free energy barrier increases when 1 > 0. The kinetic rate   is related to the rate at separation . as   =  . exp −'0 1.

(3)

Similarly, the backward rate   =  . exp '1 − 0 1. Details of the derivation are included in Section A of the Supporting Information. For sake of brevity, the pair of states has not been explicitly mentioned in the notation for  ,  ,  ,  , Δ& , Δ& , 0 , 0 and 1. Generally, 0 and 0 are constants (see Supporting Information). When the trap separation  is a function of time, the mechanical disposition and rate matrix become time-dependent. The

ACS Paragon Plus Environment

7

Journal of Chemical Theory and Computation

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 26

intrinsic rate  . , the symmetry parameter 0 and 1 are estimated from MSMs constructed at constant , as described next.



At equilibrium, detailed balance requires    =  - . Combined with *# 

*$ 



= exp −'Δ&- , we write Δ&-  = ' ( ln  /- . Thus, the free energy

difference is calculated on the basis of the equilibrium state occupations from the MSMs at different values of . The main contribution to the free energy change arises from the potential energy associated with the state as the molecule is pulled. As described below, the potential energy for a state is a quadratic function in . Excellent quadratic polynomial fit for the free energy of state S calculated as −' ( ln   as a function of  is obtained in Figure 2a, which in turn yields 1 from Equation 2. Next, the kinetic parameters  .  and 0 are estimated by fitting Equation 3 to rates in the MSMs at different values of d in Figure 2b. Section C of Supporting Information provides a values for  .  and 0 for pairs of states. Choosing . = 16Å, a detailed MSM is constructed that contains 10 states of Figure 1b, alongside the kinetic parameters of Eq. 3. We term this master-MSM as MSM(. ) or simply MSM-0. Note that MSM-0 contains a list of states and kinetic pathways that would be relevant for certain/entire range of extensions to be sampled with the model, the associated kinetic parameters include  .  and 0 , and the thermodynamic parameter 1. MSM-0 is a precursor for TD-MSMs at a variety of stretching conditions. In addition, a force-separation relation is required for each state to compute the average forces acting on the molecules.

ACS Paragon Plus Environment

8

Page 9 of 26

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

Journal of Chemical Theory and Computation

The ensemble-averaged force can be calculated as : = ∑  :  where :  is the average force experienced when the molecule is in state . This expression is valid even when the molecule is stretched quickly. Inherent is the assumption that the Markovian approximation is valid at the pulling rate and the molecule rapidly enters mechanical equilibrium when a transition to a new state occurs. When deca-alanine is modelled as a harmonic spring with the spring constant for state S being  , the state-specific force from SCMD calculation is given by   :  =   −