Hydration of biologically relevant tetramethylammonium cation by neutron scattering and molecular dynamics

Abstract

Neutron scattering and molecular dynamics studies were performed on a concentrated aqueous tetramethylammonium (TMA) chloride solution to gain insight into the hydration shell structure of TMA, which is relevant for understanding its behavior in biological contexts of, e.g., properties of phospholipid membrane headgroups or interactions between DNA and histones. Specifically, neutron diffraction with isotopic substitution experiments were performed on TMA and water hydrogens to extract the specific correlation between hydrogens in TMA (H_TMA) and hydrogens in water (H_W). Classical molecular dynamics simulations were performed to help interpret the experimental neutron scattering data. Comparison of the hydration structure and simulated neutron signals obtained with various force field flavors (e.g. overall charge, charge distribution, polarity of the CH bonds and geometry) allowed us to gain insight into how sensitive the TMA hydration structure is to such changes and how much the neutron signal can capture them. We show that certain aspects of the hydration, such as the correlation of the hydrogen on TMA to hydrogen on water, showed little dependence on the force field. In contrast, other correlations, such as the ion–ion interactions, showed more marked changes. Strikingly, the neutron scattering signal cannot discriminate between different hydration patterns. Finally, ab initio molecular dynamics was used to examine the three-dimensional hydration structure and thus to benchmark force field simulations. Overall, while neutron scattering has been previously successfully used to improve force fields, in the particular case of TMA we show that it has only limited value to fully determine the hydration structure, with other techniques such as ab initio MD being of a significant help.

---

1 Introduction

Tetramethylammonium (TMA) represents an important and ubiquitous motif in biological systems. It is found in the cellular membranes of almost all cells as the phosphorylcholine group in many phospholipids. Also, successive methylation of the amino acid lysine eventually results in a TMA functional group. This methylation is vital in histone–DNA binding and the epigenetic expression of DNA.¹ Tetraalkylammonium salts are also powerful denaturants² and are widely used as phase transfer catalysts³ and are also seen frequently in ionic liquids.⁴

TMA is one of the simplest and most spherical representatives of the so-called “hydrophobic ions”, as compared to more complex variants, such as tetraphenylphosphonium or tetrabutylammonium. TMA has short hydrophobic chains, so an intriguing question arises. How much does its hydration properties differ from those of a hypothetical perfectly spherical large ion?^5–7 Small spherical alkali cations have a fairly generic hydration structure, with the first hydration shell water oxygen atom pointing towards the cation and the hydrogen atoms away. This orientation is the strongest for lithium, getting weaker upon moving down the periodic table. How does the solvation shell structure change by the time we get to a large cation like TMA? Also, given the hydrophobic character of the methyl groups, does hydrophobicity-driven cation–cation aggregation take place in aqueous TMA solutions?⁸

Neutron scattering experiments have provided insight into TMA hydration structure in TMA chloride or bromide solutions.^5,6,9–11 In aqueous solutions, the neutron scattering signal is dominated by water–water correlations. The structural correlations of interest, such as those between TMA and water nuclei, are thus hidden among the dominant water–water correlations. A way to deal with this problem is to employ neutron diffraction with isotopic substitution (NDIS).^12–15 This method takes advantage of the fact that different isotopes have different neutron scattering properties and relies on the assumption that the mass of a nucleus does not affect the structure of the solution. This assumption of neglecting nuclear quantum effects has proven fairly robust even for materials with the largest isotope effects. Even H₂O and D₂O vary in number density by less than half a percent despite substituting two-thirds of the nuclei in the system. Still, hydrogen and deuterium are excellent nuclei to use in NDIS. Both are easily experimentally available and have among the largest neutron scattering contrasts for any element. The technique requires two identical solutions to be prepared, which differ only in isotopic concentration of one nucleus. As the structure of these solutions is assumed to be identical, the subtraction of their diffraction patterns cancels out for all components unrelated to the substituted nucleus. With its 12 identical hydrogen nuclei and the large contrast between H and D isotopes, the TMA ion is an ideal candidate to be studied by neutron scattering.^5,9,10 However, even if NDIS makes isolating a few structural correlations of interest possible, not all pairwise distributions are equally easy to understand. The intuitive way of examining TMA hydration is to look at the spatial distribution of water oxygen atoms around the central TMA nitrogen atom. Unfortunately, NDIS cannot provide this specific structure factor due to lack of suitable oxygen isotopes. The relatively easily experimentally accessible H_TMA–H_W structure factor corresponds neither to the center of the TMA nor that of the water molecule, which makes intuitive interpretations of the data significantly more complicated.

Despite these limitations, early NDIS experiments^5,9,10,16 were interpreted as evidence of apolar hydration around TMA, with an edge-on orientation of water molecules. More recently, Monte–Carlo-based empirical potential structure refinement (ESPR) simulations were used to assist the interpretation of the neutron scattering signal, confirming the apolar character of TMA hydration but suggesting that water arranges tetrahedrally around TMA.¹¹ In the past years, molecular dynamics simulations have proven to be highly valuable to help further the interpretation of neutron scattering signals,^17,18 which is also the strategy adopted in the present work. In particular, we performed NDIS experiments on a 2 m TMACl solution, allowing us to extract a single structural correlation between hydrogens of the cation H_TMA and those of water H_W. We employed DFT-based ab initio molecular dynamics (AIMD) and force field molecular dynamics (FFMD) simulations to assist in interpreting the experimental signal. Simulations with various force fields were performed to test the sensitivity of the measured signal to various changes in the hydration structure, which turned out to be the key to understanding the exact amount of information contained in the neutron scattering signal and avoiding over-interpretations. In addition, we obtained new insights into optimizing force fields to better represent the hydration shells around solutes.

2 Methods

2.1 Neutron scattering measurements

NDIS measurements were performed 23 °C using the D4C diffractometer at the nuclear reactor at the Institut Laue–Langevin in Grenoble, France.¹⁹ All the samples were loaded into the same cylindrical null scattering titanium–zirconium cell, loaded in an identical geometry in the diffractometer. The sample diameter was 5.0 mm, the wall thickness 0.75 mm and the beam height 24 mm. The neutron wavelength was 0.4985 Å. Four chemically identical solutions of 2 m TMACl in water were prepared, which differed only in H/D substitution on the TMA (h12-TMA and d12-TMA) and H/D substitution on water. The four diffraction patterns (Fig. 1a and b) were recorded for about 2 h for each D₂O solution and for 4 h for each H₂O solution.²⁰ The results were then corrected for multiple scattering and absorption and normalized against a standard vanadium scatterer.²¹

Fig. 1 (a) Total diffraction patterns for H₂O solutions of d12-TMACl (black) and h12-TMACl (grey). (b) Total diffraction patterns for D₂O solutions of d12-TMACl (black) and h12-TMACl (grey). (c) First order differences (red line, obtained by the difference of the two diffraction patterns shown in a) and (blue line, obtained by the difference of the two diffraction patterns shown in b). (d) Second order difference 140.8·(S_HTMAHW(Q) − 1), obtained through the difference of the two first-order differences shown in (c).

Taking the difference between the diffraction patterns associated with solutions that differ only by the H/D substitution on TMA (both in D₂O and H₂O solutions) yields the first-order differences and (Fig. 1c), that report on the correlation between non-exchangeable H on TMA and every other atom (X) in the system. They are respectively defined as (in units of mbarns):

These prefactors were calculated from the atomic concentration and neutron scattering lengths of the various elements in the system by standard literature methods.²²

The difference between eqn (1) and (2) yields the second order difference ΔΔS_Hnon(Q) (Fig. 1d), which reports on the single correlation between the TMA non-exchangeable H atoms and the water H atoms.

This function provides a useful internal consistency check for the accuracy of the solutions and the multiple scattering and absorption corrections performed on the data. Due to the large inelastic scattering of ¹H and the Placzek effect,²³ hydrogen-containing samples always present a dominant background. The higher the atomic concentration of ¹H, the larger the effect, primarily visible in light water samples such as in Fig. 1a. The effect is greatly diminished for the heavy water samples, as shown in Fig. 1b. The amount of inelastic scattering is largely determined by the number of ¹H nuclei per unit volume, so the first-order differences (Fig. 1c) should have exactly the same Placzek background, which should vanish completely in the second order difference (Fig. 1d). The largely constant value of the second-order difference in Fig. 1d indicates the absence of a visible background and proves that the four solutions were prepared with the same chemical composition.

2.2 Force field molecular dynamics (FFMD) simulations

We performed classical FFMD simulations of seven systems comprising 50 TMA cations or, for comparison, 50 neutral neopentane molecules matching the experimental concentration of 2 m. Force field details for each system are available in Table 1. Systems containing TMA were neutralized with chloride counterions of the CHARMM Cl_o type²⁴ (−1 charge) or CL_2s type 13 (−0.75 scaled charge). Each system contained 1388 TIP3P water molecules with Lennard-Jones parameters also on the hydrogens.²⁵ All systems were assembled using GROMACS 2021.2 and 2021.5 tools.²⁶

Table 1 TMA and neopentane simulation models used in this study. In parenthesis, the used CHARMM36 atoms types³¹ are listed

Model	Geom.	Atom (atom type) partial charge			Overall charge
		N (NTL)	C (CTL5)	H (HL)
CHARMM	TMA	−0.60	−0.35	0.25	+1.00
prosECCo	TMA	−0.61	−0.35	0.23	+0.75
Low CH dipole	TMA	−0.05	−0.10	0.10	+0.75
Center-N	TMA	+0.75	0.00	0.00	+0.75
Surface-H	TMA	0.00	0.00	0.0625	+0.75

---

	Center-bead (σ = 0.550 nm, ε = 0.83680 kJ mol^−1)
Center-bead	Sphere	+0.75	+0.75

---

		C (CT)	C (CT3)	H (HA3)
Neopentane	TMA	0.00	−0.30	0.10	0.0

---

For the TMA moiety, several parametrizations used in the head group of lipids are available in the literature, which differ mainly in their partial charges on the three atom types of the TMA group. First, we employed the CHARMM36 parameters with the full +1 charge (denoted CHARMM). Next, we used a scaled variant, prosECCo, where charges were scaled by 0.75. This physically based Electronic Continuum Correction (ECC) approach accounts for electronic polarization in a mean-field way by scaling charges.²⁷ ECC has been shown to improve significantly the description of ion pairing behavior in solution.^28,29 Based on this scaled charge model with an overall charge of +0.75, we designed three additional variants with different charge distributions but the same overall +0.75 charge: (1) low CH dipole, where partial charges are reduced; (2) central-N, where the whole charge is concentrated on the nitrogen; (3) surface-H, where all charge is redistributed on the hydrogens (see Table 1). All these TMA models have the same Lennard-Jones parameters, as prescribed in CHARMM36.

For comparison, we also performed simulations with the neutral neopentane with the same geometry as TMA using the CHARMM force field. Note that this compound is highly hydrophobic and aggregates quickly when mixed in water. Finally, we modeled a single charged sphere similar to the coarse-grained model of TMA, center-bead. In this coarse-grained-like simulations, water still has an atomistic resolution while the entire tetramethylammonium group is replaced by a single “extended atom”. The size of this so-called bead is set to match the one of TMA in the first peak of its radial distribution with the surrounding water.³⁰

All systems were simulated with GROMACS 2021,²⁶ using simulation parameters as provided by CHARMM-GUI.³² The LINCS and Settle algorithms were used to constrain the geometry of TMA hydrogens³³ and water molecules, respectively, allowing the use of 2 fs time step in the simulations.³⁴ The coordinates were saved every 10 ps for each 2 μs simulation. For neopentane that aggregates due to its hydrophobic nature, we have simulated 200 ns biased simulations placing a lower wall between neopentanes at 0.75 nm distance (KAPPA = 20000.0 EXP = 2 EPS = 1 OFFSET = 0) using the lowest collective variable³⁵ to gather enough water orientation statistics around neopentane by excluding neopentane–neopentane contacts. Simulations were performed in the isothermal–isobaric (NpT) ensemble. Temperature was controlled using a Nosé–Hoover thermostat³⁶ with a time constant of 1 ps and a reference value of 310 K, and a constant pressure of 1 bar was maintained by an isotropically coupled Parrinello–Rahman³⁷ barostat with a time constant of 5 ps. van der Waals interactions were treated using a cutoff of 1.2 nm with a force-switch at 1.0 nm using Verlet cutoff scheme³⁸ for neighbors. Long-range Coulomb interactions were accounted for using particle mesh Ewald (PME) with a 1.2 nm cutoff as implemented in GROMACS.²⁶ All simulation data can be found at https://zenodo.org (DOI: https://doi.org/10.5281/zenodo.10406618).

2.3 Ab initio molecular dynamics (AIMD) simulations

As a benchmark for the TMA solvation structure, we complemented the neutron diffraction experiments and the force field-based simulations with Born–Oppenheimer ab initio molecular dynamics simulations (AIMD) of a single TMA cation with 64 water molecules under periodic boundary conditions. The large computational cost of AIMD simulations precludes using larger systems, such as those in FFMD. The present system was not neutralized by any counterion. Force field molecular dynamics was used to preequilibrate the system and to prepare the initial configurations for the subsequent AIMD simulation, as follows, with a constant pressure simulation used to estimate the average size of the cubic simulation cell of 12.552 Å. A constant volume simulation was then used to prepare 10 initial configurations separated by 2 ns and equilibrated for 1 ns using FFMD. From these 10 equilibrated structures, AIMD simulations were performed using the generalized gradient approximation revPBE DFT functional^39–41 with the D3 dispersion correction.^42–44 Core electrons were replaced by GTH pseudopotentials,^45,46 and the triple-ζ basis set TZV2P with polarization functions was used for valence electrons.⁴⁷ A cutoff of 400 Ry was used for the auxiliary plane wave basis set in the GPW method.⁴⁸ The system was first equilibrated using AIMD Langevin dynamics for at least 16 ps with a damping constant of γ = 0.02 ps⁻¹.⁴⁹ During the production simulation, the temperature was set to 300 K via a global CSVR thermostat with a time constant of 1 ps.⁵⁰ To enhance sampling, ten 50 ps parallel AIMD simulations totaling 500 ps of trajectory were used for further structural analysis. All AIMD simulations were performed at constant volume, using the CP2K program package, version 7.1.^51–53 All simulation data can be found at https://zenodo.org (DOI: https://doi.org/10.5281/zenodo.10406618).

2.4 Density maps

Data from FFMD and AIMD simulations were processed with the help of an in-house developed software for unbiased alignment and density analysis. The analysis included atoms found within a 10 Å of the nitrogen atom in TMA, covering all TMA molecules in the system and all simulation frames. The TMA “neighborhoods” were aligned in two steps. First, a representative conformation was found for each simulation, and then the neighborhoods were aligned to that conformation using the positions of carbon atoms. In the case of center-bead approximation of TMA, the neighboring oxygen/chlorine atoms were used for the alignment instead. The alignment was performed similarly as in ref. 54, that is, using a permutation-based unconstrained alignment, such that the reference atoms used for the alignment (TMA carbons or oxygens and chlorines within its first solvation shell) do not need to be labeled and sorted. This way, the thermal noise is uniformly distributed in the density maps of all systems.

3 Results and discussion

Neutron scattering patterns were obtained for four identical TMACl solutions that differed only by their isotopic composition on the nonexchangeable H of TMA and the exchangeable water H (see Methods section). While the diffraction patterns of the solution are dominated by the signal coming from the water signal, taking differences between pairs of solutions leads to the cancellation of the signal part that does not depend on the isotopic composition. Hence, first-order differences (see Methods) report the TMA(H/D) correlation to every other atom in the system. While rich in information, this signal still contains all the intramolecular correlation peaks, which hides the information about hydration (see ESI†). Hence, we proceeded to obtain the double difference signal ΔΔS_Hnon(Q) (see Methods and Fig. 2), which reports on the single correlation between H_TMA and H_W. It thus directly probes the TMA hydration structure and is much easier to interpret than the total diffraction patterns. The obtained ΔΔS_Hnon(Q) exhibits neatly resolved features below 10 Å which characterize TMA hydration. Interpretation of the signal tends to be more intuitive in direct space. Hence, we computed the inverse Fourier-transform ΔΔG_Hnon(r) = [scr F, script letter F]⁻¹[ΔΔS_Hnon(Q)], which is directly related to the single pair-correlation function g_HTMAHW:

Fig. 2 Upper the function ΔΔS_Hnon(Q) (grey) and lower the direct Fourier transform of this function ΔΔG_Hnon(r) (grey). Shown in red is the function ΔΔS_Hnon(Q) used for the rest of this paper, terminating the data using a windows function up to ≃10 Å⁻¹. The lower red function is the real space version of the upper red function.

ΔΔG_Hnon(r) shown in Fig. 2 presents a characteristic shoulder around 3–4 Å, and a peak at 6 Å. While direct molecular interpretation of these features is not straightforward, the same pair correlation function can be easily computed from molecular dynamics simulations with two objectives in mind. First, we want to assist in interpreting the neutron diffraction signal and obtain a molecular-level picture of the TMA hydration. Second, we want to investigate how sensitive TMA hydration is to variations in the intermolecular interactions, and whether the resulting patterns in hydration structure would be distinguishable in the neutron signal. To this aim, we performed FFMD simulations of a TMACl solution, at the same 2 m concentration as used in the neutron scattering experiment, using different force field variants (see Methods, Table 1), characterizing in each case the hydration structure and the associated neutron signal.

3.1 What is the effect of overall charge on TMA solvation?

TMA is a large cation with a relatively low charge density. Hence, a question arises on how different its hydration structure is from that of neutral solutes, such as neopentane (beyond the absence of counterions and possible aggregation, as discussed below). In addition, within standard force fields, TMA, as intuitively expected, is assigned a global charge of +1. However, electronic polarization, which further screens interactions between ions, is missing in FFMD simulations using non-polarizable force fields, which may lead to artefacts such as excessive ion pairing.^27,28 A mean-field strategy to implicitly account for electronic polarization in FFMD simulations is the electronic continuum correction (ECC) approach, which is mathematically equivalent to scaling partial charges by a factor , where ε_∞ is the high-frequency dielectric constant of water. Hence, we designed an ECC version of the TMA force field (denoted as prosECCo, see Table 1), with an overall charge of +0.75.^28,55 Here, we compared the structure of the solution simulated using this scaled-charge force field with that obtained using the standard full charge CHARMM force field, as well as to that of a solution containing electrically neutral neopentane molecules for comparison. Note that neopentane is insoluble and starts precipitating fast in the simulations. Consequently, we used a biased simulation to avoid their aggregation.

For each simulation, we obtained radial distribution functions from the central nitrogen atom to the surrounding O, H_W, and Cl⁻, as well as density maps of chloride, water oxygen, and water hydrogen around TMA (see Fig. 3). Fig. 3 clearly shows that changing the charge from +1 to +0.75 has only a minor effect on the hydration structure, resulting in nearly identical radial distribution functions. If we view the TMA ion as a tetrahedron with four faces, six edges, and four corners, the water oxygen atoms (as well as the chloride counterion) are located more at the faces, less at the edges, and not at all at the corners. The amount of TMA⁺–Cl⁻ ion pairing is also very similar for the two systems, exhibiting only a small excess of solvent-shared ion pairs in the full charge force field compared to that with scaled charges (see second peak at 6–8 Å in g_NCl(r)). This is contrary to what was previously observed for small monovalent cations such as Li⁺,⁵⁶ and with divalent ions such as Ca²⁺,⁵⁷ where charge scaling changed qualitatively the number of ion pairs. In contrast, due to the low charge density of TMA, only minor polarization effects are observed for ion hydration and pairing.

Fig. 3 Comparison of charge effects in first solvation shell of TMA geometries. The first column shows charge +1.0 TMA (CHARMM), the second column +0.75 TMA (prosECCo), and the third column (neopentane). The first row shows chloride ions and the lower oxygen and hydrogen atoms. In all density maps, the density of chloride ion (green) is six times bulk density, oxygen (red) is three times bulk density, and hydrogen (white) is two times bulk density. The lower row shows the RDFs from the central atom for each model with oxygen (red), hydrogen (blue), and chloride (green). The fourth column compares the N–Cl⁻ and N/C–O RDFs for the three models.

If we now compare TMA hydration with that of neutral neopentane, the main difference is due to a size effect with the hydration layer shifted slightly further away from the central atom. The orientation of the water molecules is also different – in neopentane, the water OH bonds point very slightly towards the central atom; in contrast, in TMA, the water dipole orients towards the ion, and the OH points slightly outwards. Hence, while much less strongly oriented than around smaller ions such as lithium, the hydration of TMA still appears in these simulations significantly different from that of a hydrophobic solute. Note that the observations discussed are robust, despite neopentane statistics being worse due to less sampling due to their strong clustering propensity.

3.2 Effect of charge distribution on TMA solvation

For a given overall charge, we investigated how the charge distribution within the TMA impacted its hydration structure and its propensity to form TMA–Cl ion pairs. For water as strong hydrogen bonding moiety, the partial charge on the central oxygen is about −0.8 with +0.4 on the hydrogen, while for TMA, the partial charges in most force fields are about −0.3 on the carbon and +0.2 on the hydrogen. Starting from the scaled charge prosECCo model with an overall TMA charge of +0.75, we compared four different charge distributions – the prosECCo force field (with results presented above), a low CH dipole force field variant where the charges on the methyl C and H atoms are respectively −0.1 and +0.1, a center-N force field where all the charge is placed on the central nitrogen, and a surface-H force field where the charge is equally distributed over the surface hydrogens of the TMA only. Note that while the charge distribution is rather different between the center-N and surface-H models, they both result in a very low polarity of the C–H bonds and zero or very small charge on hydrogens. In each of the four cases, we characterized the hydration structure and ion pairing via the radial distribution functions g_NH(r), g_NO(r) and g_NCl(r) and the density maps of H_W, O and Cl around TMA (Fig. 4).

Fig. 4 For all density maps, the density of chloride ion (green) is 6× bulk density, oxygen (red) is 3× bulk density and hydrogen (white) is 2× bulk density. All the TMA cations in this figure have a +0.75 charge. Left column, prosECCo force field, center left, the polarity of the CH bond is comparable to that of neopentane, with the remaining charge on the central nitrogen. Center right, all of the charge of the ion is on the central nitrogen, and far right, all the charge is spread evenly on the hydrogen atoms. Lower are shown the RDFs from the central atom of the TMA to oxygen (red), hydrogen (blue), and chloride (green). Right are shown the same RDFs but grouped for each four FF. Upper for the chloride ion and lower for the oxygen atom.

Our calculations show that changes in the charge distribution within TMA (see Fig. 4) have a much larger effect on hydration structure than the reduction of the overall charge of the ion from +1.0 to +0.75 (see Fig. 3). Namely, lowering the charge on the surface hydrogens qualitatively changes the hydration pattern. The density maps (Fig. 4) clearly show that the water oxygens (as well as the chloride counterions, which follow the same trends) then move to the center of the faces of the TMA tetrahedron (i.e., away from the H atoms), forming bridges over the tetrahedron edges. This is manifested in the radial distribution functions g_NH(r), g_NO(r) and g_NCl(r) as a subtle increase in the bimodality of the first peak for g_NO(r) and g_NCl(r) in center-N. At the same time, changing the charge distribution does not significantly change the orientation of the hydration water OH bonds (see Fig. 4).

3.3 Effect of molecular geometry on TMA solvation

Coarse-grained models of TMA reduce the molecular geometry to a single spherical bead. Is this loss of molecular structure important or, in other words, how much does TMA behave as a simple charged sphere? To examine this issue we performed a set of three simulations with the all-atom prosECCo force field, the center-N force field variant where the charge is localized on the central N atom, and the coarse-grained force field where all the charge is at the center of a single spherical bead. Again, the density maps for O, H_W and Cl were calculated around TMA, as well as the radial distribution functions g_NH(r), g_NO(r) and g_NCl(r), see Fig. 5.

Fig. 5 For all density maps, density of chloride ion (green) is 6× bulk density, oxygen (red) is 3× bulk density, and hydrogen (white) is 2× bulk density. All the ions in this figure have +0.75 charge. Left column, prosECCo force field, center all of the charge of the ion is on the central nitrogen, and right, as center but all the tetrahedral structure of the ion has been replaced by a single VDW sphere. Lower are shown the RDFs from the central atom of the TMA to oxygen (red), hydrogen (blue), and chloride (green). Right are shown the same RDFs but group for each three FFs. Upper for the chloride ion and lower for the oxygen atom.

Interestingly, the all-atom prosECCo model is more similar in terms of the NO and N–H_W radial distribution functions to the coarse-grained charged sphere than to the center-N model (Fig. 5). Notably, the radial distribution functions of the center-N model are much more bimodal than the two others due to the combination of structural arrangements around the TMA tetrahedral features of faces, edges, and corners. This means that much of the TMA hydration structure reflects the constraints imposed by the charged sphere of a given size on the water H-bond network. At the same time, the atomistic TMA model is strikingly different from a simple charged sphere when looking at its interaction with the chloride counterion. In the former, the counterions adopt a tetrahedral geometry at the center of the faces similar to the oxygen atoms, while the latter tends to form linear Cl-TMA-Cl structures (Fig. 5).

3.4 Comparison of neutron scattering data to MD simulations

The double difference signal, ΔΔG_Hnon(r), obtained from neutron scattering experiments after Fourier transform is composed of a single radial distribution function of g_HTMAHW(r), which can be directly compared to the same radial distribution function computed from FFMD simulations with different force fields (Fig. 6). All investigated force fields (CHARMM: full charge, prosECCo: scaled, and center-N: +0.75 on N) capture the location of the main peak at about 6 Å, even if the CHARMM force field seems to provide a somewhat worse fit than the two others. Small differences are visible in the low- r range, but the comparison does not allow us to decide on the best force field. The steep rise of the experimental signal is slightly shifted for all force fields, and the shoulders around 4–5 Å while all slightly different, are never the same as in the experimental signal.

Fig. 6 The experimental function ΔΔG_Hnon(r) is shown in red, with the same function calculated from FFMD simulations in black for the CHARMM (upper), prosECCo (middle) and center-N force fields (lower).

As neutron scattering data cannot differentiate between the different 3D hydration arrangements found by FFMD, we performed DFT-based AIMD simulations and used them as a benchmark. Since the obtained FFMD density maps differ qualitatively from each other, we aimed to employ the AIMD results to determine which force field is more accurate. The very high computing cost of such simulations (see Methods for more details) necessarily limits both the size of the simulated system (a single TMA ion in a small box of 64 water molecules) and the length of the simulations (500 ps). Due to the limited box size, the calculation of radial distribution functions is thus limited to the small r range (Fig. S2 in ESI†). We also obtained 3D density maps for water oxygen atoms around TMA. Despite limited statistics, these plots tend to show that the highest density of water oxygen atoms is located at the center of the faces of the TMA tetrahedron, with bridges across the sides. The hydration density maps from AIMD (see Fig. 7) display both similarities and differences with respect to the FFMD density maps. A feature common to all density maps is that the oxygen clouds found over the faces and edges are closer to TMA than those of the H_W clouds. There is, however, a difference in how these clouds are arranged for the prosECCo and low CH dipole force fields. In each case, the H_W clouds are similar in shape but vary such that the low CH dipole force field has a greater tendency to spill over the edges of the TMA tetrahedron. For the low CH dipole force field the O clouds match the orientation of the H_W cloud, while for the prosECCo force field, an opposite pattern is observed with the triangles of the H_W and O cloud being anticorrelated (Fig. 5 and 7). Hydration of the AIMD TMA is more similar to that of low CH dipole than that of prosECCo, where prosECCo and CHARMM models have similar hydrations. Interestingly, at lower atomic densities, AIMD density maps show overlapping H_W and O_w density clouds at the corners of TMA, which is not replicated in any force field-based FFMD simulation. Taken all together, these results suggest that the orientation of water molecules in the TMA hydration shell is better captured by low CH dipole than by the other force field variants tested here.

Fig. 7 O_W and H_W density maps around TMA in various simulations. Upper O_W density is shown in red at 3.3× bulk density (approximated for AIMD simulations) and lower 2.8× bulk density for O_W (red) and 2.1× bulk density for H_W (white). The orientation of the O_W around TMA shows an inverse relationship between prosECCo and low CH dipole FFs (highlighted by blue triangles). The H_W densities are very similar for prosECCo, low CH dipole, and AIMD (highlighted by green triangles). The O_W densities for AIMD and low CH dipole are similar, however, none of the FFMDs replicate the H_W and O_W clouds off the corner of the TMA tetrahedron (highlighted in purple) at any density contour level.

To further investigate the sensitivity of the neutron scattering signal to different aspects of TMA hydration, we compared g_HTMAHW(r) (which is directly related to the experimentally measured quantity) among all the different variants of the TMA force fields. While we previously showed (Fig. 3) that different charge distributions lead to strikingly different hydration patterns around TMA and different ion-pairing behavior, resulting in different patterns visible at the H_W density map around TMA (Fig. 8), these differences do not significantly modify the g_HTMAHW(r) (see Fig. 8, where the g_HTMAHW(r) computed with different force fields are compared). Neutron scattering experiments examining the g_HTMAHW(r) correlation are thus unable to distinguish the different hydration patterns and ion pairing propensities that we have shown to exist for these different force fields. It is rather unexpected that these different three-dimensional water orientations around the different TMA force fields give such similar RDFs for g_HTMAHW(r). The origin of the similarity of the RDFs despite the significant differences in the density maps is as follows. The hydrogen density clouds have symmetric ordering around the center of the TMA molecule. However, the substituted nuclei are not at this symmetry center, and there is a correlation between each of the clouds of H_W density and each of the 12 substituted hydrogens, which makes the function g_HTMAHW(r) rather broad. The g_HTMAHW(r) component of the structural data, despite being a large fraction of the total scattering, is thus not very informative of the relevant ion hydration structure. Even if we had examined the correlations between the central nitrogen of the TMA (N_TMA) and H_W, the corresponding RDFs would still be all remarkably similar to each other for all of these force fields (see Fig. 6), and would have produced a far smaller NDIS signal. Only, if the correlation between the central nitrogen and the oxygen in water could be measured, this would enable differentiating between these force fields. Unfortunately, isolating this component of the total scattering signal with isotopic substitution is not possible, as no suitable oxygen isotopes exist for such an experiment.

Fig. 8 Density maps of water hydrogen atom around TMA, at twice bulk density. Shown to the right of each density map is the corresponding H_TMA–H_W radial distribution function, compared to the same function for the prosECCo TMA force field to help comparison.

4 Conclusions

Neutron scattering experiments with double isotopic substitution on both H_TMA and H_water were performed allowing us to single out the correlation between H_TMA and H_water. Thanks to the very high contrast of the H/D substitution and the large number of H atoms in TMA, the signal is well above the noise level. It should thus allow for a detailed characterization of the hydration structure. However, we show here that molecular interpretation of the experimental signal proves to be very challenging as all experimentally measurable RDFs (g_HTMAHW(r) and g_NH(r)) are not very sensitive to different hydration patterns caused by changes in the employed force fields. We thus can neither fully infer the most probable hydration structure from such a comparison nor validate the preferential choice of the tested force field variants. Nevertheless, the simulation results provide important insights into TMA hydration and its sensitivity to force field parameters. Shifting from CHARMM (+1.0) to prosECCo (+0.75) charges had a relatively small effect on the ion hydration. Changing the polarity of the C–H bond from −0.35 C and +0.23 H (prosECCo) to −0.1 and +0.1 (low CH dipole) largely inverts the hydration structure of oxygens around the TMA. AIMD simulations show a hydration structure very similar to low CH dipole at the faces and edges of the TMA tetrahedron. The AIMD hydration structure seen at the corners of the TMA tetrahedron is, however, not replicated by any FFMD. Thus, we suggest that the C–H bond polarity in standard CHARMM and its variants is too high to capture the TMA hydration properly. Removal of the tetrahedral structure by employing a center-bead force field for TMA and converting it to a single large bead has a relatively minor effect on the radial hydration structure of the ion. Still, it strongly changes the form of the counterion interaction. Namely, our findings imply that longer-range ion–ion ordered structures may not be accurately replicated as the structure of TMA is simplified to the level of a single bead within a coarse-grained force field. All these results suggest that caution should be taken when simulating moieties with TMA groups where hydration may play a relevant role, such as for common phospholipids and methylated lysines extensively present in biological systems.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was performed using neutron beam time allocated at the Institut Laue Langevin (proposal 8–3–844). The authors thank Henry Fischer for his help and support. H. M.-S. acknowledges the support of the Czech Science Foundation (project 19-19561S). O. T. acknowledges the Faculty of Mathematics and Physics of the Charles University (Prague, Czech Republic) where he is enrolled as a PhD student, and the International Max Planck Research School for “Many-Particle Systems in Structured Environments” (Dresden, Germany) for support. M. V. acknowledges support by the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90254), Project OPEN-28-18. E. D.-D. acknowledges support from the “Initiative dExcellence” program from the French State (Grant “DYNAMO”, ANR11-LABX-0011). P. J. acknowledges support from the European Research Council via an ERC Advanced Grant no. 101095957.

References

1. P. Cheung and P. Lau, Mol. Endocrinol., 2005, 19, 563–573.
2. P. H. von Hippel and K.-Y. Wong, J. Biol. Chem., 1965, 240, 3909–3923.
3. A. Brändström, Adv. Phys. Org. Chem., 1977, 267–330.
4. A. Mustain, B. S. Gupta, M. Taha and M.-J. Lee, New J. Chem., 2023, 47, 12304–12313.
5. J. Turner, A. K. Soper and J. L. Finney, Mol. Phys., 1990, 70, 679–700.
6. N. G. Polydorou, J. D. Wicks and J. Z. Turner, J. Chem. Phys., 1997, 107, 197–204.
7. E. W. Lang, S. Bradl, W. Fink, H. Radkowitsch and D. Girlich, J. Phys.: Condens. Matter, 1990, 2, SA195–SA200.
8. D. Bhowmik, N. Malikova, G. Mériguet, O. Bernard, J. Teixeira and P. Turq, Phys. Chem. Chem. Phys., 2014, 16, 13447–13457.
9. J. Turner, A. Soper and J. Finney, Mol. Phys., 1992, 77, 411–429.
10. J. Z. Turner, A. K. Soper and J. L. Finney, J. Chem. Phys., 1995, 102, 5438–5443.
11. E. J. Nilsson, V. Alfredsson, D. T. Bowron and K. J. Edler, Phys. Chem. Chem. Phys., 2016, 18, 11193–11201.
12. P. H. K. D. Jong and G. W. Neilson, J. Phys.: Condens. Matter, 1996, 8, 9275–9279.
13. P. E. Mason, G. W. Neilson, C. E. Dempsey and J. W. Brady, J. Am. Chem. Soc., 2006, 128, 15136–15144.
14. P. E. Mason, S. Ansell and G. W. Neilson, J. Phys.: Condens. Matter, 2006, 18, 8437–8447.
15. G. W. Neilson, P. E. Mason, S. Ramos and D. Sullivan, Philos. Trans. R. Soc., A, 2001, 359, 1575–1591.
16. J. Turner and A. K. Soper, J. Chem. Phys., 1994, 101, 6116–6125.
17. E. Pluhařová, H. E. Fischer, P. E. Mason and P. Jungwirth, Mol. Phys., 2014, 112, 1230–1240.
18. T. Martinek, E. Duboué-Dijon, Š. Timr, P. E. Mason, K. Baxová, H. E. Fischer, B. Schmidt, E. Pluhařová and P. Jungwirth, J. Chem. Phys., 2018, 148, 222813.
19. H. E. Fischer, G. J. Cuello, P. Palleau, D. Feltin, A. C. Barnes, Y. S. Badyal and J. M. Simonson, Appl. Phys. A: Mater. Sci. Process., 2002, 74, s160–s162.
20. P. E. Mason, H. E. Fischer, P. Jungwirth and S. Timr, Towards a fuller understanding of protein lipid interactions, Institut Laue-Langevin (ILL) 2015 DOI:10.5291/ILL-DATA.8-03-844.
21. G. J. Herdmants and G. W. Neilsont, J. Phys.: Condens. Matter, 1996, 4, 627–638.
22. J. E. Enderby, D. E. Williams and J. Randall, Proc. R. Soc. London, Ser. A, 1975, 345, 107–117.
23. P. E. Mason, G. W. Neilson, C. E. Dempsey and J. W. Brady, J. Am. Chem. Soc., 2006, 128, 15136–15144.
24. J. Chandrasekhar, D. C. Spellmeyer and W. L. Jorgensen, J. Am. Chem. Soc., 1984, 106, 903–910.
25. 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. K. 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, J. Phys. Chem. B, 1998, 102, 3586–3616.
26. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindah, SoftwareX, 2015, 1–2, 19–25.
27. I. Leontyev and A. Stuchebrukhov, Phys. Chem. Chem. Phys., 2011, 13, 2613.
28. E. Duboué-Dijon, M. Javanainen, P. Delcroix, P. Jungwirth and H. Martinez-Seara, J. Chem. Phys., 2020, 153, 050901.
29. V. Kostal, P. E. Mason, H. Martinez-Seara and P. Jungwirth, J. Phys. Chem. Lett., 2023, 4403–4408.
30. S. V. Bennun, M. I. Hoopes, C. Xing and R. Faller, Chem. Phys. Lipids, 2009, 159, 59–66.
31. J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O'Connor, D. J. Tobias, C. Mondragon-Ramirez, I. Vorobyov, A. D. MacKerell and R. W. Pastor, J. Phys. Chem. B, 2010, 114, 7830–7843.
32. J. Lee, X. Cheng, J. M. Swails, M. S. Yeom, P. K. Eastman, J. A. Lemkul, S. Wei, J. Buckner, J. C. Jeong, Y. Qi, S. Jo, V. S. Pande, D. A. Case, C. L. Brooks, A. D. MacKerell, J. B. Klauda and W. Im, J. Chem. Theory Comput., 2016, 12, 405–413.
33. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
34. S. Miyamoto and P. A. Kollman, J. Comput. Chem., 1992, 13, 952–962.
35. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
36. S. Nosé, Mol. Phys., 1984, 52, 255–268.
37. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
38. S. Páll and B. Hess, Comput. Phys. Commun., 2013, 184, 2641–2650.
39. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
40. Y. Zhang and W. Yang, Phys. Rev. Lett., 1998, 80, 890.
41. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1998, 80, 891.
42. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
43. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
44. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
45. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
46. M. Krack, Theor. Chem. Acc., 2005, 114, 145–152.
47. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
48. G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–488.
49. W. F. V. Gunsteren and H. J. C. Berendsen, Mol. Phys., 1982, 45, 637–647.
50. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
51. J. Hutter, M. Iannuzzi, F. Schiffmann and J. Vandevondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 15–25.
52. T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöβ, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103.
53. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
54. M. T. H. Nguyen, O. Tichacek, H. Martinez-Seara, P. E. Mason and P. Jungwirth, J. Phys. Chem. B, 2021, 125, 3153–3162.
55. R. Nencini, C. Tempra, D. Biriukov, J. Polák, D. Ondo, J. Heyda, S. O. Ollila, M. Javanainen and H. Martinez-Seara, Biophys. J., 2022, 121, 157a.
56. E. Pluhařová, P. E. Mason and P. Jungwirth, J. Phys. Chem. A, 2013, 117, 11766–11773.
57. M. Kohagen, P. E. Mason and P. Jungwirth, J. Phys. Chem. B, 2014, 118, 7902–7909.

---

Footnotes

† Electronic supplementary information (ESI) available: Additional neutron scattering and AIMD analysis. See DOI: https://doi.org/10.1039/d3cp05449g

‡ Present address: Department of Theoretical Biophysics, Max Planck Institute of Biophysics, Max-von-Laue Straße 3, 60438 Frankfurt am Main, Germany.

---