Fixman compensating potential for general branched molecules

Abstract

The technique of constraining high frequency modes of molecular motion is an effective way to increase simulation time scale and improve conformational sampling in molecular dynamics simulations. However, it has been shown that constraints on higher frequency modes such as bond lengths and bond angles stiffen the molecular model, thereby introducing systematic biases in the statistical behavior of the simulations. Fixman proposed a compensating potential to remove such biases in the thermodynamic and kinetic properties calculated from dynamics simulations. Previous implementations of the Fixman potential have been limited to only short serial chain systems. In this paper, we present a spatial operator algebra based algorithm to calculate the Fixman potential and its gradient within constrained dynamics simulations for branched topology molecules of any size. Our numerical studies on molecules of increasing complexity validate our algorithm by demonstrating recovery of the dihedral angle probability distribution function for systems that range in complexity from serial chains to protein molecules. We observe that the Fixman compensating potential recovers the free energy surface of a serial chain polymer, thus annulling the biases caused by constraining the bond lengths and bond angles. The inclusion of Fixman potential entails only a modest increase in the computational cost in these simulations. We believe that this work represents the first instance where the Fixman potential has been used for general branched systems, and establishes the viability for its use in constrained dynamics simulations of proteins and other macromolecules.

Figure 1

Fixman potential and torque plots for a three-bond (C4) chain with 90° bond angle. (a) Comparison of Fixman potential values from the analytical expression in Eq. 24 and the SOA based expression in Eq. 17. (b) Comparison of Fixman torque values from the SOA based expression in Eq. 23 and those from numerically differentiating the Fixman potential from Eq. 17.

Figure 2

Contour plot of the Fixman potential for the pair of torsion angles for a C5 chain with 109° bond angle.

Figure 3

Torsion angle distributions for the three-bond C4 chain, for (a) FLEXIBLE, (b) TORSIONAL, and (c) FIXMAN simulations. The dashed lines represent the predicted distributions, with the ones in (a) and (c) being uniform pdfs from Eq. 8, and the one in (b) from Eq. 25. The bin size for the histograms is dα = 7.2°.

Figure 4

Free energy distribution in a three-bond chain C4 system, plotted as a function of the torsion angle. The negative of the Fixman potential, $-{\mathcal{U}}_{\mathrm{f}}$, is very close to the free energy in the TORSIONAL simulation, demonstrating that the application of the Fixman potential corrects the bias in the free energy landscape of the TORSIONAL model. The free energies were calculated using bins of size 7.2° for the torsion angles.

Figure 5

Probability distribution functions for torsion angles in a fourteen-bond C15 chain. The application of the Fixman potential in the TORSIONAL model recovers the uniform pdf for all the torsion angles. The bin size for the histograms is dα = 7.2°.

Figure 6

The global minimum and maximum configurations of the Fixman potential for (a) the C5 and (b) the C11 molecule correspond to the self-intersecting planar and self-intersecting spatial conformations, respectively. In the self-intersecting planar conformation, each torsion angle is 0°. In the self-intersecting spatial conformation, each torsion angle is ≈±83°.

Figure 7

Structures of (a) alanine dipeptide and (b) valine dipeptide. The atoms highlighted with the green circles in (a) and the red circles in (b) constitute the ${\mathrm{C}}^1-{\mathrm{N}}^2-{\mathrm{C}}_{\alpha }^2-{\mathrm{C}}^2$ torsion in the respective systems. Probability distributions for these torsions for alanine dipeptide and valine dipeptide, respectively, in (c) and (d). The application of the Fixman potential in the TORSIONAL model recovers the uniform pdf observed for the FLEXIBLE model. The bin size for the histograms is dα = 7.2°.

Figure 8

Joint probability distribution function of the ${\mathrm{C}}_{\alpha }^1-{\mathrm{C}}^1-{\mathrm{N}}^2-{\mathrm{C}}_{\alpha }^2$ (α₁) and ${\mathrm{C}}^1-{\mathrm{N}}^2-{\mathrm{C}}_{\alpha }^2-{\mathrm{C}}^2$ (α₂) torsions for alanine dipeptide. (a) FLEXIBLE simulations, (b) TORSIONAL simulations, and (c) FIXMAN simulations. The bin size for each axis is dα = 18°.

Figure 9

The structure of the ten-residue β-hairpin protein, chignolin in (a). The green, red, and blue sphere atoms constitute the ${\mathrm{C}}^1-{\mathrm{N}}^2-{\mathrm{C}}_{\alpha }^2-{\mathrm{C}}^2$ torsion, the ${\mathrm{C}}^4-{\mathrm{N}}^5-{\mathrm{C}}_{\alpha }^5-{\mathrm{C}}^5$ torsion, and the ${\mathrm{C}}^7-{\mathrm{N}}^8-{\mathrm{C}}_{\alpha }^8-{\mathrm{C}}^8$ torsion, respectively. The application of the Fixman potential in the TORSIONAL model recovers the expected uniform pdf as see in (b) for the ${\mathrm{C}}^1-{\mathrm{N}}^2-{\mathrm{C}}_{\alpha }^2-{\mathrm{C}}^2$ torsion, (c) for the ${\mathrm{C}}^4-{\mathrm{N}}^5-{\mathrm{C}}_{\alpha }^5-{\mathrm{C}}^5$ torsion, and (d) for the ${\mathrm{C}}^7-{\mathrm{N}}^8-{\mathrm{C}}_{\alpha }^8-{\mathrm{C}}^8$ torsion. The bin size for the histograms is dα = 7.2°.

Figure 10

The computational time for the dynamics evaluation versus the number of clusters for the GNEIMO method with and without the GNEIMO-Fixman method. Including GNEIMO-Fixman increases GNEIMO costs by 24%.