Unfolding Dynamics of Ubiquitin from Constant Force MD Simulation

Subscriber access provided by EKU Libraries

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Unfolding Dynamics of Ubiquitin from Constant Force MD Simulation: Entropy-Enthalpy Interplay Shapes the Free Energy Landscape Anil Kumar Sahoo, Biman Bagchi, and Prabal K. Maiti J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b09318 • Publication Date (Web): 21 Jan 2019 Downloaded from http://pubs.acs.org on January 26, 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 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

Unfolding Dynamics of Ubiquitin from Constant Force MD Simulation: Entropy-Enthalpy Interplay Shapes the Free Energy Landscape Anil Kumar Sahoo,

†

∗,‡

Biman Bagchi,

∗,†

and Prabal K. Maiti

Center for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore-560012, India, and Solid State and Structural Chemistry Unit, Indian Institute of Science, Bangalore-560012, India

E-mail: [email protected]; [email protected]

∗ †

To whom correspondence should be addressed Center for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore-

560012, India

‡

Solid State and Structural Chemistry Unit, Indian Institute of Science, Bangalore-560012, India

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 27

Abstract Force probe methods are routinely used to study conformational transitions of biomolecules at single-molecule level. In contrast to simple kinetics, some proteins show complex response to mechanical perturbations that is manifested in terms of unusual force-dependent kinetics. Here, we study, via fully atomistic molecular dynamics (MD) simulations, constant force-induced unfolding of ubiquitin protein. Our simulations reveal a crossover at an intermediate force (about 400 pN) in the unfolding rate versus force curve. We nd by calculation of multidimensional free energy landscape (FEL) of the protein that the complex unfolding kinetics is intimately related to the force dependent modications in the FEL. Pearson correlation coecient analysis allowed us to identify two appropriate order parameters describing the unfolding transition. The crossover in the rate can be explained in terms of an interplay between entropy and enthalpy with relative importance changing from low force to high force. We rationalize the results by using multidimensional transition state theory.

Introduction A protein has to undergo large scale conformational transitions (folding/unfolding) once or multiple times during any cellular cycle. Understanding protein folding/unfolding pathways, its response to external perturbation and dynamics are key to in vitro protein designability, and its in vivo functions. iments.

This kind of information is often not available in bulk exper-

Single-molecule experiments

1

hold promises to detect these rare conformational

transitions, unique transition pathways and intermediate states

2

during protein folding or

unfolding. Especially, single-molecule force spectroscopy (SMFS) experiments

3,4

have been

used to measure many aspects of protein folding/unfolding process, such as (i) the free energy barrier and FELs along the pulling direction, kinetics of conformation transitions, times,

10

9

5,6

(ii) energy landscape roughness,

7,8

(iii)

(iv) more recently, transition paths and transition path

and (v) thus allows us to unravel multi-pathways folding/unfolding transitions

2


1113

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


as anticipated in the folding funnel paradigm.

14,15

Although SMFS techniques have evolved over years to the improved resolution of sub-

µs

nanometer in distance and

in time, these experiments can access only the molecular

extension or the pulling reaction coordinate (RC) during any conformational transition. If the pulling RC is the slowest among all other possible structural RCs and it is aligned along the direction of actual reaction, then the kinetics of the molecular rupture can be described by one dimensional (1D) diusion in the FEL along the pulling RC using Kramers' theory.

16,17

The above scenario is typically characterized by a straight line or a slight down turn in the force versus logarithm of the transition rate (f vs logk(f )) plot, where the force directly decreases the activation barrier for a reaction.

17

Since any macromolecule such as a protein contains thousands of atoms, its potential energy surface is intrinsically multidimensional. The energy landscapes may thus encompass multiple pathways, where pathways switch may happen by changing force; also contain intermediates.

20

11,12,18,19

it can

The reaction path can proceed in orthogonal to the actual

pulling direction for a range of forces, and the reaction direction may get aligned to the pulling direction as force is increased further.

21

Thus, many interesting scenarios for the

force dependent reaction rate can arise due to the multidimesional nature of the FEL, such as (i) increase in

k(f )

as force is increased with a upturn in

non-monotonic force dependence of of forces,

22

k(f ):

a dip in

f

k(f )

vs logk(f ) plot,

19

and (ii)

vs logk(f ) plot for intermediate range

which the 1D model completely fails to capture.

The above scenarios for

f

17

can be captured by multidimensional rate theory.

23

There

are several approaches to calculate the barrier-crossing rate in a multidimensional FEL. The earliest is transition state theory (TST) that provides the rate in terms of activation energy, attempt frequency and orthogonal vibrational frequencies. TST rate expression is given by:

k

T ST

QN

(f ) = QNi=1 −1 i=1

ωi (f ) exp[−∆F † (f )/kB T ], ωiB (f )

3


24

(1)


where

∆F † (f )

Page 4 of 27

is the free energy barrier for activation at force

f , ωi (f )

and

ωiB (f )

are,

respectively, the frequencies of vibrations in the folding well and at the saddle point, and is the total number of vibrational modes. Landauer and Swanson

25

N

provided an expression

for diusive barrier crossing in multidimensional free energy surface.

Langer

26

provided

an elegant method to nd the decay of a metastable state across a saddle by nding the imaginary frequencies.

A further relevant advance was made by Pollack who derived a

general expression for the model of a bistable potential connected linearly to a system of harmonic oscillators.

27,28

The method of Pollack employs the same technique as Langer in

expressing rate in terms of an imaginary frequency. Also, for describing molecular transitions multidimensional diusive and multiple pathways models exist in the literature that can capture many nontrivial features of the force dependent reaction rate.

21,22

However, the

interpretation of SMFS results become doubtful because dierent mechanism can give rise to similar behavior of the force dependent reaction rate. So, nding the correct mechanism or model for a given SMFS result remains a key challenge.

29

The mechanical stability and folding of ubiquitin protein under force have been well explored by both experimental SMFS studies

9,30,31

and using coarse-grained simulations.

3234

Recently, the folding of ubiquitin has been studied by all-atom molecular dynamics (MD) simulations.

35

In another all-atom MD study,

36

the structure, elastic properties and dy-

namics of unfolded ubiquitin under force have been inquired. from simulation using coarse-grained Go ¯ model Best et al.

37

Interestingly, for ubiquitin reported a turnover: nearly

force independent unfolding time at low forces, whereas force dependent unfolding time at higher forces. They suggested a plausible scenario for observing the turnover by performing Brownian dynamics simulation using a model potential. By varying pulling directions they came up with a novel reaction coordinatethe pulling direction which shows the smallest turnoverfor describing the intrinsic folding pathway. They showed that this reaction coordinate can be used to reliably extract equilibrium kinetic informations from nonequilibrium measurements.

37

However, no report of such turnover exists for ubiquitin either in experi-

4


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


ment or in all-atom simulation. Little is known about how force modies the underlying free energy surface that gives rise to such complex unfolding kinetics. Also, the thermodynamic driving forces (entropy and enthalpy) and the microscopic origin for the turnover are not known from previous studies. Here, we use fully atomistic MD simulation with explicit solvent to study force induced unfolding dynamics of ubiquitin across wide range of perturbations. and analysis details are described in Methods section.

The MD simulation

Our simulations reveal a crossover

for the force dependent unfolding rate similar to the kinetic turnover reported in the coarsegrained Go ¯ model simulation.

37

We discuss the unfolding pathways at dierent forces in

terms of the multidimensional FELs, microscopic understanding for the complex kinetics in terms of interplay between entropy and enthalpy contributions to the free energy, and atomic structural details of the unfolding process.

Methods MD Simulation Details.

The crystal structure of the ubiquitin protein (PDB ID: 1UBQ)

was taken as the initial structure.

CHARMM 36 force eld

model the protein. The protein was solvated with TIP3P

39

38

parameters were used to

water molecules. A large enough

rectangular simulation box was chosen, considering the protein remains completely solvated in its unfolded state and does not interact with its periodic images. of 6.1×6.3×21.0 nm

3

Precisely, a box size

consisting of 77242 atoms was taken for the force of 20 pN, while for

all the other higher forces a box size of 5.1×5.3×30.0 nm

3

consisting of 77143 atoms was

considered. Periodic boundary condition (PBC) was used for all the simulations. NAMD-2.10 package

40

was used to perform all the simulations. Energy minimization of

the solvated system was done for 10000 steps. The temperature of the system was slowly increased to 300 K during heating. Then, the equilibration of the system was done in NPT ensemble (at 300 K temperature and 1 atm pressure). The temperature of the system was

5



maintained using Langevin thermostat

41

Page 6 of 27

with collision frequency of 1.0 ps

−1

. The pressure

was maintained using Nose-Hoover Langevin pressure control [see section 7.5.2 of NAMD manual, Ref: http://www.ks.uiuc.edu/Research/namd/2.10b1/ug.pdf ] with a piston period of 0.1 ps and damping time constant of 0.05 ps. The production run was performed in NVT ensemble, for which the pressure control was removed. SHAKE algorithm

42

was used to con-

strain all the bonds involving a hydrogen atom during the simulations (heating, equilibration and production runs), which facilitates use of 2 fs time step for integration of the equations of motion. The long range electrostatic interactions were calculated using the Particle Mesh Ewald (PME)

43

sum technique with a real space cuto of 12 Å. A switch function (with

switchdist = 10 Å) was used to truncate the van der Walls interaction at 12 Å. The Steered Molecular Dynamics (SMD) simulations were performed by restarting the simulation after the production run, where the N-terminal force was applied to the C-terminal

Cα

Cα atom was xed, and a constant

atom to unfold the protein completely. Many SMD

unfolding simulations of the protein were performed at each force, such as 20 pN (11 runs), 100 pN (10 runs), 200 pN (10 runs), 300 pN (5 runs), 400 pN (5 runs), 500 pN (5 runs), 600 pN (7 runs), 700 pN (7 runs) and 800 pN (7 runs). The total simulation time including all equilibrium MD and SMD runs is 22

Analysis.

µs.

All the analyses were performed using in-house codes. Images were processed

using Visual Molecular Dynamics (VMD) software.

44

The radius of gyration of a molecule is dened as

v u N u1 X mi |ri (t) − Rcom (t)|2 , Rg (t) = t

N

where

Rcom

is position of the center of mass and

is the position of the

ith

atom at time

(2)

i=1

t, mi

M

is the total mass of the molecule,

is the mass of the

number of atoms in the molecule.

6


ith

atom, and

N

ri (t)

is the total

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


The root mean square deviation in distances is dened by the following equation

v u N u1 X RM SD(t) = t |ri (t) − ri0 |2 ,

N

where

N

is the total number of heavy atoms in the molecule,

atom at time

t and ri0

(3)

i=1

ri (t)

is the position of the

ith

is the position of the same atom in the native structure of the protein.

The solvent accessible surface area (SASA) of the protein was calculated by rolling a sphere of radius 1.5 Å on the surface of the protein by using VMD software. The fraction of native contacts (Q) is dened as follows

and

j

at time t,

0 rij

Nc

pairs of native contacts

is the distance between

smoothening parameter taken to be 5 Å during contact formation.

−1

i and j , and

in the native structure of the protein,

Here, the list of native contact pairs

i

and

j

Φ(t) =

where the sum runs over the

t, φ0i

Nφ

µ is a

ν , taken to be 1.8, takes care the uctuations

Ri

and

(i, j)

Rj

is constructed by

such that

|Ri − Rj | >

is less than 4.5 Å.

The fraction of native dihedral angles (Φ) is dened as

time

(4)

(i, j), rij (t) is the distance between i

considering all pairs of heavy atoms belonging to residues 3 and the distance between

45

1 X 1 , 0 Nc (i,j) 1 + exp[µ(rij (t) − νrij )]

Q(t) =

where the sum runs over the

44

46

1 X (φi (t) − φ0i )2 exp[− ], Nφ i σ2

backbone dihedral angles

(5)

φi , φi (t)

is a dihedral angle at

is the same dihedral angle in the native structure of the protein, and

The joint entropy (SXY ) of two random variables,

SXY = −

X

X

and

Y,

P (x, y)ln[P (x, y)],

x∈X,y∈Y

7


σ = 60◦ .

is dened as

(6)


where

x

and

y

are specic values of

X

and

Y,

Page 8 of 27

respectively.

The secondary structures (SS) of the protein were calculated by using VMD software. The SS contents of a particular kind (either

α

or

β)

44

is dened as the ratio of number of

residues of the protein involved in the particular kind of SS formation at time

t to that

of in

the native state of the protein.

Results and discussion Force Dependent Unfolding Kinetics.

The unfolding simulation protocol for the protein

is depicted in Figure 1a. Unfolded structures of the protein for some pulling forces are shown in Figure 1c. We have taken forces in the range of 10-1000 pN that are accessible to typical atomic force microscopy (AFM) or optical tweezer pulling experiments.

3

The end-to-end

distance (L), dened as the distance between the N-terminal atom and the C-terminal atom, is monitored in course of a simulation to calculate the unfolding time at a particular force. The time evolution of

L

at dierent forces is plotted in Figure S1 in the supporting

information (SI). A typical unfolding event is dened by a sudden increase in value in the native state of the protein (L0

≈ 3.8 nm).

After the sudden increase,

L

from its

L saturates

at a higher value depending on the force. The equilibrium extension or end-to-end distance, and the corresponding unfolded structure of the protein as a function of force are shown We dene the unfolding time (tU ) of a

in Figures S2a and S2b in the SI, respectively. protein as the time taken for

L

to cross a critical value of 10 nm for the rst time.

The

force dependent unfolding rate (kU (f )) is inverse of the mean unfolding time (τU (f )) and is given by

kU (f ) = 1/τU (f ),

with

independent simulations at a force

τU (f ) = f

and

1 N

PN

τUi (f )

i i=1 τU (f ), where

N

is the total number of

is the unfolding time for

ith

simulation.

47

The natural logarithm of the unfolding rate (lnkU (f )) is plotted against force in Figure 1b. There are two distinct force regimes, namely the low-force regime (20-300 pN) and the high-force regime (500-800 pN), with a crossover near

8

f =


400 pN. In the low-force regime

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


kU

barely increases with force, whereas in the high-force regime

kU

increases exponentially

with force. Note the similarity of the observed complex kinetics with the earlier reported kinetic turnover.

37

Figure 1: Unfolding simulation protocol and the force dependent unfolding rates. (a) The native structure of protein ubiquitin is depicted in the secondary structure (SS) representation, where its SSs are named. The N-terminal force,

f,

is applied to the C-terminal

Cα

Cα

atom (red sphere) is xed, and a constant

atom (green sphere) to unfold the protein com-

pletely. (b) Natural logarithm of the ufolding rate (kU ) as a function of the pulling force is shown. Note that the force magnitudes (10-1000 pN) used here are accessible to typical AFM experiments.

3

(c) The unfolded structures of the protein for dierent forces are shown.

See Figure S2b in the SI for more such structures at other forces.

The insensitivity of

kU

to increasing forces in the low-force regime and the overall upturn

in the lnkU (f )−f plot cannot be explained by the Kramers' theory of diusive barrier crossing in a 1D FEL along the pulling direction (F (L)) containing either a single barrier multiple barriers.

48

Any upturn or upward curvature in the lnkU (f )

to either (i) switch in the unfolding pathways dynamics.

21

11,19,22,49

−f

17

or even

plot can be related

or (ii) multidimensional unfolding

In the former case, the unfolding process takes place by switching from one

9



Page 10 of 27

pathway to another pathway after a critical force, while in the latter case the unfolding process initially proceeds in some orthogonal direction that gets continuously aligned to the pulling direction by increasing force.

Force Dependent Modications of the Projected 2D FELs. crossover in the force dependence of

kU (f ),

To understand the

we calculate the multidimensional FELs by

monitoring along the protein unfolding trajectories the pulling RC (L) and many other structural RCs that are typically considered for describing protein folding/unfolding, such as the radius of gyration (Rg ), the root mean square deviation in distances (RM SD ), the solvent accessible surface area (SASA), the fraction of native contacts (Q) and fraction of native dihedral angles (Φ). Notably, it has been revealed in earlier studies that as

Φ 46

are good RCs for many proteins.

L) and Y

as well

All the structural RCs are dened in Methods

section. We construct the 2D joint probability density, RC

Q 45

P (X, Y ), between X

(the mechanical

(any one of the structural RCs), which is averaged over many (≈10) unfolding

trajectories at each force (Figure S1 in the SI). Hence, the obtained 2D FEL, using the relation

F (X, Y ) = −kB T lnP (X, Y )

F (X, Y ),

by

is reliable. The 2D FELs for dierent pair of

RCs at forces 20, 100 and 200 pN are shown in Figure 2 and in Figure S3 in the SI. To select among the structural RCs a better one to aid interpretation along with calculate Pearson correlation coecient variables

X

and

Y

ρX,Y

ρX,Y =

are completely correlated with quite independent of

L

L

we

(see the denition and plot in Figure 3a). Two

are completely uncorrelated or independent of each other if

are fully correlated when

L,

1. Among the RCs in the position space, over the entire force regime, while

Q

ρX,Y = 0, and

Rg

and

as well as

RM SD

SASA

are

in the low-force regime (20-400 pN). Among all the structural RCs,

Φ, dened in the dihedral space, is least correlated with L.

So, any one of the

L − Q or L − Φ

free energy surfaces can be considered to elucidate the multidimensional nature; though we have limited the discussion here to the former one, and pointed out the similarities and dierences for other free energy surfaces. The FEL at 20 pN force projected in the

L−Q

10

plane reveals many interesting features


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


Figure 2: Force dependent modications in the 2D FELs for dierent pairs of RCs.

The

projected FELs for unfolding of the protein at forces: 20, 100 and 200 pN are shown in the plane of (a)

L − Q and (b) L − SASA.

The folding basin (FB) and the unfolding basin (UB)

are pointed in each plot.

11



Figure 3: (a) Pearson correlation coecient,

ρX,Y ,

between X (L) and Y (anyone of the

structural RCs) as a function of force. It is given by the covariance,

σX

and

σY

Page 12 of 27

ρX,Y = cov(X, Y )/σX σY ,

where

cov

is

are the standard deviations of X and Y, respectively. (b) Scatter

plot for the joint probability distributions in the

L − Q plane at dierent forces.

The folding

basin (FB) and the unfolding basin (UB) are pointed in the plot. At 20 pN force starting from a highly curved minimum-energy unfolding pathway joining the FB and the UB, the pathway moves gradually to become more like a line by increasing the pulling force (see direction of the arrow).

The dashed black line joining the FB and the UB is for guiding

the viewer's eye. Notice that the minimum-energy pathways are almost linear starting from force of 500 pN.

(Figure 2a). The free energy basin corresponding to the folded state of the protein shows larger uctuations in are large in

Q

compared to that in the

L

direction. In contrast, the uctuations

L compared to Q direction in the unfolding basin.

As the pulling force increases,

the folding and the unfolding basins are squeezed both in the

L

and

Q

directions (see the

free energy plots at 100 pN and 200 pN in Figure 2a). At higher forces, the suppression of uctuations in both

L

and

Q

directions for the unfolded state is due to the much extended

conformation of the protein, where all the native contacts are broken; the protein does not have much conformational freedom and hence shows entropic elasticity. the folding basin at higher forces both in the

L

and

Q

36

The squeezing of

directions is rather subtle and may

be associated with the restricted freedom of the protein to move in any other direction apart from the pulling direction. Notice that in the

L−Q

free energy surface starting form a highly curved minimum-

energy unfolding pathway joining the folding and the unfolding basins at very low force (20

12


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


pN), the pathway moves to become more linear like as the pulling force increases (Figure 2a). The movement of minimum-energy pathways as a function of force can be clearly noticed in Figure 3b.

L − SASA

We see very similar multidimensional features for the free energy surfaces

(Figure 2b) and

L−Φ

(Figure S3a in the SI), i.e., the shape of the minimum-

energy pathway changes from a curved to more linear like as the force increases. Also, the force dependent modications of both the FELs discussed

L−Q

L−Φ and L−SASA are similar to the above

FEL. In contrast, the FELs projected in the

L − RM SD

and the

L − Rg

planes do not show much multidimensional features even at low force of 20 pN (as evident from

ρX,Y ≈ 1

for

Rg

and

RM SD

in Figure 3); though the force dependent modications of

the FELs are similar to the other cases (see Figures S3b and S3c in the SI). The Pearson correlation coecient plot in Figure 3a shows that there is sharp increase in

ρX,Y

value in the high force regime (500-800 pN). Evidently, there is a crossover from

multidimensional to uni-dimensional behaviour that is intimately related to the observed crossover in

kU (f )

(see Figure 1).

Hence, from Figure 3 we see that as pulling force in-

creases the structural RCs start aligning towards the pulling RC (L).

This alignment

is almost complete for forces greater than 450 pN. These observations support the case of multidimensional unfolding dynamics.

21

Free Energy Barrier for Unfolding along L at Dierent Forces and the Role of Entropy.

Although it is known from earlier studies

low forces

L

may not be a good RC for describing protein unfolding, because force directly

couples to

L

we try to understand whether the crossover in

33,37

that for some proteins especially at

kU (f )

can be quntitatively ex-

plained by considering the force dependent modications in the projected 1D FEL Figure 4a).

50

The gure shows that the positions of the two minima of

F (L)

F (L) (see

(corresponding

to the FB and the UB) changes with force, where the latter is more aected. So, the nature of the unfolded state highly depends on force (see Figure S2b in the SI for the snapshots).

36

† The free energy barrier for unfolding (∆F ), the dierence in free energy of the transition state and the folded state, is plotted as a function of force in Figure 4b. Notice the unusual

13



Page 14 of 27

Figure 4: (a) Force dependent modications in the 1D FELs along the pulling direction

L

(F (L)) at dierent forces. The solid lines represent polynomial t to the data points. The folding basin (FB) and the unfolding basin (UB) are pointed in the plot. (b) The free energy † barrier for unfolding (∆F (f )) calculated from F (L) as a function of force. (c) Entropyenthalpy classication for the free energy change for unfolding as a function of force. The enthalpy contribution is given by: and

∆EU B−F B = ∆UU B−F B − f · ∆LU B−F B ,

where

∆UU B−F B

∆LU B−F B , respectively, are the potential energy dierence and the distance between the

unfolded and folded states. The entropy contribution (−T ∆SU B−F B ) has been evaluated by performing quasi-harmonic analysis

5153

for obtaining congurational entropies of the protein

both in the UB and in the FB. (d) The joint entropy (SXY ) between two RCs,

Y

(any one of the structural RCs, such as

Rg , RM SD, SASA, Q

force.

14


and

Φ),

X (L)

and

as a function of

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


feature of non-monotonic force dependence of

∆F † (f )

compared to the case where

∆F †

de-

creases monotonically with increasing the force. To understand this non-monotonic behavior one should separately look at the enthalpy and entropy contributions to decomposition for

∆F †

∆F † .

However, such

is a daunting task as it involves entropy (and enthalpy) calculation

at the transition state, which is very short-lived that does not allow sucient exploration. Instead, we quantify the enthalpy (∆EU B−F B ) and entropy (∆SU B−F B ) contributions for the free energy change for unfolding (∆FU B−F B ), i.e, the dierence in free energy of the unfolded and folded states (see the calculation details and plots in Figure 4c), and infer or estimate the same for

∆F † . ∆EU B−F B

with increasing the force, while

monotonically decreases, i.e, becomes more negative

∆SU B−F B

increases (becomes less negative) up to 500 pN,

and saturates after that (Figure 4c). Similarly, the enthalpy contributions to decrease with

f

according to the relation:

∆U † − f · ∆L† , where ∆U †

and

∆F †

should

∆L† , respectively,

are the potential energy dierence and the distance between the transition state and the folded state. In order to get an estimate of the entropy contributions to entropy,

SXY ,

(see Eq.

∆F † ,

we calculate the joint

6 in Methods for the denition) between two RCs,

X

a function of force for all pairs of RCs discussed before for the 2D FELs. Figure 4d, of RCs.

SXY

and

Y,

as

As shown in

decreases with force up to 400 pN, and saturates after that for all pairs

One can argue that the entropy might decrease by increasing the force due to

very fast transformation from the folded to the unfolded state at higher forces without exploring the whole conformational space, which would have been explored under very slow transformations in equilibrium folding/unfolding conditions. dependence of (∆F

†

∆F † (f )

So, the non-monotonic force

in Figure 4b can be understood as a shift from the entropy driven

increases due to entropic penalty) to enthalpy driven (∆F

†

decreases) force regimes.

This observation also perfectly correlates with the force dependence of

kU (f )

in the high-

force regime (Figure 1b). In this regime, starting from 500 pN with increasing the force decreases monotonically that facilitates exponential increase in

15


kU (f ).

∆F †


But, relative insensitivity of surprising (Figure 1b). in the rate

kU (f ).

kU (f )

Page 16 of 27

to force in the low-force regime (20-300 pN) is rather

∆F †

In this regime,

increases with force that would give a fall

To understand the insensitivity of

kU (f ),

one should also take into

account force dependence of the barrier-crossing attempt frequency (A(f )); as is known from TST

24

kU (f ) depends

that the rate

A(f )exp(−∆F † (f )/kB T ).

We nd that

on both

A(f )

A(f ) and ∆F † (f ) by

makes

kU (f )

A(f )

kU (f ) =

increases with the pulling force (see Figure 5

and the description of theoretical model for the calculation of the increase in

the relation

A(f )

in the next section). So,

compensates for the contribution to the rate coming from

∆F † .

This

relatively insensitive to force in the low-force regime.

Heuristic Model for Calculation of

A(f ).

From transition state theory,

24

attempt

frequency for unfolding of a protein is proportional to frequency in the reactant/folding basin (ω ); precisely

A = ω/2π .

We consider the folding basin to be a harmonic trap (

1 2 2 ω L ). And 2

we solve for Brownian motion of the end-to-end distance (L) of protein in that harmonic trap to calculate given by,

ω.

The mean square displacement (MSD) of

D

t→∞

kB

in the long-time limit is

54

D

E

E

lim ∆L2 (t) = lim L2 (t) − hL(t)i2 =

where

L

t→∞

is the Boltzmann constant and

T

kB T , ω2

(7)

is the absolute temperature. We calculate MSD in

the folding basin at dierent pulling force. As expected at each force MSD saturates at some nite value (see Figure 5a), and that value is equated to the right hand side of Eq. 7 to get Since the saturation value of MSD decreases as the pulling force increases,

A(f ) increases with the pulling force (Figure 5b). directly from simulation using the relation

A(f ) = kU (f )exp(∆F † (f )/kB T ), †

A(f )

where both the

(f )) for unfolding are obtained from

simulation results. The theoretical and the simulation results for

A(f )

and hence

We also calculate force dependence of

unfolding rate (kU (f )) and the free energy barrier (∆F

5c for comparison.

ω(f )

ω.

A(f )

are shown in Figure

Although the theoretical model captures qualitatively the increase in

by increasing force, quantitatively the model prediction doesn't match quite well with

the simulation result across the whole range of forces.

16


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


Figure 5: (a) Mean square displacement (MSD) of the end-to-end distance, in the folding basin is plotted as a function of time lag

ω,

∆t.

L, of the protein

(b) The frequency of oscillation,

in the folding basin calculated by the theoretical model is shown as a function of pulling

force. (c) The attempt frequency (A(f )) as a function of force divided by

A

at the lowest

pulling force (20 pN) calculated from simulation results, and the same obtained from the theoretical model (ω(f )/ω(20)) are plotted for comparison.

Mechanism for the Complex Unfolding Kinetics. discuss Is the crossover in

kU (f )

Now, we are in a position to

due to a discrete discernible switch in the pathways, i.e,

the protein unfolds in one pathway in the low-force regime, and after a critical force it switches to another parallel unfolding pathway?. Note that protein folding and unfolding are cooperative phenomena. increase in

L

Before the unfolding transition, characterized by a sudden

(see Figure S1 in the SI), depending on the force many local structures of the

protein melts that eventually determines nature of the unfolding pathways. So, to decipher the mechanism of force dependent unfolding pathways it is appropriate to monitor both melting of local structures and breaking of tertiary contacts along the unfolding trajectory. We have calculated secondary structures (SS) of the protein to detect the former. The SS

17



Page 18 of 27

contents along the unfolding trajectory at dierent forces are plotted in Figure 6. At low forces (20-200 pN), before the unfolding transition the of the protein) decrease signicantly compared to the

α-helix

β -sheets

contents (the exible parts

contents that forms the rigid

parts of the protein, and vice versa at high forces (600-800 pN). While at intermediate force (300-500 pN), both

α

and

β

contents reduce before the unfolding transition.

The above observations also correlate with the analysis for lifetimes of all the tertiary native contact pairs of the protein as a function of force (see Figure S4 in the SI). The alpha helix

αH1

(see Figure 1a and the rectangular mark in Figure S4 in the SI) especially in the

low-force regime starts melting much earlier than the unfolding transition. In contrast, due to the presence of ve hydrogen bonds between

β 1 and β 5 (trapezoid mark) the mechanically

resilient contacts break at the onset of the unfolding transition over the entire force regime. The contacts between the loop region after

β4

β2

with some part of

αH1

and the loop after

along with H3 (see Figure 1a and the rhombus mark in Figure S4 in the SI) decrease in

the low-force regime, while

β 3-β 5

contacts (ellipse mark) decrease in the high-force regime

much before the unfolding transition. From the microscopic structural analyses and the 2D FELs (see Figure 2 and notice the continuous movements of transirion paths as a function of force in Figure 3b), we do not nd any discrete switch in the unfolding pathways going from lower to higher pulling forces. The complex behavior for the unfolding kinetics (Figure 1b) is due to the force dependent modications in the multidimensional FEL. The unfolding process initially proceeds in some orthogonal direction to the pulling direction at low forces. As force increases, the unfolding direction gradually turns towards the pulling direction in a continuous fashion.

The

unfolding direction gets completely aligned in the direction of pulling at a force of 500 pN (notice the sharp increase in

ρX,Y

value in Figure 3a, and the saturation of

Figure 4d close to this force). And from this force onwards (Figure 1b).

18


kU (f )

SXY

value in

increases exponentially

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


Figure 6: The time evolution of secondary structures (α-helix and the end-to-end distance,

L,

β -sheets)

contents and

(see right axis) of the protein along a representative unfolding

trajectory at each pulling forces are shown. The pulling force is mentioned in each plot.

19



Page 20 of 27

Conclusions Our all-atom MD simulations for ubiquitin show an intriguing crossover in the unfolding rate from almost complete independence of pulling force

f

to an exponential dependence.

The results can be rationalized by using multidimensional nature of the FEL (see Figure 2). Importantly, the FELs show continuous movement of minimum-energy unfolding pathways as a function of force (Figure 3b). The crossover can be understood quantitatively by taking

† into account force dependence of both the free energy barrier for unfolding (∆F ) and the barrier-crossing attempt frequency. Overall, our work provides valuable insights on the origin of the crossover in case of ubiquitin,

37

and will also help in interpreating the force-dependent

complex kinetics reported for other proteins.

1113,18,19,22,49

It is expected that protein containing both

α

(exible) and

β

(rigid) structures might

show similar force dependent complex unfolding kinetics. Because of inherent heterogeneity in the architecture of a protein the dynamics may be heterogeneous even in equilibrium conditions, and external force can aect the dynamics of dierent parts very dierently. Thus, pulling experiments of this kind can provide valuable information on protein response and dynamics, and complexity of the underlying free energy landscape.

Acknowledgement We thank SERC, IISc, for providing computational facility @SAHASRAT machine and DAE, India for nancial support.

A.K.S. thanks MHRD for a fellowship.

B.B. thanks

DST (nanomission) and Sir JC Bose Fellowship (DST) for partial support.

Supporting Information Available Additional results and gures (PDF). This material is available free of charge via the Internet at

http://pubs.acs.org/.

20


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


References (1) Chung, H. S. Transition path times measured by single-molecule spectroscopy.

Biol.

J. Mol.

2018, 430, 409423.

(2) Marszalek, P. E.; Lu, H.; Li, H.; Carrion-Vazquez, M.; Oberhauser, A. F.; Schulten, K.; Fernandez, J. M. Mechanical unfolding intermediates in titin modules.

Nature

1999,

402, 100103. (3) Neuman, K. C.; Nagy, A. Single-molecule force spectroscopy: optical tweezers, magnetic tweezers and atomic force microscopy.

Nat. Methods

2008, 5, 491505.

(4) oldák, G.; Rief, M. Force as a single molecule probe of multidimensional protein energy landscapes.

Curr. Opin. Struct. Biol.

2013, 23, 4857.

(5) Hummer, G.; Szabo, A. Free energy surfaces from single-molecule force spectroscopy.

Acc. Chem. Res.

2005, 38, 504513.

(6) Woodside, M. T.; Block, S. M. Reconstructing folding energy landscapes by singlemolecule force spectroscopy.

Annu. Rev. Biophys.

2014, 43, 1939.

(7) Hyeon, C.; Thirumalai, D. Can energy landscape roughness of proteins and RNA be measured by using mechanical unfolding experiments?

Proc. Natl. Acad. Sci. U.S.A.

2003, 100, 1024910253. (8) Chavez, L. L.; Onuchic, J. N.; Clementi, C. Quantifying the roughness on the free energy landscape: entropic bottlenecks and protein folding rates.

J. Am. Chem. Soc.

2004, 126, 84268432. (9) Fernandez, J. M.; Li, H. Force-clamp spectroscopy monitors the folding trajectory of a single protein.

Science

2004, 303, 16741678.

21



Page 22 of 27

(10) Neupane, K.; Foster, D. A.; Dee, D. R.; Yu, H.; Wang, F.; Woodside, M. T. Direct observation of transition paths during the folding of proteins and nucleic acids.

Science

2016, 352, 239242. (11) Jagannathan, B.; Elms, P. J.; Bustamante, C.; Marqusee, S. Direct observation of a force-induced switch in the anisotropic mechanical unfolding pathway of a protein.

Natl. Acad. Sci. U.S.A.

Proc.

2012, 109, 1782017825.

(12) Guinn, E. J.; Jagannathan, B.; Marqusee, S. Single-molecule chemo-mechanical unfolding reveals multiple transition state barriers in a small single-domain protein.

Commun.

Nat.

2015, 6, 6861.

(13) Schönfelder, J.; Perez-Jimenez, R.; Muñoz, V. A simple two-state protein unfolds mechanically via multiple heterogeneous pathways at single-molecule resolution.

Commun.

Nat.

2016, 7, 11777.

(14) Onuchic, J. N.; Luthey-Schulten, Z.; Wolynes, P. G. Theory of protein folding: the energy landscape perspective.

Annu. Rev. Phys. Chem.

1997, 48, 545600.

(15) Dill, K. A.; Chan, H. S. From Levinthal to pathways to funnels.

Nat. Struct. Mol. Biol.

1997, 4, 1019. (16) Kramers, H. A. Brownian motion in a eld of force and the diusion model of chemical reactions.

Physica

1940, 7, 284304.

(17) Dudko, O. K.; Hummer, G.; Szabo, A. Intrinsic rates and activation free energies from single-molecule pulling experiments.

Phys. Rev. Lett.

2006, 96, 108101.

(18) Yew, Z. T.; Schlierf, M.; Rief, M.; Paci, E. Direct evidence of the multidimensionality of the free-energy landscapes of proteins revealed by mechanical probes.

2010, 81, 031923.

22


Phys. Rev. E

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


(19) Zhuravlev, P. I.; Hinczewski, M.; Chakrabarti, S.; Marqusee, S.; Thirumalai, D. Forcedependent switch in protein unfolding pathways and transition-state movements.

Natl. Acad. Sci. U.S.A.

Proc.

2016, 113, E715E724.

(20) Hyeon, C.; Thirumalai, D. Multiple barriers in forced rupture of protein complexes.

Chem. Phys.

J.

2012, 137, 055103.

(21) Suzuki, Y.; Dudko, O. K. Single-molecule rupture dynamics on multidimensional landscapes.

Phys. Rev. Lett.

2010, 104, 048101.

(22) Pierse, C. A.; Dudko, O. K. Distinguishing signatures of multipathway conformational transitions.

Phys. Rev. Lett.

2017, 118, 088101.

(23) Hänggi, P.; Talkner, P.; Borkovec, M. Reaction-rate theory: fty years after Kramers.

Rev. Mod. Phys. (24) Bagchi, B.

1990, 62, 251.

Molecular relaxation in liquids ; OUP USA, 2012.

(25) Landauer, R.; Swanson, J. Frequency factors in the thermally activated process.

Rev.

Phys.

1961, 121, 1668.

(26) Langer, J. S. Statistical theory of the decay of metastable states.

Ann. Phys. (N.Y)

1969, 54, 258275. (27) Grote, R. F.; Hynes, J. T. The stable states picture of chemical reactions. II. Rate constants for condensed and gas phase reaction models.

J. Chem. Phys.

1980,

73,

27152732.

(28) Pollak, E. Theory of activated rate processes: A new derivation of Kramers' expression.

J. Chem. Phys.

1986, 85, 865867.

(29) Makarov, D. E. Perspective: Mechanochemistry of biological and synthetic molecules.

J. Chem. Phys.

2016, 144, 030901. 23



Page 24 of 27

(30) Carrion-Vazquez, M.; Li, H.; Lu, H.; Marszalek, P. E.; Oberhauser, A. F.; Fernandez, J. M. The mechanical stability of ubiquitin is linkage dependent.

Biol.

Nat. Struct. Mol.

2003, 10, 738743.

(31) Schlierf, M.; Li, H.; Fernandez, J. M. The unfolding kinetics of ubiquitin captured with single-molecule force-clamp techniques.


2004, 101,

72997304.

(32) Li, P.-C.; Makarov, D. E. Ubiquitin-like protein domains show high resistance to mechanical unfolding similar to that of the I27 domain in titin: evidence from simulations.

J. Phys. Chem. B

2004, 108, 745749.

(33) Li, P.-C.; Makarov, D. E. Simulation of the mechanical unfolding of ubiquitin: probing dierent unfolding reaction coordinates by changing the pulling geometry.

Phys.

J. Chem.

2004, 121, 48264832.

(34) Kirmizialtin, S.; Huang, L.; Makarov, D. E. Topography of the free-energy landscape probed via mechanical unfolding of proteins.

J. Chem. Phys.

2005, 122, 234915.

(35) Piana, S.; Lindor-Larsen, K.; Shaw, D. E. Atomic-level description of ubiquitin folding.


2013, 110, 59155920.

(36) Stirnemann, G.; Giganti, D.; Fernandez, J. M.; Berne, B. Elasticity, structure, and relaxation of extended proteins under force.


2013, 110,

38473852.

(37) Best, R. B.; Paci, E.; Hummer, G.; Dudko, O. K. Pulling direction as a reaction coordinate for the mechanical unfolding of single molecules.

J. Phys. Chem. B

2008, 112,

59685976.

(38) Huang, J.; MacKerell, A. D. CHARMM36 all-atom additive protein force eld: Validation based on comparison to NMR data.

24

J. Comput. Chem.


2013, 34, 21352145.

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


(39) 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, 926935. (40) Phillips, J. C.;

Braun, R.;

Wang, W.;

Gumbart, J.;

Tajkhorshid, E.;

Villa, E.;

Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD.

J. Comput. Chem.

2005, 26, 17811802.

(41) Van Gunsteren, W. F.; Berendsen, H. A leap-frog algorithm for stochastic dynamics.

Mol. Simul.

1988, 1, 173185.

(42) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.

Comput. Phys.

1977, 23, 327341.

(43) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: Ewald sums in large systems.

J. Chem. Phys.

An N.log(N) method for

1993, 98, 1008910092.

(44) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics.

Graph.

J.

J. Mol.

1996, 14, 3338.

(45) Best, R. B.; Hummer, G.; Eaton, W. A. Native contacts determine protein folding mechanisms in atomistic simulations.


2013, 110, 17874

17879.

(46) Best, R. B.; Hummer, G. Coordinate-dependent diusion in protein folding.

Acad. Sci. U.S.A.

Proc. Natl.

2010, 107, 10881093.

(47) Reimann, P.; Schmid, G.; Hänggi, P. Universal equivalence of mean rst-passage time and Kramers rate.

Phys. Rev. E

1999, 60, R1.

(48) Garai, A.; Zhang, Y.; Dudko, O. K. Conformational dynamics through an intermediate.

J. Chem. Phys.

2014, 140, 135101. 25



Page 26 of 27

(49) Graham, T. G.; Best, R. B. Force-induced change in protein unfolding mechanism: discrete or continuous switch?

J. Phys. Chem. B

2011, 115, 15461561.

(50) Berezhkovskii, A.; Szabo, A. One-dimensional reaction coordinates for diusive activated rate processes in many dimensions.

J. Chem. Phys.

2005, 122, 014503.

(51) Karplus, M.; Kushick, J. N. Method for estimating the congurational entropy of macromolecules.

Macromolecules

1981, 14, 325332.

(52) Schlitter, J. Estimation of absolute and relative entropies of macromolecules using the covariance matrix.

Chem. Phys. Lett.

1993, 215, 617621.

(53) Andricioaei, I.; Karplus, M. On the calculation of entropy from covariance matrices of the atomic uctuations.

(54) Balakrishnan, V.

J. Chem. Phys.

2001, 115, 62896292.

Elements of nonequilibrium statistical mechanics ; Ane Books, 2008.

26


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


Graphical TOC Entry

27


Unfolding Dynamics of Ubiquitin from Constant Force MD Simulation

Recommend Documents