Charting molecular free-energy landscapes with an atlas of collective variables

Charting molecular free-energy landscapes with an atlas of collective variables Behrooz Hashemian, Daniel Mill´an, and Marino Arroyo∗ LaC`aN, Universitat Polit`ecnica de Catalunya–BarcelonaTech, Barcelona, Spain E-mail: [email protected] Abstract Collective variables (CVs) are a fundamental tool to understand molecular flexibility, to compute free energy landscapes, and to enhance sampling in molecular dynamics simulations. However, identifying suitable CVs is challenging, and increasingly addressed with systematic data-driven manifold learning techniques. Here, we provide a flexible framework to model molecular systems in terms of a collection of locally valid and partially overlapping CVs: an atlas of CVs. The specific motivation for such a framework is to enhance the applicability and robustness of CVs based on manifold learning methods, which fail in the presence of periodicities in the underlying conformational manifold. More generally, using an atlas of CVs rather than a single chart may help us better describe different regions of conformational space. We develop the statistical mechanics foundation for our multi-chart description and propose an algorithmic implementation. The resulting atlas of data-based CVs are then used to enhance sampling and compute free energy surfaces in two model systems, alanine dipeptide and β-D-Glucopyranose, whose conformational manifolds have toroidal and spherical topology. ∗To whom correspondence should be addressed 1 1 Introduction Collective variables (CVs) provide a coarse-grained, low-dimensional description of molecular conformations, and thus help us rationalize many molecular systems including atom clusters1–3, small molecules,4–6polypeptides7or proteins.8–10 If CVs are able to separate the metastable states of the system, then they can be used to enhance sampling in molecular simulations and to efficiently compute free energy surfaces.11 Furthermore, since CVs often represent the slowly relaxing degrees of freedom of the system, they may be the basis of reduced models for conformational dynamics.12,13 A poor selection of the CVs, however, mixes metastable states and results in hysteresis and lack of convergence of enhanced sampling algorithms.14,15 Because the identification of suitable CVs for a given system is a major challenge, there has been an effort over the last decade to develop strategies based on machine learning, more specifically dimensionality reduction methods,16 to find systematically such low-dimensional representations. The most widespread dimensionality reduction method is principal component analysis (PCA), a linear method that selects optimal projection directions to maximize variance. PCA is the basis of the essential dynamics approach to analyze collective motions of biomolecules,8,17–20 and can be used to enhance sampling in molecular dynamics.21,22 Linear methods, such as PCA, however, fail to identify nonlinear correlations in datasets of molecular conformation, e.g. as a result of bond rotations or steric interactions.23–25 To better capture the nonlinear low-dimensional geometry of the accessible region in conformational space, the more recently introduced techniques for nonlinear dimensionality reduction (NLDR)26–28 were rapidly applied to understand molecular trajectories, with the underlying hypothesis that these evolve close to a nonlinear manifold often called intrinsic manifold.9,29–31 The idea behind an important subset of NLDR methods is to find a low-dimensional representation or embedding of the set of conformations such that highdimensional distances are preserved as much as possible. The low-dimensional coordinates of the embedding then become the CVs, provided that a differentiable map is available to rep- 2 resent in low-dimensions general conformations (out-of-sample conformations), which are not necessarily those used in the process of identifying the low-dimensional coordinates.32,33 The high-dimensional distances can be Euclidean in classic multidimensional scaling,16 geodesic distances along the intrinsic manifold,9,26 diffusion distances,34,35 or distances transformed by a sigmoid function to emphasize intermediate length-scales.7Beyond the nonlinear manifold model, a NLDR method called sketch map was put forth that views the accessible region of configuration space as a network of basins connected by a spiderweb of transition pathways.1,7This view is not incompatible with the intrinsic manifold model because in principle this network of basins and transition pathways can be embedded in a nonlinear manifold, albeit leaving high-energy regions rarely visited.33 Besides using these methods to analyze well-sampled ensembles, a number of methods have been developed to enhance sampling along nonlinear CVs identified by machine learning.6,32,33,36 In these approaches, the manifold learning algorithm needs to start with a set of conformations representative of molecular flexibility to identify the underlying CVs, a training set in the language of statistical learning, which in itself may be significant challenge. This problem, however, is less severe than that of obtaining thermodynamically meaningful ensembles. Specialized methods of conformational exploration2have been proposed. Furthermore, we showed in previous work that this initial dataset can be transferred from a simpler and easier to sample model, or can be iteratively constructed in combination with enhanced simulation passes.33 Despite promising results, these dimensionality reduction approaches face a fundamental obstacle when dealing with conformational landscapes of non-trivial topology, which are very common, e.g. as a result of cyclic motions around a dihedral angle. Periodicities prevent NLDR algorithms from untangling and displaying coherently in low dimensions the relevant conformational space. Instead, dimensionality reduction methods applied to such systems can mix conformations belonging to distant basins, collapse transition pathways, and severely distort parts of the landscape. This can ruin the value of the resulting CVs. This issue was 3 identified and discussed in Ref. 36, and has plagued a number of applications of NLDR methods to molecular simulation.6,31,35,37,38 In Ref. 39, we further examined this issue, and showed that it manifests itself irrespective of the NLDR algorithm; we considered a nonmetric NLDR method (Locally linear embedding27), and two distance-preserving methods that use different notions of distance (Diffusion map34 and Isomap26). While acknowledging the issue, the authors in Ref. 36 did not provide a specific remedy to these topological obstructions in NLDR, but rather proposed an enhanced sampling algorithm not based on the classical notion of CV to alleviate the consequences of a deficient mapping to lowdimensions. Here, our goal is to develop a general method capable of dealing with systems of arbitrary topology to realize the potential of manifold learning methods for systematic identification of CVs. The idea behind our method comes from the realization that analogous issues are confronted in cartography. The globe is a two-dimensional manifold without boundary. Furthermore, as a result of Gauss’s Egregium Theorem, it is not possible to map even a portion of the surface of a sphere onto a plane without distorting it. These are fundamental topological and geometric hurdles to cartography, which are dealt with by mapping the globe using a collection (atlas) of charts, each describing with appropriate detail a region of limited extent. Each chart has boundary and is planar, just like the output of a NLDR algorithm. By reducing the lateral extent of the region covered by a chart, the metric distortion can be reduced arbitrarily. Here, we examine whether this idea to describe geographical landscapes can be transposed to conformational landscapes. Beyond facilitating the application of manifold learning methods to construct CVs, this approach could have broader applicability. After all, it is natural to expect that a complex molecule may be best described by different CVs (different charts) in different regions of its conformational space, where different steric constraints may determine molecular flexibility. With a framework to model molecular systems with an atlas of CVs, each individual chart could be suitably identified with manifold learning methods, and still describe a coherent conformational landscape. Fur- 4 thermore, the partitioning of configuration space into different regions naturally introduces a length-scale, in that dimensionality reduction methods can only provide correlations within the extent of the corresponding partition. It has been argued that it is important to observe an intermediate length-scale in the NLDR of complex molecular systems.7 However, to our knowledge free energy formalisms and enhanced sampling methods have only been developed considering a single set of CVs. We devote Section 2to revisiting the classical free energy formalisms for a system described by an atlas of multiple partially overlapping CVs. This theory shows how to compute seamlessly thermodynamical observables over a multiple-chart description, and also imposes constraints on the connection between adjacent CVs for a meaningful statistical mechanics description. In Section 3, we describe a specific computational method to implement the atlas of CVs, although the formalism can be applied to other data-driven approaches to build CVs and to other enhanced sampling strategies. Our method is based on the intrinsic manifold model. We identify the nonlinear manifold around which conformations cluster with the Isomap NLDR algorithm.26 The multiple chart description is then built using systematic graph-partitioning algorithms combined with the SandCV methodology presented earlier,33 which provides a differentiable mapping for out-of-sample conformations onto the reduced space of CVs. In section 3, it is described how to extend naturally an enhanced sampling algorithm from a single set of CVs to an atlas of CVs, focusing on the adaptive biasing force method (ABF).40 As a first step to test the applicability of the atlas of CVs framework to molecular modeling, in Section 4 we obtain atlases and perform enhanced sampling simulations to compute the free energy of two model systems: alanine dipeptide and β-D-Glucopyranose. The first of these molecules is a standard benchmark in the field of identification of CVs by statistical/manifold learning. Good CVs in terms of dihedral angles are known for these systems, which provide a basis for comparison with data-driven CVs. Furthermore, the intrinsic manifold on which the conformations evolve has non-trivial topology in both cases (that of a torus and of a sphere), which if improperly dealt with, has led in the past to failures in manifold learning 5 algorithms, see 39 and references therein. 2 Theory Before introducing an atlas description of a molecular system, we review the standard formalism based on a single CV. Consider a N-atom system with separable Hamiltonian and potential energy V(r). Mathematically, a CV is a differentiable mapping taking a conformation rin high-dimensional configurational space D=R3Nand mapping it to ξ, a point in low-dimensional space Rd, where d3N. We denote this mapping by ξ=C1(r). The free energy along this CV, up to an additive constant, can be defined as A1(ξ) = −1 βln ZD1 e−βV(r)δ(C1(r)−ξ)dr,(1) where δ(·) is the Dirac delta distribution,11,41 and 1/β =kT is the Boltzmann constant times temperature. For any given value of the collective variable ξ, the level set L1(ξ) is the set of all molecular conformations rthat are mapped onto ξ, i.e. C1(r) = ξ. Applying the co-area formula,42 the free energy can be written as an integral over the level set L1(ξ) as A1(ξ) = −1 βln ZL1(ξ) e−βVvol(DC1)−1dσξ,(2) where DC1is the Jacobian matrix of the CV, dσξis the volume element of L1(ξ) and vol(DC1) = p|DC1DCT 1|.| · | stands for the determinant. If the CV has been chosen appropriately, then the probability distribution along the level sets L1(ξ) should be not too different from a gaussian, i.e. the system should not exhibit transversal metastability. Figure 1illustrates some of these concepts, with the high-dimensional configurational space represented by a box, and the level set of the scalar CV, ξ=C1(r), represented in red. In black, we represent the intrinsic manifold for this system, which passes close to the energy basin within this red level set. It is clear from the figure that, if the transverse energy 6 D=R3N D1 D2 C1(r)=cst C2(r)=cst C1(r)=⇠ C2(r)='(⇠) Intrinsic manifold Figure 1: Graphical illustration of high-dimensional configuration space D, described in the left by collective variable C1, in the right by C2, and in the central overlapping region by both. The surfaces illustrate the level sets. The level set in the middle can be described either as C1(r) = ξor as C2(r) = ϕ(ξ). landscape was complex with multiple basins, exhibiting transverse metastability, then the notion of intrinsic manifold would become less apparent. The intrinsic manifold is not strictly required in the free energy formalism, but will become useful later. In this simple illustration, configuration space in three-dimensional, there is a single collective variable (d= 1), hence the intrinsic manifold is one-dimensional, and the level sets are two-dimensional (3N−d= 2). For the molecular systems examined later, the intrinsic manifold is two-dimensional (d= 2), whereas the level sets and configuration space are higher-dimensional and cannot be easily visualized. We are now going to use the collective variable C1to describe the system only in a subregion D1⊂Dof the full configurational space, and use a different CV, η=C2(r), in another region of configuration space D2⊂D, which partially overlaps with D1(D1∩D26=∅). Examining Eq. (2), it is clear that the free energies along these two CVs can only be related in the overlapping region if their respective level sets can be mapped to each other. For this to be the case, there must exist a one-to-one transition mapping ϕ:D1⊂Rd→D2⊂Rd such that C2(r) = ϕ◦C1(r) (3) 7 for r∈D1∩D2. If this is the case, then the hyper-surface given by C1(r) = ξcan also be characterized using the second CV as C2(r) = ϕ(ξ). We assume for now that such a transition mapping η=ϕ(ξ) exists and leave for later an explicit construction. This mapping provides a connection between different CVs and a means to build a meaningful global statistical mechanics description of the system. Indeed, we can write Eq. (2) for the second CV as A2(η) = −1 βln ZL2(η) e−βVvol(DC2)−1dση.(4) Using the chain rule on Eq. (3), we have the relation between Jacobian matrices DC2(η) = Dϕ(ξ)DC1(ξ). Furthermore, we note that vol(DϕDC1) = |Dϕ|vol(DC1).43 Thus, recalling that L2(η) = L1(ξ), we can perform a change of variables in Eq. (4) to obtain A2(η) = −1 βln ZL1(ξ) e−βV[|Dϕ|vol(DC1)]−1dσξ =−1 βln ZL1(ξ) e−βVvol(DC1)−1dσξ+1 βln |Dϕ(ξ)| =A1(ξ) + 1 βln |Dϕ(ξ)|.(5) where we have used the fact that Dϕonly depends on ξand hence can be factored outside of the integral. Equation (5) relates the free energies along the two CVs in the overlapping region. This expression highlights the fact that the free energy is subjective in that it depends on the specific parametrization along a collective variable, although observables computed from it are objective.41,44 The last term in Eq. (5) is reminiscent of a Fixman potential.42,45 Recalling that both A1and A2are computed up to an additive constant, Eq. (5) provides a compatibility relation between these constants. Now, let us compute thermodynamic observables, such as relative probabilities of states, in a multiple CV framework. Consider two conformations characterized by regions Aand B 8 belonging to D1and D2respectively. Consider also an auxiliary state C⊂D1∩D2. The relative probabilities between states Aand C, and between Cand Bcan be computed as p(A) p(C)=RC1(A)e−βA1(ξ)dξ RC1(C)e−βA1(ξ)dξ,p(C) p(B)=RC2(C)e−βA2(η)dη RC2(B)e−βA2(η)dη. Consequently, the relative probability between A⊂D1and B⊂D2takes the form p(A) p(B)=RC1(A)e−βA1(ξ)dξ RC2(B)e−βA2(η)dηRC2(C)e−βA2(η)dη RC1(C)e−βA1(ξ)dξ | {z } =1 .(6) Invoking Eq. (5) and using the change of variable formula, it follows that the second term in this equation is 1. Consequently p(A)/p(B) can be computed without reference to state C. This shows that the statistical mechanics of the system can be seamlessly formulated across multiple CVs. The arguments drawn here extend directly to an atlas of multiple partially overlapping CVs. Next, we examine the impact of having multiple CVs on accelerated free energy calculations. In many enhanced sampling methods, such as metadynamics46 or the adaptive biasing force method (ABF),40 an approximation to the thermodynamic force at ξ=C(r) along the CV, denoted by f(ξ)≈ −DA(ξ) (or fi≈ −∂A/∂ξiwith i= 1, . . . , d), is mapped onto a force on the atoms that biases the dynamics given by Fj(r) = d X i=1 ∂Ci ∂rj(r)fi(C(r)) .(7) with j= 1,...,3N, or in matrix notation F(r) = DCT(r) [f◦C(r)] .(8) In the ideal situation, fis close to −DA(ξ), and the enhanced sampling trajectory undergoes nearly free diffusion along the CVs irrespective of free energy barriers. Consider now a two- 9 where f1(ξ) is the biasing thermodynamic force conjugate to ξin the master chart and given by the ABF method. Next, we let all other active charts know about the exerted biasing force, and project this force along the corresponding CVs. Denoting C2a slave chart, we have f2(η) = Dϕ−T(η)f1(ξ).(10) The biasing force f1(ξ) and its representation in slave charts f2(η) are binned to generate a histogram of the mean thermodynamic force. Upon convergence of the algorithm, the free energy is reconstructed in each chart from the histograms as the potential of mean force. A visual summary of the procedure is shown in Figure 5. Thus, in this method there is a single trajectory but one independent free energy surface per chart. Yet, as shown in Section 2, this multiple chart description is not an obstacle to computing thermodynamical observables. 4 Results We demonstrate next the concept of atlas of SandCVs, by implementing it in a standard MD code to perform enhanced sampling simulations and free-energy calculations. MD simulations were performed in version 2.8 of NAMD50 with Langevin thermostat at 300 K. The atlas of SandCVs and the ABF algorithm for enhanced sampling are implemented in C++ and communicate with NAMD through a TCL interface. 4.1 Alanine dipeptide We focus first on alanine dipeptide (N–acetyl–N’–methyl–L–alanylamide) in vacuum. We generate a training set using multiple short high-temperature MD runs33 with the CHARM22 force field.51 The actual partitioning and low-dimensional embedding used to build an atlas of CVs, as described in the previous section, is shown in Figure 4. In Figure 6(a), we show the free energy surfaces (FES) along each of the four CVs systematically identified by the algorithm. The enhanced sampling MD trajectory seamlessly transitioned between charts 16  1 2 3 4 (a) (d) (b) 1 2 3 4 0 10 20 kcal/mol (c) ⇠ ⌘ Figure 6: Analysis with an atlas of CVs of alanine dipeptide. (a) Free energy surfaces computed for each of the four CVs resulting from the systematic partitioning in Figure 4, and obtained by a single trajectory seamlessly transiting between CVs. (b) Coherent juxtaposition of these free energy profiles, possible here because of the special flat geometry of conformational space of this molecule and the isometric character of Isomap. (c) Sampling of two states, C7eq and C7ax, of the molecule. (d) Projection of these states in dihedral space, where we also computed the free energy surface using ABF along the two dihedral angles shown in Figure 4(a). and resulted in converged FES. The length of the trajectory required for convergence was comparable to that needed with a conventional ABF simulation based on backbone dihedrals. A direct comparison of the FES obtained with the atlas of CVs and with dihedral angles [Figure 6(d)] is not strictly meaningful because they differ by Fixman-like terms. However, we can compare later thermodynamical observables obtained from these two different FES. Our approach is obviously more expensive than a conventional ABF simulation based on backbone dihedral angles. On the one hand, and leaving aside the multi-patch aspect of the method, evaluating data-based CVs, like SandCV, is more costly than evaluating an explicit CV. In SandCV, the main computation costs comes from the high-dimensional closest-point projection of a configuration onto the intrinsic manifold (and the derivatives of 17 this operation). However, we previously showed that even this cost becomes insignificant in simulations involving larger systems, e.g. in explicit water.33 On the other hand, there is a computational overhead linked to the multi-patch nature of the method, as compared to a conventional single-chart CV. We found that this specific overhead due to the atlas description is negligible. This can be understood from the fact that the main extra computation is in the evaluation of Eq. (10), a low-dimensional matrix multiplication, and of the partition of unity functions, Eq. (12), which involves calculating Euclidean distances to a few landmark configurations. Thus, the computational overhead of an enhanced simulation along an atlas of data-based CVs will not be significant for a large enough system. Next, we discuss the role of the Fixman correction [Eq. (5)] in this example. Since the statistical error in the evaluation of the free energy is in the order of kT, Eq. (5) shows that the Fixman contribution is significant if |Dϕ|> e, which for d= 2 would require for instance stretching adjacent CV spaces by a factor of 1.6 along each coordinate. Interestingly, because of the nature of Isomap and the fact that the underlying intrinsic manifold is essentially flat (a flat torus),39 the four embeddings for this molecule are nearly isometric. This means that distances between low-dimensional points in Figure 4(f) are similar to distances of the corresponding conformations in high-dimensions [Figure 4(e)] and therefore close to the lowdimensional distances in an adjacent chart. As a result, |Dϕ| ≈ 1 and in this example the Fixman correction is negligible. Thus, we can graphically “glue” the different patches to visualize a globally coherent FES, as illustrated in Figure 6(b). This observation favors building the atlas of CVs using the Isomap method for the NLDR, as compared to methods preserving other notions of high-dimensional distance. We note however that such joint representation is not possible in general, nor is it necessary to compute thermodynamic observables. To show this, we consider two states A and B, computationally described by two ensembles obtained by restrained MD simulations around conformations C7eq and C7ax of the molecule [Figure 6(c)]. For reference, these two states are graphically represented in dihedral space, where we also compute the FES [Figure 6(d)]. 18 When represented in the atlas of CVs [Figure 6(a)], we note that C7eq spreads over each of the four charts. Using Eq. (6) sampled at 1 million configurations for each state, we find using the atlas of CVs that the probability of C7eq relative to C7ax is Patlas = 56.406. To sample integrals spanning multiple charts, we counted each sample only at its master chart. The excellent agreement with the same quantity computed using a separate free-energy surface defined over dihedral space, Pdihed = 56.415, illustrates the accuracy of framework proposed here across multiple CVs, and its potential to study more complex systems with nontrivial topology for which suitable CVs are not readily available. 4.2 Six-membered rings We apply the proposed methodology to a different benchmark molecule, β-D-Glucopyranose,4,5 known to have the topology of the sphere.52 For this six-membered ring, the Cremer-Pople puckering coordinates53 reduce to spherical coordinates (Q, θ, φ), and adequately describe its molecular flexibility with small fluctuations about the “radius” Q. See Appendix Bfor details. To obtain a good training set, we performed a 10 ns metadynamics simulation with φand θas collective variables with Gaussian height of 0.1 and sigma of 0.1, starting from chair conformation (4C1), and using Glycam force field54 at 300 K using Langevin thermostat. We set time step to 0.2 fs and we sampled at each 500 steps (100 fs), obtaining 105 configurations for training set. The systematic partitioning of this training set led here to six different charts, as shown in Figure 7(a) in (θ, φ) and (Q, θ, φ) spaces. Note that two partitions would be enough to remove the topological obstruction. However, the algorithm further subdivided the training set until the quality criteria given in Appendix Awas met. Unlike the previous example, now the underlying two-dimensional conformational manifold is not flat, but rather intrinsically curved. With only two partitions, the low-dimensional embeddings in the plane significantly distort the actual distances between conformations, which in turn introduces significant Fixman terms as in Eq. (5). Instead, applying enhanced sampling along this atlas of six 19 CVs, we find that the transition mappings have very little effect on the free-energy. The resulting FES are shown in Figure 7(b), where the dashed lines delimit the overlap region. Three selected conformations are placed on the FES to illustrate the connections between CVs. Because of the intrinsic curvature of the underlying manifold, it is not possible for this molecule to stitch these six FES along six local CVs without gaps or overlaps. However, as illustrated in the previous example, the atlas of CVs provides a self-standing statistical mechanics description of the system, and it is not necessary to glue these six FES to compute thermodynamics observables. A −1 0 1 −1 0 1 1 −1 0 1 −1 0 1 6 A −1 0 1 −1 0 1 5 A −1 0 1 −1 0 1 4 B −1 0 1 −1 0 1 3 C −1 0 1 −1 0 1 2 C 1 4 5 3 2 6 0 ⇡ 2⇡ ⇡ 2 ⇡ 0 10 20 0 kcal/mol C B (a) (b) Figure 7: Atlas of collective variables for β-D-Glucopyranose, a six-membered ring carbohydrate. The systematic partitioning algorithm for this molecule results in six patches to overcome topological obstructions and alleviate geometric distortion (a), leading to six free-energy profiles (b). 20 5 Conclusion We have shown that molecular conformations and free energy landscapes can be described by an atlas of collective variables, in the same way that geographical landscapes are described by atlas of charts. Furthermore, we have provided a practical algorithm based on SandCVs33 to systematically implement this concept. This method builds a data-driven atlas of CVs identified from machine learning, which can be used to enhance sampling in MD simulations. The proposed framework generalizes the way we model molecular systems, by providing a globally meaningful thermodynamic description from locally appropriate CVs. The proposed framework could be attractive in systems that require different resolution in different regions of conformational space. More generally, the atlas of CVs provides a flexible framework for molecular modeling, which may help to better adapt the description of conformational landscapes to local features. We have shown that this method overcomes some of the main limitations of CVs obtained from machine learning, such as their inability to deal with conformational manifolds of general topology such as those containing periodicities, very common in molecular systems, or the unavoidable geometric distortion introduced when trying to embed globally a curved manifold into a planar CV space. Thus, this method may also expand the applicability of data-driven CVs to systems with arbitrarily complex conformational manifold. We note, however, that an atlas of CVs could be constructed without resorting to machine learning, provided that the conditions identified in Section 2are satisfied. The application of the proposed method to simple yet nontrivial model systems suggests that the method could be applied to more complex molecules for which globally valid collective variables are not available. Such further studies will provide a clearer picture of the potential of the proposed methods. An important question not addressed here is whether a smooth atlas of CVs obtained from machine learning will provide a correct description of kinetic properties, i.e. whether it can be used to identify reaction coordinates. This can be tested using committor histograms.55 Interestingly, even for a simple reaction such as the 21 isomerization of alanine dipeptide, the two dihedral angles (Φ,Ψ) cannot describe correctly the transition and at least three are required.56,57 However, while the two-dimensional description given by our method bears similarities with the two-dihedral description, it takes into account all atomic positions, and thus implicitly all dihedral angles in the molecule. Whether this non-locality inherent to SandCV is useful to correctly capture kinetics will be assessed in future studies. Acknowledgement We acknowledge the support of the European Research Council under the European Community 7th Framework Programme (FP7/2007-2013)/ERC Grant agreement No 240487. We also thank Alejandro Torres for helpful discussions. A Algorithm to geometrically partition the conformational manifold The atlas of collective variables approach relies on a systematic partitioning of the conformational manifold into pieces or patches that are tractable with dimensionality reduction techniques (open sets). Consider a smooth d-manifold Membedded in RDand sampled by a set of points R={r1,r2,...,rN} ⊂ RD, the conformations in the training set. The goal is to represent numerically Mfrom the data in Rthrough a collection of overlapping parametrizations and make computations on it.49 Below are the steps involved in this algorithm: 1. Partition the set of points Rinto Lgroups on the basis of proximity using METIS domain decomposition with a k-nearest neighbor graph.48 METIS tries to partition a graph in equal size subdomains with minimal shared boundary lengths. These Lgroups of points can be represented with index sets Iα, α = 1, . . . , L satisfying ∪L α=1Iα= 22 {1,2, . . . , N}and Iα∩Iκ=∅when α6=κ. 2. For each partition we create an enlarged index set Jαby combining the indices in Iα and the indices from their nearest neighbors defined by a cutoff distance dcut: Jα={a∈ {1, . . . , N}such that |ra−rb|< dcut for some b∈ Iα}. Here and elsewhere, | · | denotes the Euclidean norm. The enlarged index set obeys ∪L α=1Jα={1,2, . . . , N}, but now Jα∩ Jκ6=∅. The intersection of these index sets defines the overlapping between different patches. 3. A dimensionality reduction technique is applied to each one of the enlarged sets {ra}a∈Jα⊂RDto find their low-dimensional embedding {ξa}a∈Jα⊂Rd. In our calculations we use Isomap, a NLDR method designed to find a nearly isometric lowdimensional representation by approximately preserving the geodesic distance on the manifold. 4. The quality of the resulting embedding is measured through a reconstruction error computed as ea=ra−Pb∈Kawab rb |ra|∀a= 1, . . . , N, where Kais the index set with the first k-nearest neighbors of the a-th configuration in the low-dimensional space, and wab are the weights that best linearly reconstruct ξa from its k-nearest neighbors,27 obtained by min wab  ξa−X b∈Ka wab ξb subject to X b∈Ka wab = 1. 5. We monitor in each patch the reconstruction error. We require that max{ea}a∈Jα< Tolefor a patch to be acceptable, where Toleis a numerical tolerance typically below 0.1. If the reconstruction error exceeds the tolerance, the patch is recursively 23 number of patches 1 2 3 4 5 6 relative reconstruction error (%) 0 5 10 15 20 25 30 Average Maximum Figure 8: Average (green circles) and maximum (red squares) relative reconstruction error in the partitioning of alanine dipeptide as a function of the number of patches in the partition. It is clear from this plot that the maximum norm best discriminates the number of required patches given by the algorithm, four in this case. subdivided until the reconstruction requirement is met. In Figure 8, we show the reconstruction error for different number of patches of a training set of alanine dipeptide. 6. After the partitioning process is finished and the low-dimensional embedding for each partition is available, we build a smooth parametrization for each partition using smooth basis functions pa(ξ) associated to nodes from a uniform grid of landmarks, Mα: Ωα⊂Rd−→ RD ξ7−→ X a∈Jα pa(ξ)ya.(11) The high-dimensional control points ya∈RDare chosen such that the reconstruction error is minimized in a least-squares sense. Full details are presented in 33. Thus, we end up with a collection of partially overlapping parametrizations of the intrinsic manifold, Figure 9. We consider here the local maximum-entropy (LME) basis functions.58–60 A partition of unity subordinates to the geometric partition of the conformational mani- 24 (a) (b) (c) Figure 9: Partitioning a training set of alanine dipeptide. The partitioning procedure recursively proceeds until all the patches in the partition are tractable with nonlinear dimensionality reduction (NLDR) methods, which leads to four patches in this example. Although the patches have overlap, using a partition of unity we can unambiguously assign a single master patch to any conformation (a). These patches are also depicted over a torus of dihedral angles to highlight the topology of the intrinsic manifold (b). The low-dimensional embedding of each of the four patches of the partition, calculated separately with NLDR methods, is shown in (c). The darker colors mark the low-dimensional representation of conformations into their master patch. fold is a set of non-negative functions ψα(r) that add up to 1 everywhere, and whose support is contained in the respective patch. Given a set of non-negative reals {βa}a=1,2,...,N , we consider the Shepard partition of unity with Gaussian weight associated to patch αas ψα(r) = Pa∈Iαexp(−βa|r−ra|2) PN b=1 exp(−βb|r−rb|2).(12) These non-negative functions form a partition of unity in RD. For very large βa, these functions tend to the characteristic functions of the Voronoi cells in high-dimension associated to the group of points given by Iα. They can thus be viewed as a smooth regularization of these characteristic functions.49 The four partition of unity functions in the torus representation of alanine dipeptide is illustrated in Figure 10. These partition of unity functions also 25 (45) Fixman, M. Proceedings of the National Academy of Sciences of the United States of America 1974,71, 3050–3053. (46) Laio, A.; Parrinello, M. Proceedings of the National Academy of Sciences of the United States of America 2002,99, 12562–6. (47) Spiwok, V.; Oborsk´y, P.; Paz´urikov´a, J.; Kˇrenek, A.; Kr´alov´a, B. The Journal of Chemical Physics 2015,142, 115101. (48) Karypis, G.; Kumar, V. SIAM Journal on Scientific Computing 1998,20, 359–392. (49) Mill´an, D.; Rosolen, A.; Arroyo, M. International Journal for Numerical Methods in Engineering 2013,93, 685–713. (50) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kal´e, L.; Schulten, K. Journal of Computational Chemistry 2005,26, 1781–802. (51) MacKerell, A. D.; Bashford, D.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wi´orkiewicz-Kuczera, J.; Yin, D.; Karplus, M. The Journal of Physical Chemistry B 1998,102, 3586–3616. (52) Biarn´es, X.; Ard`evol, A.; Planas, A.; Rovira, C.; Laio, A.; Parrinello, M. Journal of the American Chemical Society 2007,129, 10686–93. (53) Cremer, D.; Pople, J. A. Journal of the American Chemical Society 1975,97, 1354– 1358. (54) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; Gonz´alez-Outeiri˜no, J.; 32 Daniels, C. R.; Foley, B. L.; Woods, R. J. Journal of Computational Chemistry 2008, 29, 622–55. (55) Li, W.; Ma, A. Molecular Simulation 2014,40, 784–793. (56) Bolhuis, P. G.; Dellago, C.; Chandler, D. Proceedings of the National Academy of Sciences of the United States of America 2000,97, 5877–82. (57) Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. The Journal of Chemical Physics 2006,125, 24106. (58) Arroyo, M.; Ortiz, M. International Journal for Numerical Methods in Engineering 2006,65, 2167–2202. (59) Rosolen, A.; Mill´an, D.; Arroyo, M. International Journal for Numerical Methods in Engineering 2010,82, 868–895. (60) Mill´an, D.; Rosolen, A.; Arroyo, M. International Journal for Numerical Methods in Engineering 2011,85, 723–751. (61) Sega, M.; Autieri, E.; Pederiva, F. The Journal of Chemical Physics 2009,130, 225102. 33 View publication statsView publication stats