# Influence of water models on water movement through AQP1

Miguel A. Gonzalez<EMAIL_ADDRESS>Universidad Rey Juan Carlos,
ESCET, 28933, Mostoles-Madrid, Spain  Alberto Zaragoza  Department of
Chemistry, The University of Utah, Salt Lake City, Utah 84112-0850, USA
Charlotte I. Lynch  University of Oxford, Department of Biochemistry, South
Parks Road, OX1 3QU, Oxford, UK  Mark S.P. Sansom  University of Oxford,
Department of Biochemistry, South Parks Road, OX1 3QU, Oxford, UK  Chantal
Valeriani  Universidad Complutense de Madrid, Facultad de Ciencias Fícias,
Departamento de Estructura de la Materia, Física Térmica y Electrónica, 28040,
Madrid, Spain

###### Abstract

Water diffusion through membrane proteins is a key aspect of cellular
function. Essential processes of cellular metabolism are driven by osmotic
pressure, which depends on water channels. Membrane proteins such as
aquaporins (AQPs) are responsible for enabling water permeation through the
cell membrane. AQPs are highly selective, allowing only water and relatively
small polar molecules to cross the membrane. Experimentally, estimation of
water flux through membrane proteins is still a challenge, and hence accurate
simulations of water permeation are of particular importance. We present a
numerical study of water diffusion through AQP1 comparing three water models:
TIP3P, OPC and TIP4P/2005. Bulk diffusion, diffusion permeability and osmotic
permeability are computed and compared among all models. The results show that
there are significant differences between TIP3P (a particularly widespread
model for simulations of biological systems), and the more recently developed
TIP4P/2005 and OPC models. We demonstrate that OPC and TIP4P/2005 reproduce
protein-water interactions and dynamics in very good agreement with
experimental data. From this study, we find that the choice of the water model
has a significant effect on the computed water dynamics as well as its
molecular behaviour within a biological nanopore.

##  I Introduction

Water is the most abundant molecule in cells. Even though water is able to
diffuse slowly through a lipid bilayer, for essential physiological processes
a high water flux is required. This can be achieved making use of protein
membranes such as aquaporins (AQPs) which allow for a high selectivity for
water permeation across biological membranes [1, 2]. AQPs are located in
several cells found in the brain, kidneys, skin, blood vessels, liver and
connective tissue. [3, 4] The lack of functionality of such cells might induce
several diseases, such as glaucoma, brain edema or congestive heart failure.
[3, 5, 6] Although evidence of the existence of water channels had been
reported several years earlier, the aquaporin structure was not described
until 1992 [7]. All members of the AQP family are small ($\sim$ 30 kDa)
intrinsic membrane proteins with a strongly hydrophobic character, and consist
of four monomers. Each monomer is a functional water channel having a
hourglass shape [8] with an amphiphilic interior and a selectivity filter. The
primary sequence of those proteins shows six trans-membrane helices (I-VI)
linked by five loops. The water flux inside the channel depends on an
asparagine-proline-alanine (NPA) motif, found in two loops located
approximately in the middle of the membrane protein channel. The interaction
between these motifs not only confers a particular shape to the channel but it
is also responsible for the orientation of water molecules inside the channel.
In combination with the NPA motif, water permeation is also modulated by the
aromatic/arginine selectivity filter which is situated next to the
extracellular exit/entrance of the channel. [9, 10] The topology of the pore
is extremely important as it enables selected small polar molecules to cross
the membrane (e.g. H2O, NH3, glycerol) whilst preventing the passage of others
(e.g. ions and larger molecules). Based on its permeability, AQPs can be
classified in two subfamilies: classical aquaporins (permeable to water) and
aquaglyceroporins (permeable to both water and glycerol).

AQP1 is a classical aquaporin with a high water permeability as described in
detail in Ref. [11] and confirmed in several numerical studies [12, 9, 10, 13,
14, 15]. There are three different approaches presented in the literature to
study water through aquaporins: 1) a continuum hydrodynamic approach [16] that
performs surprisingly well for a nanopore of molecular dimensions, capturing
key steric effects; 2) a potential of mean force approach, to study the
energetics of water transport [17, 14]; 3) “a continuous-time random-walk
model” approach of water transport across the channel [18].

As suggested by De Groot et al. in Ref. [19], approaches 2 and 3) are in a
quantitative agreement. Thus, we will follow approach 3 throughout our study
(extracting equations from Ref. [20]). Numerous molecular dynamics studies
have been conducted on the water permeation mechanism through aquaporins [12,
21, 22, 23, 24, 25, 26, 27, 11, 28, 29]. As discussed by Ozu et al. [30], many
AQP simulations have been carried out using different force fields (GROMOS or
CHARMM) for the protein and the standard older TIP3P or SPC water models. De
Groot and Grubmüller [21] demonstrated that proton translocation across AQPs
was prevented by an electrostatic barrier. However Jensen et al. [27], using a
different force field to simulate AQP and water (TIP3P), concluded that the
proton induced an interruption of the hydrogen bond chain network built into
the channel, preventing the proton from crossing the pore. Thus, it is
reasonable to suggest that models of water permeation might be affected by the
force field used to simulate water molecules.

When dealing with pure water, both Vega et al. [31] and Onufriev et al. [32]
carried out a detailed comparison between experimental and numerical values of
thermodynamic and dynamic properties, simulating water by means of several
interaction potentials. Even though they clearly demonstrated that TIP3P was
not one of the best performing potentials for pure water, TIP3P remains a
widely used water model in biomolecular simulations, in combination with
CHARMM (the most widely used protein force field). With this in mind, we
propose to study the molecular mechanism of water permeation of AQP1 (1H6I
code PDB) [12] simulated with CHARMM36m [33], in combination with the TIP3P,
OPC and TIP4P/2005 water models. Even though the latter two water models are
less frequently used in biomolecular simulations, they are known to give very
reliable thermodynamic and dynamic properties for bulk water. To our
knowledge, TIP4P/2005 has been already implemented with the AMBER ff03w [34]
force field as a default water model but has not been employed in simulations
of AQP, where the microscopic behaviour of water clearly plays a vital role.
It is important to note that we have chosen to focus our study on rigid non-
polarizable water models. There is a range of potentially more accurate
polarizable and flexible water models, however the use of these comes at
significant computational cost [35]. It is therefore beneficial to first
assess the consistency and accuracy of a range of rigid non-polarizable models
in order to establish whether the use of polarizable water models is needed
for the study of water behaviour in aquaporins. Similar comparisons of non-
polarizable water models have been conducted for the 5-HT3 receptor channel
[36] and Cucurbit[7]uril-guest systems [37]. In this article, we have studied
several properties of human AQP1 (water diffusion across the pore, diffusion
permeability and osmotic permeability) taking into account the molecular
mechanism characteristic of AQP. In this way we have selected to simulate a
range of representative (non-polarizable) models of water, comparing two more
recent water potentials to what is still the “standard model”. This allows us
to capture the likely range of behaviours for this class of model.

##  II Methods

###  II.1 Simulation details

We have carried out molecular dynamics simulations for AQP1 embedded in a
palmitoyloleoylphosphatidylcholine (POPC) bilayer membrane using the GROMACS
2016.4 software package [38]. Human AQP1 (1H6I code PDB) [12] is simulated via
the CHARMM36m force field [33], whilst the water molecules are simulated using
TIP3P [39], OPC [40], and TIP4P/2005 [41] potentials. In order to maintain the
protein fold, the simulations were carried out with backbone restraints.
Figure 1 a) shows an image of the system under study, representing the
different regions where water molecules can be found: bulk water molecules (in
brown), interfacial water molecules (in red), lipids (in blue) and AQP1 (in
purple). Panel b) represents a rendered image of the AQP1 channel, showing
water molecules (in red/white) passing through, while changing their
orientation when crossing the center.

a) b)

Figure 1: a) Rendered image of the simulated system (from top/bottom towards
the middle): bulk water molecules (in brown), interfacial water molecules (in
red), lipids (in blue) and AQP1 (in purple). b) Rendered image of the AQP1
channel, showing water molecules passing through, while changing their dipolar
moment when crossing the center. The water molecules are plotted in red/white.

Considering the three models of water, we have set up three initial
configurations with the AQP1 tetramer embedded in a lipid bilayer of 209 POPC
molecules and solvated with approx. 20000 water molecules. In order to prepare
a neutral configuration, we have added 4 Cl- ions. Each AQP1-water system was
simulated for 100 ns, and simulations were repeated four times with different
random initial velocities. We set a time step of 1 fs keeping the temperature
constant at 298 K using a velocity rescale thermostat [42] with a relaxation
time set to 0.67 ps. The Parrinello-Rahman barostat [43] was applied to
maintain the pressure constant at 1 bar using a 2.57 ps relaxation time. When
dealing with water, the Lennard-Jones (LJ) potential was truncated at 1.2 nm,
adding standard long-range corrections to the LJ energy. Ewald sums [44] were
considered to compute long-range electrostatic forces with a real space cut-
off at 1.2 nm. In order to compare the behaviour of the three water models and
their interactions with AQP1, we have selected the most relevant properties
that allow us to establish the influence of water potentials on water dynamics
through protein channels.

###  II.2 Diffusion, diffusion permeability, rate of water molecules and
average dipole moment across the channel

Water diffusion through a cell membrane is considered as a subdiffusive
process [45, 46, 47, 48, 49]. Also, Yamamoto et al.[50] demonstrated that the
water diffusion near a lipid membrane is lower than the bulk due to two
mechanisms. On the one hand, there is a divergent mean trapping time induced
by the lipid-water interaction which can be described by a continuous-time
random walk. On the other hand, the viscoelasticity of the lipid membrane
induces a long-time correlated noise described by fractional Brownian motion.
In order to compute water diffusion, we distinguished four regions depending
on the distance ($d$) between the lipid headgroup and the position of the
water molecule. As in Ref. 48, the criteria applied to create those regions
are the following: we identify a 1) bulk region whenever $d$ $>$ 0.6 nm, an 2)
interface region whenever 0.3 nm $<d<$ 0.6 nm, a 3) contact region whenever
$d<$ 0.3 nm) and an 4) AQP region defined as the region inside the lipid
headgroups (corresponding to waters in the pore). The long time diffusion
coefficient ($D$) is then computed using Einstein’s relation based on the
water mean square displacement:

$D=\frac{\left\langle\left|r(t_{0}+t)-r(t_{0})\right|^{2}\right\rangle}{2\cdot
dim\cdot t}$ (1)

where $r(t)$ is the position of the water molecule centre of mass, $t$ is the
time, and $dim$ the dimension that depends on the geometry of the system that
the molecules are moving through, with $dim=1$ for 1D geometries such as
carbon nanotubes CNT or AQP1 and $dim=3$ for bulk [51].

Another approach to calculate water diffusion inside the protein channel (AQP
region) is based on a single-water molecule column [18]. In this case, water
diffusion ($D_{z}$) is computed according to the relation [52]:

$D_{z}=\frac{k_{0}z^{2}}{2}=\frac{z^{2}p_{f}}{2\nu_{w}}$ (2)

where $k_{0}$, $z$, $p_{f}$ and $\nu_{w}$ are the rate of water molecules
across the channel, the average distance between two water molecules inside
the channel, the osmotic permeability (see below for more details) and the
volume of a single molecule, respectively.

Water permeability can be divided into osmotic permeability ($p_{f}$) and
diffusion permeability ($p_{d}$) depending on the system’s conditions. Osmotic
pressure drives the movement of water molecules through the membrane protein
when there is a solute gradient across the membrane. Hence in these conditions
water permeability would be given by the osmotic permeability of the system.

As employed in previous studies, [53, 19, 20] osmotic permeability is defined
as the proportionality constant between the net water flux ($j_{w}$) across
the protein channel and the solute concentration difference ($\Delta C_{s}$),

$j_{w}=p_{f}\Delta C_{s}.$ (3)

Taking into account the rate of water molecules across the membrane protein,
$k_{0}$, the osmotic permeability is also defined as

$p_{f}=\nu_{w}k_{0}.$ (4)

At equilibrium and in the absence of solute or solute gradient, there is still
net water moving through the channel that could be attributed to diffusion of
water molecules due to temperature. Thus, the diffusion permeability, $p_{d}$,
allows us to study the water flux across the protein under these conditions:

$j_{tr}=p_{d}\Delta C_{tr}$ (5)

where $\Delta C_{tr}$ corresponds to the difference in the tracer’s
concentration between the two sides of the membrane and $j_{tr}$ is the net
tracer’s flux (where a tracer is defined as a labelled water molecule that can
be distinguished from similar molecules, i.e. a molecule of heavy water is a
tracer for water). At equilibrium, the rate of water molecules crossing the
channel is computed by $q_{0}$, which is related to $p_{d}$ via

$p_{d}=(V_{w}/N_{a})q_{0}=\nu_{w}q_{0}.$ (6)

Note that the centre used to decide if a water molecule crossed the channel is
the average of Z position of N atoms for ASN residue inside the channel. One
molecule has been considered that it crossed the AQP when its Z coordinate is
1.5 nm from this point. This centre is updated each single time step during
the analysis.

Molecular simulations become an essential tool for computing $p_{d}$ as they
allow us to label individual water molecules and therefore distinguish them
from one another. Both properties, $p_{d}$ and $p_{f}$, are proportional to
the number of water molecules, $N$, taking part in the translocation following
the continuous-time random-walk process [18] based on Ref. [54]. For AQP1, Zhu
et al. [20] presented a value for $p_{f}/p_{d}=N+1=11.9$. This is in good
agreement to the experimental value of $p_{f}/p_{d}=13.2$ [55]. However, these
values are higher than the average number of molecules inside the protein
channel, since there are molecules of water outside the channel that
participate in the permeation across the membrane.

In order to compute the flux in Eq. 5 we need to compute $p_{d}$ and $\Delta
C_{tr}$. The latter is computed by simple counting whereas the former is
computed via Eq. 6. To compute $q_{0}$ in Eq. 6 we count the number of water
molecules that cross the membrane through the AQP1. Thus, we define the centre
of the protein and a typical distance from which a molecule is considered to
be outside the channel. The channel centre is identified in the following way:
two nitrogen atoms inside the channel are bound to water via hydrogen bonds
and their orientation is essential to describe the AQP1 permeability. Thus, we
compute the channel centre as the average of these nitrogen atom positions.
The channel length is set at 1.5 nm from the centre.

Water molecules in bulk are characterised by an average dipole moment equal to
zero. However, in nano-confined structures such as a CNT [51] or membrane
channels [11], there is a particular orientation of the molecules due to the
hydrogen bond network. Thus, the dipole moment is an important parameter to
study the molecular orientations inside the channel, both artificial and
natural. In our work, we have computed the water dipole moment $\mu$ according
to $\vec{\mu}=\sum_{i}q_{i}\vec{r_{i}}$, where $q_{i}$ and $r_{i}$ are the
charge and the position vector for molecule $i$, respectively. We focus on the
$z$ component of $\vec{\mu}$, since in our simulations the orientation of the
water molecule inside the channels is mainly projected onto this axis. Hence,
we calculate the average dipole moment of slabs, having a thickness of 0.025
nm, across the $z$ axis to compare the molecular orientation between the
channel and molecules.

##  III Results And Discussion

###  III.1 Water molecule orientations inside AQP1

As suggested in Ref.[11], water molecules arrange inside the channel forming a
single-molecule chain with a specific orientation. In our work we demonstrate
that this mechanism is reproduced by the three water models.

As shown in Fig. 2, the water molecule orientation can be illustrated by
computing the average of the $z$ component of the dipole moment inside the
channel.

Figure 2:  $z$ component of the average dipole moment $\mu_{z}$ for the three
water models. The red line is for TIP3P, the blue line for TIP4P/2005, and the
black line for OPC. Error bars have been added taking into account the
standard deviation of $\mu_{z}$. The vertical dashed lines indicate the
approximate locations of the Asn residue of the NPA motifs. Note that the
centre of the bilayer is at z = 0 nm, and the pore runs from approximately z =
-1 to +0.8 nm.

Although there is a dynamic permeation across the channel, it is possible to
distinguish the positions of water molecules, described by the slightly ragged
shape of the water dipole moment inside the protein (see Fig. 2 ). These
results demonstrate that there are two stable orientations inside the channel
as described by Ref. [11].

There are no major differences between the dipole moment profiles for the
three water models, although the degree of water dipole orientation is less
for the OPC model in the second (i.e. $z>0$) half of the pore than for the
other two water models. This implies that, even though TIP3P is the only water
model fully “tuned” with CHARMM36, both OPC and TIP4P/2005 are capable of
capturing a qualitatively similar behaviour inside the AQP. Therefore, one
could deduce that the similarity between the dipole moments computed with the
different water models is due to water interacting in a similar way with AQP1.
Water interacts with AQP via hydrogen bonds that can be only formed when water
is close enough to the asparagines located at the middle of the pore. Our
results on the dipolar moment suggest that water interacts with asparagine in
a similar way independently of the model used to mimic its behaviour.

###  III.2 Water molecule diffusion through the protein channel

As in Ref.48, we identify interfacial, pore and contact regions, depending on
the position of the water molecule with respect to the bilayer and the AQP1
pore. To compute the water diffusion constant ($D$) we use Eq. 1 with $dim=3$
for the bulk ($D_{bulk}$), interface ($D_{inter}$), and contact
($D_{contact}$) regions, and $dim=1$ for the AQP ($D_{AQP}$) region (AQP1
being almost cylindrical).

(a)

(b)

(c)

Figure 3:  Mean square displacement (nm2) for TIP3P (left), TIP4P/2005
(centre), and OPC (right). The coloured lines denote the following regions:
bulk (red), interface (blue), lipid-water contact region (green), and AQP
(black).

Moreover, we compare the value of the diffusion computed in the bulk region at
$d>0.6$ nm ($D_{bulk}$) to the one computed for each model in bulk
($D_{model}$), since the presence of the lipid membrane reduces the value of
$D$ in nearby regions (even at distances higher than 2 nm [50]). All mean
square displacements (msd) for each region are reported in Fig. 3. Two groups
of data could be distinguished for each model (although in the case of TIP3P
it is more evident). The bulk, interface and contact regions have similar
slopes, while the slope for the AQP region is lower than the rest. The
diffusion coefficient (D) values for all water models in bulk and in the
various regions are presented in Table 1 (as calculated from the mean square
displacement in Fig. 3).

Table 1: Diffusion coefficient values computed using Eq. 1 with $dim=3$ in the bulk, interface and contact region, with $dim=1$ in the AQP region for TIP3P, TIP4P/2005 and OPC water models. The values in parenthesis are the standard deviations. All diffusion coefficients are expressed in 10-5 (cm2/s). Water model | Dbulk | Dinter | Dcontact | DPORE
---|---|---|---|---
TIP3P | 4.60 (0.01) | 3.74 (0.02) | 2.93 (0.02) | 0.10 (0.01)
TIP4P/2005 | 1.87 (0.04) | 1.24 (0.1) | 0.77 (0.04) | 0.05 (0.02)
OPC | 1.94 (0.04) | 1.26 (0.08) | 0.72 (0.09) | 0.02 (0.02)
Exp. | 2.30 [31] |  |  | 0.04

When dealing with bulk water, the literature values of $D_{literature}$ at
298K and 1 bar are $5.5$x$10^{-5}$ cm2/s [56] for TIP3P, $2.06$x$10^{-5}$
cm2/s for TIP4P/2005 [41] and $2.3$x$10^{-5}$ cm2/s [40] for OPC
$D_{literature}$. These data are higher than the $D_{bulk}$ computed from the
mean square displacement in our simulations (see Fig. 3), reflecting the
phenomena of subdiffusion of water close to the lipid membrane, explained in
Ref. [50], and in addition it might be due to differences in the simulation
conditions, such as the value of the cut-off and/or the thermostat applied
during the simulation. However, the value of water diffusion given in this
article is reliable. However, it is clear that TIP4P/2005 and OPC water models
reproduce the experimental bulk water diffusion better than TIP3P. Concerning
the hydrogen bonds formed in bulk water, a recent paper[57] has shown that the
TIP3P water model yields a higher rotational relaxation time and a lower
network complexity parameter. This results in a higher diffusion coefficient
for TIP3P water in bulk. In contrast, TIP4P/2005 (a water model very similar
to OPC) yields a lower rotational relaxation time and a higher network
complexity parameter. This results in a lower diffusion coefficient for
TIP4P/2005 (and by extension for OPC) water in bulk. TIP3P bulk water has a
higher diffusion coefficient than both OPC and TIP4P/2005, due to its higher
rotational relaxation time and lower network complexity parameter. This effect
is enhanced when water is confined within the single file pore AQP, as already
suggested in Ref[51].

On the one hand, although it is possible to consider TIP3P water subdiffusing
as $D$ is compared to $D_{literature}$, TIP3P overestimates bulk water
diffusion in comparison to the experimental data of $2.3$x$10^{-5}$ cm2/s in
all regions. On the other hand, TIP4P/2005 and OPC show a consistent
subdiffusion of water molecules close to the membrane, as expected [50] both
with respect to their $D_{literature}$ and in comparison to the experimental
value. Focusing on the AQP region, i.e. inside the membrane protein, the
experimental data of water diffusion is $0.04$x$10^{-5}$ cm2/s [52]. TIP3P
gives $0.10$x$10^{-5}$ cm2/s for $D_{PORE}$, which is dramatically higher than
the experimental value and also higher than the other water models,
$0.05$x$10^{-5}$ cm2/s and $0.02$x$10^{-5}$ cm2/s for TIP4P/2005 and OPC,
respectively. When comparing the diffusion coefficient of water inside the
pore with respect to its bulk counterpart, we get a comparable ratio
independent of the model. The reason is that TIP3P water diffuses more than
the other water models not only in bulk but also when confined in a nanopore.
In addition, note that the diffusion of water inside the channel is different
with respect to the experimental value depending on the chosen model. In the
case of TIP3P the diffusion is 2.5 times higher than the experimental one. OPC
also shows a significant discrepancy, although in this case the model presents
a diffusion that is 2 times lower than the experimental one. However, the
diffusion of the TIP4P/2005 model is in optimal agreement with the
experimental value.

In order to visualise the values of $D$ across the four regions, a diffusion
profile has been plotted for each water model, see Fig. 4. This shows the
diffusion profiles for the three models by considering only the molecules that
have crossed the membrane through the protein channel, as these molecules
allow us to compute the water diffusion through the channel.

(a)

(b)

(c)

Figure 4: Diffusion profiles for TIP3P (left), TIP4P/2005 (centre), and OPC
(right) water models. The profiles have been split into four regions: bulk
(red crosses), interface (blue crosses), lipid-water contact (green crosses),
and AQP region (black crosses).

It can be observed that $D$ decreases from the bulk region to the AQP region,
starting from a value of $D$ lower than $D_{literature}$. Each region displays
a distinct diffusion profile, although there is a slight overlap between
regions. This is due to the criteria used to define the regions, i.e. an
average of the phosphorus atoms position in $z$. As a consequence, it is
possible to find a water molecule having a lower value of $D$ than the one
expected in that region. As shown in Fig. 4, the AQP region is not symmetric
and there are higher values of $D_{PORE}$ at negative distances, $d$, from the
membrane centre than at positive ones. This is because the membrane protein is
not symmetric and the diffusion around the protein is asymmetric as a
consequence.

###  III.3 Diffusion permeability and osmotic permeability

To compute the water permeability, both the diffusion (pd) and osmotic
permeability (pf), we have calculated the rate at which molecules are able to
cross the membrane ($q_{0}$ and $k_{0}$) and the diffusion of those molecules
through the AQP ($D_{z}$). The results have been collected in Table 2, where
we compare numerical and experimental data.

Table 2: Number of molecules per second (q0, multiplied by 108 (s-1) ; diffusion permeability pd, by 10-14 (cm3/s) ; osmotic permeability pf, by 10-14 (cm3/s) ([*] via $p_{f}/p_{d}=11.9$); number of molecules per second for an osmotic process k0, by 109 (s-1); water diffusion inside the AQP $D_{z}$, by 10-5 (cm2/s) ([**] via Eq. 2) and $D_{PORE}$, by 10-5 (cm2/s)  Water model | q0 | pd | pf ∗ | k0 ∗∗ | Dz ∗∗ | DPORE
---|---|---|---|---|---|---
TIP3P | 4.0 (0.3) | 1.20 (0.09) | 14.2 (1.1) | 4.85 | 0.190 | 0.10 (0.01)
TIP4P/2005 | 0.65 (0.2) | 0.194 (0.07) | 2.31 (0.8) | 0.77 | 0.03 | 0.05 (0.02)
OPC | 0.5 (0.1) | 0.149 (0.04) | 1.78 (0.4) | 0.59 | 0.02 | 0.02 (0.02)
Exp. |  |  | 5.43 [58] |  |  | 0.04

Diffusion permeability ($p_{d}$) has been computed using Eq. 4, which is based
on the number of water molecules that crossed the membrane through the
membrane protein ($q_{0}$) and the volume of a single water molecule
($\nu_{w}=2.989$x$10^{-23}$ cm3). From molecular dynamics simulations $q_{0}$
is a parameter that is easily accessible. For TIP3P $q_{0}=4.0$x108 s-1
whereas for TIP4P/2005 ($q_{0}=0.65$x108 s-1) and OPC ($q_{0}=0.5$x108 s-1)
the number of molecules per second is lower, by almost an order of magnitude.
Therefore, we find 1) on the one hand an overestimation of the TIP3P diffusion
permeability and 2) on the other hand a good performance of both TIP4P/2005
and OPC diffusion permeability (see Table 2).

The osmotic permeability $p_{f}$ has been computed from the diffusion
permeability using the ratio $p_{f}/p_{d}=N+1=11.9$ proposed by Zhu et al.
[20] for AQP1 and the values are reported in Table 2. The experimental data
for $p_{f}$ is $p_{f}=5.43$x $10^{-14}$ cm3/s [58], which is in very good
agreement with the TIP4P/2005 (2.31x$10^{-14}$cm3/s) and OPC
($1.78$x$10^{-14}$cm3/s) values. The TIP3P value, 14.2x$10^{-14}$cm3/s, is
dramatically higher than the experimental value.

Having calculated the osmotic permeability, the rate of water crossing the AQP
($q_{0}$) can be calculated by applying Eq. 2 with a value of z = 0.28 nm
[52]. From the same Eq. 2, we calculated $D_{z}$ for the different water
models, obtaining 19.0x10-7 cm2/s, 3.03x10-7 cm2/s, and 2.33x10-7 cm2/s for
TIP3P, TIP4P/2005 and OPC, respectively. TIP4P/2005 and OPC compare well with
the experimental data, 4.0 x10-7 cm2/s. However, TIP3P shows a value higher
with respect to the experimental one than the other water models.

It is worth nothing that $D_{z}$ could be compared to $D_{PORE}$ using Eq. 1
as $dim=1$. This fact allows us to test the consistency of our diffusion
results, since $D_{PORE}$ is calculated via Eq. 1 and $D_{z}$ is computed via
an independent method (via Eq. 2) based on the number of water molecules that
crossed the membrane through AQP ($q_{0}$) and its ratio with respect to the
water permeability. As it can be observed, the results for water diffusion
inside the AQP are in excellent agreement independently of the method used to
compute it, either using Eq. 1 on the AQP region or applying Eq. 2.

In summary, our results provide compelling evidence that the TIP3P water model
shows a higher diffusion rate than the other water models studied. Inside the
channel, the TIP4P/2005 model produces results in excellent agreement with
experiment. Focusing on the diffusion through the AQP protein using two
independent methods, we demonstrate that there is an overestimation of the
diffusion, resulting in overestimated values of water permeability for TIP3P.
However, for TIP4P/2005 and OPC, which have bulk diffusion values close to the
experimental one, their calculated values of water permeability are comparable
to the experimental value. This implies that TIP4P/2005 and OPC are capable of
producing reliable estimates of water diffusion through the AQP1 channel.

##  IV Conclusion

In this work, water diffusion and permeability have been studied using
molecular dynamics for three systems, consisting of AQP1 embedded in a POPC
lipid membrane and three water models, TIP3P, TIP4P/2005 and OPC. TIP3P is the
standard water model used for CHARMM36m, whereas TIP4P/2005 and OPC are water
models not frequently used in biomolecular simulations but known to better
reproduce bulk water properties [31, 40].

Recently a number of studies have addressed the relevance of water models in
e.g. protein aggregation and host-guest recognition systems [37, 59]. Our main
purpose was to compare among the behaviour of different water models and shed
light on AQP1 water permeability. This enabled us to obtain reliable results
for the molecular mechanism of water diffusion as well as both diffusion and
permeability values, near the membrane and through the AQP1. First of all, the
well-known molecular mechanism of water permeation through AQP1 has been
tested for the three models, giving similar and excellent results. All of them
show a specific orientation of water molecules inside the AQP1 as demonstrated
in Ref. [11]. Hence, using any of these water models, the water-protein
interaction is reliable and AQP1 keeps its high-selectivity permeation of
water across the membrane.

The obvious difference among water models can be detected in the transport
properties. Diffusion inside the protein has been computed using two methods,
one based on the mean square displacement considering the geometry of the
channel and another one from the water permeability, giving consistent results
for the three models. TIP3P has a diffusion constant twice that of the
experimental value, whereas TIP4P/2005 and OPC produce values remarkably close
to experiment. The overestimation of TIP3P is consistent in the whole process
of water permeation across the membrane. The most relevant region is inside
the membrane protein,  i.e., the AQP region, where the TIP3P diffusion is
higher than the experiments. On the other hand, OPC presents a lower diffusion
rate. When comparing to experiments, we conclude that TIP4P/2005 is the model
that best resembles water diffusivity inside the channel. It is worth noting
that the ratios of TIP3P/TIP4P/2005 and TIP3P/OPC increase when moving from
bulk to the interfacial region to the AQP region. Therefore, the difference
between the models seems to be more pronounced within the nanopore
environment.

Focusing on water permeability, we have computed the number of water molecules
able to cross the membrane through the AQP1, $q_{0}$. By means of this
parameter we calculated the diffusion permeability, which allowed us to obtain
the value of osmotic permeability. As in the diffusion case, TIP3P gives a
diffusion permeability between 6 and 8 times higher than TIP4P/2005 and OPC,
respectively.

The overestimation of TIP3P is observed also in the calculation of the osmotic
permeability (being almost 3 times higher than the experimental counterpart)
and in the water diffusion inside the AQP1 (being 5 times higher than the
experimental counterpart). However the calculated values obtained with the OPC
and TIP4P/2005 water models are in much better agreement with the experiments.
It is difficult to compute the osmotic permeability in experiments and there
are several values of $p_{f}$. However, in this work reliable values of water
diffusion and water permeability are reported, which demonstrate the clear
differences between the behaviour of TIP3P and the results of TIP4P/2005 and
OPC. We believe that this observation can be relevant not only for the AQP
channel but also for numerical studies of transport properties of water and
ions within nanopores and ion channels.

Simulating polarizable models require a considerable increase in computational
effort, and to date their application has generally been limited to somewhat
simpler membrane [60] and channel systems [36, 61]. We therefore have decided
to first understand the basic features of water transport through AQP1, using
only non-polarizable models. However, in the future it will be important to
consider systematically polarizable models for water, proteins and lipids [62,
61].

This study has focused on the dynamical properties of the water model in AQP
simulations. However, it should be noted that the protein and lipid parameters
have been optimized for use with e.g. TIP3P. The question therefore remains
whether in the future one should re-optimize these parameters for improved
water models or switch to polarizable models. The latter are being extended to
lipids as well as to water and proteins and are now feasible, if
computationally expensive, for membrane channel protein systems.[63] An
interim compromise may be the use of charge scaling (via the electronic
continuum correction approach) to mimic electronic polarization, which has
been used with some success for biological membranes.[64]

## Supplementary material

See the supplementary material is available free of charge at

  * •

Dipole moment for TIP4P/2005. (Figure S1)

## DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## acknowledgements

The work has been performed under the Project HPC-EUROPA3
(INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action
under the H2020 Programme; in particular, the author gratefully acknowledges
the support of M.S.P. Sansom and the computer resources and technical support
provided by ARCHER. M.A.G. thanks for the support by Ayuda Juan de la Cierva-
Incorporación (IJCI-2016-27497) from Ministerio de Ciencia, innovación y
universidades (Spain). M.S.P. Sansom acknowledges the following fundings:
BBSRC BB/N000145/1 EPSRC EP/R004722/1 and EP/V010948/1, and C. Valeriani
acknowledges fundings from MINECO PID2019-105343GB-I00.

## References

  * Beitz [2009] E. Beitz, ed.,  _Aquaporins_, Handbook of Experimental Pharmacology, Vol. 190 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2009).
  * Yang [2017] B. Yang, ed.,  _Aquaporins_, Advances in Experimental Medicine and Biology, Vol. 969 (Springer Netherlands, Dordrecht, 2017).
  * Castle [2005] N. A. Castle, “Aquaporins as targets for drug discovery,” Drug Discov. Today 10, 485–493 (2005).
  * Day  _et al._ [2014] R. E. Day, P. Kitchen, D. S. Owen, C. Bland, L. Marshall, A. C. Conner, R. M. Bill, and M. T. Conner, “Human aquaporins: Regulators of transcellular water flow,” Biochim. Biophys. Acta - Gen. Subj. 1840, 1492–1506 (2014).
  * Kruse, Uehlein, and Kaldenhoff [2006] E. Kruse, N. Uehlein, and R. Kaldenhoff, “The aquaporins.” Genome Biol. 7, 206 (2006).
  * Agre and Kozono [2003] P. Agre and D. Kozono, “Aquaporin water channels: Molecular mechanisms for human diseases,” FEBS Lett. 555, 72–78 (2003).
  * Preston  _et al._ [1992] G. M. Preston, T. P. Carroll, W. B. Guggino, and P. Agre, “Appearance of water channels in Xenopus oocytes expressing red cell CHIP28 protein.” Science 256, 385–7 (1992).
  * Walz  _et al._ [1997] T. Walz, T. Hirai, K. Murata, J. B. Heymann, K. Mitsuoka, Y. Fujiyoshi, B. Smith, P. Agre, and A. Engel, “The Three-Dimensional Structure of Aquaporin-1,” Nature 387, 624–627 (1997).
  * Chen, Wu, and Voth [2006] H. Chen, Y. Wu, and G. A. Voth, “Origins of Proton Transport Behavior from Selectivity Domain Mutations of the Aquaporin-1 Channel,” Biophys. J. 90, L73–L75 (2006).
  * Kosinska Eriksson  _et al._ [2013] U. Kosinska Eriksson, G. Fischer, R. Friemann, G. Enkavi, E. Tajkhorshid, and R. Neutze, “Subangstrom Resolution X-Ray Structure Details Aquaporin-Water Interactions,” Science 340, 1346–1349 (2013).
  * Tajkhorshid  _et al._ [2002] E. Tajkhorshid, P. Nollert, M. Ø. Jensen, L. J. W. Miercke, J. O’Connell, R. M. Stroud, and K. Schulten, “Control of the Selectivity of the Aquaporin Water Channel Family by Global Orientational Tuning,” Science 296, 525–530 (2002).
  * de Groot, Engel, and Grubmüller [2001] B. L. de Groot, A. Engel, and H. Grubmüller, “A refined structure of human aquaporin-1,” FEBS Lett. 504, 206–211 (2001).
  * Law and Sansom [2004] R. J. Law and M. S. P. Sansom, “Homology modelling and molecular dynamics simulations: comparative studies of human aquaporin-1,” Eur. Biophys. J. 33, 477–489 (2004).
  * Ariz-Extreme and Hub [2017] I. Ariz-Extreme and J. S. Hub, “Potential of Mean Force Calculations of Solute Permeation Across UT-B and AQP1: A Comparison between Molecular Dynamics and 3D-RISM,” J. Phys. Chem. B 121, 1506–1519 (2017).
  * Decker  _et al._ [2017] K. Decker, M. Page, A. Boyd, I. MacAllister, M. Ginsberg, and A. Aksimentiev, “Selective Permeability of Truncated Aquaporin 1 in Silico,” ACS Biomater. Sci. Eng. 3, 342–348 (2017).
  * Gravelle  _et al._ [2013] S. Gravelle, L. Joly, F. Detcheverry, C. Ybert, C. Cottin-Bizonne, and L. Bocquet, “Optimizing water permeability through the hourglass shape of aquaporins,” Proc. Natl. Acad. Sci. 110, 16367–16372 (2013).
  * Aponte-Santamaría  _et al._ [2017] C. Aponte-Santamaría, G. Fischer, P. Båth, R. Neutze, and B. L. de Groot, “Temperature dependence of protein-water interactions in a gated yeast aquaporin,” Scientific Reports 7, 4016 (2017).
  * Berezhkovskii and Hummer [2002] A. Berezhkovskii and G. Hummer, “Single-file transport of water molecules through a carbon nanotube,” Phys. Rev. Lett. 89, 064503/1–064503/4 (2002).
  * de Groot and Grubmulle [2005] B. L. de Groot and H. Grubmulle, “The dynamics and energetics of water permeation and protonexclusion in aquaporins,” Curr. Opin. Struct. Biol. 15, 176–183 (2005).
  * Zhu, Tajkhorshid, and Schulten [2004] F. Q. Zhu, E. Tajkhorshid, and K. Schulten, “Theory and simulation of water permeation in aquaporin-1,” Biophys. J. 86, 50–57 (2004).
  * de Groot  _et al._ [2001] B. L. de Groot, H. Grubmüller, W. B. Guggino, P. Agre, J. O’Connell, R. M. Stroud, K. Schulten, and H. Grubmuller, “Water permeation across biological membranes: Mechanism and dynamics of aquaporin-1 and GlpF,” Science 294, 2353–2357 (2001).
  * de Groot  _et al._ [2003] B. L. de Groot, T. Frigato, V. Helms, and H. Grubmuller, “The mechanism of proton exclusion in the aquaporin-1 water channel,” J. Mol. Biol. 333, 279–293 (2003).
  * Hashido, Ikeguchi, and Kidera [2005] M. Hashido, M. Ikeguchi, and A. Kidera, “Comparative simulations of aquaporin family: AQP1, AQPZ, AQP0 and GlpF,” FEBS Lett. 579, 5549–5552 (2005).
  * Hub and de Groot [2006] J. S. Hub and B. L. de Groot, “Does CO2 Permeate through Aquaporin-1?” Biophys. J. 91, 842–848 (2006).
  * Hub  _et al._ [2010] J. S. Hub, C. Aponte-Santamaría, H. Grubmüller, and B. L. de Groot, “Voltage-Regulated Water Flux through Aquaporin Channels In Silico,” Biophys. J. 99, L97–L99 (2010), arXiv:0608246v3 [arXiv:physics] .
  * Mamonov  _et al._ [2007] A. B. Mamonov, R. D. Coalson, M. L. Zeidel, and J. C. Mathai, “Water and deuterium oxide permeability through aquaporin 1: MD predictions and experimental verification,” J. Gen. Physiol. 130, 111–116 (2007).
  * Jensen, Tajkhorshid, and Schulten [2003] M. Ø. Jensen, E. Tajkhorshid, and K. Schulten, “Electrostatic Tuning of Permeation and Selectivity in Aquaporin Water Channels,” Biophys. J. 85, 2884–2899 (2003).
  * Zhu, Tajkhorshid, and Schulten [2002] F. Q. Zhu, E. Tajkhorshid, and K. Schulten, “Pressure-induced water transport in membrane channels studied by molecular dynamics,” Biophys. J. 83, 154–160 (2002).
  * Lee  _et al._ [2016] S. C. Lee, S. Khalid, N. L. Pollock, T. J. Knowles, K. Edler, A. J. Rothnie, O. R. T. Thomas, and T. R. Dafforn, “Encapsulated membrane proteins: A simplified system for molecular simulation,” Biochim. Biophys. Acta-Biomembranes 1858, 2549–2557 (2016).
  * Ozu  _et al._ [2013] M. Ozu, H. A. Alvarez, A. N. McCarthy, J. R. Grigera, and O. Chara, “Molecular dynamics of water in the neighborhood of aquaporins,” Eur. Biophys. J. 42, 223–239 (2013).
  * Vega and Abascal [2011] C. Vega and J. L. F. Abascal, “Simulating water with rigid non-polarizable models: a general perspective,” Physical Chemistry Chemical Physics 13, 19663–19688 (2011).
  * Onufriev and Izadi [2018] A. V. Onufriev and S. Izadi, “Water models for biomolecular simulations,” Wiley Interdisciplinary Reviews: Computational Molecular Science 8, e1347 (2018).
  * Huang  _et al._ [2017] J. Huang, S. Rauscher, G. Nawrocki, T. Ran, M. Feig, B. L. de Groot, H. Grubmüller, and A. D. MacKerell, “CHARMM36m: an improved force field for folded and intrinsically disordered proteins,” Nat. Methods 14, 71–73 (2017).
  * Best and Mittal [2010] R. B. Best and J. Mittal, “Protein Simulations with an Optimized Water Model: Cooperative Helix Formation and Temperature-Induced Unfolded State Collapse,” J. Phys. Chem. B 114, 14916–14923 (2010).
  * Lynch, Rao, and Sansom [2020] C. I. Lynch, S. Rao, and M. S. P. Sansom, “Water in nanopores and biological channels: A molecular simulation perspective,” Chemical Reviews 120, 10298–10335 (2020).
  * Klesse, Tucker, and Sansom [2020] G. Klesse, S. J. Tucker, and M. S. P. Sansom, “Electric field induced wetting of a hydrophobic gate in a model nanopore based on the 5-ht3 receptor channel,” ACS Nano 14, 10480–10491 (2020).
  * Çınaroğlu and Biggin [2021] S. S. Çınaroğlu and P. C. Biggin, “Evaluating the performance of water models with host–guest force fields in binding enthalpy calculations for cucurbit[7]uril–guest systems,” J. Phys. Chem. B 125, 1558–1567 (2021).
  * Abraham  _et al._ [2015] M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, “GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers,” SoftwareX 1-2, 19–25 (2015).
  * Jorgensen [1981] W. L. Jorgensen, “Quantum and statistical mechanical studies of liquids. 10. Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water,” J. Am. Chem. Soc 103, 335–340 (1981).
  * Izadi, Anandakrishnan, and Onufriev [2014] S. Izadi, R. Anandakrishnan, and A. V. Onufriev, “Building Water Models: A Different Approach,” J. Phys. Chem. Lett. 5, 3863–3871 (2014).
  * Abascal and Vega [2005] J. L. F. Abascal and C. Vega, “A general purpose model for the condensed phases of water: TIP4P/2005,” J. Chem. Phys. 123, 234505 (2005).
  * Bussi, Donadio, and Parrinello [2007] G. Bussi, D. Donadio, and M. Parrinello, “Canonical sampling through velocity rescaling,” J. Chem. Phys. 126, 014101 (2007).
  * Parrinello and Rahman [1981] M. Parrinello and A. Rahman, “Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method,” J. Appl. Phys. 52, 7182–7190 (1981).
  * Essmann  _et al._ [1995] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, “A smooth particle mesh Ewald method,” J. Chem. Phys. 103, 8577–8593 (1995).
  * Róg, Murzyn, and Pasenkiewicz-Gierula [2002] T. Róg, K. Murzyn, and M. Pasenkiewicz-Gierula, “The dynamics of water at the phospholipid bilayer surface: a molecular dynamics simulation study,” Chem. Phys. Lett. 352, 323–327 (2002).
  * Bhide and Berkowitz [2005] S. Y. Bhide and M. L. Berkowitz, “Structure and dynamics of water at the interface with phospholipid bilayers,” J. Chem. Phys. 123, 224702 (2005).
  * Yamamoto  _et al._ [2013] E. Yamamoto, T. Akimoto, Y. Hirano, M. Yasui, and K. Yasuoka, “Supporting Information to “ Power-law trappings of water molecules on the lipid membrane surface induce the water retardation ”,” Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. 87, 1–6 (2013).
  * Yang, Calero, and Martí [2014] J. Yang, C. Calero, and J. Martí, “Diffusion and spectroscopy of water and lipids in fully hydrated dimyristoylphosphatidylcholine bilayer membranes,” J. Chem. Phys. 140, 104901 (2014).
  * Toppozini  _et al._ [2015] L. Toppozini, F. Roosen-Runge, R. I. Bewley, R. M. Dalgliesh, T. Perring, T. Seydel, H. R. Glyde, V. García Sakai, and M. C. Rheinstädter, “Anomalous and anisotropic nanoscale diffusion of hydration water molecules in fluid lipid membranes,” Soft Matter 11, 8354–8371 (2015).
  * Yamamoto  _et al._ [2015] E. Yamamoto, T. Akimoto, M. Yasui, and K. Yasuoka, “Origin of subdiffusion of water molecules on cell membrane surfaces,” Sci. Rep. 4, 4720 (2015).
  * Zaragoza  _et al._ [2019] A. Zaragoza, M. A. Gonzalez, L. Joly, I. López-Montero, M. A. Canales, A. L. Benavides, and C. Valeriani, “Molecular dynamics study of nanoconfined TIP4P/2005 water: how confinement and temperature affect diffusion and viscosity,” Phys. Chem. Chem. Phys. 21, 13653–13667 (2019).
  * Horner  _et al._ [2015] A. Horner, F. Zocher, J. Preiner, N. Ollinger, C. Siligan, S. A. Akimov, and P. Pohl, “The mobility of single-file water molecules is governed by the number of H-bonds they may form with channel-lining residues,” Sci. Adv. 1, E1400083 (2015).
  * Shahbabaei and Kim [2020] M. Shahbabaei and D. Kim, “Exploring fast water permeation through aquaporin-mimicking membranes,” Phys. Chem. Chem. Phys. 22, 1333–1348 (2020).
  * Finkelstein [1987] A. Finkelstein,  _Water movement through lipid bilayers, pores, and plasma membranes: theory and reality_ , Distinguished lecture series of the Society of General Physiologists (Wiley & Sons, New York, 1987) p. 228.
  * Mathai  _et al._ [1996] J. C. Mathai, S. Mori, B. L. Smith, G. M. Preston, N. Mohandas, M. Collins, P. C. M. van Zijl, M. L. Zeidel, and P. Agre, “Functional Analysis of Aquaporin-1 Deficient Red Cells,” J. Biol. Chem. 271, 1309–1313 (1996).
  * MacKerell  _et al._ [1998] A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, and M. Karplus, “All-atom empirical potential for molecular modeling and dynamics studies of proteins,” J. Phys. Chem. B 102, 3586–3616 (1998).
  * Martelli [2021] F. Martelli, “Topology and complexity of the hydrogen bond network in classical models of water,” J. Mol. Liq. , 115530 (2021).
  * Walz  _et al._ [1994] T. Walz, B. L. Smith, M. L. Zeidel, A. Engel, and P. Agre, “Biologically active two-dimensional crystals of aquaporin CHIP,” J. Biol. Chem. 269, 1583–1586 (1994).
  * Emperador, Crehuet, and Guàrdia [2021] A. Emperador, R. Crehuet, and E. Guàrdia, “Effect of the water model in simulations of protein–protein recognition and association,” Polymers 13, 176 (2021).
  * Chen  _et al._ [2021] P. Chen, I. Vorobyov, B. Roux, and T. W. Allen, “Molecular dynamics simulations based on polarizable models show that ion permeation interconverts between different mechanisms as a function of membrane thickness,” J. Phys. Chem. B 125, 1020–1035 (2021).
  * Ngo  _et al._ [2021] V. Ngo, H. Li, A. D. MacKerell, T. W. Allen, B. Roux, and S. Noskov, “Polarization effects in water-mediated selective cation transport across a narrow transmembrane channel,” J. Chem. Theory Comput. 17, 1726–1741 (2021), pMID: 33539082.
  * Ren and Ponder [2003] P. Ren and J. W. Ponder, “Polarizable atomic multipole water model for molecular mechanics simulation,” J. Phys. Chem. B 107, 5933–5947 (2003).
  * Jing  _et al._ [2021] Z. Jing, J. A. Rackers, L. R. Pratt, C. Liu, S. B. Rempe, and P. Ren, “Thermodynamics of ion binding and occupancy in potassium channels,” Chem. Sci. 12, 8920–8930 (2021).
  * Duboué-Dijon  _et al._ [2020] E. Duboué-Dijon, M. Javanainen, P. Delcroix, P. Jungwirth, and H. Martinez-Seara, “A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization,” J. Chem. Phys. 153, 050901 (2020).