516
T. Borowski and M. Szaleniec
open-shell transition metal systems as QM(/MM) dynamics employing a reliable QM
model and basis set is computationally too demanding for investigations of reaction
path. Note, however, that such studies are perfectly doable with semi-empirical QM
methods for systems without transition metals [70]. How such a snapshot(s) should
be selected is not obvious, and there are at least several approaches used by the community. One way is to take the last saved snapshot from the simulated trajectory, as
this can be viewed as the “final” product of MD simulation. In another approach,
one arbitrarily selects specific frames from the saved trajectory. Both of the above
approaches may lead to selection of suboptimal structures because at physiological
temperatures, simulated with MD, energy of the system fluctuates quite substantially
and the trajectory visits various regions of the conformational space, also those characterized by higher free energy values and hence lower probability. On the other
hand, reaction mechanisms of most metalloenzymes proceed through very reactive
intermediates, e.g. Fe(II)-OO-R, Fe(IV)=O, and for this reason these enzymes usually function according to a “negative catalysis” paradigm [71, 72], which means
the enzyme holds the substrate(s) in a well-defined position and conformation, so
that only one reaction pathway, i.e. the desired one, is realized. Hence, it seems
reasonable to assume that, at least at this stage of the reaction when the reactive
species is present, the free energy landscape of the enzyme–substrate (or intermediate) complex will feature a dominating basin of attraction corresponding to the
preferred orientation of the reactants. Therefore, another way of selecting snapshots
from the trajectory rests on methods capable to identify such basins, i.e. on clustering algorithms. Objective of the latter is to divide a given set of objects, in our
case MD snapshots, into groups (clusters) that gather most similar objects under
the condition that differences between the clusters are much more pronounced than
variance within them. Arguably, the major basin of attraction should correspond to
the most populated cluster. An example of cluster analysis comes from our recently
published work on acireductone dioxygenase with the key settings and major results
presented in Fig. 6 [73]. A total of 500 snapshots saved during the last 10 ns of
a stable “production” period of MD were first superimposed minimizing RMSD
for selected atoms of residues lying within 10 Å of the metal. The specific atoms
used for RMSD calculations are marked with orange spheres in the figure; note that
symmetry-equivalent atoms were not taken into account as using them would lead to
discrimination between physically indistinguishable conformations. Second, using
this RMSD metrics, clustering algorithms that are available in the cpptraj program
from the Amber 16 suite were used to cluster the frames. For K-means and hierarchical agglomerative algorithms, the required number of clusters (n) can be explicitly
specified and hence, with their use clustering was performed for n ranging from 2
to 6. For the density-based algorithms (DBSCAN and dpeaks), one needs to input
either distance cut-off between points for forming a cluster (DBSCAN) or cut-offs
determining local density and distance (dpeaks). The choice of the appropriate values
of these parameters was aided by visual inspection of a “K-dist” plot (DBSCAN)
and the “decision graph” (dpeaks).
The decision graph for the dpeaks algorithm clearly showed 3 or 4 distinctive
points characterized by high density and large distance to the point of next highest
T. Borowski and M. Szaleniec
open-shell transition metal systems as QM(/MM) dynamics employing a reliable QM
model and basis set is computationally too demanding for investigations of reaction
path. Note, however, that such studies are perfectly doable with semi-empirical QM
methods for systems without transition metals [70]. How such a snapshot(s) should
be selected is not obvious, and there are at least several approaches used by the community. One way is to take the last saved snapshot from the simulated trajectory, as
this can be viewed as the “final” product of MD simulation. In another approach,
one arbitrarily selects specific frames from the saved trajectory. Both of the above
approaches may lead to selection of suboptimal structures because at physiological
temperatures, simulated with MD, energy of the system fluctuates quite substantially
and the trajectory visits various regions of the conformational space, also those characterized by higher free energy values and hence lower probability. On the other
hand, reaction mechanisms of most metalloenzymes proceed through very reactive
intermediates, e.g. Fe(II)-OO-R, Fe(IV)=O, and for this reason these enzymes usually function according to a “negative catalysis” paradigm [71, 72], which means
the enzyme holds the substrate(s) in a well-defined position and conformation, so
that only one reaction pathway, i.e. the desired one, is realized. Hence, it seems
reasonable to assume that, at least at this stage of the reaction when the reactive
species is present, the free energy landscape of the enzyme–substrate (or intermediate) complex will feature a dominating basin of attraction corresponding to the
preferred orientation of the reactants. Therefore, another way of selecting snapshots
from the trajectory rests on methods capable to identify such basins, i.e. on clustering algorithms. Objective of the latter is to divide a given set of objects, in our
case MD snapshots, into groups (clusters) that gather most similar objects under
the condition that differences between the clusters are much more pronounced than
variance within them. Arguably, the major basin of attraction should correspond to
the most populated cluster. An example of cluster analysis comes from our recently
published work on acireductone dioxygenase with the key settings and major results
presented in Fig. 6 [73]. A total of 500 snapshots saved during the last 10 ns of
a stable “production” period of MD were first superimposed minimizing RMSD
for selected atoms of residues lying within 10 Å of the metal. The specific atoms
used for RMSD calculations are marked with orange spheres in the figure; note that
symmetry-equivalent atoms were not taken into account as using them would lead to
discrimination between physically indistinguishable conformations. Second, using
this RMSD metrics, clustering algorithms that are available in the cpptraj program
from the Amber 16 suite were used to cluster the frames. For K-means and hierarchical agglomerative algorithms, the required number of clusters (n) can be explicitly
specified and hence, with their use clustering was performed for n ranging from 2
to 6. For the density-based algorithms (DBSCAN and dpeaks), one needs to input
either distance cut-off between points for forming a cluster (DBSCAN) or cut-offs
determining local density and distance (dpeaks). The choice of the appropriate values
of these parameters was aided by visual inspection of a “K-dist” plot (DBSCAN)
and the “decision graph” (dpeaks).
The decision graph for the dpeaks algorithm clearly showed 3 or 4 distinctive
points characterized by high density and large distance to the point of next highest