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

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 (HTMA) and hydrogens in water (HW). 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.


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. 1 Tetraalkylammonium salts are also powerful denaturants 2 and are widely used as phase transfer catalysts 3 and are also seen frequently in ionic liquids. 4MA 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][6][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 pe-riodic 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? 80][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.3][14][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 2 O and D 2 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,10However, 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. 11 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.

Neutron scattering measurements
NDIS measurements were performed using the D4C diffractometer at the nuclear reactor at the Institut Laue-Langevin in Grenoble, France. 19Four chemically identical solutions of 2 m TMACl in water were prepared, which differed only in H/D substitution on the TMA and H/D substitution on water.The four diffraction patterns (Fig. 1a,b) were recorded for about 2 hrs for each D 2 O solution and for 4 hrs for each H 2 O solution.The results were then corrected for multiple scattering and absorption and normalized against a standard vanadium scatterer. 20king the difference between the diffraction patterns associated with solutions that differ only by the H/D substitution on TMA (both in D 2 O and H 2 O solutions) yields the first-order differences ∆S 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): (2) These prefactors were calculated from the atomic concentration and neutron scattering lengths of the various elements in the system by standard literature methods. 21The difference between eq. 1 and eq. 2 yields the second order difference ∆∆S H non (Q) (Fig. 1d), which reports on the single correlation between the TMA non-exchangeable H atoms and the water H atoms.
(3) This function provides a useful internal consistency check for the accuracy of the solutions and the multiple scattering and ab-sorption corrections performed on the data.Due to the large inelastic scattering of 1 H and the Placzek effect, 22 hydrogencontaining samples always present a dominant background.The higher the atomic concentration of 1 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 1 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.

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 CHARMMCl_o type 23 (−1 charge) or CL_2s type 13 (−0.75scaled charge).Each system contained 1388 TIP3P water molecules with Lennard-Jones parameters also on the hydrogens 24 .All systems were assembled using GROMACS 2021.2 and 2021.5 tools. 25r 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 26 .ECC has been shown to improve significantly the description of ion pairing behavior in solution. 27,28Based 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 CHARMM36 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. 29l systems were simulated with GROMACS 2021, 25 using sim-ulation parameters as provided by CHARMM-GUI. 31The LINCS and Settle algorithms were used to constrain the geometry of TMA hydrogens 32 and water molecules, respectively, allowing the use of 2 fs time step in the simulations. 33The coordinates were saved every 10 ps for each 2 µs simulation.For neopentane that aggregates due to its hydrophobic nature, we saved coordinates every 20 fs to gather enough water orientation statistics to study its coordination before aggregation dominates (i.e., ≈ 0.5 ns).Simulations were performed in the isothermal-isobaric (N pT ) ensemble.Temperature was controlled using a Nosé-Hoover thermostat 34 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 35 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 36 for neighbors.Long-range Coulomb interactions were accounted for using particle mesh Ewald (PME) with a 1.2 nm cutoff as implemented in GROMACS. 25

Ab Initio molecular dynamics (AIMD) simulations
As a benchmark for the TMA solvation structure, we complemented the neutron diffraction experiments and the force fieldbased 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 [37][38][39] with the D3 dispersion correction [40][41][42] .Core electrons were replaced by GTH pseudopotentials, 43,44 , and the triple-ζ basis set TZV2P with polarization functions was used for valence electrons. 45A cutoff of 400 Ry was used for the auxiliary plane wave basis set in the GPW method. 46The system was first equilibrated using AIMD Langevin dynamics for at least 16 ps with a damping constant of γ = 0.02 ps −1 47 .During the production simulation, the temperature was set to 300 K via a global CSVR thermostat with a time constant of 1 ps. 48To enhance sampling, ten 50 ps parallel AIMD simulations totaling 500 ps of trajectory were used for further structural analysis.[51]

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. 52, 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.

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 H non (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 H non (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 Fourier-transform , which is directly related to the single pair-correlation function g H TMA H W : ∆∆G H non (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 dy-namics 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.

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. 26,27A 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 1/ √ ε ∞ ≃ 0.75, where ε ∞ is the highfrequency 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. 27,53Here, we compared the structure of the solution simulated using this scaled-charge force field with that obtained using the standard full charge CHARMM36 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 only the first 30 ps (using 1500 frames, 20 fs between frames) of this simulation for analysis.
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 + 54 , and with divalent ions such as Ca 2+ 55 , 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.
If we now compare TMA hydration with that of neutral neopen-tane, 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.

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).
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).

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 Fig. 3 Comparison 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.
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.
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 in-teraction 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).

Comparison of neutron scattering data to MD simulations
The double difference signal, ∆∆G H non (r), obtained from neutron scattering experiments after Fourier transform is composed of a single radial distribution function of g H TMA H W (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.
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. 2 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 (Figure 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.
To further investigate the sensitivity of the neutron scattering signal to different aspects of TMA hydration, we compared g H TMA H W (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 H TMA H W (r) (see Fig. 8, where the g H TMA H W (r) computed with different force fields are compared).Neutron scattering experiments examining the g H TMA H W (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 H TMA H W (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 H TMA H W (r) rather broad.The g H TMA H W (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.

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 H TMA H W (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.

Fig. 1 a
Fig. 1 a) Total diffraction patterns for H 2 O solutions of d-TMACl(black) and h-TMACl(grey).b) Total diffraction patterns for D 2 O solutions of d-TMACl(black) and h-TMACl(grey).c) First order differences ∆S X H 2 O Hnon (Q) (red line, obtained by the difference of the two diffraction patterns shown in a) and ∆S X D 2 O Hnon (Q) (blue line, obtained by the difference of the two diffraction patterns shown in b).d) Second order difference 140.8 • S H TMA H W (Q) − 1 , obtained through the difference of the two firstorder differences shown in c).

Fig. 2
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 Å −1 .The lower red function is the real space version of the upper red function.

Fig. 4
Fig.4For all density maps, the density of chloride ion (green) is 6x bulk density, oxygen (red) is 3x bulk density and hydrogen (white) is 2x 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.

Fig. 5
Fig. 5 For all density maps, density of chloride ion (green) is 6x bulk density, oxygen (red) is 3x bulk density, and hydrogen (white) is 2x 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.

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

Fig. 7 O
Fig. 7 O W and H W density maps around TMA in various simulations.Upper O W density is shown in red at 3.3x bulk density (approximated for AIMD simulations) and lower 2.8x bulk density for O W (red) and 2.1x bulk density for H W (white).The orientation of the O W around TMA shows an inverse relationship between prosECCo and low CH dipole (highlighted by blue triangles).The H W densities are very similar for prosECCo and low CH dipole FFs (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 blue) at any density contour level.

Fig. 8
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.

Table 1
30A and neopentane simulation models used in this study.In parenthesis, the used CHARMM36 atoms types30are listed.