Relaxation mode analysis for molecular dynamics simulations of proteins.

Molecular dynamics simulation is a powerful method for investigating the structural stability, dynamics, and function of biopolymers at the atomic level. In recent years, it has become possible to perform simulations on time scales of the order of milliseconds using special hardware. However, it is necessary to derive the important factors contributing to structural change or function from the complicated movements of biopolymers obtained from long simulations. Although some analysis methods for protein systems have been developed using increasing simulation times, many of these methods are static in nature (i.e., no information on time). In recent years, dynamic analysis methods have been developed, such as the Markov state model and relaxation mode analysis (RMA), which was introduced based on spin and homopolymer systems. The RMA method approximately extracts slow relaxation modes and rates from trajectories and decomposes the structural fluctuations into slow relaxation modes, which characterize the slow relaxation dynamics of the system. Recently, this method has been applied to biomolecular systems. In this article, we review RMA and its improved versions for protein systems.

Introduction

Molecular dynamics simulation is widely used for protein research. In general, the focus of this research is to extract information on the physical properties of individual proteins. The results from such simulations are then often compared with experimental results. Since these experiments are generally conducted in solvents, it is necessary to simulate protein and water molecular systems, which are complicated systems. These simulations are conducted for a variety of purposes such as to analyze the stability and dynamics of the structures around crystal structures and to determine folding from an extended structure into a native structure. There are three difficulties in current approaches for protein simulations (Freddolino et al. 2010). The first is the potential function of the protein systems. In recent years, it has become possible to evaluate the molecular force field by improving the sampling, and accuracy has consequently improved. The second problem is related to the sampling. With respect to the folding mechanism, simulation at the millisecond scale is necessary. Recently, it has become possible to perform simulations at the millisecond scale by using special hardware such as Anton (Lindorff-Larsen et al. 2011, 2012; Dror et al. 2012, Lane et al. 2013), but sampling problems still exist for complex systems such as ligand-binding systems and other even more complex systems. The third issue is related to the analysis methods. It is important to extract the characteristic degrees of freedom (order parameters) from the complex protein movements obtained from simulations, which are good indicators for analyzing trajectories.

In normal mode analysis, the normal mode near the minimum point of the potential energy of the protein molecule is obtained (Go et al. 1983; Brooks and Karplus 1983; Levitt et al. 1985). Langevin mode analysis investigates modes around the native structure, including the water effect (Lamm and Szabo 1986; Kottalam and Case 1990; Kitao et al. 1991; Hayward et al. 1993). An elastic network model and Gaussian network model approximately estimate normal modes with large amplitudes by using the harmonic potential of coarse-grained models (Tirion 1996; Baher et al. 1997; Tama and Sanejouand 2001; Cui and Bahar 2005; Miyashita and Tama 2008). This method extracts collective modes with large amplitudes in the case of huge protein systems such as viruses, because huge proteins have rigid-like motions (Tama and Brooks III 2002).

Principal component analysis (PCA), also called quasiharmonic analysis or the essential dynamics method (Levy et al. 1984; Ichiye and Karplus 1991; Abagyan and Argos 1992; Garcia 1992; Hayward et al. 1993; Amadei et al. 1993; Kitao and Go 1999), is one of the most popular methods adopted for analyzing the structural fluctuations around the average structure. The modes with large structure fluctuations are extracted and are regarded as cooperative movement, and the relation of these fluctuations with function has been widely investigated. The obtained modes are also used as the axis of the free-energy surface. Moreover, various other analysis methods have been proposed, such as full correlation analysis (Lange and Grubmüller 2007), subspace joint approximate diagonalization of eigenmatrices (Sakuraba et al. 2010), and wavelet analysis (Kamada et al. 2011), among others (Moritsugu et al. 2015; Matsunaga et al. 2015).

In recent years, it has become possible to perform an extensively long simulation; thus, development of dynamic analysis methods to identify the local minimum-energy states and analyze the transitions between them is required. Accordingly, many methods to analyze the dynamics and kinetics of protein simulations have been developed (Zuckerman 2010; Komatsuzaki et al. 2011; Bowman et al. 2014). In particular, the Markov state model has been presented and applied to many protein systems (Schütte et al. 1999; Swope et al. 2004; Singhal et al. 2004; Chodera et al. 2006, 2007; Chodera and Noé 2014; Noé et al. 2007; Noé and Fischer 2008; Noé and Clementi 2017 Buchete and Hummer 2008; Prinz et al. 2011; Pérez-Hernández et al. 2013; Schwantes and Pande 2013; Schwantes et al. 2014; Bowman et al. 2014; Wu et al. 2017). The Markov state model can analyze transitions between local minimum-energy states, which are identified from clustering analysis methods. This is a powerful method for analyzing dynamics in the context of both long and short simulations of proteins.

Relaxation mode analysis (RMA) was developed to investigate the “dynamic” properties of spin systems (Takano and Miyashita 1995) and homopolymer systems for Monte Carlo (Koseki et al. 1997) and molecular dynamics (Hirao et al. 1997) analyses, and has been applied to various polymer systems (Hagita and Takano 2002; Saka and Takano 2008; Iwaoka et al. 2015; Natori and Takano 2017) to investigate their slow relaxation dynamics (de Gennes 1984; Doi and Edwards 1986). Recently, RMA has also been applied to biomolecular systems (Mitsutake et al. 2011; Mitsutake et al. 2005; Mitsutake and Takano 2015; Nagai et al. 2009, 2013). RMA approximately estimates slow relaxation modes and rates from trajectories obtained from simulations.

The relaxation modes {X} satisfy Here, 〈A(t)B(0)〉 denotes the equilibrium correlation of A at time t and B at time 0: where T(Q|Q′) is the conditional probability that the system is in state Q at time t given that it is in state Q′ at time t = 0. Further, Peq(Q′) denotes the probability that the system is in state Q′ at equilibrium. The relaxation rate of X is denoted by λ. The relaxation time is given by 1/λ. Note that the relaxation modes and rates are given as left eigenfunctions and eigenvalues of the time evolution operator of the master equation of the system, respectively, from the viewpoint of the statistical mechanics (Hirao et al. 1997; Koseki et al. 1997; Mitsutake and Takano 2015) (see the “Relaxation modes {” section). The point of RMA is that we consider the variational problem, which is equivalent to the eigenvalue problem of the time evolution operator, and choose an appropriate trial function to estimate the slow relaxation modes and rates in the system (see the “RMA” section). From these processes, we obtain the generalized eigenvalue problem of the time correlation matrices for two different times. From the eigenvectors and eigenvalues, we approximately estimate slow relaxation modes and rates.

Conventional RMA approximately estimates slow relaxation modes by solving the generalized eigenvalue problem of the time correlation matrices of coordinates for two different times, C(τ + t0) and C(t0), which are calculated from the trajectory. Recently, dynamical analysis methods for molecular simulations of biopolymer systems have been developed to investigate slow dynamics. In these techniques such as time structure-based independent component analysis (tICA) (Naritomi and Fuchigami 2011, 2013), time-lagged independent component analysis (TICA) (Pérez-Hernández et al. 2013; Schwantes and Pande 2013), and dynamic component analysis (DCA) (Mori et al. 2015, 2016), time correlation matrices of certain physical quantities or states are used. (Note that tICA is a special case of RMA with t0 = 0. See Mitsutake et al. (2011) and Naritomi and Fuchigami (2011) for more details on the differences between tICA and RMA.) In tICA, TICA, and DCA, the time correlation functions C(τ) and C(0) are used, whereas C(τ + t0) and C(t0) are used in RMA. The relaxation modes and rates are given as left eigenfunctions and eigenvalues of the time evolution operator of the master equation of the system, respectively. From this point of view, RMA is related to Markov state models. (The relationship among the Markov state model, tICA, and TICA is explained in Pérez-Hernández et al. (2013), Schwantes and Pande (2013), and Mitsutake and Takano (2015).) The combination method of tICA and a Markov state model was also proposed (Pérez-Hernández et al. 2013; Schwantes and Pande 2013). A Markov state model was constructed from clustering in the subspace determined by tICA.

In this review, we first provide a definition of relaxation modes and rates from the viewpoint of the statistical mechanics in the “Relaxation modes {” section. The “RMA” section explains the original RMA (RMA with a single evolution time) and the process of RMA using coordinates for the trial function in detail. The “Improvement of RMA” section explains the improved versions of RMA, including RMA with multiple evolution times, principal component RMA (PCRMA), two-step RMA, and Markov-state RMA (MSRMA). Finally, in the “Application of RMA to a system with large conformational changes” section, we present results from studies in which RMA was applied to a system with large conformational changes. The “Conclusions” section provides conclusions and perspectives on the state of the field.

Relaxation modes {X} and rates λ

In this section, we provide the definition of relaxation modes and rates from the viewpoint of the statistical mechanics (Risken 1989; Zwanzig 2001). The relaxation modes {X} satisfy Eq. 1. The relaxation modes and rates are given as left eigenfunctions and eigenvalues of the time evolution operator of the master equation of the system, respectively. We first explain the relation in three types of simulations satisfying the detailed balance condition.

In a Monte Carlo simulation satisfying the detailed balance condition, the time evolution of the probability P(Q;t) that the biomolecule is in a state at time t is described by a master equation: Here, Γ(Q|Q′) denotes the (Q, Q′)-component of the time evolution matrix Γ, and denotes the summation over all possible states. Γ(Q|Q′) is also chosen so that the detailed balance for the equilibrium distribution function Peq(Q) is satisfied:

In the Brownian dynamics simulation, the time evolution of coordinates ,(i = 1,⋯ ,N) is given by the Langevin equation for a biomolecule with N atoms: Here, (t) denotes the position of the i th atom at time t, and ζ is the friction constant. The interaction between atoms is described by the potential U({}) = U(1,…,). The random force (t) acting on the i th atom is a Gaussian white stochastic process and satisfies where w, k, and T denote the α-component of (α = x, y, or z), the Boltzmann constant, and the temperature of the system, respectively. The Smoluchowski equation equivalent to Eq. 5 can be written as Here, Q = {1,…,} denotes a point in the phase space of the system, and P(Q, t)dQ denotes the probability that the system is found at time t in an infinitesimal volume dQ at point Q in the phase space. The time evolution operator Γ satisfies the detailed balance condition (Risken 1989): where . Here, Γ(Q)δ(Q − Q′) and the adjoint operator Γ‡(Q)δ(Q − Q′) act only on Q in δ(Q − Q′). In the matrix representation, so that Γ(Q)δ(Q − Q′) = Γ(Q|Q′) and Γ‡(Q)δ(Q − Q′) = Γ(Q′|Q), the detailed balance condition is the same as that in Eq. 4.

In a molecular dynamics simulation with the Langevin thermostat, the time evolution of coordinates ,(i = 1,⋯ ,N) is given by the Langevin equation for a biomolecule with N atoms: with Here, (t) and (t) denote the position and the velocity of the i th atom at time t, respectively. The mass of the i th atom is denoted by m and ζ is the friction constant.

The Kramers equation, equivalent to Eqs. 9 and 10, can be written as Here, Q = {1,…,,1,…,} denotes a point in the phase space of the system. The time evolution operator Γ satisfies the detailed balance condition: where and Peq(Q) = Peq(𝜖Q). Here, 𝜖Q denotes the time-reversed state of the state Q, namely, 𝜖Q = {𝜖11,…,𝜖,.𝜖1,…,𝜖2} with In the matrix representation, the detailed balance condition is written as follows:

The time evolution equation of P(Q;t) of Eqs. 7 and 11 corresponds to Eq. 3 in the matrix representation. In Monte Carlo and Brownian dynamics, because only coordinates are the degrees of freedom in the system, 𝜖Q = Q, the detailed balance condition in all three cases is given by Eq. 14.

We now consider the eigenvalue problem of the time evolution operator Γ(Q|Q′) of the master equation: Here, ϕ(Q) and ψ(Q) are the left and right eigenfunctions of the time evolution operator Γ with eigenvalue λ, respectively. When we define a quantity through then . The eigenfunctions are chosen to satisfy the orthonormal condition: The equilibrium time-displaced correlation function of ϕ(Q) and is given by the following: where T(Q|Q′) = e−Γ(Q|Q′) is the conditional probability that the system is found at time t at Q given that the system is at Q′ at time 0.

If two quantities A(Q) and B(Q) are expanded as then the time correlation function of A and B in the equilibrium state is given by Thus, in terms of ϕ(Q) and , the correlation function 〈A(t)B(0)〉 is decomposed into a sum of exponentially relaxing contributions. Therefore, we use two sets of functions, {ϕ(Q)} and , as relaxation modes, and refer to {λ} as their relaxation rates. The relaxation modes and rates are given as left eigenfunctions and eigenvalues of the time evolution operator of the master equation of the system, respectively.

RMA

RMA with a single evolution time, t0

RMA approximately estimates slow relaxation modes and rates from trajectories obtained from simulations. Herein, we explain how to obtain the slow relaxation modes and rates. The point of this method is that we consider the variational problem, which is equivalent to the eigenvalue problem of the time evolution operator, and choose an appropriate trial function in order to estimate the slow relaxation modes and rates in the system.

We consider the equations for the conditional probability: The eigenvalue problem in Eqs. 22 and 23 is equivalent to the variational problem with and the stationary value of gives the eigenvalue exp(−λτ). RMA treats the variational problem of Eqs. 24 and 25 using trial functions instead of the eigenvalue problem of Eqs. 22 and 23. To choose the trial function given by a linear combination of important relevant quantities, we can evaluate the relaxation modes and rates from simulation data.

Herein, we consider a biopolymer composed of N atoms and only treat the coordinates, because the velocities have faster relaxations (∼ picosecond order) than coordinates in protein systems. We assume that is a 3N-dimensional column vector that consists of a set of atomic coordinates relative to their average coordinates with where is the coordinate of the i th atom of the biopolymer in the center-of-mass coordinate system, and 〈〉 is its average. Note that because we consider the coordinates only, holds.

In RMA, we use the following function as an approximate relaxation mode: with Here, R(Q) is the i th component of . The quantity R(t;Q) is the expectation value of R after a period t starting from a state Q and satisfies R(t;Q)| = R(Q). The parameter t0 is introduced in order to reduce the relative weight of the faster modes contained in , and it is expected that Eq. 28 becomes a better approximation as t0 becomes larger.

For the trial function (28), defined by Eq. 25 is given by where C(t) is a component of a 3N × 3N symmetric matrix C(t) defined by Then, the variational problem of Eq. 25 becomes a generalized eigenvalue problem The orthonormal condition of Eq. 18 for X is written as Equations 32 and 33 determine the relaxation rates λ and the corresponding relaxation modes f. We chose the indices of λ so that 0 < λ1 ≤ λ2 ≤⋯ holds. Here, the relation which is equivalent to the detailed balance condition of Eq. 14 with 𝜖Q = Q, and the Markovian property are used.

The inverse transformation of Eq. 28 is given by with

The time correlation functions of R are reproduced by for t ≥ t0. Here,

Because we are considering position coordinates only, the detailed balance condition yields the following consequences: C(t) is a symmetric matrix, C(t) = C(t); {λ} are real and positive, which corresponds to pure relaxation. We refer to this method as the “RMA method with a single evolution time,” which is t0/2.

In practice, the time correlation matrices for the two different times are calculated through simulations. Then, by solving the generalized eigenvalue problem, {λ} and {X} are obtained from the eigenvalues and eigenvectors, respectively. To examine the validity of the present analysis, the autocorrelation functions C(t) are reconstructed from the estimated eigenvalues and eigenvectors and are compared with those directly calculated via simulation.

Herein, we comment on the trial function. When RMA was first introduced to a spin system, states of spins on a lattice were used as the trial function (Takano and Miyashita (1995)). When RMA was first introduced to polymer systems, the coordinates of polymers were used as the trial function (Koseki et al. (1997); Hirao et al. (1997)). In polymer systems, the Rouse modes, which were derived from the theory of polymer physics (Doi and Edwards 1986), correspond to the relaxation modes. Rouse modes are given as linear combinations of coordinates. Thus, when RMA was applied to polymer systems, the modes obtained by RMA were compared with the Rouse modes. In protein systems, PCA using coordinates has been widely used. In PCA, the eigenvalue problem of the covariance matrix of coordinates is solved. Therefore, when we first applied RMA to a hetero polymer system (protein system), it seemed to be better to use coordinates as trial functions. The results of RMA and PCA were directly compared with each other. Recently, we have proposed to use physical quantities with slow motions as the trial functions and PCRMA and two-step RMA have been introduced (see the “Improvement of RMA” section). However, RMA using coordinates as the trial functions has an advantage that we can easily convert the information on the slow relaxation modes to the information in coordinate space.

RMA for protein systems

In homopolymer systems, relaxation of the positions of a polymer relative to the center of the mass is investigated. This means that the translational degrees of freedom are removed from the coordinates of the polymer. Because the rotational degrees of freedom remain, the rotational relaxation of the polymer is observed as slow relaxations. In protein systems, it is of interest to evaluate fluctuations of the conformations of a biomolecule around its average conformation. Thus, the translational and rotational degrees of freedom are removed from the sampled conformations of a biomolecule. In practice, treatment of the generalized-eigenvalue problem for removing the translational degrees of freedom in the homopolymer system was given by Koseki et al. (1997). Herein, we explain how to treat the generalized eigenvalue problem for removing the translational and rotational degrees of freedom when using the coordinates for the trial function (Mitsutake et al. 2011). The point of this process is that the generalized eigenvalue problem for real symmetric matrices can be easily solved numerically if the matrices are positive definite. Therefore, we shift the zero eigenvalues to finite positive values without changing the other eigenvalues and the corresponding eigenvectors.

A schematic illustration of the process for RMA using coordinates for the trial function is shown in Fig. 1. First, we remove the translational and rotational degrees of freedom as well as conduct PCA (Eckart 1935; McLachlan 1979). After the average structure converges, the origin of the coordinate system is chosen to be the center of the mass of the average positions, 〈〉 with i = 1,…,N, and the axes of the coordinate system are chosen to be the principal axes of the moment of the inertia tensor of the average positions. We calculate and C′(t): where , , and are unit vectors given by and , , and are unit vectors given by

Fig. 1

Schematic illustration of the RMA process using the coordinate for the trial function

The values of and are usually set to zero. These unit vectors satisfy the following relations: and where α, β = x, y, z and a, b = tr,rot. Then, we solve the generalized eigenvalue problem for C′(t0 + τ) and C′(t0), , with the orthonormal condition . The unit vectors are eigenvectors of this generalized eigenvalue problem with eigenvalues . We denote as the eigenvectors other than . Because , holds. Therefore, are identical with the eigenvectors = (f,f,…,f)T of the generalized-eigenvalue problem for C(t0 + τ) and C(t0) with the same eigenvalues exp(−λτ). Thus, and exp(−λτ) can be obtained by solving the generalized eigenvalue problem for C′(t0 + τ) and C′(t0), which are real symmetric positive definite matrices.

After obtaining relaxation modes and rates, we confirm whether or not the slow relaxation modes and rates obtained using τ and t0 are appropriate. For this purpose, the convergences of slow relaxation times as a function of τ are examined. The autocorrelation functions C(t) are reconstructed from the estimated eigenvalues and eigenvectors and are compared with those directly calculated via simulation (especially the slow relaxation behavior). After examining the validity, we use the obtained relaxation modes and rates for analysis.

Improvement of RMA

Selection of τ and t0 and relevant quantities for the trial function

The relaxation times {1/λ} and the {X} obtained via RMA depend on the manner in which t0 and τ are selected in practice. For simplification, we here consider the case of one physical quantity, R. From the variational problem of Eqs. 24 and 25, the relaxation time 1/λ is obtained from the gradient of the straight line connecting two points at t = t0 and t = t0 + τ in the semi-log plot of the correlation function C(t) = 〈R(t)R(0)〉−〈R〉2 versus t, as shown in Fig. 2a. If the time correlation function of the physical quantity contains several {1/λ}, and if we choose t0 = 0 (tICA case) or a small t0 and small τ, as shown in Fig. 2a (green line), the obtained 1/λ does not correspond to the slow relaxation behavior of log C(t) at long times. To investigate the slow relaxation, we wish to choose values of t0 and τ that are as large as possible, as shown in Fig. 2a (blue line). However, the choice of a longer t0 and τ is also limited, because of the decreasing accuracy of the time correlation function over long time periods. Therefore, we must choose the appropriate t0 and τ.

Fig. 2

Schematic illustration of RMA with a single evolution time t0 (a), and multiple evolution times (1) using t1 and t2 (b) and (2) using t (c)

We can improve the RMA explained above by using two different approaches: introduction of multiple evolution times and using the different relevant physical quantities obtained from coordinates (and velocities) for the trial function. For the first improvement, we describe two types of methods with multiple evolution times, as shown in Fig. 2b, c. (The detailed descriptions are given by Nagai et al. (2013), Natori and Takano (2017), and Karasawa et al. (2017).) For the second improvement, we describe the PCRMA (Nagai et al. 2013), in which the relevant physical quantities for the trial function are given by the PC modes with large structural fluctuations and the two-step RMA (Natori and Takano 2017; Karasawa et al. 2017), which are in turn given by the slowest relaxation modes roughly obtained by RMA. Moreover, the MSRMA (Mitsutake and Takano 2015) is also proposed. We will describe these two improved RMAs in detail below.

RMA with multiple evolution times

RMA with multiple evolution times t1 and t2(1)

In this method, the following trial functions are used as approximate relaxation modes: Note that two evolution times, t1/2 and t2/2, are used instead of a single evolution time, t0/2. Because the contributions of faster modes in time-evolved for t1/2 and those for t2/2 are different, the approximate relaxation modes can extract the faster modes, which cannot be extracted by the approximate relaxation modes using a single evolution time (see Fig. 2b). Using Eq. 45 as a trial function for the variational problem, the following generalized eigenvalue problem is obtained: with . Here, C(t) is a 6N × 6N matrix defined by and is an 3N × 3N matrix defined by where μ1, μ2 = 1 or 2. The orthonormal condition is written as The inverse transformation of Eq. 45 is given by with where . The time correlation functions of R are reproduced by where

RMA with multiple evolution times t (2)

When the relevant physical quantities in the trial function exhibit different relaxations, it is preferable to use different evolution times for the different physical quantities, as shown in Fig. 1c. That is, if we know the characteristic time scales of the relevant physical quantities, we can choose a specific evolution time t for each relevant physical quantity R based on its characteristic time scale. This RMA method is referred to as “RMA with multiple evolution times {t/2}.” In this method, we use the following trial function: The parameter t is introduced in order to reduce the relative weight of the faster modes contained in R. Further, it is expected that Eq. 54 would yield a superior approximation for larger t values.

The variational problem becomes a generalized-eigenvalue problem: Here, C(t) = 〈R(t)R(0)〉 and the orthonormal condition for X is expressed as Equations 54, 55, and 56 determine the relaxation rates λ and the corresponding relaxation modes. We chose the indices of λ such that 0 < λ1 ≤ λ2 ≤⋯ holds. The inverse transformation of Eq. 54 is given by with The time correlation functions of R are given by for t ≥ (t + t)/2. Here,

RMAs to automatically reduce the degrees of freedom of relevant quantities for the trial function

RMA requires relatively high statistical precision of the time correlation matrices because of treatment for the generalized eigenvalue problem; thus, it is difficult for RMA to handle a large number of degrees of freedom directly. We must therefore reduce the number of degrees of freedom automatically.

In an original RMA, the coordinates (and velocity) are used for the trial function. The results may change depending on which relevant quantities are used for the trial function because their correlation functions are fitted using t0 and τ. (For the Markov state model, the dependence of relaxation times on the selection of states is discussed in Swope et al. (2004) and Pérez-Hernández et al. (2013).) It is better to use the relevant quantities that include the slow behavior. For the second improvement, we describe the PCRMA in which the relevant quantities are given by the PC modes with large structural fluctuations, and the two-step RMA in which the quantities are given by the slowest relaxation modes roughly obtained by the first RMA. A schematic illustration of PCRMA and two-step RMA is given in Fig. 3.

Fig. 3

Schematic illustration of PCRMA (a) and two-step RMA (b)

PCRMA

To apply RMA to a protein system by reducing its degrees of freedom, we proposed an improved method, which is referred to as the PCRMA method (Nagai et al. 2013). In this method, PCA is carried out first, and then, RMA is applied to a small number of principal components with large fluctuations (Φ =(). We use the following function as an approximate relaxation mode: Because the degrees of freedom is reduced to N and the relevant quantities with large variance tend to have slow relaxations, the slow relaxation times can be estimated by setting t0 and τ as large values. Note that because the selected principal components also contain faster relaxation modes, as shown in Fig. 4, Nagai et al. (2013) also combined PCRMA with the RMA using multiple evolution times (1) explained above. Note that in PCRMA, if the Nth or more PC modes (with relatively small fluctuations) have slow relaxation, the slow behaviors may not be extracted; thus, there is a possibility that the slow relaxations would not be estimated with small structural fluctuations.

Fig. 4

Schematic illustrationfor PCRMA

Two-step RMA

Using a similar process to that of PCRMA, we proposed a two-step RMA method (Natori and Takano 2017; Karasawa et al. 2017). Based on our experience, the slow {X} obtained from the conventional RMA with small t0 and τ contains the true slow {X} (Mitsutake et al. 2011), although the {1/λ} values are underestimated. The slow relaxation modes obtained by the first RMA may contain the true slow relaxation modes. Thus, we use the slow relaxation modes roughly obtained from the first RMA as the relevant quantities for the trial function. In this technique, RMA with a single evolution time using small t0 and τ is implemented first, and {X} and {λ} are roughly estimated. We then apply the second RMA to a small number of the obtained slowest {X}. We denote the number of {X} used in the second RMA as Nm. In the second RMA, we also use the previously presented technique of RMA with multiple evolution times (2), because the characteristic time scales of the {X} obtained from the first RMA are roughly given by the relaxation times {1/λ}. In the second RMA, we use the following trial function: Here, X(Q) is the relaxation mode obtained from the first RMA and is determined from 1/λ. A detailed explanation is given by Natori and Takano (2017) and Karasawa et al. (2017).

In the second RMA, the time interval τ can be chosen to be large, because the number of degrees of freedom is reduced and the physical quantities {X} exhibit slow relaxations. Using the second RMA, the estimation accuracy of the relaxation modes and times can be improved.

Markov state RMA

As mentioned above, in RMA, the relaxation modes and rates are given as left eigenfunctions and eigenvalues of the time evolution operator of the master equation of the system, respectively. From this point of view, RMA is related to Markov state models. Herein, we consider the relation between RMA and Markov state models and propose the new method of MSRMA.

In the simplest Markov state model, the phase space of the system, where only the position coordinates are considered, is divided into clusters (subsets) S, i = 1,…,n. First, the joint probability that the state of the system Q is in the j th cluster at time 0 and is in the i th cluster at time τ > 0 is calculated in a simulation. Second, the transition probability that the state of the system is found in the i th cluster after time τ starting from a state in the j th cluster is calculated by where is the probability that the state of the system is found in the j th cluster, which is estimated in the simulation. Then, by solving the eigenvalue problem for the transition matrix , the p th eigenvector and its eigenvalue are obtained. The eigenvector corresponds to the equilibrium state and its eigenvalue . Other eigenvectors represent structural transitions and the corresponding eigenvalues give their relaxation time scales as Note that in the Markov description, it is important that the states are defined in a kinetically meaningful way (Swope et al. 2004; Pérez-Hernández et al. 2013). We need to define the states that are classified by order parameters representing the dynamics and kinetics of the system. Even with a good choice of states, in order for a Markov description of the process to be accurate, the time interval τ should also be chosen carefully. In other words, for the Markov description to work, the time interval of the transition matrix τ must be chosen appropriately so that it is as large as the slowest relaxation time of the states. When plotting as a function of τ, slowly converges to the appropriate time scale when τ is increased. In addition, when a much longer τ than the slowest relaxation time of the states is used, the Markov state model is not expected to be accurate. Thus, we usually set the time interval τ to the value when the variation of τ is sufficiently flat (Swope et al. 2004; Pérez-Hernández et al. 2013).

The abovementioned procedure of the Markov state model is related to the following procedure of RMA. We consider an approximate relaxation mode given by where δ(t;Q) is defined in the same way as R(t;Q) in Eq. 29 from δ(Q) given as a function of the state Q of the system by Then, the generalized eigenvalue problem is given by with where According to the definition of δ(Q), it follows that is the joint probability .

If we set t0 = 0, the generalized eigenvalue problem (68) becomes the eigenvalue problem (64) with or , because and . Thus, the Markov state model is a special case of MSRMA with t0 = 0.

Because δ(t0/2;Q) in Eq. 66 reduces the contributions of faster modes in δ(Q), the solutions of the generalized eigenvalue problem (68) provides better approximations to the slow relaxation modes and rates as t0 becomes larger. Therefore, the relaxation times obtained by the Markov state model are expected to be improved by solving Eq. 68 with t0 > 0 rather than Eq. 64.

Application of RMA to a system with large conformational changes

In this section, we apply RMA to a protein system simulation to show the effectiveness of RMA. The selection of order parameters in simulations is important to analyze the trajectory. PCA, which is a static analysis method, extracts large structural fluctuations from simulations, and the obtained PC mode is used to obtain the order parameters. Moreover, it has now become possible to perform long simulations such as those of unfolded and folded protein structures, and when the simulation involves large structural changes, the difference between local minimum-energy states is relatively small compared with that between the folded and unfolded states. In this case, it is difficult for PCA to extract the effective modes or order parameters to accurately identify the local minimum-energy states. By contrast, RMA extracts slow relaxation modes. It is thought that the local minimum-energy states are usually stable so that the system remains in this state for a long time during a simulation. The order parameters with slow relaxation may correspond to the directions between local minimum-energy states. Thus, slow relaxation modes may be suitable order parameters to identify local minimum-energy states and the transitions between them. To validate this concept, we applied RMA to the 10-residue peptide, chignolin in water near its folding transition temperature.

The detailed results are described in Mitsutake et al. (2011). Chignolin consists of a 10-amino acid sequence, GYDPETGTWG and adopts a β-hairpin turn structure (Honda et al. 2004). Several simulations of chignolin have been reported to date (Satoh et al. 2006; Suenaga et al. 2007; Harada and Kitao 2011; Kührova et al. 2012; Okumura 2012). Previous research has shown that chignolin has a stable (native) and a misfolded state, which are both found as hairpin-like structures (see Fig. 5c). These two states have a common turn structure from Asp3 to Glu5 but slightly different hydrogen bond patterns. RMA requires a relatively high level of statistical precision for the time correlation matrices and therefore requires a long simulation where many transitions between local minimum-energy states occur. In addition, we sought to analyze the system with large conformational changes. Thus, we performed a 750-ns molecular dynamics simulation of chignolin in aqueous solution near the transition temperature from an extended structure (Case et al. 2014). We observed many transitions among structures, including the native, misfolded, and unfolded states, by performing the simulation at 450 K. We used the coordinates of C atoms on the backbone as coordinates so that the degrees of freedom were 30. After removing the translational and rotational motions from the coordinates of C atoms, PCA and RMA were carried out on the coordinates of C atoms (see Fig. 1). For RMA, we set t0 and τ to 10.0 and 20.0 ps, respectively.

Fig. 5

The free-energy surfaces for a the first PC mode Φ1 and the second PC mode Φ2, and for b the first slowest RM and the second slowest RM in the case of t0 = 10.0 ps and τ = 20.0 ps. c Snapshots of the native, misfolded, intermediate, and unfolded states classified by RMA, and d distributions for the native (red), misfolded (green), and intermediate (blue) states on the free-energy surface of the first PC mode and the second PC mode. e Relaxation times of the second relaxation mode obtained by MSRMA as a function of the time interval τ. In e, the line of t0 ps corresponds to the results of a simple Markov state model. The figure was reproduced from Mitsutake and Takano (2015)

Figure 5 shows the free-energy surfaces obtained from PCA (a) and RMA (b). From the free-energy surface of PCA, the native and misfolded states were not distinguished because the conformational difference between them is much smaller than the conformational fluctuations of the system (the third PC mode distinguished the native and misfolded states). By contrast, in RMA, the transition between the native and misfolded structures is slow, and the slowest relaxation mode was found to be the axis distinguishing them. This analysis showed that the slow relaxation mode is a good order parameter to distinguish the native and misfolded structure. Interestingly, we could also identify the intermediate structure. By extracting the structures in the center part of the free-energy surface shown in Fig. 5b, the cluster was formed with a turn structure common to the native and misfolded structures. Because the structures at both terminals fluctuate, a cluster of intermediate structures forming a turn is also obtained, while ignoring the fast relaxing movement of both terminals. The upper part of the free-energy surface shown in Fig. 5b corresponded to the extended structure. Figure 5c shows the characteristic structures for the four states. When plotting the points for the obtained intermediate structure on the free-energy surface of PCA in Fig 5d, the points were distributed widely because both terminals fluctuate. Thus, RMA can identify the characteristic structure, even when it is only partially formed. From the free-energy surface obtained by RMA, it is clarified that chignolin folds to the native or misfolded structures through the intermediate (turn) structure from the extended structures.

Because the structures were classified into a smaller number of states using the free-energy surface obtained by RMA, we then applied the Markov state model and MSRMA to analyze these four states: native, misfolded, intermediate, and unfolded states. Figure 5e shows the relaxation time τ = 1/λ obtained by MSRMA as a function of τ when t0 = 0,10,50,100,200, and 500 ps. Because the first eigenvector corresponds to the steady state with infinite relaxation time τ1 = ∞, we show the second slowest relaxation times. The line of t0 = 0 corresponds to the results of a simple Markov state model. In the case of t0 = 0, the τ values slowly approach the appropriate time scale, i.e., the values for plateau regions or peak values of the solid lines, when τ is increased. For the lines of t0 > 0, the values of τ quickly approach the appropriate time scale, i.e., those corresponding to the values for plateau regions or peak values. Thus, the slow relaxation times can be improved when applying MSRMA with t0 > 0, which is introduced to reduce the relative weight of the faster modes.

Overall, RMA can be used to effectively analyze long simulations at room temperature and is also useful for investigating systems with large conformational changes, such as intrinsically disordered proteins and protein folding.

Conclusions

In this paper, we have reviewed the method and application of RMA, a dynamic analysis method for protein simulations. We described the definition of relaxation modes and rates, which correspond to the left eigenfunctions and eigenvalues of the time evolution operator of the master equation of the system, respectively. After providing the definition, we explained how to estimate the slow relaxation modes and rates from simulation data. We also summarized several new RMAs proposed, including RMA with multiple evolution times, PCRMA, two-step RMA, and MSRMA. Finally, to demonstrate the effectiveness of RMA, we briefly presented the analysis results of the unfolding/folding simulation of the 10-residue peptide chignolin detected near the transition temperature. The simulation results showed that the relaxation mode is a good order parameter for not only extracting the transition between the native state and misfolded state but also for identifying the intermediate state, which is partially folded. This suggests that RMA is suitable to investigate a system with large structural changes and naturally denatured protein systems. Although RMA is efficient for a longer simulation than the longest relaxation time of the system, it can also extract rare events in a finite-time simulation such as that conducted at the microsecond scale. By examining the extent to which the correlation function can be reconstructed, we can clarify the information that can be obtained on dynamics using the obtained relaxation modes and rates. Theoretical studies to compare data of the Markov state model with experimental data from nuclear magnetic resonance and neutron scattering analyses have emerged recently (Xia et al. 2013; Lindner et al. 2013; Zheng et al. 2013; Bowman et al. 2014). In the future, it will also be important to interpret the theoretical relationships in light of experimental data.