Hydration of biologically relevant tetramethylammonium cation by neutron scattering and molecular dynamics † Philip Mason * a , Tomas Martinek a , Balázs Fábián a ‡ , Mario Vazdar b , Pavel Jungwirth a , Ondrej Tichacek a , Elise Duboué-Dijon c , and Hector Martinez-Seara * a 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.

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. 4 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][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 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? 8 Neutron scattering experiments have provided insight into TMA hydration structure in TMA chloride or bromide solutions. 5,6,[9][10][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][13][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,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. 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.

Methods

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. 19 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,b) were recorded for about 2 hrs for each D 2 O solution and for 4 hrs for each H 2 O solution [START_REF] Mason Philip | Towards a fuller understanding of protein lipid interactions[END_REF] The results were then corrected for multiple scattering and absorption and normalized against a standard vanadium scatterer. [START_REF] Herdmants | [END_REF] Taking 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):

X D 2 O H non (Q) and ∆S X H 2 O H non (Q) (Fig.
∆S X D 2 O H non (Q) = 90.2 • S H TMA H W (Q) + 39.3 • S H TMA O (Q) + 6.5 • S H TMA C (Q) + 2.3 • S H TMA N (Q) + 2.3 • S H TMA Cl (Q) + 4.3 • S H TMA H TMA (Q) -144.8, (1)
∆S X H 2 O H non (Q) = -50.6 • S H TMA H W (Q) + 39.3 • S H TMA O (Q) + 6.5 • S H TMA C (Q) + 2.3 • S H TMA N (Q) + 2.3 • S H TMA Cl (Q) + 4.3 • S H TMA H TMA (Q) -4.0.
(2) These prefactors were calculated from the atomic concentration and neutron scattering lengths of the various elements in the system by standard literature methods. 22 The 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.

∆∆S H non (Q) = ∆S X D 2 O H non (Q) -∆S X H 2 O H non (Q) = 140.8 • (S H TMA H W (Q) -1) (3)
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 1 H and the Placzek effect, 23 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 24 (-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 25 . All systems were assembled using GROMACS 2021.2 and 2021.5 tools. 26 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 physi-cally based Electronic Continuum Correction (ECC) approach accounts for electronic polarization in a mean-field way by scaling charges 27 . 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 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. 30 All systems were simulated with GROMACS 2021, 26 using simulation parameters as provided by CHARMM-GUI. 32 The LINCS and Settle algorithms were used to constrain the geometry of TMA hydrogens 33 and water molecules, respectively, allowing the use of 2 fs time step in the simulations. 34 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 35 to gather enough water orientation statistics around neopentane by excluding neopentane-neopentane contacts. Simulations were performed in the isothermal-isobaric (N pT ) ensemble. Temperature was controlled using a Nosé-Hoover thermostat 36 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 37 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 38 for neighbors. Long-range Coulomb interactions were accounted for using particle mesh Ewald (PME) with a 1.2 nm cutoff as implemented in GROMACS. 26 All simulation data can be found at zenodo.org (DOI: 10.5281/zenodo.10406618).

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 neu- tralized 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][40][41] with the D3 dispersion correction [42][43][44] . Core electrons were replaced by GTH pseudopotentials, 45,46 , and the triple-ζ basis set TZV2P with polarization functions was used for valence electrons. 47 A cutoff of 400 Ry was used for the auxiliary plane wave basis set in the GPW method. 48 The system was first equilibrated using AIMD Langevin dynamics for at least 16 ps with a damping constant of γ = 0.02 ps -1 49 . During the production simulation, the temperature was set to 300 K via a global CSVR thermostat with a time constant of 1 ps. 50 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][52][53] All simulation data can be found at zenodo.org (DOI: 10.5281/zenodo.10406618).

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.

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 inverse Fourier-transform

∆∆G H non (r) = F -1 [∆∆S H non (Q)],
which is directly related to the single pair-correlation function g H TMA H W :

∆∆G H non (r) = ∆G X D 2 O H non (Q) -∆G X H 2 O H non (Q) = 140.8 • (g H TMA H W (r) -1) (4)
∆∆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 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.

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 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. 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 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 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 + 56 , and with divalent ions such as Ca 2+ 57 , 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 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.

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

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 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 and7). 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 de- spite 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 po- 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) Total diffraction patterns for H 2 O solutions of d12-TMACl(black) and h12-TMACl(grey). b) Total diffraction patterns for D 2 O solutions of d12-TMACl(black) and h12-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 first-order differences shown in c).

Fig. 2

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

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

Fig. 4

 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

 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.

  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 (Figure5 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

Fig. 6

 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

 7 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 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 blue) at any density contour level.

  J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Fig. 8

 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

 1 TMA and neopentane simulation models used in this study. In parenthesis, the used CHARMM36 atoms types 31 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)
	center-bead	sphere		+0.75		+0.75
			C (CT)	C (CT3)	H (HA3)
	neopentane	TMA	0.00	-0.30	0.10	0.0

1-12

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

| 1-12

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. EDD acknowledges support from the "Initiative d'Excellence" program from the French State (Grant "DY-NAMO", ANR11-LABX-0011). P.J. acknowledges support from the European Research Council via an ERC Advanced Grant no. 101095957.

Conflicts of interest

There are no conflicts to declare.

Notes and references


Luc Souilla
email: lucsouilla@gmail.com

Pascal Amedro

Shawnda A Morrison5

Children With Cardiac Disease and Heat Exposure: Catastrophic Converging Consequences?

Keywords: thermoregulation, child, environmental epidemiology, heatwaves, hyperthermia, physical literacy, congenital heart, disease, acquired heart disease

The detrimental impact of extreme heat exposure on the health and well-being of children is widely acknowledged. The direct and indirect effects of climate change have led to an increased risk of certain cardiovascular events which may be particularly harmful to children who are born with, or develop, heart disease. Purpose: To highlight the worrying paucity of investigative research aimed at differentiating how higher ambient temperatures further tax an already compromised cardiovascular system in children.

Methods:

This commentary describes basic thermoregulatory concepts relevant to the healthy pediatric population and summarizes common heart diseases observed in this population. Results: We describe how heat stress and exercise are important factors clinicians should more readily consider when treating children with heart disease. Countermeasures to physical inactivity are suggested for children, parents, clinicians, and policymakers to consider.

Conclusions:

As sudden, excessive heat exposures continue to impact our rapidly warming world, vulnerable populations like children with underlying heart conditions are at greater heat health risk, especially when coupled with the negative physical activity and fitness trends observed worldwide.

--

The current climate breakdown is considered unprecedented in Earth's history. Under all scenarios, global warming is expected to reach +1.5 °C in the near term (2021-2040), with a 50% probability of reaching higher levels, even considering the most optimistic scenario(s) [START_REF]Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]. By 2100, 75% of the world's population may be subject to at least 20 days of "extreme heat conditions" per year [START_REF] Mora | Global risk of deadly heat[END_REF] as the frequency, intensity, and duration of heat waves [START_REF] Helldén | Climate change and child health: a scoping review and an expanded conceptual framework[END_REF] will also increase [START_REF] Mcgregor | Vagues de chaleur et santé: guide pour l'élaboration de systèmes d'alerte[END_REF]. Presently, there is a clear U-shaped relationship between increased ambient temperatures and all-cause mortality (4), and risk of emergency department visits [START_REF] Bernstein | Warm Season and emergency department visits to U.S. children's hospitals[END_REF]. The main expected cause of morbidity in children includes heat strain, exhaustion, heat stroke [START_REF] Hifumi | Heat stroke[END_REF], and, sometimes, death [START_REF]Human Health: Impacts, Adaptation and Co-Benefits, Chapter 11[END_REF]. Since sudden increases in ambient temperature (ie, heatwaves) exert immediate, multisystem organ effects, maintaining a healthy and robust cardiovascular system is critical to adequately respond to these unexpected high heat events. The cardiovascular system plays an important dual role for humans exposed to hot conditions; it initiates increases in heart rate to facilitate the additional cardiac output needed to manage blood flow redistribution, and it delivers this blood to needed locations (eg, skin). Thermoregulatory responses to the heat have been well studied in adult humans [START_REF] Blois | The effects of climate change on cardiac health[END_REF] but much less so in the pediatric population [START_REF] Morrison | Moving in a hotter world: Maintaining adequate childhood fitness as a climate change countermeasure[END_REF][START_REF] Smith | Pediatric thermoregulation: considerations in the face of global climate change[END_REF]. Prepubertal responses to heat stress are usually driven with a greater emphasis on "dry" mechanisms of heat exchange, that is, shunting a greater proportion of cardiac output to the skin [START_REF] Smith | Pediatric thermoregulation: considerations in the face of global climate change[END_REF], and under most situations, this response is sufficient to avoid overt heat illness. Importantly, when ambient temperature is near/equivalent to body temperature, dry mechanisms become far less effective, increasing the risk of an uncompensable situation where heat syncope, heat cramps, and nonfatal heat stroke may result. Children's sweat composition differs from fully mature adults [START_REF] Meyer | Fluid and electrolyte loss during exercise[END_REF]; it can be higher in sodium which may place children at greater risk of hypovolemic heat stress. Children who are not able to increase cardiac output as efficiently as healthy children may experience a thermoregulatory disadvantage and greater heat health risk. The latest United Nations Children's Fund technical brief released in May, 2023, highlights the need to protect children from heat stress, first by attempting to widen communication surrounding pediatric thermoregulatory mechanisms, and also by setting factors that define child vulnerability [START_REF]Protecting children from heat stress: a technical note[END_REF]. Although this report considers heart disease as a risk factor to heat exposure, there is minimal evidence displayed regarding risk stratification.

Children With Cardiac Disease and Heat Exposure

Heart conditions vary widely, with some affecting the structure, electrical system, or blood vessels of the organ. Congenital heart disease (CHD), the most common type of birth defect, is defined as a significant abnormality of the heart and/or great vessels at birth, with different levels of complexity: simple, moderate, and complex CHD. Although mortality has declined to 34% in 2017 [START_REF] Gbd | Global, regional, and national burden of congenital heart disease, 1990-2017: a systematic analysis for the global burden of disease study 2017[END_REF], CHD still exerts a heavy economic burden and significant mortality, with the existence of common complications (eg, heart failure, arrhythmia, and pulmonary hypertension). Worldwide, inequities are overrepresented in diagnosis, treatment, and mortality of CHD children between low-, middle-, and high-income countries.

Acquired or genetic heart diseases, such as pediatric cardiomyopathies (22; eg, dilated or hypertrophic), alter cardiac function, leading to impaired cardiac output. Inherited arrhythmias (eg, Long QT), another genetic abnormality caused by dysfunctional ion channels in the heart, are recognized as a cause of sudden death in structurally normal hearts [START_REF] Tester | The molecular autopsy: should the evaluation continue after the funeral?[END_REF]. Acquired heart disease can also be caused by infection, for example, rheumatic fever or infective endocarditis.

Vulnerability to heat exposure is driven by individual susceptibility. Globally, there are around 9 million deaths annually attributed to greenhouse gases, with 62% of them due to cardiovascular disease [START_REF] Gbd | Risk Factors Collaborators. Global burden of 87 risk factors in 204 countries and territories, 1990-2019: a systematic analysis for the global burden of disease study 2019[END_REF]. Thus, the combined effects of heat, air pollution, socioeconomic, and health status are (each) responsible for avoidable acute cardiovascular disease events. While children with common CHD such as small ventricular or atrial defects may not be at greater risk of heat-related cardiac events, children with severe cardiac complications like pulmonary hypertension or heart failure will not be able to compensate for the increased circulatory demand needed during high heat events. Children with inherited arrhythmia are at higher risk of experiencing heat-related cardiac events since body temperature is an independent determinant of heart rate, and heatstroke can impact electrocardiographic conductance which prolongs QT interval and ventricular tachycardia conduction defect [START_REF] Akhtar | Electrocardiographic abnormalities in patients with heat stroke[END_REF]. Notably, children with heart conditions are often prescribed drugs such as beta-blockers, diuretics, or angiotensin-converting enzyme inhibitors. Beta-blockers decrease blood flow to the skin and cardiac output, and diuretics and angiotensin converting enzyme inhibitors can also increase the risk of hypovolemia and heat stroke [START_REF] Layton | Heatwaves, medications, and heat-related hospitalization in older Medicare beneficiaries with chronic conditions[END_REF].

Vulnerability to heat exposure is also driven by sociocultural factors (eg, poverty) and geographical factors (Figure 1). Children with heart disease living in hotter and more humid regions are more likely to have inadequate access to health care [START_REF]Climate Change 2023: Synthesis Report[END_REF]. According to the World Health Organization, child mortality related to heat exposure is predicted to reach over 100,000 deaths per year by 2050, with the greatest climate change-related mortality occurring in South Asia [START_REF]Quantitative Risk Assessment of the Effects of Climate Change on Selected Causes of Death[END_REF]. This, combined with the fact that mortality rate for CHD declines according to a country's increasing sociodemographic index, forecasts that the majority of death will disproportionally occur in low and lower-middle sociodemographic index countries [START_REF] Gbd | Global, regional, and national burden of congenital heart disease, 1990-2017: a systematic analysis for the global burden of disease study 2017[END_REF]. To summarize, higher ambient temperature leads to heat-related illness diagnoses, severity, and prevalence of heart diseases in children, although direct data on this relationship are lacking. Hot environments are associated with poorer air quality, which can lead to more CHD cases, and medications commonly prescribed to CHD patients are independently associated with greater risk of heat-related illnesses. The technical brief from UNICEF did not address these complex interactions.

Practical Implications of Exercising in the Heat

Moving one's body in hot conditions places further burden on the cardiovascular system, since a by-product of muscle contraction is heat production, and exercise (depending on the intensity and duration of the activity) can be a pyrogenic activity itself [START_REF] Bradford | Exercise can be pyrogenic in humans[END_REF]. Physical activity (PA) represents a very wide range of pursuits, from sport participation to outdoor activities, or even exercise rehabilitation, which is recognized as a central component of cardiac rehabilitation. For instance, children and adolescents playing American football are at higher risk for heat stress since the sport requires protective gear which reduces heat dissipation and climate change acts as both a cause or consequence to greater physical inactivity ( 21) which currently affects >80% of adolescents globally [START_REF] Guthold | Global trends in insufficient physical activity among adolescents: a pooled analysis of 298 population-based surveys with 1•6 million participants[END_REF]. Children with heart conditions already present with lower cardiorespiratory fitness compared to healthy controls [START_REF] Amedro | Cardiopulmonary fitness in children with congenital heart diseases versus healthy children[END_REF][START_REF] Souilla | Cardiorespiratory fitness, muscle fitness, and physical activity in children with long QT syndrome: a prospective controlled study[END_REF], and this decline remains so over time [START_REF] Gavotto | Aerobic fitness change with time in children with congenital heart disease: A retrospective controlled cohort study[END_REF]. Higher fitness elicits efficient heat balance since fitter people sweat more and can tolerate higher core temperatures (24) than unfit people. There are multiple structural and functional cardiovascular adaptations which lessen external system demands when people move in the heat (28) (Figure 2). It is therefore critical to adapt public health policies in favor of achieving adequate PA to be better equipped and meet the challenges of global warming and other environmental challenges.

The UNICEF technical brief provides essential clues for public health policymakers advocating for a multisectoral approach to minimize the detrimental health outcomes for children exposed to the heat. The specific approaches can be viewed in the report, but briefly, UNICEF describes a risk communication strategy called "B.E.A.T. the Heat," where B.E.A.T. is an acronym for B=BE AWARE of heat stress and protect yourself, E = EASILY IDENTIFY the symptoms of heat stress, A = ACT IMMEDIATELY to protect the individuals affected, and T = TAKE to a health facility if the person is presenting with severe symptoms. For children with heart conditions, in addition to the approaches specified in the policy brief, cardiologists and clinicians should be mindful of the risk stratification regarding diagnosis and clinical status of each child.

Countermeasures to Physical Inactivity in Children With Heart Disease Robust and Diverse Physical Education Programs in Schools

Globally, schools are considered an important avenue to reach children for health promotion, physical education, and after-school sport [START_REF] Juric | Effectiveness of a population-scaled, school-based physical activity intervention for the prevention of childhood obesity[END_REF]. For children with cardiac conditions, individualized adaptation strategies must be implemented to ensure safe participation in PA. It is important for physical education teachers to work closely with cardiologists to adapt exercise sessions to each student's unique needs.

Coordinated Kinesiology and Cardiology Expertise in Hospital Settings

Children with cardiac conditions meet many health professionals in the hospital system, but direct involvement with exercise specialists remains rare worldwide. PA promotion and education should be incorporated in every health care setting. Exercise rehabilitation programs (8) are an excellent opportunity to expand in scope since they are already considered a "nondrug" care management approach. Proper medication should match exercise prescription to avoid increasing the child's risk of hypovolemia, and it is important to stress that the cardiovascular prevention community has a critical role in global health [START_REF] Gulati | The role of the preventive cardiologist in addressing climate change[END_REF]. Physical Fitness Inclusion in Holistic Heat Health Policies Health policies need to prioritize physical fitness in children with cardiac conditions by promoting greater active transportation and nature-based infrastructure. Public health policies may also address significant disparities in access to care and mortality rates for children with heart disease in diverse socioeconomic countries. There remain significant gaps in this research area. Thus, the "Top 4" suggestions for managing increased heat risk projections for children with cardiac conditions are: (1) health professionals, including cardiologists, need to be knowledgeable on the risks that a hot environment poses to children with cardiac conditions. (2) Physical education teachers and kinesiologists must consider the type, intensity, and environmental conditions to accommodate PA in this population. (3) Health policies must create greater health equity with children at the fore (eg, indoor activities, water access, and food for refueling). ( 4) Conduct more experimental research on how children respond to passive and active heat exposures, especially those with cardiac diseases, and across all income-strata countries.

Conclusions

The expected increase in the frequency, intensity, and duration of heatwaves worldwide will pose additional heat health risks to humans. Children with cardiac conditions already present with lower cardiorespiratory fitness, so it will be necessary to consider the specific heat health impacts of a warming world on this vulnerable population.