The Eigenproblem Translated for Alignment of Molecules

Abstract

Molecular conformation as a subproblem of the geometrical shaping of the molecules is essential for the expression of biological activity. It is well known that from the series of all possible sugars, those that are most naturally occurring and usable by living organisms as a source of energy—because they can be phosphorylated by hexokinase, the first enzyme in the glycolysis pathway—are D-sugars (from the Latin dextro). Furthermore, the most naturally occurring amino acids in living cells are L-sugars (from the Latin laevo). However, a problem arises in dealing with the comparison of their conformers. One alternative way to compare sugars is via their molecular alignment. Here, a solution to the eigenproblem of molecular alignment is communicated. The Cartesian system is rotated, and eventually translated and reflected until the molecule arrives in a position characterized by the highest absolute values of the eigenvalues observed on the Cartesian coordinates. The rotation alone can provide eight alternate positions relative to the reflexes of each coordinate.

1. Introduction

The topological description of a molecule requires knowledge of the adjacencies (the bonds) between the atoms as well as their identities (the atoms). If this problem is simplified to the extreme, by disregarding the bond types and atom identities, then the adjacencies are simply expressed as 0 or 1 in the vertex adjacency matrix ([Ad]) and the identities are expressed as 0 or 1 in the identity matrix ([Id]). The characteristic polynomial (ChP) is the natural construction of a polynomial in which the eigenvalues of the [Ad] are the roots of the ChP, as follows:

λ is an eigenvalue of [Ad] ↔ it follows that [v] ≠ 0 eigenvector such that λ·[v] = [Ad]·[v]

(λ·[Id] − [Ad])·[v] = 0); since v ≠ 0 → [λ·Id − Ad] is singular → det([λ·Id − Ad]) = 0.

Therefore, the characteristic polynomial is defined by:

$$ChP\stackrel{\mathrm{def}}{=}\left|\lambda ·\left[Id\right]-\left[Ad\right]\right|.$$

The characteristic polynomial is a polynomial in λ of the degree of the number of atoms. The eigenproblem (the determination of eigenvalues and eigenvectors) is applicable to any Hessian [1] matrix [A] ([Ad] → [A]). The mixed derivatives of a scalar-valued function f are the entries off the main diagonal in the Hessian. Assuming that the derivatives are continuous, the order of differentiation does not matter (a result known as Schwarz’s, Clairaut’s, or Young’s theorem), and then the Hessian of f is a symmetric matrix.

Indeed, this is the case (a symmetric matrix) for the (vertex) adjacency matrix, and for the distance matrix—both topological (by bonds) and geometrical (by the atom coordinates).

Related to this problem is the issue of determining the best rotation to relate two sets of vectors. To this issue, a solution was proposed by calculating a symmetric matrix of Lagrange multipliers which is used to minimize the residuals of the linear association between the vectors [2]. Later, different approaches were proposed, such as geometric hashing [3], clique detection [4], the embedding problem [5], Gaussian molecular representation, Gaussian overlap optimization [6], and others covered in [7]. Some of the proposed solutions go a different way, involving physical means forcing the alignment [8,9], while the formulation of similarity metrics was one of the most recently proposed computational alternatives [10]. The alignment serves as a tool for other studies, including similarity analysis [11], docking [12], and structure–activity relationships [13].

The eigenproblem in relation to geometrical alignment was stated before in the context of surface analysis and control [14], and also can go another direction into the context of the molecule. In this context, the molecule is seen as more than a simple unweighted undirected molecular graph with undistinguishable atoms [15].

The eigenproblem of molecular alignment is analyzed in this paper.

2. Materials and Methods

The alignment of molecules can be stated in many ways, as listed in the introduction. For instance, one approach is to search for topological alignment, and another is to search for geometrical alignment. To anticipate the type of molecular alignment, it is necessary to employ the latter method—to search for geometrical alignment.

A molecule is taken here as an example from PubChem CID 444173 ((2R,3S,4R,5R)-oxane-2,3,4,5-tetrol), as shown in Figure 1.

For convenience, hydrogen atoms are excluded from the data and the analysis. The next table (

Table 1. 3D structural data for CID 444173 (heavy atoms, geometric coordinates, and atom symbols).

x (Å)	y (Å)	z (Å)	Atom and Label
0.7428	−1.4498	−0.0709	O	1
−1.1425	1.1688	1.3882	O	2
1.1461	1.0581	−1.4377	O	3
−2.754	−0.3648	−0.3408	O	4
2.7344	−0.2934	0.3835	O	5
−0.7774	1.0064	0.0187	C	6
0.7504	0.9905	−0.0675	C	7
−1.3475	−0.306	−0.532	C	8
1.3187	−0.2968	0.5474	C	9
−0.671	−1.513	0.1111	C	10

) contains the relevant information for the heavy atoms in the reference molecule.

The general way of constructing a characteristic polynomial is to provide an identity matrix [Id] and a Hessian matrix (herein labeled as [A]). If considering the topology of the molecule, then it is necessary to have the information regarding the connections between the atoms (e.g., bonds). Since all bonds are single bonds for the selected molecule, listing the atoms pairs of the bonds is enough (

Table 2. Topology data for CID 444173 (list of bonds between heavy atoms).

(1, 9)

(1, 10)

(2, 6)

(3, 7)

(4, 8)

(5, 9)

(6, 7)

(6, 8)

(7, 9)

(8, 10)

).

A deeper look into the eigenproblem (|λ∙I − A| = 0) is performed in the next section, with a specific focus on changing of the mathematical properties of the eigenproblem when the adjacencies in [A] change from symmetric to anti-symmetric.

3. Results and Discussion

From Table 2, the adjacency matrix [Ad] is immediate—zeros represent the entries without a bond between the labeled atoms, while ones appear otherwise. The adjacency matrix is Hessian. For convenience, its characteristic polynomial is:

ChP(λ; CID444173, “[Ad]”) = 1∙λ¹⁰ − 10∙λ⁸ + 31∙λ⁶ − 35∙λ⁴ + 11∙λ² − 1∙λ⁰.

As can be seen, the degree of the polynomial is 10, which is equal to the number of the (connected) atoms in the molecule. The general rule is that a characteristic polynomial is always of a degree equal to the size of the square matrices [I] and [A] (see the before given equation), from which it was derived.

The same strategy can be applied if the adjacency matrix [Ad] is replaced by the distance matrix [Di]. For convenience,

Table 3. Adjacency and distance matrices for CID 444173 (heavy atoms).

Ad	1	2	3	4	5	6	7	8	9	10	Di	1	2	3	4	5	6	7	8	9	10
1	0	0	0	0	0	0	0	0	1	1	1	0	4	3	3	2	3	2	2	1	1
2	0	0	0	0	0	1	0	0	0	0	2	4	0	3	3	4	1	2	2	3	3
3	0	0	0	0	0	0	1	0	0	0	3	3	3	0	4	3	2	1	3	2	4
4	0	0	0	0	0	0	0	1	0	0	4	3	3	4	0	5	2	3	1	4	2
5	0	0	0	0	0	0	0	0	1	0	5	2	4	3	5	0	3	2	4	1	3
6	0	1	0	0	0	0	1	1	0	0	6	3	1	2	2	3	0	1	1	2	2
7	0	0	1	0	0	1	0	0	1	0	7	2	2	1	3	2	1	0	2	1	3
8	0	0	0	1	0	1	0	0	0	1	8	2	2	3	1	4	1	2	0	3	1
9	1	0	0	0	1	0	1	0	0	0	9	1	3	2	4	1	2	1	3	0	2
10	1	0	0	0	0	0	0	1	0	0	10	1	3	4	2	3	2	3	1	2	0

lists these two matrices.

It can be checked (but is also true for the general case) that the distance matrix is Hessian, and therefore a characteristic polynomial can be computed for it as well. The next equation lists the ChP computed for the distance matrix [Di]:

ChP(λ; CID444173, ”[Di]”) = 1∙λ¹⁰ − 313∙λ⁸ + 3488∙λ⁷ − 15456∙λ⁶ − 34720∙λ⁵ − 40832∙λ⁴ − 23808∙λ³ − 5376∙λ²

The natural extension of this matrix is to employ 3D distances instead of topological distances. Of course, one consequence of this is that the characteristic polynomial would no longer have integer coefficients.

The next table lists the 3D distance matrix (distances were cut to four significant digits) and the roots of the associated characteristic polynomial.

One interesting remark to the data listed in

Table 4. 3D distance matrix and its eigenvalues for CID 444173 (heavy atoms).

3D	3D Distances										Eigenvalues
	1	2	3	4	5	6	7	8	9	10
1	0	3.541	2.885	3.671	2.347	2.890	2.440	2.427	1.429	1.427	−8.429
2	3.541	0	3.638	2.817	4.264	1.427	2.395	2.430	2.985	3.008	−6.218
3	2.885	3.638	0	4.294	2.769	2.413	1.428	2.983	2.410	3.509	−2.922
4	3.671	2.817	4.294	0	5.536	2.432	3.767	1.421	4.169	2.421	−1.893
5	2.347	4.264	2.769	5.536	0	3.762	2.406	4.183	1.425	3.627	−1.275
6	2.890	1.427	2.413	2.432	3.762	0	1.530	1.533	2.524	2.523	−1
7	2.440	2.395	1.428	3.767	2.406	1.530	0	2.510	1.536	2.884	−0.65
8	2.427	2.430	2.983	1.421	4.183	1.533	2.510	0	2.876	1.526	0
9	1.429	2.985	2.410	4.169	1.425	2.524	1.536	2.876	0	2.372	3.60 × 10−15
10	1.427	3.008	3.509	2.421	3.627	2.523	2.884	1.526	2.372	0	22.386

is that all roots are real (this is the general behavior for the roots of a characteristic polynomial).

Up until this point, the ideas presented in this paper have been reported before. Herein follows the extension to the extant knowledge. What if the same formula is applied to define the ChP for Cartesian coordinate distance matrices instead of for the Euclidian distance matrix?

Next three tables (

Table 5. First Cartesian coordinate (“x”) distances matrix for CID 444173 (heavy atoms).

Dx	1	2	3	4	5	6	7	8	9	10
1	0	1.8853	−0.4033	3.4968	−1.9916	1.5202	−0.0076	2.0903	−0.5759	1.4138
2	−1.8853	0	−2.2886	1.6115	−3.8769	−0.3651	−1.8929	0.2050	−2.4612	−0.4715
3	0.4033	2.2886	0	3.9001	−1.5883	1.9235	0.3957	2.4936	−0.1726	1.8171
4	−3.4968	−1.6115	−3.9001	0	−5.4884	−1.9766	−3.5044	−1.4065	−4.0727	−2.0830
5	1.9916	3.8769	1.5883	5.4884	0	3.5118	1.9840	4.0819	1.4157	3.4054
6	−1.5202	0.3651	−1.9235	1.9766	−3.5118	0	−1.5278	0.5701	−2.0961	−0.1064
7	0.0076	1.8929	−0.3957	3.5044	−1.9840	1.5278	0	2.0979	−0.5683	1.4214
8	−2.0903	−0.2050	−2.4936	1.4065	−4.0819	−0.5701	−2.0979	0	−2.6662	−0.6765
9	0.5759	2.4612	0.1726	4.0727	−1.4157	2.0961	0.5683	2.6662	0	1.9897
10	−1.4138	0.4715	−1.8171	2.0830	−3.4054	0.1064	−1.4214	0.6765	−1.9897	0

,

Table 6. Second Cartesian coordinate (“y”) distances matrix for CID 444173 (heavy atoms).

Dy	1	2	3	4	5	6	7	8	9	10
1	0	−2.6186	−2.5079	−1.0850	−1.1564	−2.4562	−2.4403	−1.1438	−1.1530	0.0632
2	2.6186	0	0.1107	1.5336	1.4622	0.1624	0.1783	1.4748	1.4656	2.6818
3	2.5079	−0.1107	0	1.4229	1.3515	0.0517	0.0676	1.3641	1.3549	2.5711
4	1.0850	−1.5336	−1.4229	0	−0.0714	−1.3712	−1.3553	−0.0588	−0.0680	1.1482
5	1.1564	−1.4622	−1.3515	0.0714	0	−1.2998	−1.2839	0.0126	0.0034	1.2196
6	2.4562	−0.1624	−0.0517	1.3712	1.2998	0	0.0159	1.3124	1.3032	2.5194
7	2.4403	−0.1783	−0.0676	1.3553	1.2839	−0.0159	0	1.2965	1.2873	2.5035
8	1.1438	−1.4748	−1.3641	0.0588	−0.0126	−1.3124	−1.2965	0	−0.0092	1.2070
9	1.1530	−1.4656	−1.3549	0.0680	−0.0034	−1.3032	−1.2873	0.0092	0	1.2162
10	−0.0632	−2.6818	−2.5711	−1.1482	−1.2196	−2.5194	−2.5035	−1.2070	−1.2162	0

and

Table 7. Third Cartesian coordinate (“z”) distances matrix for CID 444173 (heavy atoms).

Dz	1	2	3	4	5	6	7	8	9	10
1	0	−1.4591	1.3668	0.2699	−0.4544	−0.0896	−0.0034	0.4611	−0.6183	−0.1820
2	1.4591	0	2.8259	1.7290	1.0047	1.3695	1.4557	1.9202	0.8408	1.2771
3	−1.3668	−2.8259	0	−1.0969	−1.8212	−1.4564	−1.3702	−0.9057	−1.9851	−1.5488
4	−0.2699	−1.7290	1.0969	0	−0.7243	−0.3595	−0.2733	0.1912	−0.8882	−0.4519
5	0.4544	−1.0047	1.8212	0.7243	0	0.3648	0.4510	0.9155	−0.1639	0.2724
6	0.0896	−1.3695	1.4564	0.3595	−0.3648	0	0.0862	0.5507	−0.5287	−0.0924
7	0.0034	−1.4557	1.3702	0.2733	−0.4510	−0.0862	0	0.4645	−0.6149	−0.1786
8	−0.4611	−1.9202	0.9057	−0.1912	−0.9155	−0.5507	−0.4645	0	−1.0794	−0.6431
9	0.6183	−0.8408	1.9851	0.8882	0.1639	0.5287	0.6149	1.0794	0	0.4363
10	0.1820	−1.2771	1.5488	0.4519	−0.2724	0.0924	0.1786	0.6431	−0.4363	0

) list those results (the number of digits is displayed according to the input data—see Table 1).

It can be observed that the Cartesian coordinates distance matrices are no longer symmetric matrices, but are in fact anti-symmetric, meaning that M_i,j = −M_j,i.

The beauty of the result shown by taking a look at the eigenvalues. The next table (

Table 8. Eigenvalues for CID 444173 (heavy atoms).

	x1	x2	x3	x4	x5	x6	x7	x8	x9	x10
[Ad]	−2.318	−1.556	−1.334	−0.506	−0.41	0.41	0.506	1.334	1.556	2.318
[Di]	−8.429	−6.218	−2.922	−1.893	−1.275	−1	−0.65	0	3.60 × 10−15	22.386
[3D]	−8.722	−5.093	−3.145	−2.521	−1.474	−1.257	−1.05	−0.911	−0.784	25.007
[Dx]	15.299∙i	−15.299∙i	0	0	0	0	0	0	0	0
[Dy]	9.629∙i	−9.629∙i	0	0	0	0	0	0	0	0
[Dz]	6.973∙i	−6.973∙i	0	0	0	0	0	0	0	0

) lists the eigenvalues for all matrices.

It should be noted that the values listed in Table 7 reveal some computational errors. It is obvious (and it is so) that 3.6 × 10⁻¹⁵ is actually a “0” and it is necessary to be aware of this type of error coming from “machine epsilon” [16] which is about 10⁻⁷ for “single” precision, 10⁻¹⁶ for “double” precision, and about 10⁻¹⁹ for “extended” precision. Most floating-point implementations use “double” precision and thus the listed value (3 × 10⁻¹⁵) “fits in range”.

More important, as can be observed (see Table 7), the eigenvalues of [Dx], [Dy], [Dz] are all 0 excepting (always) two—which are (always) paired and (always) imaginary (i = √−1, see Table 8). This is the opposite of the traditional case of symmetric matrices, when the values are (always) real. This is the beauty of the result.

Moreover, it should be noted that the polynomial can be expressed with real-value coefficients as a product of a polynomial of degree 2 and a monomial of degree (n − 2), as listed in

Table 9. The polynomials of [Dx], [Dy], and [Dz] for CID 444173 (heavy atoms).

Matrix (A)	|λ·I − A| Polynomial
[Dx]	λ8∙(λ2 + 234.0448052)
[Dy]	λ8∙(λ2 + 92.7157814)
[Dz]	λ8∙(λ2 + 48.6224414)

.

A consequence is hidden behind this result—to obtain those two coefficients (which are actually the first and third coefficients, independent of how many atoms are in the molecule) it is necessary to obtain their roots. Therefore, it is not necessary to run an “eigenvalues” routine to obtain them; it is enough to run only two steps of a coefficient determination program (such as that described in [17]), which will produce a result much more quickly.

So, what if we conduct a rotation of the molecule?

For example, by rotating the molecule by 15° (15/180 radians; coordinates are given in Table 1), the values for the polynomials are changed—see

Table 10. The polynomials of [Dx], [Dy], and [Dz] rotated (15°, 15°, 15°) for CID 444173 (heavy atoms).

Matrix (A)	|λ·I – A| Polynomial
[Dx]	λ8∙(λ2 + 162.836846)
[Dy]	λ8∙(λ2 + 90.150945)
[Dz]	λ8∙(λ2 + 47.28921)

. First it should be pointed out that the polynomial is no longer invariant due to the choice of the system of coordinates. If invariants are sought, this is not a good situation—but for the purpose of addressing the alignment problem, this setup is very useable.

In the general case, with a₁, a₂, and a₃ as rotation angles defining the rotation matrices (given below), it is necessary to maximize the variance along the axes of coordinates.

$$\left[\begin{array}{lll}1\hfill & 0\hfill & 0\hfill \\ 0\hfill & \cos (a_2)\hfill & \sin (a_2)\hfill \\ 0\hfill & -\sin (a_2)\hfill & \cos (a_2)\hfill \end{array}\right],\left[\begin{array}{lll}\cos (a_1)\hfill & 0\hfill & -\sin (a_1)\hfill \\ 0\hfill & 1\hfill & 0\hfill \\ \sin (a_1)\hfill & 0\hfill & \cos (a_1)\hfill \end{array}\right],\left[\begin{array}{lll}\cos (a_0)\hfill & \sin (a_0)\hfill & 0\hfill \\ -\sin (a_0)\hfill & \cos (a_0)\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

This results a two-step algorithm, described below:

• Since rotation by a₀ leaves untouched the “z” coordinate, the first problem is to find a value of a₀ such that the squared sum of the eigenvalue(s) for the [Dx] matrix is minimized (or its coefficient from Table 9, which is x₁·x₂ = x₁·${\overline{x}}_1$ = −x₁² = −x₂², is maximized);
• Next, we need to leave untouched the “x” coordinate—which was already fitted in the first step. For this, we may want to employ rotation by a₂, such that the squared sum of the eigenvalue(s) for the [Dy] matrix is minimized (or its coefficient from Table 9 is maximized);
• There is no third step involving the third rotation matrix, because by maximizing (or minimizing) the first two coordinates, we have already employed all coordinates (x and y in the first step; y and z in the second).

Therefore, at this point we have the alignment of the molecule.

The problem of molecular 3D alignment involving the modified characteristic polynomial (eigenproblem) becomes a combinatorial problem since, after eigenvector minimization by each (two out of three) Cartesian coordinate, we obtain the molecules in their proper alignment or in the mirror of the proper alignment, when “x_i ← −x_i” and/or “y_i ← −y_i” and/or “z_i ← −z_i” transformation will align it.

Of course, a question may arise: what is the meaning of such alignment? This research is ongoing, but so far it has been found that this alignment corresponds to the minimization of the rotation inertia of the coordinates. In other words, the thinnest part of the molecule aligns with one coordinate, and then the thinnest part of what remains (so the molecule can be rotated around that axis) aligns with the second coordinate.

Revising the results communicated here, it should be noted that the classical eigenproblem is addressed to symmetric matrices—such as are the topological adjacency and topological distance matrices (shown in Table 3) and the geometrical distance matrix (Table 4). The peculiarity of the Cartesian distance matrices (shown in Table 5, Table 6 and Table 7) is the fact that they are anti-symmetric, sometimes called skew-symmetric matrices. This is, in mathematical terms, a strong property—as strong as the property of symmetry (please note that here the symmetry describes the matrices—namely, matrix A is symmetric if A = A^T and it is anti-symmetric if A = −A^T). On the other hand, the elements of the Cartesian coordinate matrices are mirrored relative to the main diagonal—this property is called reflection symmetry, line symmetry, or mirror symmetry—which makes these matrices very suitable for the same set of operations that are typically employed for symmetric matrices. Further, among the known properties of skew-symmetric matrices is the fact illustrated in Table 8—If A is a real skew-symmetric matrix and λ is a real eigenvalue, then λ = 0, i.e., the nonzero eigenvalues of a skew-symmetric matrix are purely imaginary”. Since a skew-symmetric matrix is similar to its own transposition, they must have the same eigenvalues. It follows that the eigenvalues (λ) of a skew-symmetric matrix always come in pairs (±λ), a property which is also illustrated in Table 8.

It should be noted that the generation of Cartesian coordinates from the diagonalization of adjacency or distance-related matrices is quite standard in mathematical chemistry. For instance, the methods to generate fullerene cages from Schlegel diagrams are normally embedded in fullerene sw packages (see for example [18]). Thus, the results communicated here may have useful applications in this regard.

4. Conclusions

The change from symmetry to anti-symmetry in the adjacency matrix of the eigenproblem moves the eigenvalues from real space into imaginary space. When the eigenequation is applied to the Cartesian space of the molecule instead of the topological or Euclidean spaces, the resultant roots (corresponding to the eigenvalues) are all 0 (multiple roots) excepting two, which are always imaginary (and complementary). The rotation of a molecule induces into the Cartesian space a way of aligning the molecule by maximizing the magnitude of the roots in a preselected order of the Cartesian axes. This property can be further exploited for the alignment of multiple molecules, when for highly symmetric molecules the alignment problem is turned into the (S₂)³ conformational problem.

Though the programs provided in the Supplementary Materials can be used to align any molecule, they are not communicated as a novel tool. Aligning a molecule by its Cartesian coordinates via the simultaneous alignment of many molecules—such as for molecular docking purposes—will require further study.