The Symmetric Difference Distance: A New Way to Evaluate the Evolution of Interfaces along Molecular Dynamics Trajectories; Application to Influenza Hemagglutinin

Abstract

We propose a new and easy approach to evaluate structural dissimilarities between frames issued from molecular dynamics, and we test this methodology on human hemagglutinin. This protein is responsible for the entry of the influenza virus into the host cell by endocytosis, and this virus causes seasonal epidemics of infectious disease, which can be estimated to result in hundreds of thousands of deaths each year around the world. We computed the three interfaces between the three protomers of the hemagglutinin H1 homotrimer (PDB code: 1RU7) for each of its conformations generated from molecular dynamics simulation. For each conformation, we considered the set of residues involved in the union of these three interfaces. The dissimilarity between each pair of conformations was measured with our new methodology, the symmetric difference distance between the associated set of residues. The main advantages of the full procedure are: (i) it is parameter free; (ii) no spatial alignment is needed and (iii) it is simple enough so that it can be implemented by a beginner in programming. It is shown to be a relevant tool to follow the evolution of the conformation along the molecular dynamics trajectories.

1. Introduction

Influenza virus causes seasonal epidemics of infectious disease and can even trigger pandemics. Inferring from U.S. data [1], these epidemics may cause hundreds of thousands of deaths each year around the world, mainly concerning the elderly, children, and immunosuppressed people, resulting in a real public health issue. Moreover, there were several severe influenza pandemics, such as the 1968 Hong Kong one [2] and the 1918 one, which resulted in about 20 millions deaths [3]. More recently, modeling studies attributed about 200,000 deaths to the 2009 influenza pandemic [4]. Initiation of virus infection involves multiple influenza Hemagglutinin (HA) binding [5] to sialic acids typically by α2,3 or even α2,6 linkages at the cell surface. The virus is then endocytosed through clathrin-dependent pathways, and the fusion of the viral particle with the cell’s membrane happens due to the acidification (about pH = 5) in the late endosomes. Indeed, this pH variation induces a folding change of HA homotrimers, resulting in the three fusion peptides’ exposure to the cell’s membrane [6]. HA is thus one of the two neutralizing antibody targets with neuraminidase, which is the second surface glycoprotein, allowing the virus to spread out from the cells [7].

In a recent modeling study to design new inhibitors against HA, two Molecular Dynamics (MD) simulations were performed, one at pH = 7 and one at pH = 5, this acidic pH being the one inducing a conformational change [8]. A motivation of this previous study was to compare the simulations at both pH values. Each protomer has two domains: HA1, responsible for the binding, and HA2 for the fusion, containing respectively 327 and 160 amino acid residues. The tridimensional structural data of HA were found in the Protein Data Bank (PDB), Code 1RU7, taken from [3] (see Figure 1). This structure is the one of A/Puerto Rico/8/1934 H1N1.

Starting from an initial conformation, 40 conformations were extracted every 0.5 ns for two performed simulations (see [8] for technical details about the protocol). The problem considered here is to evaluate, for each set of 41 conformations, the dissimilarities between pairs of conformations. For this purpose, we needed a simple procedure, easy to implement, and without spatial alignment because these latter generally involve unwanted contributions of meaningless parts of macromolecules, these parts being difficult to identify.

Usually, the dissimilarity between two conformers of the same macromolecule is measured through the computation of the RMSD (Root Mean Squared Deviation). For a macromolecule containing n heavy atoms, RMSD is the quadratic mean of the lengths of the n pairs of atoms (n atoms for the first conformer, n atoms for the second conformer), minimized for all rotations and translations of either of these two conformers. An analytical solution for this optimization problem is known: the optimal translation is such that the barycenter of the two conformers should coincide, and the optimal rotation is expressed with quaternions (see Appendix in [10] and Appendix A.5 in [11]). While the calculus of an RMSD is easy to compute, this quantity offers two drawbacks. The first one is the large numerical impact on the RMSD value from non-relevant parts of the macromolecule, such as the flexible ends of the backbone, or such as the irrelevant parts of the macromolecule far from the domain of interest for the experimentalist. For instance, in a protein–protein complex, the most relevant part is the protein–protein interface, because it is the area that can stand the weak bonds responsible for the stability of the complex. This is also the case for a protein–ligand complex. In the case of HA, the relevant part of the trimer is the three interfaces between the pairs of protomers (each of the three pairs of protomers defines an interface), because they are the location of interprotomers’ binding interactions. The second drawback of the RMSD is that, when it is computed over all heavy atoms, the side chains have an important numerical contribution to the RMSD value, although the location of the side chains is often considered to be meaningless due to their flexibility. For the reasons mentioned above, we restricted the evaluation of the dissimilarities between HA conformers generated by MD simulations to the C_α atoms of the three polypeptidic chains’ interfaces.

Thus, for each conformer, we have a set of C_α atoms. The simplest way to evaluate the dissimilarity between two sets is to compute their Symmetric Difference Distance (SDD). The SDD between two sets having respectively $n_1$ and $n_2$ elements is $n_1+n_2-2n_{12}$, $n_{12}$ being the number of elements of the intersection of the two sets [12]. The properties of a distance are recalled in Appendix A, including why they are useful. The computation of the SDD between the two sets of C_α atoms associated with two conformations of our HA trimers (i.e., there is one set of C_α atoms per trimer) is explained in Section 3. Here, the unit of this distance is the number of C_α atoms. It is recalled that, opposite what is encountered during RMSD calculations, the sets of C_α atoms defined by the interfaces are never pairwise associated, and in general, they have different cardinalities.

2. Results and Discussion

Starting from an initial conformation at time $t_0$, the MD simulation at pH = 7 generated 40 conformations, at times $t_1$, $t_2$, …, $t_{40}$, and similarly, the MD simulation at pH = 5 generated 40 conformations, the time step between two conformations being 0.5 ns [8]. Thus, there was a total of 41 conformations for each of these two MD simulations. We computed for the two pH values the 40 SDDs between the 40 generated conformations and their respective initial conformations at $t_0$. These distances values are plotted in Figure 2. The range of the distance values was 39–72 at pH = 7 (the unit of distances is the number of C_α atoms). It was slightly smaller at pH = 5: 25–52. The maximal values were reached for the frames generated respectively at $t_{26}$ and $t_{24}$. The maximal ratios of the SDDs values $n_1+n_2-2n_{12}$ to the total number of interface residues $n_1+n_2$ were 16.3% at pH = 7 and 17.2% at pH = 5. These moderately low values indicate that the generated conformations did not much differ from their respective initial ones at $t_0$. The correlation coefficient between the two sets of distances was 0.507. There was no reason to expect a high correlation coefficient since the SDD between the two initial conformations was 262, which was much higher than the range values at pH = 7 and at pH = 5.

We also computed for the two pH values the 40 RMSD between these 40 conformations and their respective initial conformations at $t_0$: see Figure 3. The range of RMSD values was 1.45–3.09 Å at pH = 7. It was slightly smaller at pH = 5: 1.65–2.58 Å. The RMSD values were small relatively to the size of the HA: the radius of the smallest sphere enclosing the tridimensional structure of the 11,510 atoms of HA was 68.55 Å for 1RU7. These small RMSD values indicate that the generated conformations did not significantly differ from their respective initial ones at $t_0$. More importantly, the trend and the main peaks visible in Figure 3 are visible in Figure 2. This means that working on C_α atoms at the interfaces rather than computing RMSD on all heavy atoms did not seem to cause any major information loss. Thus, our approach is pertinent.

Then, we computed for the two pH values the SDDs between the 40 successive pairs of conformations. These distance values are plotted in Figure 4. The range of distance values was similar for both pH values: 26–62 at pH = 7 and 24–52 at pH = 5. The maximal values were reached between the frames generated at $t_{22}$ for both simulations. The maximal ratios of the SDDs values $n_1+n_2-2n_{12}$ to the total number of interface residues $n_1+n_2$ were 13.9% at pH = 7 and 14.4% at pH = 5. These values were slightly lower than the ones of Figure 2 for pH = 7, and they were of the same magnitude as the ones of Figure 2 for pH = 5. For both simulations, the variations of the SDDs along time did not induce a trend to deviate more and more from the initial conformation, as shown in Figure 2. Thus, we can assume that the conformation is stable for the simulations at both pH values, at least for the interval of time considered (20 ns). The correlation coefficient between the two sets of distances was 0.468. There was no particular reason to expect a high correlation coefficient.

We also computed for the two pH values the RMSD values between the 40 successive pairs of conformations: see Figure 5. The range of RMSD values were similar for both pH values: 1.40–2.60 Å at pH = 7 and 1.42–2.54 Å at pH = 5. The RMSD values were small: we recall that the radius of the smallest sphere enclosing the 11.510 atoms of HA was 120.5 Å for 1RU7. Such a similarity between the range values at both pH is observed in Figure 4. The jumps observed in Figure 5 are also visible in Figure 4. Due to the relatively small ranges of RMSD values, the jumps seen Figure 5 are difficult to interpret. However, the maximal ratios of the SDDs values to the total number of interface residues that we already mentioned (around 14%) indicate that the interfaces may contribute significantly to the conformational changes along time. Again, working on C_α atoms at the interfaces rather than computing RMSD values for all heavy atoms did not seem to cause any major information loss.

After having observed that the HA conformation remains stable over the 40-ns simulations, we needed to define a mean conformation for docking purposes. Indeed, considering such an average structure in the course of a docking simulation, rather than a structure determined from X-ray diffraction experiments, enables us to take into account the relaxation of the protein in the solvent, at T = 310 K [8]. We defined this mean conformation with the algorithm described in [14]: it is the one that minimizes the sum of the squared distances to all other conformers. This mean conformation is the one of the frame generated at $t_{30}$ for pH = 7, and it is the one of the frame generated at $t_{15}$ for pH = 5. The residues at the interfaces of the mean conformers are given in Appendix B. The SDD between these two mean conformations was equal to 70 C_α atoms. This was few compared to the 1448 residues of each mean conformer, but it was larger compared to the respective 224 and 214 residues of the two interfaces, which had 184 residues in common. The RMSD between the 1448 pairs of C_α atoms of these mean conformers was 3.54 Å, while the largest atom-pair length was 13.7 Å (obtained for the first serine of the HA2 domain). The optimal superposition of the two interfaces is shown in Figure 6. It was realized with the CSR freeware (http://petitjeanmichel.free.fr/itoweb.petitjean.html), which implements a non-parametric algorithm performing a spatial alignment without a cutoff distance and without any input pairwise correspondence [15].

We have also performed an optimal superposition of the mean frames at $t_{30}$ for pH = 7 and at $t_{15}$ for pH = 5, on the basis of all pairs of 11.510 heavy atoms. The display of this optimal superposition in Figure 7 has been restricted to the two sets of 1448 C_α atoms because displaying both sets of 11.510 heavy atoms would have been much more confusing. Even the display of the two sets of 1448 pairs of C_α atoms in Figure 7 remains confusing. This is a general problem: the display of optimally-superposed macromolecules produced overloaded images, and viewing the result remained confusing even with interactive tools on the screen of a workstation. Thus, comparing the contents of Figure 6 and Figure 7 shows another reason why our approach based on interfaces is useful: it produces lightened images of optimally-superposed structural data, which would have been cumbersome to generate with existing visualization tools.

3. Methods

3.1. Computation of Interfaces within Macromolecular Complexes

The full procedure to compute a SDD between two macromolecular complexes relies on the computation of a protein–protein interface within each complex. According to [16], there are three main families of methods to compute interfaces:

• The cutoff method.
• The loss of accessible surface area upon binding.
• The Voronoi tessellation method.

Machine learning methods [17,18] are not considered here because they are irrelevant in our context. The cutoff method requires an arbitrary cutoff distance value, which may have a strong impact on the resulting computed interface [19]. Computing accessible surface areas requires fixing the values of atomic radii. Significantly different sets of radii values are found in the literature, depending on how they are defined [20,21,22]. It was shown that these radii have a considerable impact on surface areas [23,24]. The existence of parameters external to the input data and having a strong numerical impact on the results are drawbacks of these two families of methods. Thus, a non-parametric method is desirable. The third method, based on the Voronoi tessellation, in its original variant was parameter-free [25,26]. It was implemented in PROVAT software [27]. The full mathematical description of the Voronoi tessellation can be found in [28]. It is out of the scope of this paper. To summarize, each atom lies inside a convex polyhedral cell having its polygonal faces located at mid-distance from its neighboring atoms. Thus, two atoms are neighbors if their Voronoi cells share a common face. However, computing interfaces with the Voronoi method generates large cells, which induce the existence of meaningless long distances between neighboring atoms.

This is why we needed to compute interfaces with another free parameter method: we retained the PPIC software. PPIC is publicly available with its documentation on a repository located at http://petitjeanmichel.free.fr/itoweb.petitjean.html. The method implemented in PPIC was introduced only in a very recent preprint [29]. It extends the approach of [30] used to compute the interfaces in protein–ligand complexes. Its input is the tridimensional structure of a protein–protein complex or of a protein–ligand complex. This complex has two parts (molecule or macromolecule), named A and B. The algorithm has two steps:

• Generate the first part of the interface, constituted by the non-redundant set of all nearest neighbors of the atoms of A among the atoms of B.
• Generate the second part of the interface, constituted by the non-redundant set of all nearest neighbors of the atoms of B among the atoms of A.

Thus, the interface has two parts, i.e., two half interfaces, one in A and one in B. The roles of A and B are symmetric in the algorithm. From [31], it is known that each of these two parts is a subset of the half interface that would be computed by the Voronoi tessellation method. A nice consequence is that we do not observe anymore meaningless long distance pairs of neighboring atoms (i.e., one atom in A and one atom in B).

As a by-product, PPIC outputs the list of the RNNs (Reciprocal Nearest Neighbors). It is recalled that, in the Euclidean space, the nearest neighbor y of some point x is not always such that the nearest neighbor of y is x, e.g., consider three aligned points with abscissas $x=0$, $y=2$, and $z=3$, for which the nearest neighbor of x is y, while the nearest neighbor of y is z, not x.

When all atoms (or all heavy atoms) of the complex are used as inputs in the calculation of the interface, the pairs of atoms defined by the list of the RNN are a rough estimate of the location of potential interacting atom pairs. Moreover, to estimate the location of PPIs (Protein–Protein Interactions), looking at interacting atom pairs among the nearest neighbors makes more sense than looking at interacting atom pairs at farther distances than the nearest neighbors.

PPIC was shown to be effective to compute interfaces on a database of 1050 protein homo- and hetero-dimers ([29]; data from [32]). It was used to compute HIV2 protease–ligand interfaces ([29]; data from [33]). For both datasets, it was observed that the best agreement between the interface calculations produced on heavy atoms by PPIC and those produced by the cutoff method occurred for a cutoff near 3.6–3.7 Å (see Figure 1 in [29]). This is in agreement with the maximal donor–acceptor distance established in [34] at 3.9 Å. It is significantly smaller than the 4.5 Å cutoff distance between non-H atoms used by other authors [33,35,36].

3.2. Computation of Interfaces in Macromolecular Polymers

For a macromolecule containing two partners, A and B, the interface is constituted by two sets of atoms, one in A and one in B. In the case of a complex containing k polypeptidic chains, there are $\left(k\right(k-1\left)\right)/2$ interfaces to compute. This is the case of HA, for which $k=3$: there are three subsets, A, B, and C. Since the atoms of all chains are in the same macromolecular unit, we defined the global interface as the intersection of these $\left(k\right(k-1\left)\right)/2$ interfaces as the non-redundant list of atoms was extracted. Therefore, in the case of a macromolecular polymer, the global interface is constituted by only one list of atoms, not two. This is a difference from the case of macromolecular complexes.

3.3. Evaluation of the Dissimilarity between Two Interfaces

An interface is constituted by one or two sets of atoms. In the case of protein–ligand complexes containing several ligands, the interface may contain more than two sets of atoms. Considering two interfaces containing one set of atoms each, their dissimilarity can be evaluated by their SDD [12].

When each of the two interface contains K sets, $K\ge 1$, and when these two K-tuples of sets are pairwise associated, there are K SDDs to compute. A distance between these two K-tuples can be defined to be the sum of these K SDDs. It can be checked from the definition in Appendix A that it defines indeed a distance on the set of these K-tuples.

When the two K-tuples of sets are not pairwise associated, which in fact means that we consider unordered K-tuples, a distance between these two K-tuples can be defined to be the minimal value of the sum of the K SDDs, this minimum being taken over all $K!$ possible pairwise correspondences between the two K-tuples. It can be checked from the definition in Appendix A that this defines indeed a distance on the set of such K-tuples.

3.4. Comparison with Other Dissimilarity Measures

Several other structural dissimilarity criteria are encountered in the literature, which are qualified as metrics [37] (see also [38,39,40]), such as hydrogen bonds, distance from surface, the number of residues, or the number of heavy or polar atoms, or the number of waters in the vicinity of a specific region, RMSF (Root Mean Square Fluctuation), SASA (Solvent Accessible Surface Area), and gyration radius.

We outline here an ambiguity about the word “metric”: in many papers, it is used to name a structural dissimilarity criterion, which does not have the mathematical meaning mentioned in Appendix A. Even the dissimilarity criteria based on surfaces (e.g., SASA) cannot lead to defining a metric, except in the case of convex sets [41]. As mentioned at the end of Section 1, the SDD can be qualified as a metric [12]. The RMSD can be qualified as a metric when it is viewed as a distance induced by the Frobenius matrix norm, and this is why we compared the SDD with the RMSD. This happens in structural biology when we consider two matrices $A_1$ and $A_2$ of n lines and three columns, each matrix containing the spatial coordinates of a set of n points (the two sets of n points are thus pairwise associated). Setting $A=A_1-A_2$, $\sqrt{n}\mathrm{RMSD}=\sqrt{trace\left(A^{t}A\right)}$, where $A^t$ is the transpose of A.

3.5. Steps of the Methodology

The steps of the methodology to evaluate the evolution of interfaces along MD trajectories are summarized below. We assumed that the macromolecule of interest contained at least two partners, such as a protein and a ligand, two proteins or two protomers, etc.

• Generate the frames of the MD simulation.
• For each frame, generate the global interface with the procedure described in Section 3.2.
• For each couple of successive frames, evaluate the dissimilarity between the interfaces with the SDD, as described in Section 3.3.
• Follow the evolution of the interface along the trajectory using the SDD as a coordinate varying as a function of the time.

4. Conclusions

Our approach to evaluate structural dissimilarities between frames issued from MD calculations has several advantages over the traditional ones involving either RMSD calculations or interface calculations:

• It is parameter free.
• No spatial alignment is needed, thus no non-trivial numerical solver is needed.
• The problem of molecular graph symmetries occurring in some contexts for residues Val, Leu, Arg, Phe, Tyr, Glu, and Asp, which is almost always neglected when computing RMSD values, does not exist in our approach.
• All the steps of our algorithm can be coded by a beginner in programming.
• The dissimilarity between interfaces is measured with a distance (see Appendix A).
• Unwanted contributions of meaningless parts of macromolecules can be discarded (e.g., disordered parts in macromolecules, etc.).
• Images of optimal superpositions of full macromolecules are too overloaded compared to those of optimal superpositions of interfaces.

The software PPIC implementing our non-parametric algorithm of interfaces calculations is publicly available for free at http://petitjeanmichel.free.fr/itoweb.petitjean.html. The user can optionally run the cutoff method or the Voronoi tessellation method.

As an example, we presented results about two MD simulations of influenza HA (PDB Code 1RU7, 1448 residues). Knowing from MD simulations at pH = 7 and pH = 5 that HA is stable, the magnitude of numerical values that we observed can serve as a first reference basis to discuss the results of MD simulations of other macromolecules or complexes. Our approach is neither claimed to overcome the previous ones, nor is it devoted to replacing them: it is just an additional tool devoted to helping structural biologists, which is simple to use and which can be easily reprogrammed by beginners.