Population reversal driven by unrestrained interactions in molecular dynamics simulations: A dialanine model

The widely used “Ewald summation”,¹ a method first introduced to take into account electrostatic interactions in a crystal, and much later adapted in different versions to the computational study of molecules,² has proved to be very successful. In a few words, this method transforms a sum with very slow decay in distance into a sum with fast decay in the reciprocal space, allowing to take into account the electrostatic contributions of all charges, including those very far away.³ An underlying assumption of the method is that we are dealing with periodic boundary conditions (PBC). Specifically, the system we are dealing with is replicated indefinitely in each dimension, forming an infinite lattice (Figure 1). But this assumption has a downside when simulating systems that don’t form a lattice—like it is done in most simulations in computational biology: the distance among each element in a given lattice cell and its corresponding element in another cell is a constant. The same restriction applies to their relative orientations.

Consequently, studying a protein system, any protein-protein interactions between proteins in different cells are strongly restricted. If the system is made of highly charged proteins, or proteins with a high dipolar moment, there will be a strong electrostatic force between them, which could potentially lead to their aggregation, if protein-protein interactions are considered in full. This is not allowed when doing the usual one-protein simulation with Ewald, since a change in protein-protein distance is strictly forbidden. Ewald artifacts (inaccuracies created by the use of artificial periodicity) were studied before,^4–6 but it was Weber and McCammon⁷ who specifically studied the problem presented here: due to the computational power available at the time, it was difficult to clearly distinguish between the different minima in the dialanine free energy space, hence leading to the neglect of the effects of this particular issue.

The vast majority of biological activities, which are often concentration dependent, are affected by protein-protein interactions. In physiological functions like protein cooperativity, in physiological/pathological aggregation and in the more extreme case of protein crowding,^8–10 the presence of other proteins clearly affects protein’s average conformation and function. Induction of even minor conformational changes can cause key variations in protein activity and function. What is currently lacking is that in typical Molecular Dynamics simulations (MD) with PBC there is only one protein in the simulation box (Figure 1, left panel). This arrangement doesn’t allow, as we discuss above, full interactions among different proteins. All the system’s features driven by protein-protein interactions are missing with this simulation arrangement.

A step ahead is to test if having more proteins in the same simulation box (Figure 1, right panel) allows capturing those features. In this work we will show the effect of protein-protein interaction on the thermodynamics of the system, in MD simulations. As a model system we will use a peptide, N-acetyl-N’-methyl-L-analylamide, usually known as dialanine. In particular, we will show how the free energy landscape of an individual molecule is affected by the presence of other freely interacting molecules of its same type. In order to explore these effects, simulations were performed with many copies of the same molecule (dialanine) in the same box, so as to allow full peptide-peptide (p-p) interactions, and not be restricted by the fixed interactions typical of Ewald type PBC. In order to enhance p-p interactions and make its effects stronger, we also designed dialanine molecules with increased artificial dipole moments.

We selected the well-known dialanine^11–14 as our model system because it features a free energy diagram with well separated minima and a simple geometry with only two dihedral angles: ψ and φ (see Figure 2).

We prepared 6 different systems, differing in the number of dialanines in the central cell, and the artificial dipole moment of those peptides, as shown in Table I. We built dialanine with different dipole moments changing the charges of the two carbons at the ends of the molecule. We added ±0.3 e, ±0.5 e and ±1.0 e, where e is the charge of an electron. The notation we use for these systems is αA-βD, where α is the number of dialanines in the central cell, and β is the charge added to the system in order to increase its dipole moment.

The intrinsic dialanine dipole moment ranges from 0.3 to 8.5 D, depending on conformation. The dipole moment values of dialanine prepared with additional charges span the intervals (4,14), (9,21) and (19,38) D for the cases with ±0.3 e, ±0.5 e and ±1.0 e, respectively. Note that a dipole moment of 9 D (the average value for the case with ±0.3 e) or, equivalently, a moment of 0.9 D per heavy atom, is a value not far from the average for the structures in the PDB.¹⁵ Hence, the values here used, although artificial, are representative of real proteins.

All the systems were built using the LEAP program available in AMBERTOOLS 1.5¹⁶ and solvated by locating them at the center of a rectangular parallelepiped, which also defined the size of the central cell of the periodic box. For the 1-dialanine cases, the minimum initial distance between each atom of the peptide and the side of the box was 15 Å; for the 8-dialanine cases, that distance was increased to 18 Å.

For each of the systems, we performed a minimization in order to remove bad contacts, using 5,000 steps of steepest descent, followed by another 5,000 steps of conjugate gradient. Then, we performed a MD simulation where the temperature was raised from 0 to 298 K, during 1 ns, using a Langevin thermostat with a collision frequency of 2 ps⁻¹. This was followed by a 5 ns pressure equilibration, via a weak coupling isotropic barostat with a relaxation time of 2 ps, while keeping the temperature constant at 298 K, again with the Langevin thermostat and same parameters as before. Finally, a 200 ns constant number of atoms, volume and temperature (NVT) production run was performed. All this procedure was carried out using the AMBER12¹⁶ package, with the AMBER ff10 force field and the TIP3P¹⁷ water model. The use of the SHAKE¹⁸ algorithm for hydrogen-containing bonds allowed us to set the time step to 2 fs. PBC were applied using the particle mesh Ewald method.² Snapshots were saved every 10 ps for later analysis.

In order to calculate the precise free energy landscape of the ψ dihedral, we employed the umbrella sampling method. We used 36 sampling windows, one every 10 degrees, where the ψ angle was forced toward the center of the window with the help of a harmonic constraint with an elastic constant of 20 kcal/(mol Ų). In the cases with multiple dialanines, the sampling was taken only on one of them. A MD simulation was performed for each window, for 5 ns, using the same parameters that were used in the production run. The results were analyzed with the Grossfield’s lab implementation of the Weighted Histogram Analysis Method, WHAM.¹⁹ 

The value of the dihedral angles of dialanine, φ and ψ, were calculated for each frame in the different trajectories. Since φ is more than 95% of the time in the same region, between -150°and -50°, we focus our study in ψ, who shows a much stronger variation as it will be shown below. In terms of convergence of the MD, the number of transitions between the two minima of ψ, are in the range (900, 1200), for the different systems we employed. These large numbers, consistent with the simulation of a small molecule (2 residues), strongly suggest that the length of the simulation is enough to correctly cover the landscape.

Figure 3 shows the distribution of ψ, where we see that:

1. Without an enhanced dipole, the behavior of ψ is the same for the cases with single or multiple dialanines, i.e., the p-p interaction existing in the multiple case is not sufficient to change the behavior.

2. The addition of a dipole to the case with only one dialanine is still not enough to generate a significant change.

3. The addition of a dipole to the case with several dialanines drastically changes the behavior of the system; the higher p-p interaction generated by the enhanced dipoles, together with the possibility of close interactions since now there are several (8) freely moving dialanines in the same cell, creates forces strong enough to disrupt the equilibrium of the previous cases.

Figure 4 shows that this change in equilibrium is not caused by an increase in concentration, as one might initially think. The change, as mentioned above, is due to the strong p-p interactions, in addition to the availability of freely moving dialanines, which allow close inter-dialanine interactions.

It is interesting to note that the distribution of dialanine’s end-to-end distance, shown in Figure 5 and defined as the distance between the two outermost carbon atoms, closely follows the distribution of ψ. The global maximum changes its location from ∼7 Å to ∼5.5 Å when the p-p interaction is fully allowed (by way of introducing more dialanines into the simulation box).

Using the well-known relation between frequency and free energy, F_i = -k_BT ln(N_i/N_max), Figure 6 shows the free energy landscape of the ψ vs φ space for the 1A-0.0D and 8A-1.0D cases. There, we see the above mentioned drastic change in the distribution of ψ, a population reversal (global minimum switches from 155° to -10°), while there is no significant change in the distribution of φ, which has a constant global minimum around -90°.

In order to obtain a highly precise free-energy profile of the ψ coordinate, we applied the umbrella sampling method to this angle. It is important to note that this calculation does not depend on the thorough coverage of the entire configurational space available to the molecule, as it was required for the previous evaluation represented in Figure 6. The reason for such difference is that umbrella sampling imposes relatively strong restraining potentials on the coordinate studied, as explained later, which has the benefit of drastically reducing the amount of configurational space to be covered by the MD.

Using 36 windows to cover the 360° range of the dihedral, and checking that there is a good overlap of the windows, we employed WHAM to get the free energy profile showed in Figure 7 for two different cases. As expected from the frequency distribution results, the profiles show that there is a global minimum located at 155° for the 1A-1.0D case. On the other hand, when p-p interactions are fully allowed, the minimum moves to -8°. ΔF (ψ:150° → 0°) changes from +0.4 to -0.9 kcal/mol, which amounts to a total change of 1.3 kcal/mol, a very significant figure for a small peptide like dialanine. The profiles of 1A-0.0D, 8A-0.0D, not shown in the figure for clarity purposes, are comparable to the one of 1A-1.0D—red curve in figure. Likewise, the profiles for 8A-0.3D and 8A-0.5D are comparable to the one of 8A-1.0D—blue curve.

An error estimation for the free energy calculations was done, by checking for convergence of the results in terms of 1) the equilibration length of the windows (obtaining a free energy when 0, 0.5, 1, 1.5, 2 and 2.5 ns of the MD for each window is discarded), and 2) the production length of the window (obtaining a free energy when using a production MD length of 3, 3.5, 4, 4.5 and 5 ns). The worst case error thus obtained is around 0.1 kcal/mol, well below the 1.3 kcal/mol figure encountered for the difference between the cases in Figure 7.

This study shows that the cases with one dialanine, with or without dipole, and the cases with eight dialanines and no dipole (1A-0.0D, 1A-1.0D and 8A-0.0D; from now on, these cases will be called I) have the same behavior in terms of the equilibrium configurations of dihedral angles. In other words, the free energy landscape of dialanine systems with non—or weakly—interacting proteins, are indistinguishable. Conversely, in the case of 8 dialanines with an artificially enhanced dipole moment (8A-1.0D, from now on it will be called case II) there is a strong modification in the free energy landscape of a single dialanine, induced by inter-molecule interactions. Therefore, to affect the behavior of a single protein it is not sufficient to either introduce a stronger interaction, like a high artificial dipole moment, or to allow the dialanine to interact with other peptides in the same cell. We need to have both, other peptides and a stronger dipole-dipole interaction. And this can, and in the present study does, drastically perturb the equilibrium of the system, as shown in Figure 7.

A qualitative difference between the 2 mentioned cases, I and II, is that the conformation given by a dipole-dipole coupling (anti-parallel alignment of two dialanines) is frequent only in II. And the corresponding distribution of charges favors a contracted conformation. And this, in turn, is associated with structures corresponding to ψ values around 0° (see Figure 5). This intuitive observation is confirmed by the numbers in Table II, extracted from the simulations.

This is, of course, a qualitative explanation of the change in the equilibrium conformation. In fact, a proper discussion on equilibrium conformation changes should consider the free energy by adding an entropic component to the above numbers. But we expect minor variations by doing that, not relevant to this qualitative argument.

Thus, the ψ = 150° conformation is energetically favored when there is no dipole coupling, or it is very weak (cases I), whereas ψ = 0° , is energetically favored when there’s a strong dipolar coupling (case II). In this way, the change in the location of the global minimum of ψ is explained by the strong dipolar moments (that can generate strong p-p interactions) plus the presence of many freely moving dialanine molecules (permitting the close coupling of dialanines). As shown, the removal of any of those two elements results in the original location of the minimum, at ψ = 0° .

In summary, this work shows the importance of carefully evaluating the appropriateness of using Ewald type PBC when only one peptide is considered in the simulation cell, as is commonly done, since this could alter the equilibrium of the system to a strong degree. Such simplified approach will be appropriate only if the system does not produce strong p-p interactions. Otherwise, a simulation with several peptides in the periodic cell, allowing peptides to get closer to each other, will be necessary.

We thank Prof. Ivet Bahar for stimulating discussions and for her support. F. P. acknowledges the RiMED Foundation, Italy, for their financial support.