Factorization in molecular modeling and belief propagation algorithms

1.   Introduction

It is well acknowledged that high dimensional problems are difficult to analyze and predict. Multimedia data (e.g., pictures, videos) and molecular systems (e.g., biomolecules, materials, nano-systems) are typical examples. Since analytical solutions are either extremely challenging to achieve or non-existant for accurate and detailed description of probability distributions in high dimensional spaces, various numerical approximations are essential. The two most important and distinct paths are direct utility of dimensionality reduction algorithms ^{[1,2,3,4,5,6]} and factorizations ^[7,8,9], and both are extensive fields. The former path is to transform a $n$-dimensional probability distribution $P_n(x_1, x_2, \ldots, x_n)$, usually approximately, into a $m$-dimensional $P_m(x_1, x_2, \ldots, x_m)$ with $n \gg m$, while the latter path is to factorize $P_n(x_1, x_2, \ldots, x_n)$, either exactly or approximately, into a product of $K$ low dimensional terms $\prod _{i = 1}^K P_{m_i}(x_{i1}, x_{i2}, \ldots, x_{im_i})$, with the union of variables in all $K$ terms being $(x_1, x_2, \ldots, x_n)$. In this perspective, we focus on different factorizations utilized in molecular simulations and belief propagation algorithms, and hope to intrigue interests in development of more accurate and efficient approximations in relevant applications. In a probabilistic description of multivariate systems, there are three important types of problems. The first type is to find the most probable (set of) configuration(s) at a given resolution of representation, and more generally the global landscape of wells with significant probability weights, which is termed free energy landscapes in molecular systems ^[10,11,12] and multimodal learning in AI ^[13]. A highly related type of problems are to locate transition paths connecting various wells on the global landscape. Understanding such transition paths is of critical importance in accelerating sampling of other wells when starting from an arbitrary well of a landscape. Significant efforts have been invested in finding transition paths in molecular simulations ^[14,15,16]. The second type is to focus only on a subset of variables (or a single one in the extreme case), and the process of computing distributions of a subset of variables from the global joint distribution is termed marginalization. For example, marginalizing a two dimensional distribution $P_2(x, y)$ with both $x, y \in (0, 1)$ to the variable $x$ is accomplished by the following two summations: $P_1 (x = 0) = P_2 (x = 0, y = 0) + P_2 (x = 0, y = 1)$ and $P_1 (x = 1) = P_2 (x = 1, y = 0) + P_2 (x = 1, y = 1)$. There are $2^2$ terms to be calculated, and for a joint distribution with $n$ variables each has $|X|$ discrete values, $|X|^n$ terms need to be calculated, indicating the exponential computational complexity of marginalization. The third type is the parameterization of model parameters, including parameterizing force fields for molecular modeling and factor/node potentials for belief propagation. After introduction of typical factorization schemes in molecular modeling and belief propagation, we discuss strengths of different factorization in various type of problems. The organization of this perspective is summarized below. Following this introduction, Section 2 gives a brief introduction to exact belief propagation algorithms. Section 3 is on junction trees and approximate loopy belief propagations. Sections 4 and 5 describe factorization in regular molecular simulations and the local distribution theory. Section 6 provides discussions on the potential synergy and distinctions among these factorizations, and deep learning. Finally conclusions are presented.

2.   Exact belief propagation algorithms for trees

As mentioned in the introduction, marginalization of high dimensional joint distributions is difficult, and in many cases intractable. Fortunately, factorization significantly reduces this cost whenever it is possible to be performed exactly, or approximately with sufficient accuracy. A $n$-variable system $\boldsymbol{x} = (x_1, x_2, \ldots, x_n)$ forms an undirected graph with each variable being a node and an edge between all pairs of interacting variables. In graph theory ^[17], a clique is a set of nodes with an edge between each pair of nodes in this set, and a maximal clique is one that will be destroyed by an addition of a new node into this set. For a graph $\mathcal{G}$ with a set of maximal cliques $\mathcal{C}$ ^[17]:

with $f(x_c)$ being the potential over the variable set $c$ and $Z$ being the partition function $\sum_{\boldsymbol{x}}\prod_{c\in \mathcal{C}}f(x_c)$, which is summation of the product of all clique potentials over all possible configurations. In a tree, the maximal clique set is the edge set $\mathcal{E}$. For the specific example with 6 nodes shown in Figure 1, Eq (2.1) may be written as the following:

Figure 1.  A schematic illustration of the marginalization in Eqs (2.3), (2.5) and (2.8). Blue arrows indicate direction of message passing, and corresponding messages are shown beside each arrow. (a) Illustration of the process presented in Eq (2.3) for the marginal $p_1(x_1)$; (b) Illustration of the process presented in Eq (2.5) for the marginal $p_3(x_3)$ and (c) illustration of the belief propagation presented in Eq (2.8).

DownLoad: Full-Size Img PowerPoint

Given this factorization, the computational cost bound of marginalization goes from $|x|^6$ to $5\times |x|^2$. Specifically, the marginal probability $p_1(x_1)$ may be calculated in the marginalization process shown below (and is illustrated in Figure 1(a)) :

and $p_1(x_1)$ may be subsequently obtained with the normalization shown below:

Similarly, the marginal probability $p_3(x_3)$ may be obtained as shown below and is illustrated in Figure 1(b):

and $p_1(x_3)$ may be subsequently obtained with the normalization shown below:

In the calculation of both marginals $p_1(x_1)$ and $p_1(x_3)$, the following messages are computed:

Computing the marginal for each node in a tree takes $n-1$ messages, and consequently $n(n-1)$ messages are necessary for $n$ nodes when computed respectively (see Figure 1(a), (b)). However, instead of computing marginals for each variable separately, $2(n-1)$ messages over $n-1$ edges of an $n$-node tree maybe computed, such that all marginals maybe obtained without repetitive message calculation as illustrated in Figure 1(c) and demonstrated in the equation below.

This message passing is the sum-product, or the belief propagation algorithm for a general tree structure with pair potentials ^[7,8,9].

the message from $i$ to $j$ ($m_{i \rightarrow j}$) is sent once $i$ receives messages from all of its neighbors $N(i)$ except the node $j$. After two sweeps of message passing, the marginal for variable $i$ is obtained by:

When a unary potential is included for each node, we have:

One may easily compute the maximum probability configuration $\boldsymbol{x}^*$ by replacing the sum operations in the Eq (2.8) or (2.10) with max operations, and by bookkeeping the corresponding state of each variable in each maximum operation for backtracking. This results in the max-product algorithm for calculating the maximum probability configuration.

3.   Junction tree algorithm, loopy and gaussian belief propagation

In graphs where one or more loop(s) exist(s), the belief propagation algorithm ensures neither convergence nor correct calculation of marginals, and the same is true regarding the max product algorithm. Exact belief propagation maybe achieved for a general graph by constructing a chordal graph through triangulation ^[9], which involves adding edges to a graph such that no chordless loop of the length greater than 4 exists (see Figure 2 for a schematic illustration).

Figure 2.  A graph with a five-edge loop (left) is triangulated with addition of two edges (dashed blue lines) and become chordal (right).

DownLoad: Full-Size Img PowerPoint

Maximal cliques in a triangulated graph can be joined to form a junction tree. The most important property of such a tree is the junction tree property, which ensures that if a node $x_i$ is in two maximal cliques $c_j$ and $c_k$, then it must be a member of every maximal clique that is on the path from $c_j$ to $c_k$ ^[18,19]. The computational complexity of a junction tree is determined by the maximum tree width, which is the size of the largest maximal clique. Many different chordal graphs with different tree widths may be generated from a given graph with loops, and a number of algorithms are available to perform triangulation and construct chordal graphs ^{[20,21,22,23]}. However, no existing algorithm is available to ensure finding the best triangulation scheme, which is expected to generate one or more chordal graphs with the smallest possible tree width.

Even without assurance that exact marginals would be obtained by the belief propagation algorithm for graphs with loops, we may carry out such message passing computations anyway, and the algorithm is accordingly termed loopy belief propagation (LBP). Turbo code ^[24] is an outstanding example of utilizing LBP for telecommunication coding. From the parallel version ^[25] of the belief propagation algorithm (see Eq (3.1)) shown below:

Messages are usually initialized to be 1 and normalized at each time step. If a loopy belief propagation did converge, then Eq (3.1) suggests that the converged message would be a fixed point of the messaging updating function. For a general graph with pair potentials, the fixed point of belief propagation message updates are some local extrema of the Bethe variational problem ^[26]. The details of both the Bethe approximation and other treatments of loopy belief propagation are well beyond the scope this perspective, and interested readers are encouraged to find relevant information elsewhere ^[27,28].

Belief propagation was initially developed for discrete variables ^[7]. The summation in message computation becomes integration for continuous variables. Such integration is usually intractable for arbitrary distributions. Fortunately, gaussian distributions may be both analytically integrated and resulting in gaussians. This property may be utilized in belief propagation to treat continuous variables as gaussians. In cases where such direct approximation with gaussians is not sufficient, linearization with the first order Taylor expansion may be utilized to improve the accuracy ^[29].

4.   Factorization in conventional molecular modeling

Conventional molecular simulation is a mature methodology with a wide variety of tools. These methods are routinely utilized in many research fields, including but not limited to chemistry, biology and materials science with more than 30,000 relevant publications each year ^[30]. An essential component of the foundation for molecular simulation is a parameter set termed force fields, which represents interaction energies between/among various molecular degrees of freedom as a function of relevant coordinates. Conventional molecular simulation is deemed by many of its practitioners as an independent research field from, while being in fact a well defined form of, statistical machine learning. The force fields parameterization is the learning stage and the simulation and sampling constitute the inference process.

The total energy of a given molecular configuration $\boldsymbol{x} = (x_1, x_2, \ldots, x_n)$ for a simple Lennard-Jones (LJ) system is defined as ^[31]:

Such LJ terms are usually utilized to describe various Van de Waals interactions between neutral atoms/particles with respective parameters $\epsilon$ and $\sigma$. Apparently, this equation is essentially a factorization:

with $\boldsymbol{x}$ representing the set of coordinates for all comprising atoms/particles of a specific configuration. Let us denote each particle as a node, and the corresponding interaction between a pair of particles as an edge. Each factor represents an interacting pair of particles. This factorization has a very similar form as that in the belief propagation in a tree (Eq (2.1)). However, there are three major differences. Firstly, in Eq (2.1), the graph is static with fixed number of edges. However, for a graph formed by LJ particles, edges are dynamic in two aspects. One is that as the system propagates with time, edges break as some pairs of particles drift from within to beyond cutoff and new ones form as others drift from beyond to within. The other is that for pairs within cutoff, their interactions change as a function of distance. Secondly, nodes in Eq (2.1) have discrete states, while properties of LJ particles (e.g., position and momentum) are continuous. Thirdly, in Eq (2.1), each factor is over a maximal clique comprising a pair of connected nodes. For LJ particles, when the cutoff for direct interaction is taken as $3\sigma$, all particles within a sphere with the cutoff as its diameter are certainly within a single clique (since they all interacts) as shown in Figure 3(a). The number of particles ranges from a few to more than a dozen. LJ potential is factorized over pairs of particles instead of maximal cliques. Such factorization is apparently an approximation, and the resulting error is not easy to quantify for arbitrary LJ interactions. Apparently, the exact number of particles in each maximal clique can be different as local packing of LJ particles continuously changes over space and time. In a graph, two maximal cliques share one or more nodes are effectively connected. For densely packed molecular systems that may be represented as LJ particles, each maximal clique usually overlaps with multiple neighboring maximal cliques. Therefore, it is essentially a sure thing that most maximal cliques participate in loops. Constructing junction trees on such systems certainly results in huge tree width. Exact inference is consequently prohibitive, and the utility of LBP is necessary.

Figure 3.  (a) An example of 13 Lennard-Jones particles of diameter $\sigma$ enclosed within a sphere of diameter $3\sigma$ ($r_{cutoff}$). Three particles partially covered by other particles are shown with dashed circles. These 13 particles form a clique with 78 edges, all should be represented by a single factor/potential for exact inference while in conventional molecular simulation, 78 edges correspond to 78 factors. (b) A schematic illustration of LDT theory in which each spherical space of radius $r_{cutoff}$ is represented as a single factor, a central particle experiences mediated interactions from particles beyond $r_{cutoff}$ but within $2r_{cutoff}$, only interactions further than $2r_{cutoff}$ are treated as long-range interactions.

DownLoad: Full-Size Img PowerPoint

For more general molecular systems, energy contributions are more diverse ^[32,33,34] with addition of various bonded (including bonds, bends and dihedrals) (see Figure 4) and electrostatic interaction terms between ions or partially charged polar atoms, and are usually calculated as shown below:

Here $i, j, k, l$ and $esi$ are indices for bonds, bends, dihedrals, Van de Waals and electrostatic interactions respectively. Taking exponentiation of this equation results in:

Figure 4.  Schematic representation of bonded interactions. Left: A bond between the atom $1$ and the atom $2$, Center: A bend formed by the bond $1–2$ and the bond $2–3$, Right: A dihedral formed by the plane $1–2–3$ and the plane $2–3–4$.

DownLoad: Full-Size Img PowerPoint

Bonded terms are fundamentally local in space. Van der Waals interactions are usually modeled by LJ potentials. Similar to simple LJ systems, interactions in maximal cliques are approximately represented by assuming independent local interactions involving single (VdW and bonds), double (bends) or triple (dihedrals) edge(s). There is also challenge of addressing continuous variables with non-gaussian distributions and dynamic edges. More importantly, electrostatic interactions are long-range. Should such interactions to be treated exactly, then factorization is not possible and a molecular system needs to be treated as one single maximal clique. The exponential complexity with respect to the size of the maximal clique renders exact treatment intractable in many realistic molecular systems, so treating long range interactions as being independent from both local ones and each other is an essential approximation in computational molecular science. In Eq (4.5), all interaction terms are treated as being independent regardless of interaction range in space. Fortunately, such approximations have been widely tested in molecular simulations with great successes ^[35]. In reality, naive pairwise calculation of electrostatic interaction diverges and usually some form of Ewald summation is utilized to ensure convergence ^[36,37].

5.   Factorization in the local distribution theory

Just like the repetitive local computation in the calculation of marginals as demonstrated in Eqs (2.3), (2.5) and (2.7), ubiquitous repetitive local sampling (RLS) exists in molecular simulations ^[30,38,39]. From a graph perspective, such RLS is fundamentally associated with the intrinsic dynamic property of graphs describing molecular systems. The local distribution theory (LDT) was developed to address RLS and tremendously increases computational efficiency without reducing accuracy/resolution ^[30,39]. As shown in Eq (5.1), LDT is another approximate factorization scheme for facilitating statistical machine learning of dynamic graphs in high dimensional space.

Here, $\Phi = (\phi_1, \phi_2, \ldots, \phi_k)$ is relevant thermodynamic and environmental variables, $\boldsymbol{x} = (x_1, x_2, \ldots, x_n)$ are molecular coordinates. $R = (r_1, r_2, \ldots, r_m) \; (m \le n, \; m = n \text{ is preferred})$ are dynamic local regions, each of which represent a dynamic collection of coordinates ($r_i = (x_{i1}, x_{i2}, \ldots, x_{il})$) for molecular degrees of freedom. A translation, rotation and permutation invariant coordinate transformation is necessary to feed $r_i$ into neural networks or any other powerful nonlinear fitting machinery for local free energy computation. Effectively, taking a local region $r_i$ and calculating its local free energy (high dimensional potential of mean force) is similar to taking two particles in conventional molecular simulation and calculating their two-body potential of mean force. $\prod_{i = 1}^mq(\Phi, r_i)$ is the product of local distributions and the global correlation factor $\frac{Q(\Phi, R)}{\prod_{i = 1}^mq(\Phi, r_i)}$ is approximately factorized into mediated ($\sum F_{MED}(\Phi, R)$) and long range ($\sum F_{LR}(\Phi, R)$) contributions. The learning stage of LDT is to train local distributions, which cache statistical significance of various molecular configurations within local spaces defined by a preset cutoff from the origin of a local coordinates originating from a given molecular DOF. Such local distributions, defined by trainable neural network parameters with selected neural network architectures, store both energetic and entropic contributions from extensive local sampling data of either experimental or computational origin. By caching and repetitively utilizing such local distributions, LDT eliminates RLS to a great extent, and therefore saves computational resources in sampling dynamic graphs. In the inference stage, local distributions are assembled under sampling constraints due to mediated interactions. Direct long range interactions are treated as independent from local ones (Eq (5.1)), and selection of specific method does not impact local computation. Searching for the best methodology to match the LDT scheme is an important issue to be tackled in the future. Hierarchical definition of "local" with different spatial coverage might be an efficient path to be explored for large systems with hierarchical dynamics spanning multiple time scales.

In conventional molecular simulations, pairs of interactions (or edges in a graph representing a molecular system) within cutoff are included as independent factors (Figure 3(a)). As illustrated in Figure 3(b), both local interactions confined within cutoff spherical space and mediated interactions are jointly considered in LDT. Independence assumptions are therefore significantly weaker and accordingly, improvement of accuracy is expected with tremendously increased efficiency due to substitution of RLS by local distributions. Therefore, LDT has potential to address the tradeoff dilemma between accuracy and efficiency in molecular simulations. However, this does come with a significant one-time price to pay for accurate fitting of local distributions, which is fundamentally a density estimation task in high dimensional system, an open problem with many feasible approaches and respective limitation ^[40,41].

6.   Distinction and potential synergy among different factorization schemes

One fundamental difference between factorization in the belief propagation algorithm and that in the molecular modeling of physical systems is that the former deals with a static graph $(V, E)$ while the latter addresses a dynamic one $(V, E_{dynamic})$. The variation of edges in molecular systems is continuous, but may be simplified as binary in the case of a protein contact map ^[42]. Specifically, as the distance between two atoms/particles changes over time, so the corresponding interaction potential. In the belief propagation, an isolated maximal clique has a fixed factor potential, the computation is to figure out how interactions with neighboring cliques influence the distribution of internal states of nodes. However, we assume that forces exerted on each atom/particle does not change within a time step of molecular dynamics simulations, we may similarly deem a molecular system in a sufficiently short time period as a static graph to apply belief propagation algorithms. When the exact belief propagation formulation is possible, it provides two great advantages. One is the rapid determination of the most probable configuration, and the other is the explicit normalization of local probability distributions, which is a mechanism missing in the molecular modeling methodology tool box, and is potentially helpful in keeping errors under control if successfully applied in molecular systems. To realize such a goal, the first step is to define nodes with internal states. For example, in a full quantum mechanical treatment, one may define each atom as a node with rich internal electronic states. For proteins, residues may be defined as distinct nodes with internal rotameric states, which may be approximately represented either discretely or continuously. Apparently, definition of nodes is flexible for a given type of molecular systems and requires significant effort to explore. For a given molecular system with specific definition of nodes, we may perform loopy belief propagation to search for the most probable global configuration (or a number of configurations with the largest statistical weights) of internal node states, followed by sampling of different edge configuration with LDT or some other properly selected sampling strategy. Three important factors for the application of the belief propagation in molecular systems are definition of nodes, selection/development of approximate loopy belief propagation algorithms and selection of sampling strategy. It is likely that each may impact and influenced by the other two factors. Exploring different combinations of these algorithms are likely fertile grounds for developing next generation molecular modeling methodologies. An additional advantage of formulating a molecular system with the belief propagation is that marginalization, which is of critical importance in docking and design of key spots in biomolecular interactions, becomes straight forward.

While exact belief propagation need only one full sweep (to and from each node) of message passing, LBP requires iteration, and therefore the speed of convergence becomes a practical challenge. For the first type of problems stated in the introduction (to find maximum probability configuration in particular and to map out full landscape in general), the process is essentially some form of walking over the probability landscape and settling down at the most easy-to-reach local well from a start configuration. To map out the global landscape, one has to design ways to escape from any given local well and reach other potentially deeper or shallower wells. Unfortunately, no systematic theory is available for belief propagation algorithms to achieve this goal. The rich history of enhanced sampling of free energy landscapes ^{[11,43,44,45,46,47,48]}, coarse graining ^[49,50] (which is effectively smoothing of free energy landscape), and transition path sampling ^[14,15,16] (which focus on finding transition paths between/among different free energy wells) in molecular systems may provide hints to inventing new algorithms for the accelerated convergence of LBP.

Belief propagation algorithms reduces repetitive local computation in the redistribution of marginals for internal node states by surrounding connected nodes. LDT reduces the repetitive local sampling of dynamic edges by utilizing pre-trained neural networks which implicitly holds the probability density of a given configuration. Deep learning has found extensive applications in molecular simulations ^[51,52], especially in the construction of deep potentials (parameterization of force fields) ^{[53,54,55,56]}. It is likely that there are more potential ways of combining different factorization schemes in molecular simulations, belief propagation and deep learning to be invented. For example, when compared with other conventional dimensionality reduction algorithms, deep learning has the advantage of being much more flexible, and such dimensionality reduction might be combined with molecular simulation and belief propagation to construct novel approximate multi-scale solutions.

Matrix factorization is widely utilized in transforming input for deep learning ^[57,58,59]. It is important to note that such factorization is a global behavior in the probability distribution space. This is in strong contrast to decomposing a global distribution into local ones in the belief propagation and molecular simulations.

7.   Conclusions and prospects

Belief propagation is a family of algorithms in a large class of factorization methodology called probabilistic graphical models with focus on static graphs. These algorithms are mainly developed by informational scientists and are applied in coding, communication and machine learning. Variants of factorization in molecular simulations, including the local distribution theory, are fundamentally factorization focusing on dynamic graphs representing complex molecular systems, and are mainly developed by computational chemists and physicists. Both are extensive research fields with very rich contents and numerous important contributions from many scientists. The most outstanding difference is that factorization in belief propagation has an exact solutions for simple and junction trees but lacks methodology development in enhanced sampling, while factorization in molecular simulations has rich history in accelerating sampling of free energy landscapes but no exact solution. Comprehensive survey of these two fields is apparently well beyond the scope of present work. This perspective is intended to emphasize this discipline gap, and intrigue interest for developing next generation algorithms that assimilates the best of both fields.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

We appreciate the Interdisciplinary Integration and Innovation Project of JLU with grant number JLUXKJC2021ZZ05. We thank Professor Zhonghan Hu for comments on the local distribution theory when the manuscript was in preparation.

Conflict of interest

The authors declare there is no conflict of interest.