Computing the Dissociation Constant from Molecular Dynamics Simulations with Corrections for the Large Pressure Fluctuations—Aquaglyceroporins Have High Affinity for Their Substrate Glycerol

TL;DR

This paper shows that correcting pressure fluctuations in simulations reveals glycerol has high affinity for aquaglyceroporin proteins.

Contribution

A novel correction factor for pressure fluctuations in MD simulations is derived and applied to compute accurate ligand–protein binding affinities.

Findings

Without correction, glycerol's affinity for aquaglyceroporins appears very low (~1/M).

With correction, glycerol's affinity is in the 1/mM to 1/µM range, indicating high binding affinity.

The correction factor depends on volume changes between bound and unbound states of the ligand.

Abstract

In this paper, we consider the inevitable large fluctuations of pressure in typical molecular dynamics (MD) simulations of ligand–protein binding problems. In simulations under the constant pressure of one bar, the pressure artifactually fluctuates over the range of ±100 bars or more. This artifact can cause gross inaccuracy in the apparent binding affinity computed as the ratio of the probability for the ligand to be bound inside the protein and the probability for the ligand to be outside the protein. Based on statistical thermodynamics, we derive a correction factor for the ligand–protein binding affinity to compensate for the artifactual pressure fluctuations. The correction factor depends on the change in the system volume between the bound and the unbound states of the ligand. We conducted four sets of MD simulations for glycerol affinities with four aquaglyceroporins: AQP10, AQP3, AQP7, and GlpF. Without the correction factor, the apparent affinity of glycerol with each of these four aquaglyceroporins is computed directly from the simulations to be very low (~1/M). With the correction factor applied, glycerol’s affinity is computed to be 1/mM to 1/µM. In conclusion, glycerol has high affinity for its native facilitator aquaglyceroporins, which is in contrast to the current literature not correcting the artifactual consequences of the large pressure fluctuations in typical in silico experiments.

Figures

Full text

1. Introduction

The determination of the binding affinity (or its inverse, the dissociation constant) of a protein for its substrate molecule is fundamental to the understanding of biophysical/biochemical processes. In principle, binding affinity can simply be computed from molecular dynamics (MD) simulations by counting the protein–substrate binding and unbinding events from the simulated dynamics. However, typical simulations of a system consisting of a few hundred thousand atoms are not suitable for such a simple task, for the following reason: A typical model system carries fluctuations in pressure around a few hundred bars for an NPT simulation (N indicates constant number of molecules in a model system, P means a constant pressure of one bar, and T represents a constant temperature). See Figure 1 for an example. Based on statistical thermodynamics [1], this undesirable artifact is inevitable due to the smallness of the system size (as explicitly pointed out in NAMD User’s Guide https://www.ks.uiuc.edu/Research/namd/2.14/ug/node39.html (accessed on 28 December 2025)), the pressure fluctuations are inversely proportional to the volume of a model system). This huge pressure fluctuation naturally favors unbinding protein from its substrate and thus leads to the computed affinity (apparent affinity) being significantly lower than the true affinity [2]. In other words, the computed apparent dissociation $[eqn]$ can be significantly larger than the true dissociation constant $[eqn]$ .

In this paper, we take a theoretical consideration of the artifactual pressure fluctuations and derive a correction factor for the computed apparent dissociation constant to estimate the true $[eqn]$ . We run four sets of MD simulations for the binding affinities of four aquaglyceroporins (AQGPs) (human AQP10, AQP3, AQP7, and E. coli GlpF) for their substrate glycerol. In each case, the apparent affinity computed from the simulation is very low ( $[eqn]$ ~100 mM) but the corrected estimation for the true $[eqn]$ is in the µM range. This confirms that aquaglyceroporins have reasonably high affinities for their substrate glycerol.

The significance of this study stems from the following fact: AQGPs, constituting a subfamily of aquaporin (AQP) proteins [3,4,5], play essential roles in human physiological processes. They facilitate the transport of glycerol and some other nutrient molecules across the cell membrane [6], which are fundamental to many physiological/pathophysiological processes such as fat metabolism and metabolic diseases [7]. This research produces insights into two controversies about whether or not an aquaglyceroporin has high affinity for its substrate glycerol and why typical simulations suggest very low aquaglyceroporin–glycerol affinity. Here is a brief summary of the current state: In 1994, in vitro experiments showed that E. coli aquaglyceroporin GlpF facilitates the unsaturable uptake of glycerol up to 200 mM into Xenopus oocytes [8]. These experiments suggested that GlpF has very low affinity for its substrate glycerol. From 2008 to 2014, however, multiple experiments [9,10,11] demonstrated that human aquaglyceroporins AQP7, AQP9, and AQP10 conduct the saturated transport of glycerol. In these experiments, the Michaelis constants were measured to be around 10 µM, indicating that human AQGPs have high affinities for glycerol. Examining the crystal structures of AQGPs (GlpF in 2000 [12], P. falciparum PfAQP in 2008 [13], AQP10 in 2018 [14], and AQP7 in 2020 [15,16,17]), they all have glycerol molecules inside the AQGP channel and near the channel openings on both the intracellular (IC) and the extracellular (EC) sides. These structures tell us that all four AQGPs have affinities for glycerol. If we insisted that unsaturated transport precludes high affinity, the previously stated data would suggest inconsistency among the experiments in the current literature. However, in an in silico in vitro study [18] of glycerol uptake into human erythrocytes through AQP3 [19], we showed that “an AQGP (having high affinity for its substrate glycerol) can conduct glycerol transport that is unsaturated up to 400 mM. The transport pathway for unsaturated transport through a high affinity facilitator protein was shown to involve two glycerol molecules next to each other, both bound inside an AQP3 channel (one at the high affinity site and one at a low affinity site) for the transport of one glycerol molecule across the cell membrane [18]. It is the substrate–substrate interactions (mostly repulsion due to steric exclusion) inside a single-file channel that make it easy for two glycerol molecules cooperatively to move one substrate molecule across the AQGP channel via the high affinity site.” Therefore, in vitro experiments exhibiting unsaturated glycerol transport do not preclude high affinity between an AQGP and its substrate glycerol. All in vitro experiments either allow or show the AQGP–glycerol affinity being high. This study demonstrates why the simulations of the current literature fail to predict such high affinities and how to correct for the artifactual large fluctuations in typical MD studies.

2. Methods

2.1. Model System Setup and Simulation Parameters

We followed the well-tested steps in the literature to set up our model systems and to run MD simulations. Therefore, this subsection is mostly direct quotes from the current literature with minor adaptations: “We employed CHARMM-GUI [20,21,22] to build an all-atom model of an AQGP tetramer embedded in a 120 Å × 120 Å patch of membrane (lipid bilayer consisting of 193 phosphatidylethanolamine/POPE, 119 phosphatidylcholine/POPC, and 80 cholesterol/CHL1 molecules). The coordinates of AQP10 (PDB: 6F7H), AQP3, AQP7 (PDB: 6QZI), and GlpF (PDB: 1FX8) were taken from Ref. [14], Ref. [23], Ref. [15] and Ref. [12], respectively. The positioning of the AQGP tetramer was determined by matching the hydrophobic side surface with the lipid tails and aligning the channel axes perpendicular to the membrane. The AQGP–membrane complex was sandwiched between two layers of TIP3P waters, each of which was approximately 30 Å thick. The system was then neutralized and salinated with Na^+^ and Cl^−^ ions to a salt concentration of 150 mM. Finally, glycerol was added to the system to a concentration of $[eqn]$ . In this way, four all-atom model systems were built for the four AQGPs: System AQP10, System AQP3, System AQP7, and System GlpF.

We employed NAMD 2.13 (for initial stage of a simulation when constraints were needed) and NAMD 3.0 (for all equilibrium runs when no constraints were used) [24,25] as the MD engines. We used CHARMM36 parameters [26,27,28] for inter- and intra-molecular interactions. We followed the literature’s standard steps to equilibrate the system [17,29,30,31]. Then, we ran unbiased MD for a few thousand ns with constant pressure at 1.0 bar (Nose–Hoover barostat) and constant temperature at 303.15 K (Langevin thermostat). The Langevin damping coefficient was chosen to be 1/ps. The Langevin piston period was 50 fs and the Langevin piston decay was 25 fs. The periodic boundary conditions were applied to all three dimensions. The particle mesh Ewald (PME) was used for the long-range electrostatic interactions (grid level: 128 × 128 × 128). The time step was 2.0 fs. The cut-off for long-range interactions was set to 10 Å with a switching distance of 9 Å.”

2.2. Direct Computation of Apparent Affinity from MD Simulations

We ran the MD simulation for each system for a sufficiently long time. After a system was fully equilibrated, we continued the MD run for hundreds of ns and used that part of the MD trajectory to compute the probability $[eqn]$ for an AQGP channel being occupied with a glycerol molecule. We considered a channel as being occupied when a glycerol molecule was found in the single-file region of the channel, 7.1 Å to the IC/EC side from the NAA/NPS motifs. We illustrate such a case of occupancy in the left panel of Figure 2. When the system is in its equilibrium state, the probability $[eqn]$ for a given glycerol concentration $[eqn]$ , from which we computed the dissociation constant from the binding probability: $[eqn]$ .

2.3. Theoretical Formulation of Correction for the Pressure Fluctuation

The dissociation constant of a substrate–protein complex can be computed from two partial partitions:

$$[eqn]$$

Here, N is the number of molecules in the system; P, the pressure; and T, the absolute temperature. The partial partition $[eqn]$ is the NPT ensemble partition of the system under the condition that the substrate molecule is away from the protein, freely sampling the volume of $[eqn]$ where $[eqn]$ is the standard concentration of 1 mol/L. The partial partition $[eqn]$ is the NPT ensemble partition of the system under the condition that the substrate and the protein are bound together as a complex. In this, the two partial partitions have identical dimensionalities and units. In an in silico implementation of the NPT ensemble, the preservation of N constant is trivial and the maintenance of T constant is well accomplished with fluctuations much smaller than the temperature T but the fluctuations in P are much greater than the constant pressure of $[eqn]$ bar in a typical MD study of the literature (See, e.g., Figure 1). The root mean square of the pressure fluctuations $[eqn]$ is much greater than the mean pressure $[eqn]$ . (The brackets represent statistical average over P.) Therefore, the computed value, the apparent dissociation constant,

$$[eqn]$$

which is not equal to the true dissociation constant

$$[eqn]$$

Consider approximating the Gibbs free energy $[eqn]$ as follows ( $[eqn]$ denoting the Boltzmann constant):

$$[eqn]$$

Here, $[eqn]$ is the mean volume of the system. Assuming the simplest form of the pressure fluctuations (with a Gaussian distribution), we have the following approximations for the partial partitions:

$$[eqn]$$ $$[eqn]$$

Therefore, we arrive at the following relationship:

$$[eqn]$$

Solving for $[eqn]$ and carrying out the Gaussian averages, we obtain the following:

$$[eqn]$$

Here, $[eqn]$ and $[eqn]$ are, respectively, the volume of system when the substrate and the protein are bound together and the volume of the system when the substrate is away from the protein. Further, we can rewrite the above equation as

$$[eqn]$$

It is a fundamental law of statistical thermodynamics [34] that $[eqn]$ where $[eqn]$ is the bulk modulus of the system. Then, we observe that the multiplier

$$[eqn]$$

can be easily evaluated from in silico simulations. This multiplier is the central result of this research, the correction factor for an accurate estimation of the true dissociation constant from simulations with large pressure fluctuations.

Considering $[eqn]$ , $[eqn]$ and the small volume difference, e.g., $[eqn]$ , leads to $[eqn]$ . Indeed, the apparent dissociation constant computed directly from MD simulations is generally very far from the true dissociation constant, depending on the volume change when a substrate molecule is dissociated from the protein.

3. Results

We ran 2000 ns of MD simulation of System AQP10 (illustrated in SI, Figure S1) after the initial runs for system preparation/setup. The root mean square deviation (RMSD) curves of AQP10 tetramers (shown in SI, Figure S2) show that System AQP10 reached equilibrium during last 1000 ns of the simulation. All atoms of a protomer were included in the RMSD calculations. Therefore, we conducted equilibrium statistical analyses of data from this part of the MD for the estimation of glycerol–AQP10 affinity. In Figure 3, we plotted the probability for an AQP10 channel to be occupied by one or more glycerol molecules. This probability gives the apparent dissociation constant $[eqn]$ , suggesting low affinity for glycerol–AQP10, in line with the current literature. Also, for the last 1000 ns of the simulation, we analyzed the correlation between the volume of System AQP10 and number of glycerol molecules located inside the four channels of AQP10, which is shown in Figure 4. This analysis is important because a glycerol molecule displaces different numbers of water molecules when it is inside an AQP10 channel and when it is outside the protein (in the bulk of the saline). Consequently, the volume of System AQP10 fluctuates around different mean values. Linear regression of the data gives rise to the volume difference $[eqn]$ (which is equal to the negative of the linear regression slope). Namely, the system volume averages higher by 5.98 mL/mol when a glycerol molecule is relocated from the bulk outside AQP10 to inside an AQP10 channel and, correspondingly, a number of water molecules move out of the channel to the bulk outside AQP10. The positive volume difference reflects on the fact that glycerol displaces fewer water molecules inside a channel than when it is outside the channel, which is illustrated in Figure 2. Using this volume difference in the correction factor, we obtain the estimated true dissociation constant $[eqn]$ , demonstrating that AQP10 has reasonably high affinity for its substrate glycerol.

In studies similar to System AQP10, we conducted a 2600 ns MD run for System AQP3, a 1500 ns MD run for System AQP7, and a 1000 ns run for System GlpF. The simulation data are shown in SI, Figures S3–S14. In each case, the MD was run for at least 500 ns after the system reached full equilibrium, as indicated by the RMSD curves for the four monomers of an AQGP tetramer (SI, Figure S4 for System AQP3, Figure S8 for System AQP7, and Figure S12 for System GlpF). The last equilibrium parts of the trajectories (the last 600 ns for System AQP3, the last 500 ns for System AQP7, and the last 500 ns for System GlpF) were used to extract the probability for a glycerol molecule to be inside an AQGP channel (shown in SI, Figure S5 for System AQP3, Figure S9 for System AQP7, and Figure S13 for System GlpF), from which the apparent dissociation constants were computed. The results are tabulated in Table 1. The values of the apparent dissociation constant $[eqn]$ all suggest very weak affinity between an AQGP and its substrate glycerol, which, again, is in line with the current literature. These equilibrium parts were further analyzed for the correlations between the system volume and the number of glycerol molecules inside an AQGP tetramer, which are shown in SI, Figure S6 for System AQP3, Figure S10 for System AQP7, and Figure S14 for System GlpF. Linear regression produces a negative slope, giving rise to a positive volume difference $[eqn]$ in each of the three cases (tabulated in Table 1) and producing an estimated true $[eqn]$ in micro molars.

4. Conclusions

Pressure fluctuations in a typical in silico study of the current literature can lead to predictions of very low affinities for protein–substrate binding problems. The apparent dissociation constant can be many times greater than the true dissociation constant if dissociation of the substrate from the protein corresponds to an increase in the system’s volume. The in silico studies of four aquaglyceroporins (human AQP10, AQP3, AQP7, and E. coli GlpF) all suggest low affinities without the correction factor to count for the large artifactual fluctuations in pressure. Once the correction factor is included, the estimated values of the true dissociation constant all lead to the conclusion that an aquaglyceroporin has reasonably high to very high affinity for its substrate glycerol.