Structural Bioinformatics Library
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 $ \sblELCL $ and $ \sblELCE $.

Using sbl-energy-landscape-comparison-euclid.exe

Pre-requisites

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 $ B $ of a local minimum. Abusing notations, the footprint of the basin in the conformational space $ \calC $ is also denoted $ B $. Denoting $ k_B $ is the Boltzmann constant and $ Z $ the partition function.

  • For small systems, the mass of the system may be obtained by integrating the Boltzmann factor over the basin region, namely:

    $ w(B) = \frac{1}{Z} \int_B \exp{ \frac{-V(c)}{k_B T}} dc. $

  • If minima are found by optimization (e.g., basin hopping), the $ w(B) $ can be estimated using the eigenvalues of the Hessian (curvature) matrix evaluated at those points [113] . 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 $ n_B $ falling in the basin region:

    $ w(B) = \frac{n_B}{n}. $

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 $ \PELsource $ and the demand landscape $ \PELdemand $. We also denote $ n_s $ and $ n_d $, respectively, their numbers of basins.

The local minima and the associated basins of the source landscape are denoted $ s_i $ and $ \basinof{s_i} $, respectively. For the demand landscape, local minima and basins are denoted $ d_j $ and $ \basinof{d_j} $.

We also assume that the weights of the basins have been computed. Denoting these $ \sbasinw $ and $ \dbasinw $ for a source and demand basin respectively, we define the sum of weights $ W_s = \sum_i \sbasinw $ and $ W_d = \sum_j \dbasinw $. Finally, so that the distance between the two aforementioned local minima is $ \dCalC(\PELsourcemin{i}, \PELdemandmin{j}) $.

To compare the two landscapes, we use the earth mover distance [101] , which is a particular case of mass transportation [110] . Intuitively, the technique fills basins in the target (aka demand) landscape using mass from basins of the source landscape. A basin from $ \PELsource $ can be split into several parts, and equivalently, a basin from $ \PELdemand $ can be filled from several basins from $ \PELsource $ (Fig. fig-comparison-of-basins-energy-landscapes). Denote $ \flowij $ the mass from $ \basinof{s_i}\in \PELsource $ moved into $ \basinof{d_j}\in \PELdemand $. The cost of moving $ \flowij $ units of mass depends linearly on the distances between the minima $ \PELsourcemin{i} $ associated with $ \basinof{s_i} $ and $ \PELdemandmin{j} $ associated with $ \basinof{d_j} $. A transport plan from $ \PELsource $ to $ \PELdemand $ is defined by triples $ (\PELsourcemin{i}, \PELdemandmin{j}, \flowij) $, with $ \flowij>0 $. Note that there are at most $ n_s\times n_d $ such triples.

Finding the optimal i.e. least cost transport plan amounts to solving the following linear program (LP):

$ LP \begin{cases} \mbox{Min} \sum_{i=1,\dots,n_s, j=1,\dots,n_d} \flowij \times \dCalC(\PELsourcemin{i}, \PELdemandmin{j})\\ \sum_{i =1,\dots,n_s} \flowij = \dbasinw & \forall j \in 1,\dots,n_d, \\ \sum_{j = 1,\dots,n_d} \flowij \leq \sbasinw & \forall i \in 1,\dots,n_s, \\ \flowij \geq 0 & \forall i \in 1,\dots,n_s, \forall j \in 1,\dots,n_d. \end{cases} $

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 $ \PELdemand $ need to be filled, while the third one indicates that a basin from $ \PELsource $ cannot provide more than it contains. (To simplify matters, we have assumed that $ W_s \geq W_d $. Handling the case $ W_s < W_d $ poses no difficulty, and the reader is referred to [101] .)

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 [101] :

$ \begin{cases} \numedgesPemd = \sum_{i,j \mid \flowij > 0} 1,\\ \flowPemd = \sum_{i,j} \flowij,\\ \costPemd = \sum_{i,j} \flowij \dCalC(\PELsourcemin{i}, \PELdemandmin{j}),\\ \distPemd = \costPemd / \flowPemd. \end{cases} $

Comparing two energy landscapes $ \PELsource $ and $ \PELdemand $
The landscape $ \PELsource $ is partitioned into the basins associated to its local minima (three of them on this example), and likewise for $ \PELdemand $ (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 $ S $ from $ \PELsource $ – that is $ S $ connects selected local minima in $ \PELsource $. Let $ D $ be the set of vertices from $ \PELdemand $ such that for each vertex $ d_j $ in $ D $, there is at least one edge emanating from a vertex $ s_i $ of the subgraph $ S $ with $ \flowij>0 $. The transport plan is called valid iff the subgraph $ D $ 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 $ \PELsource $ induces, through edges carrying strictly positive flow, a connected set of basins from $ \PELdemand $.


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 $ \text{\algoemdcccg} $: computes a transport plan given a upper bound $ C $ on the transport cost.

  • algorithm $ \text{\algoemdcccgi} $: 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 $ [0,C_{max}] $ of possible costs is given. The two modes are:

Recursion mode: refined. Given two costs $ C_{inf} $ and $ C_{sup} $, $ \text{\algoemdcccg} $ is called three times for the different costs $ C_{inf} $, $ C=(C_{inf}+C_{sup})/2 $, and $ C_{sup} $. Then:

  • a) If the two total volumes of flow with cost $ C_{inf} $ and $ C $ are different, then $ \text{\algoemdcccgi} $ is called with the interval $ [C_{inf},C] $.

  • b) If the two total volumes of flow with cost $ C $ and $ C_{sup} $ are different, then $ \text{\algoemdcccgi} $ is called with the interval $ [C,C_{sup}] $.

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 $ {s_1,s_2} $ of the source graph export towards the vertices $ d_1 $ and $ d_3 $ of the demand graph. The subgraph of the demand graph induced by these vertices is not connected – there is no edge linking $ d_1 $ to $ d_3 $.

Input: Specifications and File Types

The input are transition graphs of energy landscapes of the same molecule. They are produced by programs of the application Transition_graph_of_energy_landscape_builders, and stored in XML archives. A calculation is launched as follows:

> sbl-energy-landscape-comparison-euclid.exe --transition-graph data/himmelbleau_transition_graph_noisy.xml --transition-graph data/himmelbleau_transition_graph_simplified.xml --with-connectivity-constraints --symmetric-mode --directory results --verbose --output-prefix --log 
The main options of the program sbl-energy-landscape-comparison-euclid.exe are:
–transition-graph string: transition graph XML archives (used twice for target and source transition graphs)
–with-connectivity-constraints string: run Earth Mover Distance with connectivity constraints
–symmetric-mode string(= when using connectivity constraints): moving also from target to source


File Name

Description

XML transition graph Transition graph of noisy energy landscape of Himmelbleau function
XML transition graph

Transition graph of energy landscape of Himmelbleau function

Input files for the run described in section Input: Specifications and File Types .

Output: Specifications and File Types

PreviewFile Name

Description

General: log file

Log file

Log file containing high level information on the run of sbl-energy-landscape-comparison-euclid.exe

Module EMD comparisons: Earth Mover distance with and without connectivity constraints

EMD XML archive

Low level information on the run of EMD from source to target

EMD XML archive

Low level information on the run of EMD from target to source

Click it EMD graph graphviz file

Graph representation of EMD from source to target (see Graphviz)

Click it EMD graph graphviz file

Graph representation of EMD from target to source (see Graphviz)

Output files for the runs described in section Input: Specifications and File Types, classified by modules – see Fig. fig-elc-workflow .

Using sbl-energy-landscape-comparison-lrmsd.exe

Input: Specifications and File Types

The input are transition graphs of energy landscapes of the same BLN 69 molecule. They are produced by programs of the application Transition_graph_of_energy_landscape_builders, and stored in XML archives. A calculation is launched as follows:

> sbl-energy-landscape-comparison-lrmsd.exe --transition-graph data/bln69_transition_graph.xml --transition-graph data/bln69_transition_graph.xml --with-connectivity-constraints --directory results --verbose --output-prefix --log 
The main options of the program sbl-energy-landscape-comparison-lrmsd.exe are:
–transition-graph string: transition graph XML archives (used twice for target and source transition graphs)
–with-connectivity-constraints string: run Earth Mover Distance with connectivity constraints


File Name

Description

XML transition graph Transition graph of BLN69 energy landscape
Input files for the run described in section Input: Specifications and File Types .

Output: Specifications and File Types

PreviewFile Name

Description

General: log file

Log file

Log file containing high level information on the run of sbl-energy-landscape-comparison-lrmsd.exe

Module EMD comparisons: Earth Mover distance with and without connectivity constraints

EMD XML archive

Low level information on the run of EMD from source to target

Click it EMD graph graphviz file

Graph representation of EMD from source to target (see Graphviz)

Output files for the runs described in section Input: Specifications and File Types, classified by modules – see Fig. fig-elc-workflow .

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 [69] . Practically, various solvers can be used, e.g. the one from the Computational Geometry Algorithms Library [32] , lp_solve, the CPLEX solver from IBM, etc. In the following, the algorithm solving the linear program of Eq. (eq-emd-lp) is called $ \text{\algoEMDLP} $ .

Earth mover distance with connectivity constraints

Finding transport plans respecting connectivity constraints turns out to be a hard combinatorial problem [26] . The problem is not in APX, which means that if $ \Pcc \neq \NPcc $ holds, then, no polynomial algorithm with constant approximation factor exist.

Following definition def-tp-cc, we provide two algorithms:

  • $ \text{\algoemdcccg} $: computes a transport plan respecting connectivity constraints, given a bound C.

  • $ \text{\algoemdcccgi} $: computes different transport plans for different bounds not given by the user. More precisely, since the maximum cost is upper-bounded by $ C_{max} $, $ \text{\algoemdcccgi} $ computes different flow solutions by calling $ \text{\algoemdcccg} $ 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 $ \text{\algoemdcccg} $ for earth mover distance with cost and connectivity constraints}, has been reported in [26] .


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:

$ \begin{cases} \numedgesAemdccc = \sum_{i,j \mid \flowij > 0} 1,\\ \flowAemdccc = \sum_{i,j} \flowij,\\ \costAemdccc = \sum_{i,j} \flowij \dCalC(\PELsourcemin{i}, \PELdemandmin{j}),\\ \distAemdccc = \costAemdccc / \flowAemdccc. \end{cases} $
It can be shown that the Earth Mover Distance as defined by Eq. (eq-emd-lp}) yields a metric, provided that $ \dCalC $ is itself a metric, and that the sum of weights for the source and the demand graphs are equal [101] . On the other hand, the EMD with connectivity constraints fails to satisfy the triangle inequality [26] . 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:
$ \begin{cases} \ratiosymFlow = \frac{ \min(\flowAemdccc(A, B), \flowAemdccc(B, A)) } { \max(\flowAemdccc(A, B), \flowAemdccc(B, A))} ,\\ \\ \ratiosymCost = \frac{ \min(\costAemdccc(A, B), \costAemdccc(B, A)) } {\max(\costAemdccc(A, B), \costAemdccc(B, A))}. \end{cases} $


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.

The Traits Class

T_Energy_landscape_comparison_traits:

This class defines the main types used in the modules of the workflow. It is templated by the classes of the concepts required by these modules. This design makes it possible to use the same workflow within different(biophysical) contexts to make new programs. To use the workflow T_Energy_landscape_comparison_workflow , one needs to define:

  • what is the representation of a conformation (e.g number type of coordinates).
Template Parameters
GeometricKernelTraits class defining various geometric objects, in particular the number type used for representing the coordinates of a conformation, and the base representation of a point in dimension D – see class CGAL::Cartesian_d from the CGAL library.
ELSampleRepresentation of a conformation enriched by a height that is typically its energy.
DistanceFunctor returning the distance between two conformations. It has to define the types Point representing the conformation, and the type FT representing the returned value type.

The Workflow Class

T_Energy_landscape_comparison_interface_workflow: