Limitations of non-polarizable force fields in describing anion binding poses in non-polar synthetic hosts

Transmembrane anion transport by synthetic ionophores has received increasing interest not only because of its relevance for understanding endogenous anion transport, but also because of potential implications for therapeutic routes in disease states where chloride transport is impaired. Computational studies can shed light on the binding recognition process and can deepen our mechanistic understanding of them. However, the ability of molecular mechanics methods to properly capture solvation and binding properties of anions is known to be challenging. Consequently, polarizable models have been suggested to improve the accuracy of such calculations. In this study, we calculate binding free energies for different anions to the synthetic ionophore, biotin[6]uril hexamethyl ester in acetonitrile and to biotin[6]uril hexaacid in water by employing non-polarizable and polarizable force fields. Anion binding shows strong solvent dependency consistent with experimental studies. In water, the binding strengths are iodide > bromide > chloride, and reversed in acetonitrile. These trends are well captured by both classes of force fields. However, the free energy profiles obtained from potential of mean force calculations and preferred binding positions of anions depend on the treatment of electrostatics. Results from simulations using the AMOEBA force-field, which recapitulate the observed binding positions, suggest strong effects from multipoles dominate with a smaller contribution from polarization. The oxidation status of the macrocycle was also found to influence anion recognition in water. Overall, these results have implications for the understanding of anion host interactions not just in synthetic ionophores, but also in narrow cavities of biological ion channels.


SI Text 1 Force field Parametrization for Biotin Macrocycles
For non-polarizable, pairwise additive force fields, parameters for the methylated and unmethylated biotin macrocycles were generated with the second generation of General AMBER Force Field (GAFF2) 1 and the Parsley force field from the Open Force Field Initiative. 2 To assign force field parameters to biotin [6]uril hexamethyl ester, we used antechamber and tleap for the GAFF2 force field. 3 For the nonbonded interactions, van der Waals parameters are inherited from the AMBER force fields and partial charges are assigned according to the so-called Austin Model 1 bond charge corrections (AM1-BCC) model. Mulliken charges are obtained by a semi-empirical AM1 calculation, followed by a bond charge correction (BCC) scheme. 1 The Parsley force field relies on an AM1-BCC model for charge assignment as well 2 . Simulations with pairwise additive force fields (GAFF2 and Parsley) were carried out using the three sites transferable intermolecular potential (TIP3P) water model 4 . We used the parameters for halide ions from Li et al. and other parameters for non-bonded interactions were taken from the AMBER force field. 5,6 The AMOEBA force field (Atomic Multipole Optimized Energetics for Biomolecular Applications) incorporates a polarizable multipole framework up to quadrupole moments. 7 We used the AMOEBA03 parameters for water. 8 For the parameters for monovalent ions (halide anions, and potassium) and acetonitrile, we used the AMOEBA09 parameter set. 7 The parametrization of biotin [6]uril hexamethyl ester for the AMOEBA force field followed the POLTYPE protocol. 9 We used density functional theory calculations to derive electrostatic parameters. The defaults in the POLTYPE protocol are ab initio calculations which were also used by Laury et al. for the parametrization of the cucurbit [8]uril molecule. 10 The initial structure was optimized using the semi-empirical PM3 method and then optimized with B3LYP level of theory with a 6-311G(1d,1p) basis set and empirical dispersion gdBJ using the Gaussian 16 software package. The optimized structure was used as an input for a distributed multipole analysis using Stone's GDMA 11 that determined initial atomic multipole estimates such as the charge, the 3 components of the dipole vector and the quadrupole tensor. The Tinker program POLEDIT was used to choose local multipole frames and atomic polarizabilities and to define polarization groups. A polarization group in AMOEBA defines a group of atoms whose permanent multipoles do not polarize one another. In AMOEBA, the multipole moment vector is decomposed into a permanent contribution, independent of the environment and an induced component that can change according to the local electric field of other atoms. Within a polarization group, the permanent multipoles are excluded from contributing to the field at an atomic site, and hence only the other induced dipoles polarize the atom in question. 12 A single point calculation at B3LYP level of theory with cc-ptvz basis and empirical dispersion gd3BJ was used as a reference point for electrostatic potential fitting. A Cartesian grid of points, at which to calculate the electrostatic potential, was created with the tinker POTENTIAL program. 13 The quantum mechanical (QM) electrostatic potential was evaluated at each of these grid points with the Gaussian CUBEGEN program. We then evaluated the AMOEBA electrostatic potential at each grid point and fitted the dipole and quadrupole components to minimize the root mean square error between the QM and AMOEBA electrostatic potential with the tinker POTENTIAL program while keeping the partial charges (the monopoles) fixed. 13 Finally, we assigned parameters for bonds, angles, stretch-bends, outof-plane bends, torsions and van der Waals interactions by similarity to the AMOEBA09 parameter set with the tinker VALENCE program. 7,13 The AMOEBA09 parameter file contains atom classes for many functional groups such as saturated carbon and hydrogen, amine nitrogen, sulfides and carboxylic acids. 4

Ion Parameters
Within molecular mechanics force-fields, the 12-6 Lennard-Jones potential is the most commonly used potential to model van der Waals interactions. The repulsive term is proportional to !"# and the attractive term to !$ : (1) where refers to the zero of the potential and to the well depth of its minimum.
Each atom i is characterized by a parameter pair ( , ( . The pairwise interaction between atom i and j can be modeled with Lorentz-Berthelot combination rules for radius and energy: In simulations with Parsley 2 , which we selected as the non-polarizable force-field, we used the ion parameters from Li et al., which were parameterized to reproduce hydration free energies. 6 Fig. S1a shows the 12-6 Lennard-Jones potential for the four halide anions. With increasing ion size, both and increase. In simulations with ion parameters from the AMBER03 parameter set, the ions are not stably bound inside the cavity and this likely reflects the fact that the AMBER03 Lennard-Jones potential for halide anions is much shallower (Fig. S1c) compared to the other parameter sets considered here.
To test whether the difference in radius contributed to the lack of binding of the chloride anion to the centre of the cavity, we fitted the Lennard-Jones potential to the buffered 14-7 potential in Fig. S4. The Lennard-Jones-fit resembling the AMOEBA buffered 14-7 potential does not preserve the hydration free energy of chloride in water. The chloride with fitted LJ parameters is not stably bound in the cavities of recently commented on in a study by Smith et al. 15 in which they warn about overfitting the van der Waals radii of ions. The chloride parameters from Li et al. 6 relate the radius and the well depth via a socalled "noble-gas" curve which is a fit based on the VDW parameters of noble gas atoms to correlate the radius and the well depth in the LJ equation. Compared to additive force-fields, the AMOEBA force field 7, 17 uses different potentials for non-bonded interactions. Van der Waals interactions are described by a buffered 14-7 potential of the form Where ( ) denotes the reduced interatomic distance depending on the atomic polarisabilities ( . In the AMOEBA force field, all atoms have the same dimensionless width parameter = 0.39 of the smeared charge distribution.

SI Text 2. Effects of Refitting Ion Parameters
The chloride with fitted LJ potential was placed centrally in the ring, energy minimized, and carefully equilibrated. 1ns NVT, 1ns NPT, 1ns NPT with restraints of 1500 and 1000 and 500 kJmol/nm^2 restraints respectively to ensure that the ion is centrally bound in the cavity. To allow flexibility in the xyplane of the ring, we applied flat-bottom restraints similar in set-up as in the umbrella simulations. All restraints were lifted for production runs of 100 ns.
We set up three repeats with the set-up described above and in all three repeats, the chloride with modified VDW parameters left the cavity almost immediately. In the combined simulation time of 300 ns, we did not observe any re-binding event (Fig. S4b).
We also calculated the potential of mean force profile along the reaction coordinate for the chloride with modified VDW parameters and for fluoride, which has an even smaller VDW radius (Fig. S4c). Both ions cannot be stably bound by the biotin [6]uril host. Even though the interaction energies between host and ions with smaller VDW radius are favorable (Fig. S4d, e), the ions do not bind centrally to the cavity.