4.10.1. Diffusion map — MDAnalysis.analysis.diffusionmap

Authors:

Eugen Hruska, John Detlefs

Year:

2016

Copyright:

Lesser GNU Public License v2.1+

This module contains the non-linear dimension reduction method diffusion map. The eigenvectors of a diffusion matrix represent the ‘collective coordinates’ of a molecule; the largest eigenvalues are the more dominant collective coordinates. Assigning physical meaning to the ‘collective coordinates’ is a fundamentally difficult problem. The time complexity of the diffusion map is \(O(N^3)\), where N is the number of frames in the trajectory, and the in-memory storage complexity is \(O(N^2)\). Instead of a single trajectory a sample of protein structures can be used. The sample should be equilibrated, at least locally. The order of the sampled structures in the trajectory is irrelevant.

The Diffusion Map tutorial shows how to use diffusion map for dimension reduction.

More details about diffusion maps are in [1][2][3][4].

4.10.1.1. Diffusion Map tutorial

The example uses files provided as part of the MDAnalysis test suite (in the variables PSF and DCD). This tutorial shows how to use the Diffusion Map class.

First load all modules and test data

import MDAnalysis as mda
import MDAnalysis.analysis.diffusionmap as diffusionmap
from MDAnalysis.tests.datafiles import PSF, DCD

Given a universe or atom group, we can create and eigenvalue decompose the Diffusion Matrix from that trajectory using DiffusionMap and get the corresponding eigenvalues and eigenvectors.

u = mda.Universe(PSF,DCD)

We leave determination of the appropriate scale parameter epsilon to the user, [4] uses a complex method involving the k-nearest-neighbors of a trajectory frame, whereas others simple use a trial-and-error approach with a constant epsilon. Currently, the constant epsilon method is implemented by MDAnalysis.

dmap = diffusionmap.DiffusionMap(u, select='backbone', epsilon=2)
dmap.run()

From here we can perform an embedding onto the k dominant eigenvectors. The non-linearity of the map means there is no explicit relationship between the lower dimensional space and our original trajectory. However, this is an isometry (distance preserving map), which means that points close in the lower dimensional space are close in the higher-dimensional space and vice versa. In order to embed into the most relevant low-dimensional space, there should exist some number of dominant eigenvectors, whose corresponding eigenvalues diminish at a constant rate until falling off, this is referred to as a spectral gap and should be somewhat apparent for a system at equilibrium with a high number of frames.

import matplotlib.pyplot as plt
f, ax = plt.subplots()
upper_limit = # some reasonably high number less than the n_eigenvectors
ax.plot(dmap.eigenvalues[:upper_limit])
ax.set(xlabel ='eigenvalue index', ylabel='eigenvalue')
plt.tight_layout()

From here we can transform into the diffusion space

num_eigenvectors = # some number less than the number of frames after
# inspecting for the spectral gap
fit = dmap.transform(num_eigenvectors, time=1)

It can be difficult to interpret the data, and is left as a task for the user. The diffusion distance between frames i and j is best approximated by the euclidean distance between rows i and j of self.diffusion_space.

4.10.1.2. Classes

class MDAnalysis.analysis.diffusionmap.DiffusionMap(u, epsilon=1, **kwargs)

Non-linear dimension reduction method

Dimension reduction with diffusion mapping of selected structures in a trajectory.

eigenvalues

Eigenvalues of the diffusion map

Type:

array (n_frames,)

run()

Constructs an anisotropic diffusion kernel and performs eigenvalue decomposition on it.

transform(n_eigenvectors, time)

Perform an embedding of a frame into the eigenvectors representing the collective coordinates.

Changed in version 2.2.0: DiffusionMap now also accepts AtomGroup.

Parameters:

• u (MDAnalysis Universe or AtomGroup or DistanceMatrix object.) – Can be a Universe or AtomGroup, in which case one must supply kwargs for the initialization of a DistanceMatrix. Otherwise, this can be a DistanceMatrix already initialized. Either way, this will be made into a diffusion kernel.

• epsilon (Float) – Specifies the method used for the choice of scale parameter in the diffusion map. More information in [2][3][4], Default: 1.

• **kwargs – Parameters to be passed for the initialization of a DistanceMatrix.

class MDAnalysis.analysis.diffusionmap.DistanceMatrix(universe, select='all', metric=<function rmsd>, cutoff=-4.0, weights=None, **kwargs)

Calculate the pairwise distance between each frame in a trajectory using a given metric

A distance matrix can be initialized on its own and used as an initialization argument in DiffusionMap.

Parameters:

• universe (~MDAnalysis.core.universe.Universe or ~MDAnalysis.core.groups.AtomGroup) – The MD Trajectory for dimension reduction, remember that computational cost of eigenvalue decomposition scales at O(N^3) where N is the number of frames. Cost can be reduced by increasing step interval or specifying a start and stop value when calling DistanceMatrix.run().

• select (str, optional) – Any valid selection string for select_atoms() This selection of atoms is used to calculate the RMSD between different frames. Water should be excluded.

• metric (function, optional) – Maps two numpy arrays to a float, is positive definite and symmetric. The API for a metric requires that the arrays must have equal length, and that the function should have weights as an optional argument. Weights give each index value its own weight for the metric calculation over the entire arrays. Default: metric is set to rms.rmsd().

• cutoff (float, optional) – Specify a given cutoff for metric values to be considered equal, Default: 1EO-5

• weights (array, optional) – Weights to be given to coordinates for metric calculation

• verbose (bool, optional) – Show detailed progress of the calculation if set to True; the default is False.

atoms

Selected atoms in trajectory subject to dimension reduction

Type:

~MDAnalysis.core.groups.AtomGroup

results.dist_matrix

Array of all possible ij metric distances between frames in trajectory. This matrix is symmetric with zeros on the diagonal.

Added in version 2.0.0.

Type:

numpy.ndarray, (n_frames, n_frames)

dist_matrix

Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please use results.dist_matrix instead.

Type:

numpy.ndarray, (n_frames, n_frames)

Example

Often, a custom distance matrix could be useful for local epsilon determination or other manipulations on the diffusion map method. The DistanceMatrix exists in diffusionmap and can be passed as an initialization argument for DiffusionMap.

import MDAnalysis as mda
import MDAnalysis.analysis.diffusionmap as diffusionmap
from MDAnalysis.tests.datafiles import PSF, DCD

Now create the distance matrix and pass it as an argument to DiffusionMap.

u = mda.Universe(PSF,DCD) dist_matrix = diffusionmap.DistanceMatrix(u, select=’all’) dist_matrix.run() dmap = diffusionmap.DiffusionMap(dist_matrix) dmap.run()

Changed in version 1.0.0: save() method has been removed. You can use np.save() on DistanceMatrix.results.dist_matrix instead.

Changed in version 2.0.0: dist_matrix is now stored in a MDAnalysis.analysis.base.Results instance.

Changed in version 2.2.0: DistanceMatrix now also accepts AtomGroup.

Changed in version 2.8.0: DistanceMatrix is now correctly works with frames=… parameter (#4432) by iterating over self._sliced_trajectory

References

If you use this Dimension Reduction method in a publication, please cite [2].

If you choose the default metric, this module uses the fast QCP algorithm [5] to calculate the root mean square distance (RMSD) between two coordinate sets (as implemented in MDAnalysis.lib.qcprot.CalcRMSDRotationalMatrix()). When using this module in published work please [5].

[1]

J Porte, Ben Herbst, Willy Hereman, and Stéfan van der Walt. An introduction to diffusion maps. In 11 2008.

[2] (1,2,3)

Ronald R. Coifman and Stéphane Lafon. Diffusion maps. Applied and Computational Harmonic Analysis, 21(1):5–30, 2006. Special Issue: Diffusion Maps and Wavelets. URL: https://www.sciencedirect.com/science/article/pii/S1063520306000546, doi:https://doi.org/10.1016/j.acha.2006.04.006.

[3] (1,2)

Andrew L. Ferguson, Athanassios Z. Panagiotopoulos, Ioannis G. Kevrekidis, and Pablo G. Debenedetti. Nonlinear dimensionality reduction in molecular simulation: the diffusion map approach. Chemical Physics Letters, 509(1):1–11, 2011. URL: https://www.sciencedirect.com/science/article/pii/S0009261411004957, doi:https://doi.org/10.1016/j.cplett.2011.04.066.

[4] (1,2,3)

Mary A. Rohrdanz, Wenwei Zheng, Mauro Maggioni, and Cecilia Clementi. Determination of reaction coordinates via locally scaled diffusion map. The Journal of Chemical Physics, 134(12):124116, 2011. doi:10.1063/1.3569857.

[5] (1,2)

Douglas L. Theobald. Rapid calculation of RMSDs using a quaternion-based characteristic polynomial. Acta Crystallographica Section A, 61(4):478–480, Jul 2005. doi:10.1107/S0108767305015266.