Slow dynamics of a protein backbone in molecular dynamics simulation revealed by time-structure based independent component analysis Yusuke Naritomi and Sotaro Fuchigami Citation: The Journal of Chemical Physics 139, 215102 (2013); doi: 10.1063/1.4834695 View online: http://dx.doi.org/10.1063/1.4834695 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/21?ver=pdfcov Published by the AIP Publishing Articles you may be interested in On the influence of hydrated imidazolium-based ionic liquid on protein structure stability: A molecular dynamics simulation study J. Chem. Phys. 139, 115102 (2013); 10.1063/1.4821588 BDflex: A method for efficient treatment of molecular flexibility in calculating protein-ligand binding rate constants from Brownian dynamics simulations J. Chem. Phys. 137, 135105 (2012); 10.1063/1.4756913 Slow dynamics in protein fluctuations revealed by time-structure based independent component analysis: The case of domain motions J. Chem. Phys. 134, 065101 (2011); 10.1063/1.3554380 Application of principal component analysis in protein unfolding: An all-atom molecular dynamics simulation study J. Chem. Phys. 127, 165103 (2007); 10.1063/1.2796165 Application of time series analysis on molecular dynamics simulations of proteins: A study of different conformational spaces by principal component analysis J. Chem. Phys. 121, 4759 (2004); 10.1063/1.1778377

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

THE JOURNAL OF CHEMICAL PHYSICS 139, 215102 (2013)

Slow dynamics of a protein backbone in molecular dynamics simulation revealed by time-structure based independent component analysis Yusuke Naritomi1 and Sotaro Fuchigami2,a) 1

Department of Supramolecular Biology, Graduate School of Nanobioscience, Yokohama City University, 1-7-29 Suehiro-cho, Tsurumi-ku, Yokohama 230-0045, Japan 2 Department of Medical Life Science, Graduate School of Medical Life Science, Yokohama City University, 1-7-29 Suehiro-cho, Tsurumi-ku, Yokohama 230-0045, Japan

(Received 17 September 2013; accepted 12 November 2013; published online 3 December 2013) We recently proposed the method of time-structure based independent component analysis (tICA) to examine the slow dynamics involved in conformational fluctuations of a protein as estimated by molecular dynamics (MD) simulation [Y. Naritomi and S. Fuchigami, J. Chem. Phys. 134, 065101 (2011)]. Our previous study focused on domain motions of the protein and examined its dynamics by using rigid-body domain analysis and tICA. However, the protein changes its conformation not only through domain motions but also by various types of motions involving its backbone and side chains. Some of these motions might occur on a slow time scale: we hypothesize that if so, we could effectively detect and characterize them using tICA. In the present study, we investigated slow dynamics of the protein backbone using MD simulation and tICA. The selected target protein was lysine-, arginine-, ornithine-binding protein (LAO), which comprises two domains and undergoes large domain motions. MD simulation of LAO in explicit water was performed for 1 μs, and the obtained trajectory of Cα atoms in the backbone was analyzed by tICA. This analysis successfully provided us with slow modes for LAO that represented either domain motions or local movements of the backbone. Further analysis elucidated the atomic details of the suggested local motions and confirmed that these motions truly occurred on the expected slow time scale. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4834695] I. INTRODUCTION

Molecular dynamics (MD) simulation is a powerful tool that is widely used to elucidate dynamic behavior of proteins such as equilibrium fluctuations, large conformational changes, and protein folding, and to reveal molecular mechanisms of their functions at an atomic resolution.1–3 Recent progress in computer hardware and software technology has enabled us to simulate protein motions in a realistic virtual environment on a time scale long enough to compare them directly with experimental results, thereby affording us detailed insight that could not be acquired through physical experiments alone. In particular, the supercomputer Anton, which specializes in MD simulations, is capable of performing them on a scale of millisecond,4 and has been utilized in many successful applications.5–10 One of the greatest advantages of MD simulation is that we can trace the motion of all atoms involved in a simulation system, which is basic and critical information for understanding protein dynamics. On the other hand, not all motions are important or necessary for this understanding, and most of them could be considered to be just insignificant noise. Thus, the extraction and selection of information is essential. However, this process is not always easy because of the massive size of the obtained data arising from the large system sizes and long time scales involved. a) [email protected]

0021-9606/2013/139(21)/215102/10/$30.00

Many analytical methods have been developed to identify significant motions of a protein from its simulation results. The most widely used one is principal component analysis (PCA), by which collective motions of a protein can be detected as low-frequency modes.11–16 It has been found that a small number of these modes are sufficient to describe large fluctuations of a protein and significantly contribute to its conformational change. However, PCA is not necessarily appropriate for representing and examining the dynamic behavior of proteins. To overcome this limitation of the PCA approach, various other methods have been proposed and applied to clarify the dynamics of proteins and peptides.17–24 We also have proposed time-structure based independent component analysis (tICA) as a method for identifying and revealing the slow dynamics of a protein from simulation results, and have validated its usefulness.25 In this previous study, we focused on domain motions of a target protein by using rigid-body domain analysis; this approach reduced the number of simulated protein motions tremendously, from 706 to only 6 degrees of freedom, thereby neglecting many kinds of motions occurring in the protein. However, a protein shows not only domain motions, but also a variety of motions involving main and side chains. If the time scales of these motions are sufficiently slow, then tICA could be used to identify and extract them. Indeed, Schwantes and Pande used tICA to successfully project high-dimensional data of a folding simulation of NTL9 onto lower-dimensional space of slow degrees of freedom, leading to improvement of its Markov state

139, 215102-1

© 2013 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

215102-2

Y. Naritomi and S. Fuchigami

J. Chem. Phys. 139, 215102 (2013)

model.26 However, the usefulness and limitations of the tICA approach are not fully understood. In the present study, we investigate protein backbone dynamics on a long time scale by using an all-atom MD simulation and tICA. For the target protein, we select lysine-, arginine-, ornithine-binding protein (LAO), which undergoes a large conformational change upon ligand binding, and subsequently performed a 1-μs MD simulation of apo-LAO in explicit water. Applying tICA to the trajectory of Cα atoms yields slow modes of the LAO backbone. Protein motions described by these modes are analyzed and characterized to reveal their conformational dynamics in atomic detail. The relevance of these motions to protein function is also discussed.

II. MATERIALS AND METHODS A. Molecular dynamics simulation

We performed an all-atom MD simulation of LAO in explicit water. LAO consists of 238 amino acid residues and 3649 atoms, and contains two domains, the large (residues 1– 90, 191–238) and small (residues 91–190) domains, as shown in Fig. 1(a) (blue and red, respectively). The crystal structure of LAO in the ligand-free open form (PDB ID: 2LAO27 ) was used as the initial structure. After adding missing atoms (including hydrogen atoms), the initial structure of LAO was solvated in a cubic box filled with 25 392 water molecules and three sodium ions to neutralize the net charge of the system.

(a)

The obtained simulation system contained 79 828 atoms in total. MD simulation of LAO in explicit water was performed using the MD program MARBLE28 with the CHARMM22/CMAP force field for the protein and ions29, 30 and the TIP3P model for water molecules.31 Under periodic boundary conditions, van der Waals interactions were switched smoothly to zero over a range of 8–10 Å, and Coulomb interactions were evaluated through the particle mesh Ewald method.32 A symplectic integrator for rigid bodies was used for the system integration to constrain the bond lengths and angles involving hydrogen atoms so that the functional groups –CHx , –NHx (x = 1, 2, 3), –SH, and –OH in the protein and water molecules were treated as rigid bodies.28 A time step of 2 fs was used. After energy minimization by steepest descent, the system was equilibrated for 100 ps in the NPT ensemble (at 1 atm pressure and 300 K temperature) with harmonic restraints (force constants of 1 kcal/mol/Å2 ) on the positions of nonhydrogen atoms of the initial structure. The restraints were then gradually decreased to zero over 100 ps, and further equilibration for 500 ps without restraints was subsequently carried out. A production run was subsequently performed for 1 μs in the NVE ensemble.

B. tICA

Here, we briefly summarize the approach of tICA, which we previously proposed for the characterization of protein slow dynamics.25 When tICA is applied to the real-valued ndimensional time series x(t) = t (x1 (t), · · · , xn (t)), where the superscript t denotes transpose, the following two matrices are required: the covariance matrix C = (x(t) − x(t)) t (x(t) − x(t))

(1)

and the time-lagged covariance matrix with a certain lag time t = t0 ¯ = (x(t) − x(t)) t (x(t + t0 ) − x(t)), C

RMSD [Å]

(b)

(2)

3

where · · · denotes the time average. The value of the lag time t0 should be carefully decided according to the time scales of the dynamics involved. tICA is performed by solving the following generalized eigenvalue problem:

2

¯ = CFK, CF

4

1

0

0

200

400

600

800

1000

Time [ns] FIG. 1. MD simulation of LAO. (a) The initial conformation of LAO based on its crystal structure (PDB ID: 2LAO), which is composed of two domains: the small domain (residues 91–190, red) and large domain (residues 1–90 and 191–238, blue). (b) Time evolutions of RMSD of Cα atoms from the crystal structure for the whole (black), small domain (red), and large domain (blue) of LAO. The 10 ns moving average of RMSD for the whole protein is also shown in brown.

(3)

where K = diag(k1 ,· · ·, kn ) and F = ( f 1 · · · f n ) are eigenvalue and eigenvector matrices, respectively. The resultant eigenvalues and eigenvectors are generally expressed in ¯ is usually asymmetric. For a complex numbers because C ¯ should be symmetric; time-reversible system, however, C hence, we used a symmetrized time-lagged covariance matrix   ¯ + tC ¯ /2, instead of C ¯ to obtain the real-valued ¯ sym = C C eigenvalues and eigenvectors. The obtained eigenvectors are not orthogonal to each other, but satisfy the relation t f i C f j = δij . Thus, { f i } forms a non-orthonormal basis, and the corresponding dual basis {g i } can be defined by g i = C f i to yield the following

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

215102-3

Y. Naritomi and S. Fuchigami

J. Chem. Phys. 139, 215102 (2013)

relation: t

f i g j = δij .

(4)

In this non-orthonormal basis, the time series x(t) is decomposed as   x(t) = g i t f i x(t) = ai (t)g i , (5) i

i

where ai (t) = f i x(t) describes the trajectory along the direction of g i , that is, not f i but g i represents the ith independent component (ICi) of the time series as a mode vector. By using this non-orthogonal basis, the total fluctuation of a protein, which represented by the trace of its covariance matrix TrC of the trajectory, can be decomposed into contributions from each IC mode as follows:    t t gi C f i = gi gi = g i 2 . (6) TrC = t

i

i

i

Thus, the contribution ratio from the ith IC mode is described as g i 2 . (7) TrC The corresponding eigenvalue of the ith IC mode provides information on its time scale, which is roughly estimated by ci =

τi = −

t , ln ki

(8)

where the eigenvalue ki must be a positive number. When the autocorrelation function of ICi exhibits a single exponential decay, the above estimated time scale τ i corresponds exactly to the mean lifetime of the decay.

PCA. PCA is performed by solving an eigenvalue problem of the covariance matrix as follows: C ei = λi ei ,

where λi and ei (i = 1, · · · , 3N ) are eigenvalues and eigenvectors, respectively. Usually, the first 3N − 6 eigenvalues in descending order are non-zero positive, and their corresponding eigenvectors, ei (i = 1, · · · , 3N − 6), form an orthonormal basis for the subspace of the internal degrees of freedom of the protein. The other six eigenvectors, ei (i = 3N − 5, · · · , 3N ), describe the external motions of the protein, and their eigenvalues are practically zero. Thus, the internal dynamics of the protein can be represented as q int (t) = t (e1 , · · · , e3N−6 ) (q  (t) − q  (t)),

Protein conformation is generally represented as a vector in 3N-dimensional Cartesian coordinates q = t (q1 , · · · , q3N ), where N is the number of atoms comprising the protein. Thus, protein dynamics can be described as its time series, q(t). In order to analyze the dynamics of the protein, its external motions are typically eliminated from the time series by superimposing each conformation in the time series on its average structure, yielding a time series of the protein’s internal motions, q  (t). The covariance matrix of q  (t), C = (q  (t) − q  (t)) t (q  (t) − q  (t)),

(9)

has six zero-eigenvalues in principle, which correspond to the external degrees of freedom describing the translation and rotation motions of the protein. This singularity of the covariance matrix causes a practical difficulty in solving the generalized eigenvalue problem in Eq. (3): its standard solver, which is available in software packages such as LAPACK, requires the matrix corresponding to the covariance matrix to be symmetric positive-definite. One of the ways to overcome this difficulty is to focus on only the internal degrees of freedom of the protein and describe protein dynamics in the corresponding (3N − 6)dimensional subspace. This can be accomplished by using

(11)

which is a time series of the (3N − 6)-dimensional vector. The covariance matrix calculated from q int (t) is easily confirmed to be a diagonal matrix composed of the 3N − 6 nonzero positive eigenvalues, i.e., diag(λ1 , · · · , λ3N−6 ), which is obviously symmetric positive-definite. Hence, by analyzing the modified time series of the internal dynamics of the protein, q int (t), instead of the original, q  (t), the aforementioned practical difficulty in solving the generalized eigenvalue problem can be avoided. When tICA of q int (t) is performed, the resultant eigenint vectors f int i and their corresponding vectors g i are represented not in 3N-dimensional Cartesian coordinates but in the (3N − 6)-dimensional space determined by the internal degrees of freedom of the protein, which is defined by ei (i = 1, · · · , 3N − 6). However, these vectors can be easily converted to vectors in 3N-dimensional Cartesian coordinates by basis transformation as follows: h = (e1 , · · · , e3N−6 ) hint ,

C. Practical procedure for tICA

(10)

(12)

int where h and hint are f i and f int i , respectively, or g i and g i , respectively.

D. Domain-motion analysis of IC

When tICA is applied to a trajectory of a protein as described in Sec. II C, an obtained mode vector of ICi, g i , represents protein movement from its average structure. In the case that the protein consists of two domains, A and B, the protein motion described by ICi generally involves both domain motion and conformational change within each domain. To analyze the extent of domain motions in the protein motion described by ICi, we consider a decomposition of g i into , for interdomain motions described by extertwo parts: g inter i , for intradomain nal motions of two rigid domains, and g intra i motions that deform the domains. The decomposition of g i can be performed as follows. First, we focus on a part of g i corresponding to domain A, gA i , which describes the motion of domain A. To detect rigid domain motion of domain A involved in g A i , we then prepare six mode vectors of the external motions for domain A, t A A ˆ j eˆ k = δj k . These vectors can be eˆ A j (j = 1, · · · , 6), where e obtained easily through PCA of domain A, and also by analytical expression of translational and rotational vectors of

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

215102-4

Y. Naritomi and S. Fuchigami

J. Chem. Phys. 139, 215102 (2013)

domain A. By using these mode vectors, rigid domain motion involved in g A i can be estimated as   A t A ˆj . eˆ j · g A rA (13) i = i e j

Thus, conformational change of domain A itself involved in A A A gA i is expressed as s i = g i − r i . For domain B, the same procedure is performed to yield its rigid domain motion r Bi and its conformational change sBi . Finally, compiling the rigid inter B . domain motions of the two domains r A i and r i yields g i intra The intradomain motions deforming the domains, g i , can B be obtained by compiling sA i and s i , and the results satisfy inter intra the relation g i = g i + g i . Accordingly, the ratio of interdomain motions of g i can be calculated by g inter 2 /g i 2 . i III. RESULTS A. Protein backbone motion analyzed by tICA

We performed a 1 μs MD simulation of LAO in its ligand-free open conformation [see Fig. 1(a)] to investigate the slow dynamics of its backbone. Figure 1(b) shows the time evolutions of the root-mean-square displacement (RMSD) of Cα atoms from the crystal structure for whole LAO, its small domain, and its large domain. It indicates that whole LAO shows large and slow fluctuations in the simulation, whereas the two domains are quite stable and rigid. This clearly demonstrates that LAO is not fully stable in the conformation of the crystal structure, and undergoes significant conformational changes primarily by rigid domain motions. In addition, the large domain was found to change its conformation occasionally and temporarily, even though the magnitude of the change is small. These conformational changes of LAO obviously have slow time scales; thus, tICA is expected to reveal their details. To characterize the slow dynamics of the LAO backbone that is hidden in its random fluctuations, we focused on only Cα atoms of LAO and analyzed their trajectory using tICA with a lag time of 1.0 ns. Since LAO consists of 238 amino acid residues, the dimension of the trajectory analyzed is 714 = 238 × 3. After eliminating the six external degrees of freedom as described in Sec. II C, tICA yields 708 ICs as mode vectors for the internal degrees of freedom of LAO in descending order of their time scales as defined by Eq. (8). The time scales of the five slowest modes, IC1 to IC5, are 28.0, 13.4, 10.7, 6.6, and 4.5 ns, respectively, as summarized in Table I. All of these time scales are longer than 1.0 ns, which is the lag time used in the tICA; therefore, their modes are sufficiently robust (as discussed in our previous study25 ), and describe significant motions of the backbone. On the other

hand, the contribution ratios [given by Eq. (7)] of the five ICs are not in descending order, as shown in Table I. Moreover, their values are not as high as those of the collective modes for domain motions, indicating that these ICs might describe not domain motions but instead some local motions. Figures 2(a)–2(e) illustrate the slow time-scale motions of LAO described by IC1 to IC5, respectively. Strikingly, four of the five ICs (all except IC2) were found to have a few Cα atoms that possessed much larger movements compared with those of the other atoms, which were better explained by ICinduced displacements of Cα atoms [Figs. 2(f)–2(j)]. The observed non-uniformity of displacements suggests that these four ICs primarily describe local motions of the backbone around the Cα atoms with remarkable displacements. For clarity, we defined a “remarkably mobile Cα atom” in an IC mode as an Cα atom whose squared displacement contributes more than 10% to the summed squared displacement over all Cα atoms. According to this definition, IC1 has three remarkably mobile Cα atoms: those of residues R218, Q219, and D220, which are located at the end of the α-helix of the large domain, as shown in Fig. 2(a). Large movements of these Cα atoms would be induced by partial unwinding of the α-helix. For IC3, the two Cα atoms of residues A15 and P16 were identified as remarkably mobile: these residues form a turn in the large domain [see Fig. 2(c)]. The local motion described by IC3 might break the hydrogen bond at this turn, leading to unfolding of the local region including these residues. The remarkably mobile Cα atoms for IC4 and IC5 are single Cα atoms belonging to residues G24 and K186, respectively. As shown in Figs. 2(d) and 2(e), the former residue is located at the protruding loop of the large domain, and the latter at the nonsecondary structure region of the small domain, close to the domain boundary. These large displacements of single Cα atoms are expected to be made possible by compensating motions of the backbone, such as the well-known crankshaft motion. However, no remarkably mobile Cα atoms in IC2 were detected, indicating that IC2 describes not a local motion but a collective one. Indeed, the mode vector of IC2 shown in Fig. 2(b) clearly shows a twisting motion of the domains. We further examined the motions described by the ICs by quantifying the extent of rigid domain motions. After calculating the ratios of interdomain motions as described in Sec. II D, IC2 gave the highest ratio among the five slowest modes (Table I), a fact which is consistent with the aforementioned observations. For the other ICs, the ratios of interdomain motions are significantly lower than that of IC2, but are not vanishingly small. This raises the suspicion that the remarkable mobility of the Cα atoms in those ICs does not originate from local motions of the backbone, but is rather artifacts caused by domain motions. If so, then remarkably

TABLE I. Summary of the properties of IC modes. Properties Time scale (ns) Contribution ratio (%) Remarkably mobile Cα atoms Ratio of interdomain motion (%)

IC1

IC2

IC3

IC4

IC5

28.0 5.6 R218, Q219, D220 60.6

13.4 11.9 — 92.1

10.7 2.6 A15, P16 65.7

6.6 1.6 G24 76.7

4.5 0.6 K186 49.0

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

Y. Naritomi and S. Fuchigami

J. Chem. Phys. 139, 215102 (2013)

(a)

(f)

180°

Displacement [Å]

215102-5

1.5

1

0.5

0

50

100

150

200

Residue number

(b)

180°

Displacement [Å]

(g)

1.5

1

0.5

0

50

100

150

200

Residue number

(c)

180°

Displacement [Å]

(h)

1.5

1

0.5

0

50

100

150

200

Residue number

(d)

180°

Displacement [Å]

(i)

1.5

1

0.5

0

50

100

150

200

Residue number

(e)

180°

Displacement [Å]

(j)

1.5

1

0.5

0

50

100

150

200

Residue number FIG. 2. The five slowest modes obtained by tICA: (a)–(e) Mode vectors of IC1 to IC5 depicted by arrows on Cα atoms of the average structure. For clarity, each mode vector is magnified by a factor of 10, and arrows are not shown for atoms whose displacements are smaller than 0.1 Å. Arrows for remarkably mobile Cα atoms are highlighted in magenta. (f)–(j) Displacements of Cα atoms described by IC1 to IC5 (black line) and domain-motion-removed ICs (red and blue lines for the small and large domains, respectively). Magenta-filled diamonds are plotted for remarkably mobile Cα atoms.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

215102-6

Y. Naritomi and S. Fuchigami

J. Chem. Phys. 139, 215102 (2013)

IC2 [Å]

2 0 -2 -4

0

200

400

800

600

1000

Time [ns] FIG. 3. Trajectory of IC2. Its 10 ns moving average is also shown in brown. The corresponding value of the crystal structure of LAO used as the initial conformation is shown in a red broken line.

mobile Cα atoms would not exhibit large displacements in , which are obtained by redomain-motion-removed ICs, g intra i moving rigid domain motions from ICs, g i (see Sec. II D). However, Figs. 2(f)–2(j) clearly indicate that these Cα atoms maintain their large displacements even after removing domain motions, whereas displacements of many other Cα atoms become smaller than the originals. Thus, IC1, IC3, IC4, and IC5 were confirmed to describe significant local motions. The backbone motions detected by tICA could explain the slow time-scale behavior observed in the RMSD trajectories shown in Fig. 1(b), which is considered to arise from do-

IC1 [Å]

(f)

(a)

2

main motion of LAO and deformation of its large domain. As expected, IC2 would best describe domain motions, inducing large and slow fluctuations of LAO, because it has the highest ratio of interdomain motion and is the only IC primarily describing domain motion among the five slowest ICs. This is confirmed by the fact that the trajectory of the IC2 (Fig. 3) clearly exhibits slow dynamical behavior similar to that observed in the RMSD trajectory of whole LAO. This similarity is quantitatively supported by the high correlation coefficient (0.71) between 10 ns moving averages of these trajectories. On the other hand, the large domain can be deformed by the backbone motions induced by IC1, IC3, and IC4. As can be seen in Figs. 4(a), 5(a), and 6(a), the trajectories of these ICs illuminate their distinct dynamical behaviors, and comparisons with the RMSD trajectory of the large domain give strong evidence that IC3 describes the observed domain deformation. As demonstrated above, tICA can adequately detect and extract protein backbone motions on a slow time scale, regardless of whether they are domain motions or local ones. The domain motion described by an IC is easily understood simply by analyzing its mode vector, whereas local motions remain unclear, requiring further investigation to reveal their atomic details. In Secs. III B–III E, we scrutinized the local motions indicated by the four ICs.

D220 G221

0

Q219

Dihedral angle [deg]

-2

T222

(b)

100

R218

0

Y223

-100

(c)

300 200 100

(g)

8

(d)

G221

D220 Q219

Distance [Å]

6 4

T222

2 8

(e)

Y223

6

R218

4 2 0

200

400

600

800

1000

Time [ns] FIG. 4. Local motion of the backbone of LAO detected by IC1: Trajectories of (a) IC1, (b) D220ψ, (c) G221φ, (d) R218O–G221N, and (e) D220O–Y223N. The corresponding values of the crystal structure of LAO used as the initial conformation are shown in red broken lines. Backbone conformations of LAO close to the remarkably mobile Cα atoms of IC1 (f) in the crystal structure and (g) in the final snapshot of the 1 μs simulation. Backbone carbon atoms of the residues having a remarkably mobile Cα atom are indicated in magenta. Related backbone–backbone hydrogen bonds that were stable and unstable in the simulation are shown as cyan and magenta broken lines, respectively.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

Dihedral angle [deg]

IC3 [Å]

215102-7

Y. Naritomi and S. Fuchigami

4 2 0 -2 -4

J. Chem. Phys. 139, 215102 (2013)

F29

(f)

(a)

P16 F17 A15

I27

(b)

300

S18

250

Y14

S19

200

(c)

200

D11

150 100

F29

(g)

8

F17

(d)

P16

Distance [Å]

6 4

I27

2

A15

S19

(e)

8

S18

Y14

6 4 2 0

200

400

600

800

1000

D11

Time [ns] FIG. 5. Local motion of the backbone of LAO detected by IC3: Trajectories of (a) IC3, (b) D11φ, (c) S18ψ, (d) Y14O–F17N, and (e) P16O–F29N. The corresponding values of the crystal structure of LAO used as the initial conformation are shown in red broken lines. Backbone conformations of LAO close to the remarkably mobile Cα atoms of IC3 (f) in the crystal structure and (g) in the snapshot at 470 ns in the temporal conformation. Backbone carbon atoms of the residues having a remarkably mobile Cα atom are indicated in magenta. Related backbone–backbone hydrogen bonds that were stable and unstable in the simulation are shown as cyan and magenta broken lines, respectively.

B. Local motion detected by IC1

IC1, the slowest mode of the LAO backbone detected by tICA, was found to describe the local motion around residues R218, Q219, and D220. In the crystal structure used as the initial conformation of the simulation, the backbone amides of these three residues form hydrogen bonds (R218N–L214O, Q219N–T215O, and D220N–E216O, respectively, where O and N denote oxygen and nitrogen atoms in the backbone) at the C-terminal end of the α-helix [see Fig. 4(f)]. In addition, the backbone carbonyl oxygen atom of R218 forms a hydrogen bond with the backbone amide of G221, i.e., R218O–G221N. Accordingly, the backbone conformation of the residues in question is highly stabilized by these hydrogen bonds. Nevertheless, tICA indicated that some local conformational change occurred around these residues in the simulation. The trajectory of IC1 shown in Fig. 4(a) describes dynamical behavior of the local motion detected by IC1, suggesting that the region in question is stable in the crystal conformation for more than 500 ns, but then (at 545.2 ns) changes its conformation to another one and remains in this new conformation throughout the rest of the simulation. This transition is definitely responsible for the slowest time scale of IC1, even though the time scale estimated by Eq. (8) is

much smaller than expected from the transitional behavior because of the definition of the time scale. By comparing the final snapshot at 1 μs [Fig. 4(g)] with the crystal structure [Fig. 4(f)], we can confirm that the conformational change of the backbone occurs at the depicted region by crankshaft motion of the peptide bond between D220 and G221, in which two neighboring dihedral angles D220ψ and G221φ rotate in concert so as not to cause a large conformational change. In fact, this crankshaft motion induces little change in the position and orientation of the two α helices (residues 204–219 and 222–228) before and after the peptide bond. On the other hand, its movement breaks two hydrogen bonds observed in the crystal conformation (R218O–G221N and D217O– T222N), and induces the formation of a new hydrogen bond, D220O–Y223N, stabilizing the conformation after transition. The transitional behavior of the crankshaft motion is also verifiable by the trajectories of its relevant dihedral angles and hydrogen bonds, as can be seen in Figs. 4(b)–4(e), which show that the transitions completely coincide with that of the IC1 trajectory [Fig. 4(a)]. In addition to the aforementioned transition, the IC1 trajectory shows two other small temporal transitions at around 140 and 470 ns in the opposite direction. Although this transitional behavior is not observed in the trajectories of dihedral angles responsible for the crankshaft motion described

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

215102-8

Y. Naritomi and S. Fuchigami

J. Chem. Phys. 139, 215102 (2013)

6

IC4 [Å]

(e)

(a)

4

S19

2 0

D21

-2

(b)

200

G24

I27

100

E25

0

(c)

100

(f) S19

0 -100

Distance [Å]

Dihedral angle [deg]

300

D21 (d)

6 4

I27

2 0

200

400

600

800

1000

G24

E25

Time [ns] FIG. 6. Local motion of the backbone of LAO detected by IC4: Trajectories of (a) IC4, (b) G24φ, (c) G24ψ, and (d) D21O–G24N. The corresponding values of the crystal structure of LAO used as the initial conformation are shown in red broken lines. Backbone conformations of LAO close to the remarkably mobile Cα atom of IC4 (e) in the crystal structure and (f) in the snapshot at 810 ns in the temporal conformation. Backbone carbon atoms of the residue having a remarkably mobile Cα atom are shown in magenta. Related backbone–backbone hydrogen bonds that were stable and unstable in the simulation are shown as cyan and magenta broken lines, respectively.

by IC1, it is quite similar to behaviors exhibited by the RMSD trajectory of the large domain [blue line in Fig. 1(b)] and the IC3 trajectory [Fig. 5(a)]. Thus, these transitions are considered to be caused not by the crankshaft motion indicated by IC1, but by the local motion around residues 14 and 15 described by IC3, whose atomic details are explored in Sec. III C. This consideration is also supported by the fact that the second-largest peak in the Cα atom displacements of the domain-motion-removed IC1 is at residue 15 [see Fig. 2(f)]. In other words, the dynamical behavior of the IC1 trajectory can be explained by two distinct local motions, which are likely to occur independently but have some dependence on each other. We believe that this dependency is not a real relation, but rather an artifact caused by insufficient sampling of transition events.

C. Local motion detected by IC3

The tICA results for IC3 suggest that some local motion occurred around residues A15 and P16 in the large domain of LAO. As shown in Fig. 5(f), these residues form a β-turn in the crystal structure, the conformation of which is stabilized by four neighboring backbone–backbone hydrogen bonds: Y14N–S18O, Y14O–S18N, Y14O–F17N, and P16O– F29N. If the remarkably mobile Cα atoms of residues A15 and P16 undergo actual movements larger than those of the other Cα atoms, then most of these hydrogen bonds would break, which would lead to a local unfolding of this region.

Figure 5(a) shows the trajectory of IC3, indicating that the local conformational change occurred twice at around 140 and 470 ns. These changes are completely coincident with those observed in the RMSD trajectory of the large domain shown in Fig. 1(b). To elucidate the atomic details of the observed temporal conformational changes, a snapshot at 470 ns [shown in Fig. 5(g)] was compared with the crystal structure [Fig. 5(f)]. The comparison clearly demonstrates that the four hydrogen bonds stabilizing the crystal conformation of this region are all lost in the temporal conformation because of conformational changes in the backbone of not only those residues having remarkably mobile Cα atoms (A15 and P16) but also of residues 11–18. This range is supported by the observation that the trajectories of two dihedral angles, D11φ and S18ψ, represent the same temporal behavior as can be seen in Figs. 5(b) and 5(c). Similar behavior also found in the trajectories of two hydrogen bonds, Y14O–F17N [Fig. 5(d)] and P16O–F29N [Fig. 5(e)], cause local unfolding and refolding of this region. In the temporal conformation, the backbone amide of F29 forms a new hydrogen bond with the backbone carbonyl oxygen of F17 instead of P16 as in the crystal structure. This new hydrogen bond extends the antiparallel β sheet by one residue, thereby changing the orientation of the protruding loop of the large domain (residues 19–27) slightly and resulting in a corresponding peak of Cα atom displacements for IC3 [see Fig. 2(h)]. These local motions described by IC3 are responsible for the conformational change of the large domain observed in their RMSD trajectory.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

Y. Naritomi and S. Fuchigami

Dihedral angle [deg]

IC5 [Å]

215102-9

4 2 0 -2 -4

J. Chem. Phys. 139, 215102 (2013)

(d)

(a)

D187

V185 K186

100

F191

0 -100

(b)

S92

(e) e) D187

300

V185

200 100

(c) 0

200

400

600

800

1000

K186 F191

S92

Time [ns] FIG. 7. Local motion of the backbone of LAO detected by IC5: Trajectories of (a) IC5, (b) K186ψ, and (c) D187φ. The corresponding values of the crystal structure of LAO used as the initial conformation are shown in red broken lines. Backbone conformations of LAO close to the remarkably mobile Cα atom of IC5 (d) in the crystal structure and (e) in the snapshot at 900 ns in the temporal conformation. Backbone carbon atoms of the residue having a remarkably mobile Cα atom are shown in magenta. Related backbone–backbone hydrogen bonds that were stable in the simulation are shown as cyan broken lines.

D. Local motion detected by IC4

E. Local motion detected by IC5

IC4 indicates that the Cα atom of residue G24 could move largely without significant movements of its neighboring Cα atoms. G24 is located at the protruding loop of the large domain, which is stabilized by a β-hairpin formation with four hydrogen bonds (D21O–G24N, D21N–E25O, S19O– I27N, and S19N–I27O) in the crystal structure, as shown in Fig. 6(e). As long as these hydrogen bonds are stably formed, the large movement of the Cα atom of G24 is not possible; therefore, realizing the movement would require the breaking of at least one hydrogen bond, probably D21O–G24N. The dynamical behavior of the local motion described by IC4 can be observed in the trajectory of IC4 shown in Fig. 6(a). Many sharp peaks in the trajectory indicate that the Cα atom of G24 moves temporarily and repeatedly from its position in the crystal structure to another position. Figure 6(f) shows a snapshot at 810 ns of the temporal conformation around G24, which demonstrates large movement of the Cα atom of G24 arising from a backbone conformational change of G24 and its adjacent residues. This change is also supported by the fact that the trajectories of two dihedral angles, G24φ [Fig. 6(b)] and G24ψ [Fig. 6(c)], exhibit temporal transitional behavior that is highly correlated with the IC4 trajectory [Fig. 6(a)]. For this conformational change, the hydrogen bond D21O–G24N is broken in the temporal conformation, as clearly seen in the snapshot at 810 ns [Fig. 6(f)] and the trajectory of distance between D21O and G24N [Fig. 6(d)]. Breaking of the hydrogen bond D21O–G24N further induces temporal breaking of the adjacent hydrogen bond D21N– E25O, which enhances movement of the Cα atom of G24; the other two hydrogen bonds in the β-hairpin are stable during the simulation (data not shown). Thus, the atomic details of the local motion described by IC4 are comprehensible and unsurprising.

IC5 was found to describe the local backbone motion around K186, which is located in the vicinity of the domain boundary. As shown in Fig. 7(d), the backbone of K186 does not form any hydrogen bonds, but two adjacent residues, V185 and D187, do (forming V185N–S92O and D187O– F191N, respectively). Thus, some dihedral angles in the region between these two hydrogen bonds would change collectively, inducing a movement of the Cα atom of K186. The trajectory of IC5 shown in Fig. 7(a) appears to have some small transitions in the minus direction, but they are not as clear as those observed in the IC1, IC3, and IC4 trajectories. On the other hand, the trajectories of two dihedral angles, K186ψ and D187φ, exhibit many temporal transitions, as can be seen in Figs. 7(b) and 7(c), respectively. We found that the transitions in the dihedral trajectories corresponded well with those of the IC5 trajectory. Negatively correlated behavior of the two dihedral angles K186ψ and D187φ indicates that crankshaft motion of the peptide bond between K186 and D187 occurs and changes the backbone conformation locally and temporarily. Indeed, its occurrence is confirmable by observing the backbone conformation around K186 in the snapshot at 900 ns shown in Fig. 7(e). Since the two aforementioned hydrogen bonds in the vicinity of the peptide bond are stable (data not shown), the crankshaft motion can displace the Cα atom of K186 only slightly. Nevertheless, its displacement was detected by the sensitive tICA methodology. IV. DISCUSSION

Our tICA of trajectory of Cα atoms within LAO indicated that not only domain motion but also multiple local motions of its backbone occurred on the slow time scale in the simulation. Four local motions were further investigated to confirm

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13

215102-10

Y. Naritomi and S. Fuchigami

their occurrence and to elucidate their dynamic behavior and atomic details. However, it remains unclear whether these local motions are relevant to protein function, i.e., to specific ligand binding. The structures of LAO when bound with three ligands (lysine, arginine, and ornithine) in the closed conformation have been determined by X-ray crystallography (PDB ID: 1LST, 1LAF, and 1LAH, respectively).33 We explored these liganded structures with special focus on the regions of local motions suggested by tICA, and considered their relations to ligand binding. We found that none of the four regions in any of the three liganded structures were in the conformation induced by the corresponding local motion described by IC1, IC3, IC4, and IC5 [see Figs. 4(g), 5(g), 6(f), and 7(e)]; in fact, all were in the same conformation as those of the unliganded structure utilized for the simulation. Thus, detected local motions do not seem to be directly required for LAO to bind its ligand stably. In the process of ligand binding, however, some of the local motions might yet play a significant role. Most possible of these is the IC3-induced local unfolding, because its region involves the residues D11 and Y14, which interact with the ligand in the closed conformation. The possible role of detected local motions, including the reality of these motions and their slow time scales, could be confirmed by a further experiment on protein dynamics, for example, by employing nuclear magnetic resonance (NMR) spectroscopy. V. CONCLUSIONS

The tICA approach was applied to a 1 μs trajectory of LAO focusing on only Cα atoms to identify slow dynamics of its backbone. We examined protein motion described by the five slowest ICs obtained: only IC2 represented a twisting motion of the two domains of LAO, but the other four ICs indicated local motions of the backbone. Additional investigations confirmed that these local motions actually occurred on the expected slow time scale, and that tICA was successful. This is more evidence that tICA is powerful method for analyzing slow dynamics of a protein, regardless of whether a motion occurs globally or locally and of whether its magnitude is large or small. Some local motions of the LAO backbone detected by tICA were so slow that only few events were observed in the simulation. It is noted that time scales estimated for such local motions would not be sufficiently accurate. Furthermore, such rare events do not always happen in every trajectory. Thus, in order to examine slow dynamics of a protein comprehensively and quantitatively, it will be necessary to analyze a much longer trajectory or a variety of trajectories by the tICA approach. ACKNOWLEDGMENTS

We gratefully thank Professor Akinori Kidera and Associate Professor Mitsunori Ikeguchi for valuable discussion and help. This work was supported by a Grant-in-Aid for Scientific Research to SF (Grant No. 23770180) from the Ministry of Education, Culture, Sports, Science, and Technol-

J. Chem. Phys. 139, 215102 (2013)

ogy of Japan. The simulations were performed at the Department of Medical Life Science, Yokohama City University, and at the Research Center for Computational Science, Okazaki, Japan. 1 J.

L. Klepeis, K. Lindorff-Larsen, R. O. Dror, and D. E. Shaw, Curr. Opin. Struct. Biol. 19, 120 (2009). 2 R. O. Dror, R. M. Dirks, J. P. Grossman, H. Xu, and D. E. Shaw, Annu. Rev. Biophys. 41, 429 (2012). 3 T. J. Lane, D. Shukla, K. A. Beauchamp, and V. S. Pande, Curr. Opin. Struct. Biol. 23, 58 (2013). 4 D. E. Shaw, R. O. Dror, J. K. Salmon, J. P. Grossman, K. M. Mackenzie, J. A. Bank, C. Young, M. M. Deneroff, B. Batson, K. J. Bowers, E. Chow, M. P. Eastwood, D. J. Ierardi, J. L. Klepeis, J. S. Kuskin, R. H. Larson, K. Lindorff-Larsen, P. Maragakis, M. A. Moraes, S. Piana, Y. Shan, and B. Towles, in Proceedings of the Conference on High Performance Computing, Networking, Storage and Analysis (SC09) (ACM, New York, NY, 2009). 5 D. E. Shaw, P. Maragakis, K. Lindroff-Larsen, S. Piana, R. O. Dror, M. P. Eastwood, J. A. Bank, J. M. Jumper, J. K. Salmon, Y. Shan, and W. Wriggers, Science 330, 341 (2010). 6 R. O. Dror, D. H. Arlow, P. Maragakis, T. J. Mildorf, A. C. Pan, H. Xu, D. W. Borhani, and D. E. Shaw, Proc. Natl. Acad. Sci. U.S.A. 108, 18684 (2011). 7 K. Lindorff-Larsen, S. Piana, R. O. Dror, and D. E. Shaw, Science 334, 517 (2011). 8 M. O. Jensen, V. Jogini, D. W. Borhani, A. E. Leffler, R. O. Dror, and D. E. Shaw, Science 336, 229 (2012). 9 A. Arkhipov, Y. Shan, R. Das, N. F. Endres, M. P. Eastwood, D. E. Wemmer, J. Kuriyan, and D. E. Shaw, Cell 152, 557 (2013). 10 S. Piana, K. Lindorff-Larsen, and D. E. Shaw, Proc. Natl. Acad. Sci. U.S.A. 110, 5915 (2013). 11 A. Kitao, F. Hirata, and N. Go, Chem. Phys. 158, 447 (1991). 12 T. Ichiye and M. Karplus, Proteins 11, 205 (1991). 13 A. Amadei, A. B. M. Linssen, and H. J. C. Berendsen, Proteins 17, 412 (1993). 14 S. Hayward and N. Go, Annu. Rev. Phys. Chem. 46, 223–250 (1995). 15 A. Kitao and N. Go, Curr. Opin. Struct. Biol. 9, 164 (1999). 16 H. J. C. Berendsen and S. Hayward, Curr. Opin. Struct. Biol. 10, 165 (2000). 17 K. Moritsugu and A. Kidera, J. Phys. Chem. B 108, 3890 (2004). 18 Y. Matsunaga, S. Fuchigami, and A. Kidera, J. Chem. Phys. 130, 124104 (2009). 19 P. H. Nguyen, Proteins 67, 579 (2007). 20 O. F. Lange and H. Grubmüller, Proteins 70, 1294 (2008). 21 S. Sakuraba, Y. Joti, and A. Kitao, J. Chem. Phys. 133, 185102 (2010). 22 R. I. Cukier, J. Chem. Phys. 135, 225103 (2011). 23 A. Mitsutake, H. Iijima, and H. Takano, J. Chem. Phys. 135, 164102 (2011). 24 T. Nagai, A. Mitsutake, and H. Takano, J. Phys. Soc. Jpn. 82, 023803 (2013). 25 Y. Naritomi and S. Fuchigami, J. Chem. Phys. 134, 065101 (2011). 26 C. R. Schwantes and V. S. Pande, J. Chem. Theory Comput. 9, 2000 (2013). 27 B. H. Oh, J. Pandit, C. H. Kang, K. Nikaido, S. Gokcen, G. F. L. Ames, and S. H. Kim, J. Biol. Chem. 268, 11348 (1993). 28 M. Ikeguchi, J. Comput. Chem. 25, 529 (2004). 29 A. D. MacKerell, Jr., D. Bashford, M. Bellott, R. L. Dunbrack, Jr., J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. JosephMcCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher III, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. WiorkiewiczKuczera, D. Yin, and M. Karplus, J. Phys. Chem. B 102, 3586 (1998). 30 A. D. MacKerell, Jr., M. Feig, and C. L. Brooks III, J. Comput. Chem. 25, 1400 (2004). 31 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). 32 U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995). 33 B. H. Oh, G. F. L. Ames, and S. H. Kim, J. Biol. Chem. 269, 26323 (1994).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.240.225.44 On: Sat, 20 Dec 2014 13:30:13