Application of the Brown Dynamics Fluctuation-Dissipation Theorem to the Study of Plasmodium berghei Transporter Protein PbAQP

In this article, the Brownian dynamics fluctuation-dissipation theorem (BD-FDT) is applied to the study of transport of neutral solutes across the cellular membrane of Plasmodium berghei (Pb), a disease-causing parasite. Pb infects rodents and causes symptoms in laboratory mice that are comparable to human malaria caused by Plasmodium falciparum (Pf). Due to the relative ease of its genetic engineering, P. berghei has been exploited as a model organism for the study of human malaria. P. berghei expresses one type of aquaporin (AQP), PbAQP, and, in parallel, P. falciparum expresses PfAQP. Either PbAQP or PfAQP is a multifunctional channel protein in the plasma membrane of the rodent/human malarial parasite for homeostasis of water, uptake of glycerol, and excretion of some metabolic wastes across the cell membrane. This FDT-study of the channel protein PbAQP is to elucidate how and how strongly it interacts with water, glycerol, and erythritol. It is found that erythritol, which binds deep inside the conducting pore of PbAQP/PfAQP, inhibits the channel protein's functions of conducting water, glycerol, etc. This points to the possibility that erythritol, a sugar substitute, may inhibit the malarial parasites in rodents and in humans.


INTRODUCTION
The fluctuation-dissipation theorem (FDT) is a corner stone of statistical mechanics. Since the Einstein relation for the Brownian motion [1], it has been fully established with a great many variations and applications. See, e.g., [2][3][4][5][6][7][8][9][10]. Amazingly and exactly, FDT relates the dissipative characteristics of a non-equilibrium system (such as diffusivity for mass flux in response to a concentration gradient and conductivity for electric current in response an electric field) to the equilibrium fluctuations of the corresponding microscopic quantities (velocities and currents, respectively). The former (diffusivity/conductivity) gives how a system responds to a perturbation imposed on the system while the latter represents the statistics of the microscopic quantities (velocities/currents) that constantly fluctuate even in an equilibrium state. In theoretical-computational studies, we can calculate diffusivity via the Einstein relation from the time-correlation function of the particle velocities sampled in molecular dynamics simulations. We can also compute conductivity via the Kubo formula from the time-correlation function of the microscopic current densities. Conversely, in this paper, I aim to elucidate the equilibrium properties of a system from how the system respond to a prescribed set of perturbations (condition changes) imposed on the system. Namely, the goal of this study is to compute the Gibbs free energies under biochemically different conditions: when a relevant small molecule is bound to a protein (bound state) and when that small molecule is far away from the protein (apo state). Specifically, the protein of interest here is the membrane protein, aquaglyceroporin PbAQP, that facilitates transport of water, glycerol, ammonia, and several other small neutral molecules across the plasma membrane of Plasmodium berghei. The strategy of this study is to force/pull the relevant small molecules from the bound state to the apo state in allatom simulations and quantify how the system respond to such imposed perturbations, from which the Gibbs free energies of binding will be extracted.
Plasmodium falciparum infects humans and causes the most severe form of malaria leading to several 100 thousands of deaths per year even today [11]. Plasmodium berghei infects rodents and causes malarial symptoms that are comparable to P. falciparum malaria. Since genetic engineering in P. berghei is more achievable, P. berghei has been exploited as model organism for the study of P. falciparum malaria in humans. For example, P. berghei aquaglyceroporin (PbAQP) gene knock-out has been conducted to determine the essential roles played by this channel protein in the liver and blood stages of the parasite infection in laboratory mice [12,13].
Maintenance of hydrohomeostasis, uptake of nutrients, and excretion of metabolic wastes across the cell membrane are essential for all cells to proliferate or even to survive. Some of these essential tasks are accomplished by two subfamilies of aquaporins (AQPs) [14][15][16][17][18][19][20]: the water channel proteins (aquaporins) and the glycerol channel proteins (aquaglyceroporins) that also conduct water. In fact, the human body expresses a total of 13 types of AQPs, having multiple types of AQPs in a single cell of various tissues to redundantly fulfill the essential tasks of conducting water, glycerol, etc. Intriguingly, P. berghei only expresses a single type of AQP, the multifunctional aquaglyceroporin PbAQP for the multiple essential tasks and so does P. falciparum [21][22][23][24][25][26][27][28]. In a recent study [13], Promeneur et al. determined the role played by PbAQP in the liver stage of P. berghei, adding to their previous study in the blood stage [12]. They elucidated that PbAQP facilitates the uptake of exogenous glycerol for the rapid phospholipid syntheses in hepatic merozoites. PbAQP-null parasites are less virulent than the wild type: Delayed prepatent period, slower replication, and smaller merosomes. These findings highlight the possibility of curing malaria by inhibiting P. berghei/falciparum with PbAQP/PfAQP inhibitors.
In this article, we present an atomistic computational study of PbAQP to elucidate how and how strongly this channel protein interacts with water, glycerol, erythritol, etc. The computed results show that glycerol and erythritol both bind deep inside the conducting pore of PbAQP near the Asn-Leu-Ala (NLA) and Asn-Pro-Ser (NPS) motifs, inhibiting the channel protein's functions of conducting water and excreting ammonia etc. Analyzing the fluctuations and the interactions with the channel, erythritol is shown to fit better than glycerol to the intra-channel space near the NLA-NPS motifs. Consequently, erythritol's affinity of binding is 50 folds stronger than glycerol. Therefore, erythritol competitively inhibits glycerol uptake through PbAQP in addition to its inhibition of water conduction and ammonia excretion through this channel protein. This prediction of erythritol inhibiting PbAQP based on the homology modeling is consistent with earlier studies of two other aquaglyceroporins based on their crystal structure, PfAQP and the Escherichia coli aquaglyceroporin GlpF [21,23,[29][30][31][32][33][34]. In particular, the computational studies of PfAQP showed that the free-energy landscape of water was flat throughout this protein channel [30,[34][35][36] in accordance with its high permeability [21]. The glycerol-PfAQP affinity was computed to be in the low to sub mM range [34,37] which is sufficiently strong to interfere with water transport through the channel [38]. The erythritol-PfAQP affinity was computed to be many folds stronger than glycerol [34] and thus prediction was made that erythritol inhibits PfAQP's multiple functions from water conduction to glycerol uptake. All these insights consistently point to the possibility that erythritol, a sugar substitute, may inhibit the malarial parasite in rodents and in humans.

Theoretical Formulation
Biophysical/biochemical processes occur down the Gibbs free energy gradient. Therefore, it is necessary and effective to map out the free energy profile in terms of the potential of mean force (PMF) that is defined as the Gibbs free energy of the entire system when a small number of degrees of freedom selected to characterize a process (namely, the reaction coordinates) are located at a set of given values. The PMF gradients (or changes) along a relevant path give a quantitative characterization of the biophysical/biochemical process. The PMF is a function of state. The PMF difference between two states (State A and State B) should be independent of the paths connecting A and B. This basic property allows for non-equilibrium sampling of forced paths connecting A and B in the computation of the PMF change, PMF = PMF (B) − PMF(A). In this paper, we employ the Brownian dynamics fluctuation-dissipation theorem (BD-FDT) [39] to compute the PMF differences from the non-equilibrium work done to the system when the system is forced along a given path from A to B and the same from B to A: Here k B is the Boltzmann constant. T is the absolute temperature of the thermal bath the system is coupled to (via the stochastic random forces in the Brownian or Langevin dynamics). W A→B is the work done to the system when it goes from A to B. Likewise, W B→A for B to A. The brackets indicate statistical mean over the paths sampled between A and B. In the Supplementary Material 1, we give mathematical details of an exactly solvable model system to illustrate the validity of BD-FDT. The focus of this paper, however, is the application of BD-FDT to the biophysical processes facilitated by PbAQP. It is noted that the BD-FDT formula has been proven for quasi-static paths connecting A and B sampled by slow pulling and for arbitrary non-equilibrium paths sampled by fast pulling as well [39].

System Setup
This study was based on an all-atom model of PbAQP embedded in the cell membrane. Using the PfAQP crystal structure (PDB code: 3C02) without the ligands, a homology model of PbAQP was built using the SWISS-MODEL workspace [40] (modeling details given in the Supplementary Material 2). Following the well-tested steps of the literature [e.g., [30,32,41,42]], the PbAQP tetramer, was embedded in a patch of fully hydrated palmitoyloleylphosphatidyl-ethanolamine (POPE) bilayer. The PbAQP-POPE complex was sandwiched by two layers of water, each of which is ∼40 Å thick. The system was then neutralized and salinated with Na + and Cl − ions to a salt concentration of 150 mM. The entire system, consisting of 145,302 atoms, was 108 Å × 107 Å × 119 Å in dimension when fully equilibrated (illustrated in Figure S1). The Cartesian coordinates were chosen such that the origin is at the geometric center of the PbAQP tetramer. The xy-plane is parallel to the lipid-water interface and the z-axis is pointing from the periplasm to the cytoplasm. The PbAQP-membrane system with a glycerol in each of the four channels of the PbAQP tetramer was built by replacing the pore-filling waters in the NLA-NPS regions with glycerol molecules. Likewise, the system with an erythritol in each of the tetrameric channels.

Simulation Parameters
All the simulations of this work were performed using NAMD 2.13 [43]. The all-atom CHARMM36 parameters [44, 45] were adopted for inter-and intra-molecular interactions. The pressure and the temperature were maintained at 1 bar and 298.15 K, respectively. The Langevin damping coefficient was chosen to be 5/ps. The periodic boundary conditions were applied to all three dimensions, and the particle mesh Ewald (grid level: 128 × 128 × 128) was used for the long-range electrostatic interactions. The time step of 1 fs was used for short-range interactions and 4 fs was used for long-range forces. The cut-off for long-range interactions was set to 10 Å with a switching distance of 9 Å. In all simulations, the alpha carbons on the trans-membrane helices of PbAQP were fixed to their respective z-coordinates after the systems were fully equilibrated (their x-and y-degrees of freedom are not constrained, free to fluctuate in the xy-planes parallel to the membrane).

Free-Energy Profiles and Affinities
The free-energy profile of water/glycerol/erythritol, in terms of the PMF as a function of the center-of-mass z-coordinate of the molecule, was computed from steered molecular dynamics (SMD) runs over multiple sections as described in [33] except that, in each section, the z-coordinate is advanced by 1.0 Å at a slower SMD speed of 1.0 Å/ns. This slow pulling speed was used so that the computed results converge with known accuracy. The dissociation constant (inverse binding affinity) of glycerol/erythritol was computed from the PMF difference between dissociated state (when the molecule is in the cytoplasm away from the protein) and the bound state (when it is inside the channel near the NLA-NPS motifs), the z-fluctuations at the binding site, and the xy-fluctuations at the cytoplasmic orifice of the channel (shown in Figures S2, S3). At the binding site, the small molecule fluctuates as governed by the free stochastic Brownian/Langevin dynamics with no constraints. However, the interactions with the channel-lining residues of PbAQP confine its x-and y-coordinates and larger fluctuations along the z-axis. At the cytoplasmic orifice (z = 10 A), the z-coordinate was fixed while the x-and y-coordinates can freely fluctuate but they are confined by the interactions with the PbAQP residues. All these fluctuations give rise to normal diffusion [10] within the available domains bounded by the channel protein.

PbAQP Conducts Water and Small Polar Molecules
The top panel of Figure 1 shows how waters interact with the pore-lining residues that form a PbAQP channel. The coordinates are from the last frame of the 60 ns equilibrium molecular dynamics (MD) run of the PbAQP-membrane system without glycerol/erythritol after 20 ns initial run of setting up the systems. Movie 1 illustrates the dynamic behavior during the 40 ns MD run starting from the last frame of the 60 ns MD run. Like in all other AQPs (especially PfAQP), waters line up in a single file inside each PbAQP channel. Even though there are rooms near the NLA-NPS motifs where the two half transmembrane helices meet (1Å < z < 6Å in Figure 1), waters stay in single file due to hydrogen-bond networks with themselves and with the porelining residues of the channel. All these are in agreement with the literature of water and glycerol-conducting aquaglyceroporins [21,29,30], which indicates accuracy of our homology modeling. The bottom panel of Figure 1 shows the Gibbs free energy of the system, namely, the PMF as a function of a water molecule's z-coordinate when it traverses the PbAQP channel. This freeenergy profile correlates exactly with the single-file lining of waters inside the channel (Figure 1, top panel). The seven minima in the PMF curve are located at the single-file lining of waters (−8Å < z < 10Å). Quantitatively, the difference between the highest PMF peak and the deepest PMF well inside the channel is <2.0 kcal/mol. This very low free-energy barrier means that water can traverse the PbAQP channel as efficiently as free diffusion, in agreement with the current literature that aquaglyceroporins conduct water as well as (or better than) the water-specific AQPs [47].
It is worth noting that ammonia, a metabolic waste, is similar to water in size and polarity. It can move through PfAQP as easily as water [26,34,37,48] and thus is expected to permeate PbAQP because PbAQP is exactly parallel (in channel characteristics) to PfAQP [21,23] and GlpF [29,30]: Lining waters up in single file and disallowing them to pass over one another inside the channel despite their dynamic fluctuations [47].

Glycerol or Erythritol Clogs up the PbAQP Channel
To illustrate the glycerol-PbAQP interaction, we took the full equilibrated structure of the PbAQP-membrane system, deleted the waters near the NLA-NPS motifs, placed a glycerol inside each    of the tetrameric PbAQP channels, and ran MD for 120 ns. After the initial 20 ns, the channel characteristics of PbAQP (illustrated in Figure 2 and Movie 2) remain unaltered, in close resemblance of the PfAQP channel characteristics: Waters remain in single file and glycerol stays inside the channel fluctuating near the NLA-NPS motifs. The glycerol so dwelling inside a channel does not cause distortions to the pore-lining residues but completely severs the single-file line of waters into two. Therefore, it inhibits PbAQP's function of conducting water and other small polar solutes including ammonia, which is similar to what was computed for PfAQP [34]. This close resemblance is because PbAQP is 62% identical to PfAQP in amino-acid sequence. In functional assays, PfAQP was found to conduct water, ammonia, glycerol, etc. [21,26]. In the crystal structure of PfAQP, a glycerol was located inside each channel near the NLA-NPS motifs interrupting the single-file lining of waters [23]. Here, we observe that PbAQP acts likewise.
To quantify the glycerol-PbAQP interaction, we computed the Gibbs free energy of binding, namely, the difference in the free energies of the system when the glycerol resides inside the PbAQP channel (shown in the left panel of Figure 2) and when it is in dissociated state (shown in the right panel of Figure 2). The binding free energy involves the PMF along the dissociation path of pulling a glycerol from the bound state to the dissociated state (shown in the central panel of Figure 2) and the fluctuations of a glycerol when inside the channel (shown in Figure 3). The entire dissociation path over a z-displacement of 25 Å was sampled as 25 sections of 1.0 Å each. Over each section, four forward paths of pulling a glycerol along the z-axis at a speed of 1.0 Å/ns and four reverse paths of pulling a glycerol along the opposite direction at the same speed. At each end point of the 25 sections, the system was equilibrated for 4.0 ns while the z-coordinate of a glycerol was held constant. Altogether, we conducted 300 ns SMD runs for the PMF curve shown in Figure 2. The PMF and the fluctuations (shown in Figure 3) combined to give the Gibbs free energy of glycerol-PbAQP binding at −5.55 ± 2.07 kcal/mol.
Replacing the four glycerols inside the tetrameric PbAQP channels with four erythritols, the 120 ns MD run showed that erythritol is also stable at the same binding site near the NLA-NPS motifs (shown in the left panel of Figure 4 and Movie 3). Like glycerol, an erythritol located inside the channel does not cause distortions to the pore-lining residues but occludes the channel from conducting water or small solutes. Unlike glycerol, the fluctuations of erythritol are smaller   (compare Figure 5 with Figure 3), indicating that erythritol inhibits PbAQP's conduction of water and ammonia more effectively than glycerol.
The SMD runs for computing the erythritol PMF are parallel to the case of pulling glycerols for 300 ns. The central panel of Figure 4 shows the PMF along the path of dissociating erythritol from the binding site (Figure 4, left panel) to the intracellular space (Figure 4, right panel). The PMF (shown in Figure 4) and the fluctuations (shown in Figure 5) combined to give the Gibbs free energy of erythritol-PbAQP binding at −7.87 ± 1.55 kcal/mol.

Erythritol Competitively Inhibits Glycerol Uptake Through PbAQP
Recent experiments of Promeneur et al. [13] showed that P. berghei relies on the PbAQP channel for efficient infection and proliferation. The parasite uses PbAQP for glycerol uptake required by its high demand of lipids. Inhibiting PbAQP with an inhibitor that has a higher affinity than glycerol would be similar to knocking out the PbAQP gene to disable the parasite's uptake of glycerol and excretion of ammonia through its sole aquaglyceroporin PbAQP. Therefore, it is interesting to compare the glycerol-PbAQP affinity and the erythritol-PbAQP affinity for the possibility that erythritol competitively inhibits glycerol uptake through PbAQP as they have the same binding site (Figures 2, 4, 6).
In Table 1, we summarize the glycerol-and erythritol-PbAQP interactions. Erythritol's binding free energy is −2.3 kcal/mol stronger than glycerol, which is outside (but not far from) the margin of error. Correspondingly, the erythritol-PbAQP affinity (dissociation constant k D = 2 µM) is ∼50 times the glycerol-PbAQP affinity (k D = 0.1 mM). It is reasonable to conclude that erythritol competitively inhibits glycerol uptake through PbAQP. At a typical blood glycerol level, e.g., 0.1 mM, erythritol at 1.0 mM concentration would inhibit PbAQP's conduction of glycerol to <1%.

Atomistic Interactions and Steric Fit
Analyzing the relevant atomistic interactions, we can see how and how strongly a glycerol/erythritol interacts with the protein and its aqueous environments. We can also see that the uncertainties of this computational work are probably lower than what might be suggested by the margins of error in the computed free energies ( Table 1).  First, the hydrogen bonds. Glycerol (erythritol) is hydrophilic. Away from the protein where a glycerol (erythritol) molecule is surrounded by waters, it can form six (eight) hydrogen bonds with the surrounding waters. When it is inside the protein channel, it forms two hydrogen bonds with waters and one or two hydrogen bonds with the pore-lining residues. This seems to suggest that hydrogen-bond interactions favor the unbound state of glycerol (erythritol) where it has three (four) more hydrogen bonds to its three (four) hydroxyl groups, if we ignore the hydrogen bonds between waters. When a glycerol (erythritol) is inside the channel, it displaces two or three (three or four) waters that are bonded to each other and to the pore-lining residues by three or four (five or six) hydrogen bonds. Namely, in the bound state, a glycerol (erythritol) breaks four (six) hydrogen bonds to waters. In the dissociated state, a glycerol (erythritol) displaces four (five) waters from the volume it occupies and thus disrupts 10 (12.5) hydrogen bonds between the waters if assuming 2.5 hydrogen bonds per water in the aqueous bulk. Altogether, the system has six (6.5) more hydrogen bonds when a glycerol (erythritol) is inside the protein than in the aqueous bulk. Therefore, hydrogen-bond interactions favor the binding of a glycerol (erythritol) inside the channel.
Second, the entropic forces. There is enough room inside the channel near the NLA-NPS motifs to accommodate a glycerol (erythritol). In fact, the dihedral energy of a glycerol (erythritol), shown in Figure 7, simply fluctuates around the same level along the dissociation path leading from the NLA-NPS region to cytoplasmic bulk. Otherwise, the dihedral energy would be significant higher in the bound state inside the channel than in the dissociated state outside the protein. Not surprisingly, the glycerol fluctuations at the binding site is higher than erythritol (the root mean square fluctuation values, RMSF, are shown in Table 1). The room there is greater than the volume a glycerol occupies but is nearly a perfect fit to accommodate an erythritol. In this aspect, the entropic contribution to the free energy of binding favors glycerol over erythritol by −1.1 ± 0.3 kcal/mol (which is estimated as k B Tln RMSD erythritol RMSD glycerol 3 with values given in Table 1).
Third, the van der Waals (vdW) interactions. These forces were computed with the smallest uncertainties in this work. In fact, the vdW interactions between a glycerol (erythritol) and the protein are all attractive. The erythritol-PbAQP attraction is −18.60 ± 1.29 kcal/mol and the glycerol-PbAQP attraction is −14.80 ± 1.01 kcal/mol. The separation between these two vdW attractions are far outside the margins of error.
Overall, both glycerol and erythritol bind inside the PbAQP near the NLA-NPS motifs mostly due to the vdW interactions and the hydrogen-bond compensations. The erythritol-PbAQP vdW attraction is stronger than the glycerol-PbAQP vdW attraction by 3.8 kcal/mol, which is the dominant factor leading to the 50 folds stronger affinity of erythritol than glycerol.
Based on the computed binding free energy of erythritol, −7.87 ± 1.55 kcal/mol (Table 1), the possible values of its k D range from 0.15 to 27 µM. Therefore, at a concentration in the low or sub mM range, erythritol should strongly inhibit P. berghei. It is particularly interesting that erythritol is a sugar substitute generally regarded as safe. If proven effective to inhibit P. berghei in laboratory rodents, it is expected to inhibit P. falciparum in humans without significant adverse side effects.

DATA AVAILABILITY STATEMENT
The datasets including all the parameters, coordinates and scripts to reproduce the results of this study is available from the Harvard Dataverse: https://doi.org/10.7910/DVN/ ZOAQJW.

AUTHOR CONTRIBUTIONS
The author conducted the simulations, analyzed the data, and wrote the manuscript.

FUNDING
The author acknowledges NIH Grant #GM121275.