Free-energy landscape, principal component analysis, and structural clustering to identify representative conformations from molecular dynamics simulations: The myoglobin case

Abstract

Several molecular dynamics (MD) simulations were used to sample conformations in the neighborhood of the native structure of holo-myoglobin (holo-Mb), collecting trajectories spanning 0.22 μs at 300 K. Principal component (PCA) and free-energy landscape (FEL) analyses, integrated by cluster analysis, which was performed considering the position and structures of the individual helices of the globin fold, were carried out. The coherence between the different structural clusters and the basins of the FEL, together with the convergence of parameters derived by PCA indicates that an accurate description of the Mb conformational space around the native state was achieved by multiple MD trajectories spanning at least 0.14 μs.

The integration of FEL, PCA, and structural clustering was shown to be a very useful approach to gain an overall view of the conformational landscape accessible to a protein and to identify representative protein substates. This method could be also used to investigate the conformational and dynamical properties of Mb apo-, mutant, or delete versions, in which greater conformational variability is expected and, therefore identification of representative substates from the simulations is relevant to disclose structure–function relationship.

Introduction

Myoglobin (Mb) is a small hemo-protein involved in oxygen transport and storage, which has been widely used for functional and structural studies [1], [2], [3], [4]. Mb belongs to the globin family, which is characterized by a globular fold comprising eight α-helices (A to H, from the N- to the C-terminus) linked to each other by short loop regions [1]. Over the last decades, a number of new globins with different active-site structures has been discovered and characterized, providing interesting variants on the theme [5], [6], [7]. Moreover, the possibility of reversible heme removal and the large number of mutant and deleted versions of Mb available [8], [9], [10], [11], [12], [13], [14], [15] has extended considerably the appeal of myoglobin for structure and dynamics investigation [2], [3].

It is well known that Mb undergoes functionally relevant conformational transitions upon ligand binding, pH variation, but also in the native structure, which has been analyzed by IR spectroscopy [16], [17], NMR [18], [19], electrospray mass spectrometry [20], [21], time-resolved crystallography [22], [23] and, molecular dynamics simulations [24], [25], [26], [27], [28]. Most of these studies converge to the conclusion that full insight into conformational changes of Mb necessitates a detailed description of the conformational ensemble of the protein in the native state [29]. In this context, it is important to underline that proteins are very complex systems and that the energy surface describing the native state contains multiple minima corresponding to very similar, but slightly different conformations, or conformational substates [30], [31], [32], [33].

Proper sampling of protein energy landscape requires the generation of a large sample of molecular conformations and computational methods can be applied to this task [33], [34], [35]. However, analysis of the conformational space, even for a relatively small molecule, is very demanding due to the extremely high dimensionality of the systems [34]. Moreover, a limitation common to all computational sampling procedures is that, although they generate large conformation ensembles, it is hard to assess the extent of sampling [34], [36], [37].

Multiple molecular dynamics (MD) simulations are well suited for generating ensembles of structures [38], [39], [40], [41], [42] at a fixed temperature, and, in fact, they have been successfully used to investigate the conformational landscape near a protein native state [38], [39]. However, it has been observed that independent trajectories of the same system can sample the same local region although following distinct paths, whereas in other cases some trajectories sample more than one of the major local regions for significant time periods. Thus, extending the trajectories for a limited time period does not guarantee more extensive sampling [38], [40]. In addition, the ensemble average structure is typically taken as representative of the trajectory snapshots, even if, in the case of flexible structures, the average structure may not be adequately representative of the ensemble [36]. In fact, few studies have been devoted to a detailed analysis of the contribution of individual protein trajectories to the conformational sampling, and to the identification of representative structures from MD ensembles [34], [43], [44], [45]. Therefore, care has to be used since the average structures resulting from independent trajectories could be very different, and therefore the identification of conformations belonging to low energy basins becomes crucial.

MD simulations can provide an atomic-level picture of protein motion and a representation of the free-energy landscape (FEL) [46], [47] close to the protein native state [43], [44]. However, since it is not possible to represent the FEL as a function of 3N-6 coordinates, it is essential to choose an appropriate set of reaction coordinates that allows to distinguish among different conformational substates. Reaction coordinates are usually obtained by Principal Component Analysis (PCA), which describes the largest amplitude protein motions during a simulation [34], [35], [43], [48], [49], [50]. Once the reaction coordinates are chosen, the FEL can be derived from the probability density function [47], [51].

In fact, since the FEL is projected on a low-dimensional space, it is relevant to verify that conformations that map to the same free-energy basin are also characterized by similar three-dimensional structures. Therefore a further approach for the analysis of the conformational sampling is offered by structural family clustering, a method that collects together conformations according to geometric similarity [34], [51], [52].

In order to explore how to achieve an adequate sampling of the native conformational space through MD simulations and to derive representative structures from the MD ensemble, we performed a set of eleven independent MD simulations at 300 K to sample the conformational space of holo-myoglobin (holo-Mb) in the neighborhood of the native structure, collecting trajectories spanning 0.22 μs. In particular, principal component and FEL analyses, integrated by several cluster analyses, were carried out to gain an overall view of the conformational space accessible to this globular protein. It turns out a high coherence among the different clusters and the basins of the FEL and convergence of parameters derived by PCA, which is not achieved by lower conformational sampling. In particular, our results suggest that multiple MD trajectories spanning 140–160 ns of simulations are sufficient to ensure an adequate sampling of a protein characterized by the globin fold, as myoglobin.