The experimental techniques tend to provide ensemble averaged information and are limited to probing dynamics within narrow windows of time-scales, depending on the instrument resolution. Computational simulations allow bridging multiple time-scales and provide detailed atomistic insights into protein motions. Agarwal and co-workers performed computational studies of cyclophilin A and identified a network of protein residues whose motions influenced the reactive trajectories in the active-site. For ubiquitin, flexibility at ms time-scales have provided some insights into the conformational diversity of how ubiquitin may recognize its binding partners. Similar insights are also available for lysozyme from atomistic simulations; however, it is unclear if these motions translate into transitions between sub-states. Therefore, it would be ideal to simultaneously characterize both the flexibility of the protein and possible transitions enabled by the protein’s flexibility between sub-states that are functionally relevant. The achievable time scales of computational simulations Catharanthine sulfate continue to slowly reach towards biologically relevant time-scales. The large number of conformations sampled during single or multiple molecular dynamics simulations poses a challenge for analysis. Computational tools to analyze and identify conformational sub-states in the multi-level hierarchy that will enable to intuitively understand the biophysical basis of conformational diversity and its relevance to protein function are still limited. The conformations sampled during MD simulations correspond to a highly multi-dimensional data set due to the large number of degrees of freedom associated with the protein. Characterizing the highdimensional multi-variate data, which is embodied in these MD simulations, is a long standing problem in statistics and related fields. Indeed, descriptions of the conformational landscapes spanned by protein motions have typically relied on finding motion directions that can provide biophysically meaningful interpretations. Note, we realize that the conformational ensemble can be projected onto low dimensional representations based on a variety of methods. However, the challenge lies in identifying groups of conformations that provide new insights into the mechanism of protein function. QAA is based on pursuing higher order statistics of positional deviations associated with the conformational data sampled during the MD simulations. Using three different proteins – human ubiquitin, T4 lysozyme and enzyme cyclophilin A – we show that QAA identifies and characterizes the conformational sub-states relevant to function. Based on the inspection of the conformation populations in the sub-states using parameters such as internal energy or other biophysically relevant order parameters, we observe that the identified sub-states contain crucial structural and dynamical elements relevant to promoting the designated function of each of these proteins. A recursive application of QAA yields a multi-level motion hierarchy with global modes dominating the top level and subsequent levels revealing progressively localized motions within the proteins. Additionally, the rare-conformational transitions associated with the interconversion between the identified sub-states allows vital insights into these protein’s structure, motions and function. Individual atoms exhibit significantly anharmonic positional deviations. However, to understand coupling between different protein regions, we examine the joint positional deviations of atom pairs and measure for comparison how a well known approach in the literature, called quasi-harmonic analysis, models the underlying distributions. When the deviations are more Gaussian-like, the QHA basis vectors, which maximize variance, align well with the Ginsenoside-F5 intrinsic orientation of the data. However, when the source distributions combine Gs or Gs, the intrinsic orientations of the data can be non-orthogonal, necessitating higher-order correlations. Under these circumstances, QHA does not capture the intrinsic motions in its sole pursuit of variance.