Template C++ / Python API for developping structural bioinformatics applications.
User Manual
Authors:F. Cazals and T. Dreyfus and D. Mazauric
Energy_landscape_comparison
This package package provides method to perform various comparisons of two samples energy landscapes, taking into account the location of local minima, their occupancy probabilities, and possible their connexsions as encoded in a transition graph.
Goals
Comparing two (sampled) energy landscapes is of interest in various settings, e.g., to assess the coherence of two force fields for a given system (atomic, coarse grained), to compare two related systems (e.g. a wild type and mutant protein), or simply to compare simulations launched with different initial conditions (and check whether the same regions in conformational space have been visited).
In comparing two landscapes, two categories of criteria are of interest, namely
features of the basins, in particular the local minima and their associated occupancy probabilities, called masses for the sake of genericity in the sequel.
transitions between these basins.
This package provides methods to compare (sampled) energy landscapes, in two guises:
Earth Mover Distance: a comparison method exploiting solely the location of local minima, and the masses of the basins.
EMD with connectivity constraints: a comparison method also taking into account the connectivity of the basins.
These functionalities are available in the programs and .
Prerequisites
Landscapes modeled as transition graphs
Landscapes and vertex weighted TG. In the sequel, we assume that the energy landscape (EL) is coded in a compressed transition graphs, as defined in section Energy landscapes . We also assume that each vertex is endowed with a mass, typically the volume or the occupancy probability of its catchment basin. The reader is also referred to the package dedicated to the construction of transition graph, see Transition_graph_of_energy_landscape_builders .
In the following we provide two comparison methods: the first one deals with features of the basins only, while the second one additionally exploits the information on transitions.
Basins and their masses.Consider the basin of a local minimum. Abusing notations, the footprint of the basin in the conformational space is also denoted . Denoting is the Boltzmann constant and the partition function.
For small systems, the mass of the system may be obtained by integrating the Boltzmann factor over the basin region, namely:
If minima are found by optimization (e.g., basin hopping), the can be estimated using the eigenvalues of the Hessian (curvature) matrix evaluated at those points [188] . Such a calculation amounts to focusing on the vibrational entropy of the system.
If, on the other hand, the samples are obtained from a thermodynamic ensemble, as is typically the case for molecular dynamics and Monte Carlo procedures, Boltzmann weighting is automatically satisfied; the basin weight can correspondingly be estimated from the number of points falling in the basin region:
Earth mover distance
Consider two landscapes, for which masses of the basins have been computed.
For the sake of exposure, we call these two landscapes the source landscape and the demand landscape . We also denote and , respectively, their numbers of basins.
The local minima and the associated basins of the source landscape are denoted and , respectively. For the demand landscape, local minima and basins are denoted and .
We also assume that the weights of the basins have been computed. Denoting these and for a source and demand basin respectively, we define the sum of weights and . Finally, so that the distance between the two aforementioned local minima is .
To compare the two landscapes, we use the earth mover distance[165] , which is a particular case of mass transportation [184] . Intuitively, the technique fills basins in the target (aka demand) landscape using mass from basins of the source landscape. A basin from can be split into several parts, and equivalently, a basin from can be filled from several basins from (Fig. fig-comparison-of-basins-energy-landscapes). Denote the mass from moved into . The cost of moving units of mass depends linearly on the distances between the minima associated with and associated with . A transport plan from to is defined by triples , with . Note that there are at most such triples.
Finding the optimal i.e. least cost transport plan amounts to solving the following linear program (LP):
The first equation is the linear functional to be minimized, while the remaining ones define linear constraints. In particular, the second one expresses the fact that every basin from need to be filled, while the third one indicates that a basin from cannot provide more than it contains. (To simplify matters, we have assumed that . Handling the case poses no difficulty, and the reader is referred to [165] .)
Based on this linear program, we introduce the total number of edges, the total flow, the total cost, and their ratio, known as the earth mover distance[165] :
If the sum of basins from both landscapes match i.e. are identical, then the LP is symmetric. If not, one should consider solving the two LP.
Comparing two energy landscapes and
The landscape is partitioned into the basins associated to its local minima (three of them on this example), and likewise for (four local minima). Comparing the landscapes is phrased as a mass transportation problem on the bipartite graph defined by the two sets of minima. Note that sets of connected basins from the top landscape are mapped to connected basins of the bottom landscape.
Earth mover distance with connectivity constraints
The previous comparison ignores transitions between local minima. To take these connections into account, we modify the method by imposing connectivity constraints to transport plans. To see whether a transport plan is valid, pick any connected subgraph from – that is connects selected local minima in . Let be the set of vertices from such that for each vertex in , there is at least one edge emanating from a vertex of the subgraph with . The transport plan is called valid iff the subgraph is also connected. We summarize this discussion with the following definition–see also (Fig. fig-emd-LP-not-edge-connectivity-violated):
A transport plan is said to respect connectivity constraints iff any connected set of basins from induces, through edges carrying strictly positive flow, a connected set of basins from .
An important remark is the following: a transport plan respecting connectivity constraints may not fully satisfy the demand. (It can actually be shown that there exist instance such that no transport plan fully satisfies the demand.)
Since EMD-CC does not necessarily admit a solution fully satisfying the demand, we define the problem Earth Mover Distance with Cost and Connectivity Constraints :
The problem EMD-CCC (maximum-flow under cost and connectivity constraints problem) aims at computing the largest volume of flow that can be supported respecting the connectivity constraints, and such that the total cost is less than a given bound C.
Note that the previous definition calls for two algorithms:
algorithm : computes a transport plan given a upper bound on the transport cost.
algorithm : computes transport plans whose cost lies in a given interval.
The latter algorithm actually has 2 recursion modes. To see which, assume that an interval of possible costs is given. The two modes are:
Recursion mode: refined. Given two costs and , is called three times for the different costs , , and . Then:
a) If the two total volumes of flow with cost and are different, then is called with the interval .
b) If the two total volumes of flow with cost and are different, then is called with the interval .
Recursion mode: coarse. Similar to refined mode, except that only option b) is considered.
The solution of the linear program may not satisfy connectivity constraints
A transport plan between a source and a demand graph each consisting of a linear chain of four vertices. The vertices of the edge of the source graph export towards the vertices and of the demand graph. The subgraph of the demand graph induced by these vertices is not connected – there is no edge linking to .
Using sbl-energy-landscape-comparison-euclid.exe and sbl-energy-landscape-comparison-lrmsd.exe
In the following, we specify the two programs sbl-energy-landscape-comparison-euclid.exe and sbl-energy-landscape-comparison-lrmsd.exe, which differ by the type of metric used to compare distances between the points associated to vertices defining the graphs.
Main specifications
Main options.
The main options of the program sbl-energy-landscape-comparison-euclid.exe are: –transition-graphstring: transition graph XML archives (used twice for target and source transition graphs) –with-connectivity-constraintsstring: run Earth Mover Distance with connectivity constraints –symmetric-modestring(= when using connectivity constraints): moving also from target to source
(xml format, suffix emd_engine.xml) Transport plan in xml format, from first landscape (say A) to the second landscape (say B)
(xml format, suffix engine_symmetric.xml) As above, but for the transport plan from B to A
(dot format, suffix transportation_plan.dot) Transportation plan in dot format, to be rendered with graphviz (dot, circo)
(dot format, suffix transportation_plan_symmetric.dot) As above, but for the transport plan from B to A
Optional output. Comments.
Algorithms and Methods
Earth mover distance
Once the weights of basins have been computed, solving the linear program of Eq. (eq-emd-lp) has polynomial complexity [115] . Practically, various solvers can be used, e.g. the one from the Computational Geometry Algorithms Library [52] , lp_solve, the CPLEX solver from IBM, etc. In the following, the algorithm solving the linear program of Eq. (eq-emd-lp) is called .
Earth mover distance with connectivity constraints
Finding transport plans respecting connectivity constraints turns out to be a hard combinatorial problem [40] . The problem is not in APX, which means that if holds, then, no polynomial algorithm with constant approximation factor exist.
Following definition def-tp-cc, we provide two algorithms:
: computes a transport plan respecting connectivity constraints, given a bound C.
: computes different transport plans for different bounds not given by the user. More precisely, since the maximum cost is upper-bounded by , computes different flow solutions by calling for different costs in the range [0,C-max] in a loop. There are two different modes for determining such set of costs.
However, a greedy algorithm providing admissible solutions respecting connectivity constraints, denoted for earth mover distance with cost and connectivity constraints, has been reported in [40] .
Using the previous algorithm(s), in a manner analogous to Eq. (eq-emd-lp-dist), we define the total number of edges, the total flow, the total cost, and their ratio:
It can be shown that the Earth Mover Distance as defined by Eq. (eq-emd-lp) yields a metric, provided that is itself a metric, and that the sum of weights for the source and the demand graphs are equal [165] . On the other hand, the EMD with connectivity constraints fails to satisfy the triangle inequality [40] . The EMD with connectivity constraints is not symmetric either, but is easily made so in taking a symmetric function of the one sided quantities. To assess this lack of symmetry, we introduce the following ratios, respectively geared towards the flow and the cost:
Programmer's Workflow
The programs of Energy_landscape_comparison described above are based upon generic C++ classes, so that additional versions can easily be developed.
In order to derive other versions, there are two important ingredients, that are the workflow class, and its traits class.
Note that is this case, the transportation plan matches that observed for the plain EMD
In [10]:
sbl_exe_name="sbl-energy-landscape-comparison-euclid.exe"tg1="data/himmelbleau_tg.xml"tg2="data/himmelbleau_transition_graph_noisy.xml"lp_solver=None# since EMD is usedsymmetric=TrueconnectivityConstraints=Truecmp_landscapes(sbl_exe_name,tg1,tg2,lp_solver,symmetric,connectivityConstraints,False)
Using executable /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-energy-landscape-comparison-euclid.exe
Executing /user/fcazals/home/projects/proj-soft/sbl-install/bin/sbl-energy-landscape-comparison-euclid.exe --transition-graph data/himmelbleau_tg.xml --transition-graph data/himmelbleau_transition_graph_noisy.xml --directory tmp-results-sym-True-cc-True -o --log --symmetric-mode --with-connectivity-constraints
All output files in tmp-results-sym-True-cc-True: ['sbl-energy-landscape-comparison-euclid__f_himmelbleau_tg__f_himmelbleau_transition_graph_noisy__with_cc_constraints__algo_reg__edge_selection_min-cost__recursion_refined__symmetric_mode__emd_engine_symmetric.xml\n', 'sbl-energy-landscape-comparison-euclid__f_himmelbleau_tg__f_himmelbleau_transition_graph_noisy__with_cc_constraints__algo_reg__edge_selection_min-cost__recursion_refined__symmetric_mode__emd_engine.xml\n', 'sbl-energy-landscape-comparison-euclid__f_himmelbleau_tg__f_himmelbleau_transition_graph_noisy__with_cc_constraints__algo_reg__edge_selection_min-cost__recursion_refined__symmetric_mode__log.txt\n', 'sbl-energy-landscape-comparison-euclid__f_himmelbleau_tg__f_himmelbleau_transition_graph_noisy__with_cc_constraints__algo_reg__edge_selection_min-cost__recursion_refined__symmetric_mode__transportation_plan.dot\n', 'sbl-energy-landscape-comparison-euclid__f_himmelbleau_tg__f_himmelbleau_transition_graph_noisy__with_cc_constraints__algo_reg__edge_selection_min-cost__recursion_refined__symmetric_mode__transportation_plan_symmetric.dot\n']
XML: 1 / 1 files were loaded
EMD_comparators: summary
Total cost: 0.06592567189222925
Total flow: 0.9898
Normalized cost: 0.06660504333423849
XML: 1 / 1 files were loaded
EMD_comparators: summary
Total cost: 0.06592567189222925
Total flow: 0.9897999999999999
Normalized cost: 0.06660504333423849