The hydroxyl radical is the primary reactive oxygen species produced by the radiolysis of water and is a significant source of radiation damage to living organisms. Mobility of the hydroxyl radical at low temperatures and/or high pressures is hence a potentially important factor in determining the challenges facing psychrophilic and/or barophilic organisms in high-radiation environments (e.g., ice-interface or undersea environments in which radiative heating is a potential heat and energy source). Here, we estimate the diffusion coefficient for the hydroxyl radical in aqueous solution using a hierarchical Bayesian model based on atomistic molecular dynamics trajectories in TIP4P/2005 water over a range of temperatures and pressures.

Ionizing radiation is a feature of both terrestrial and extraterrestrial environments, presenting challenges as well as opportunities for living organisms. On the one hand, radiation can damage biological molecules either directly or via the production of reactive chemical species that modify or degrade them.1 On the other hand, ionizing radiation can act as a non-photosynthetic energy source for microbial communities, in some cases producing chemical species suitable for chemotrophy2 or in other cases maintaining a habitable environment by temperature elevation and/or maintenance of liquid water pockets via ice melting.3–5 While some terrestrial organisms are known to tolerate high levels of radiation either transiently or on an ambient basis,6–8 the role of ionizing radiation in determining habitability in a broader biological context remains largely open.

For conventional microbial organisms in aqueous environments, the predominant threat posed by ionizing radiation is the formation of reactive oxygen species due to the radiolysis of water.9 Many species are produced, including hydrogen peroxide, the superoxide radical, and the hydroperoxyl radical;10 however, the hydroxyl radical (·OH) is the dominant source of prompt radiation damage for solvated biomolecules due to its high production rate, reactivity, and unsuitability for enzymatic processing.11 In terrestrial mesophilic organisms under typical cellular conditions, ·OH survives on sufficiently long time scales to diffuse to and damage biological macromolecules,12 prominently including proteins and DNA. Although the problems associated with DNA damage are well-appreciated, proteins are the major cellular targets of ·OH-mediated radiation damage.13 Post-translational modification of proteins due to radical interactions can lead to a wide array of potentially lethal consequences, including formation of reactive peroxide species,14 formation of insoluble aggregates,15 loss of enzymatic function, and destabilization of functional complexes; at minimum, such damage increases the rate of protein expression and controlled degradation required for homeostasis, thereby raising the metabolic cost of cellular survival. Although ·OH production is a threat to irradiated organisms in any environment, some environmental conditions may partially or substantially ameliorate it. In particular, environments that favor the scavenging of ·OH by other chemical species reduce the level of chemical stress to which organisms are subject and may thus allow for greater radiation tolerance. Because many reactions involving ·OH are (or are near) diffusion-limited kinetics, the ·OH diffusion coefficient is of particular relevance to extremophile biochemistry in irradiated environments.

Changes in ·OH diffusion may have implications for radiochemistry in high pressure/low temperature environments. The direct products of γ-irradiation of water are formed in “spurs,” localized regions of high concentration for these highly reactive species (most directly, ·OH, ·H, and hydrated electrons).16,17 When a γ-ray interacts with a water molecule, it generates a highly reactive excited state,

H2O+γH2O,
(1)

which can either directly decay to form a hydrogen atom and a hydroxyl radical,

H2OH+OH,
(2)

or it can lose an electron to the surrounding solution, resulting in a radical cation that, in turn, decomposes to yield a hydroxyl radical and a proton,

H2OH2O++eaq,
(3)
H2O+H++OH.
(4)

These reactions serve as the starting point for the formation of a complex mixture of reactive oxygen species18 and/or reactions of these primary products with biomolecules.19 The recent advent of experimental techniques enabling direct detection of the long-hypothesized H2O+ intermediate20 and observation of attosecond dynamics in liquid water21 as well as the discovery of new reaction pathways22,23 have led to renewed interest in the fundamental radiochemistry of aqueous solutions. Analysis of the relevant reactions requires reasonable estimates of the diffusion coefficients of key chemical species, as small solutes can undergo anomalous diffusion in water,24–26 particularly under extreme conditions. ·OH, the focal species of this work, has been the subject of previous studies; however, these have generally focused on high temperature, e.g., Ref. 27. The different patterns of ·OH diffusion seen in the low-temperature, high-pressure regime suggest considerable value in further experimental studies of these environments.

Among the environments of particular interest for novel biochemistry are those involving low temperatures and/or high pressures. Such environments occur in deep ocean and ice/rock interfaces on Earth and in subsurface oceans in the outer solar system. Because the diffusion coefficient for ·OH at low temperature and high pressure has not been measured to date, we here estimate it using atomistic molecular dynamics (MD) simulations, employing a novel hierarchical Bayesian inference scheme to infer the diffusion coefficient while correcting for finite size effects. In the process, we also parameterize a Chemistry at Harvard Molecular Mechanics (CHARMM)-compatible model for ·OH for use with the TIP4P/2005 water model (chosen for its ease of implementation with standard MD platforms and its performance in reproducing properties of bulk water over a wide temperature and pressure range). We summarize our posterior inference in the form of a simple log–log polynomial model that can be used to reproduce our simulation-based estimates of D·OH over a range of temperatures and pressures.

The remainder of this paper is structured as follows: Sec. II describes our procedures, in particular, including the parameterization of the ·OH model (Sec. II A) and inference for the diffusion coefficient (Sec. II C). Our results are summarized in Sec. III, and Sec. IV concludes this paper.

Our interest is inferring the diffusion coefficient of the hydroxyl radical in aqueous solution, D·OH, as a function of temperature and pressure. We begin with parameterization of a model for ·OH in TIP4P/2005 water, followed by our simulation design. We then describe our approach for inferring the diffusion coefficient, D·OH, from simulated water and radical trajectories. Results are shown in Sec. III.

Unlike diffusion of H+, which involves a large anomalous contribution due to the Grotthuss mechanism,28 the diffusion of ·OH is expected to be dominated by classical (hydrodynamic) effects. A transport reaction analogous to that of H+ for ·OH would require the abstraction of a hydrogen atom (not merely a proton) from H2O, an extremely unfavorable reaction under normal conditions. Although more subtle quantum effects (like those hypothesized to lead to enhanced diffusion of OH29–31) are possible, the experimental observation that the diffusion coefficient of ·OH under ambient conditions is close to that of water32 strongly suggests that such effects are unlikely to be large. As such, we here model the diffusion of ·OH using (classical) atomistic molecular dynamics. To perform atomistic simulations of ·OH in solution, we parameterize a CHARMM-based33 model for ·OH in TIP4P/2005 water.34 We employ TIP4P/2005 because of its strong performance in reproducing the diffusion constant of water over a wide range of temperatures and pressures.35 Partial charges, bond length, mass, and force constant for ·OH are taken from Pabis, Szala-Bilnik, and Swiatla-Wojcik,36 who performed combined density functional theory (DFT) and MD studies of ·OH in Bopp, Jancsò and Heinzinger (BJH) water at physiological temperature and pressure;37 unfortunately, the non-bonded interactions employed cannot be directly adapted to the CHARMM forcefield, and hence, it is necessary to parameterize them directly. The non-bonded interactions in question are defined by a Lennard–Jones potential of the form

εLJ(RAB|εA,εB,rminA/2,rminB/2)=εAεBrminA/2+rminB/2RAB12+2rminA/2+rminB/2RAB6,

where RAB is the distance between atoms A and B, ɛA and ɛB are species-specific well-depth parameters, and rmin A/2 and rmin B/2 are “half-radii” that determine the zero-point of the interatomic force. Here, we must determine these parameters for the two respective atoms of ·OH, given the TIP4P/2005 parameters (which we take as fixed). As our interest is in D·OH, we optimize the ·OH non-bonded interaction parameters (ɛO, ɛH, rmin O/2, and rmin H/2) so as to reproduce the measured value of D·OH = 0.23 Å2/ps for ·OH in water at 298 K and atmospheric pressure,32 holding all other factors constant.

Our protocol proceeded as follows: We began with a quasi-random search of the parameter space, drawing 250 points from the intervals ɛO, ɛH ∈ (−0.5, 0) kcal/mol, rmin O/2 ∈ (0, 3) Å, and rmin H/2 ∈ (0, 2) Å using a four-dimensional Halton sequence (bases 2, 3, 5, and 7). (Halton sequences are low-discrepancy sequences38 that are useful for adaptively sampling parameter spaces more evenly than uniform sampling but without the limitations of gridding.) For each parameter vector, a 1 ns atomistic simulation of one ·OH in TIP4P/2005 water under periodic boundary conditions was performed, with frames sampled every 0.5 ps (integrator step size 2 fs). Simulations were initialized with a cubic box of 20 Å side length at 1 atm and 298 K; two adjustment phases of 100 ps each were performed (with box sizes adjusted for Particle Mesh Ewald (PME) calculations after each phase) prior to the production run, with both adjustment and production phases performed with the NpT ensemble. Langevin dynamics with an interval of 1/ps were employed for temperature control, and a Langevin–Nosé–Hoover piston with a period of 100 ps was used to maintain constant pressure.39,40 Rigid bonds were maintained for all waters, with the O–H bond of ·OH left flexible. All simulations were performed using Nanoscale Molecular Dynamics (NAMD),41 with initial conditions created using Visual Molecular Dynamics (VMD),42 psfgen,43 and Packmol.44 Each simulated trajectory was then unwrapped using the protocol of von Bülow, Bullerjahn, and Hummel45 to account for changing box sizes, and each frame was centered at its centroid to correct for net drift. Molecular positions were extracted via the resulting oxygen atom coordinates.

To obtain initial estimates for D·OH at each parameter value, the covariance-based estimator of Bullerjahn, von Bülow, and Hummel46 was applied to each processed ·OH trajectory; this estimator is computationally efficient and was found to work well in pilot runs using both ·OH and H2O under these simulation conditions. These raw estimates were then corrected for finite sample sizes using the analytical correction factor of Yeh and Hummer,47 

D=DMD+ζkBT2πηL,

where D is the diffusion constant at infinite size, DMD is the diffusion constant under periodic boundary conditions (PBCs) obtained from MD simulation, ζ ≈ 2.837 3 is a numerical constant, η is the shear viscosity of the solvent, and L is the box length. For η, the TIP4P/2005 viscosity of 8.55 × 10−4 Js/m3 from González and Abascal48 was employed, and L was taken to be the cube root of the mean box volume over the simulation. The size-corrected estimates of D·OH were then retained for further analysis.

To obtain an initial estimate of the non-bonded parameters, the absolute relative error (ARE) in the D·OH estimate (i.e., |0.23DOĤ|/0.23) was computed for each parameter vector. Inspection of the resulting estimates revealed two regions of high performance (Fig. S1). To distinguish among the competing regions, energies were calculated for the interaction of a single ·OH and H2O in both conventional hydrogen bonding and “flipped” (i.e., unfavorable H–H or O–O configurations) over a range of 1–5 Å (Figs. S2–S4). The single region for which properly oriented hydrogen bonds were more favorable and the ARE was low was selected for further analysis. (The most favorable point in this stage was ɛO ≈ −0.224, ɛH ≈ −0.291, rmin O/2 ≈ 0.802, and rmin H/2 ≈ 0.752.)

Following this first stage of calibration, we performed a subsequent stage of refinement by repeating the procedure with search using larger simulated systems over a smaller parameter range; this secondarily confirmed that estimates were robust to system size. As before, a four-dimensional Halton sequence was used to select parameters, with 25 draws taken over the range ɛO ∈ (−0.26, −0.18) kcal/mol, ɛH ∈ (−0.31, −0.27) kcal/mol, rmin O/2 ∈ (0.6, 1.1) Å, and rmin H/2 ∈ (0.72, 0.78) Å. For each parameter vector, the above simulation and analysis protocol was followed, with the exception that initial box size was increased to 50 Å. Estimated size-corrected AREs were obtained for each parameter value, and the parameter vector yielding the minimum error was selected for final use. (As above, energy calculations were employed to verify that the model favored the correct ·OH–H2O orientation.)

The final parameters for the ·OH model (including both predetermined and calibrated parameters) are rOH = 0.9751 Å, EOH = 475.6 kcal/mol/Å2, mO = 15.994 Da, mO = 1.008 Da, ɛO = −0.235 kcal/mol, ɛH = −0.306 8 kcal/mol, rmin O/2 = 0.781 85 Å, and rmin H/2 = 0.746 939 Å. The relatively small value of rmin O/2 seems to be necessary or proper interaction with the TIP4P/2005 water model, as larger values lead either to inappropriate donor/acceptor orientations for inaccurate values of D·OH (Figs. S3 and S4); the selected value leads to reasonable water/radical interactions, as shown in Fig. 1. We note in passing that the interaction energies here found to generate realistic diffusion behavior are smaller than those obtained by Du and Francisco49 based on DFT-based potential energy scanning calculations (with both differing from36), which suggests some limits on the effectiveness of ab initio calculations for use in parameterizing TIP4P/2005 solutes when targeting bulk properties. The size-corrected diffusion constant estimate for the final model in the last selection round was DOĤ0.222 Å2/ps for an ARE of ∼4% vs the target value of 0.23 Å2/ps.

FIG. 1.

Non-bonded interaction energy for H2O vs ·OH as a function of configuration and distance between proximate atoms. The final model correctly favors standard donor/acceptor orientations over “flipped” orientations, with the H2O donor interaction slightly favored over the ·OH donor interaction.

FIG. 1.

Non-bonded interaction energy for H2O vs ·OH as a function of configuration and distance between proximate atoms. The final model correctly favors standard donor/acceptor orientations over “flipped” orientations, with the H2O donor interaction slightly favored over the ·OH donor interaction.

Close modal

Given the parameterized model, we simulate the diffusion of ·OH in TIP4P/2005 water at multiple temperatures and pressures. The basic simulation strategy is equivalent to that described in Sec. II A, with the following modifications. First, in order to ensure a high level of precision in our trajectory calculations, we sample trajectories every 0.25 ps instead of every 0.5 ps; additionally, we add an extra 100 ps to the pre-production run for each trajectory, which is not included in the 1 ns employed for estimation. Rather than employing the covariance estimator for D·OH, we use a Bayesian approach as described in Sec. II C 2. Finally, because the analytical correction for system size depends on the viscosity (which is not known for TIP4P/2005 over the range of temperatures and pressures studied here), we employ statistical corrections involving variable system sizes as explained below.

We perform simulations at 263, 273, 283, and 298 K, and at 1, 10, 100, 1000, and 10  000 atm pressure (full factorial design). To allow statistical correction for finite size effects (and to reduce simulation-related error), we perform 30 replicate simulations for each condition, with initial box sizes evenly spaced from 20 to 50 Å. The unfolded trajectories from each size replicate are then used to infer the diffusion coefficient, as described below.

To infer D·OH from simulation, we must account for both transients that affect the observed within-trajectory diffusion rate, and finite size effects that lead to systematic variation across trajectories for systems of differing size. As discussed, e.g., by Bullerjahn, von Bülow, and Hummel,46 classical methods such as regression of the mean squared deviation on time have poor statistical properties, leading to inefficient inference. Relatedly, each trajectory provides relatively little information regarding the target solute (due to its low particle number), compared to the solvent. These issues motivate statistical approaches that are not only more efficient but that pool information across size replicates and between the solvent and solute to improve inference. Here, we use a two-stage Bayesian inference strategy for this purpose, first obtaining local posterior estimates of DMD using a modified Brownian motion process model and then integrating these local estimates via a hierarchical model that combines estimates of DMDOH and DMDH2O across systems of varying size to obtain final estimates of D·OH.

1. Local estimation of DMD

Although the covariance estimator of Bullerjahn, von Bülow, and Hummel46 is both computationally and statistically efficient at high temperatures, it performs less well when the diffusion coefficient becomes small relative to the background noise for which it controls (see discussion in Bullerjahn, von Bülow, and Hummel46). Here, we thus use a strategy of Bayesian estimation for the local (size uncorrected) diffusion constant DMD, which both makes more complete use of data and provides regularization of the resulting estimator. The model employed here is based on the model (somewhat tacitly) underlying the generalized least squares estimators of Bullerjahn, von Bülow, and Hummel,46 namely, a latent Brownian motion process with a Gaussian observation mechanism. Given regularly spaced observations Y = (Y1, Y2, …) at times 1, 2, …, the process may be defined in one dimension by

Xt+1=Xt+Zt,
(5)
Yt+1=Xt+1+Wt,
(6)

where Zt is iid N(0, σ2) and Wt is iid N(0, a2). Physically, X here represents a “true” or idealized Brownian motion with independent perturbations given by Z, while W reflects idiosyncratic noise factors arising from non-Brownian transients. The diffusion constant corresponds to D = σ2/2 (in the squared distance units of Y divided by the time between steps). We may observe that this is a (discretely measured) Gaussian process with the covariance function K(Yi, Yj) =  min(i, j)σ2 + a2I(i = j), and hence, the likelihood is given by (conditioning on and centering the first observation)

p(Y=y|σ2,a2)=MVN(y|0,K(y)),
(7)

where 0 is the 0-vector and K(y) is the Gram matrix of the observed sequence. This is straightforward to work with, although computationally expensive when the number of time points becomes large, and can be pooled across dimensions in the isotropic case (as is done here).

To define priors on the variance parameters, we first observe that on our physical scale of interest (Å2/ps) and over the range of conditions considered here, it is a priori unlikely that σ2 will exceed 0.5; likewise46 it is reasonable to expect a2 to be comparable to (or possibly smaller than) σ2. We thus use independent half-Gaussian priors for σ2 and a2, with a scale of 0.5, which is relatively flat over the region of interest while discouraging strongly unphysical values. (Note that this is equivalent to L2 regularization of the variance parameters.)

For parameter estimation, it is natural here due to both the Gaussian structure of the problem and the large data size to employ maximum a posteriori (MAP) estimation, invoking the Laplace approximation50 to obtain posterior standard deviations. We perform estimation using a custom R51 script, with direct optimization of the log posterior surface using Broyden-Fletcher-Goldfarb-Shanno (BFGS);52 the mclust package53 was used for efficient calculation of the multivariate Gaussian log-likelihood. The posterior standard deviation of D̂ was obtained via the Hessian of the negative log posterior about the posterior mode (exploiting the linear relationship between D and σ2). Due to the cost of computing the Gram matrix, the initial unfolded and drift-corrected trajectories were downsampled from 0.25 to 0.5 ps resolution and split into two segments of 500 ps length (i.e., 1000 observations); these were pooled in the likelihood calculation. In the case of H2O trajectories, trajectories for all water molecules were pooled and jointly analyzed. This process led to estimates of the posterior mode (assumed equal to the mean, under the Laplace approximation) and standard deviation for DMD for both ·OH and H2O under each condition and at each box size. These estimates were then integrated to estimate DH2O and D·OH in each condition as described below.

2. Estimation of D·OH

Estimation of small-molecule diffusion constants is challenging due both to the need to correct for finite size effects and a high level of idiosyncratic variation between trajectories that is difficult to account for; moreover, a single molecule trajectory provides relatively little information per simulation run (as opposed to the large number of solvent molecule trajectories obtained on each run). Here, we address both issues via a hierarchical Bayesian model that pools information between H2O and ·OH trajectories and that incorporates multiple sources of variation. The model (whose structure is described pictorially in the plate diagram of Fig. 2) is defined as follows:

FIG. 2.

Structure of the diffusion constant model. White circles indicate latent quantities, while observed random quantities are indicated by shaded circles; fixed values are uncircled. Elements within the central plate are replicated N times, while those outside are pooled.

FIG. 2.

Structure of the diffusion constant model. White circles indicate latent quantities, while observed random quantities are indicated by shaded circles; fixed values are uncircled. Elements within the central plate are replicated N times, while those outside are pooled.

Close modal

We begin with the observation that if D·OH and DH2O are the respective bulk diffusion constants for ·OH and H2O, then

DOH=DMDOHiα/Li,
(8)
DH2O=DMDH2Oiα/Li,
(9)

where DMDOHi and DMDH2Oi are the diffusion coefficients for a PBC system with length scale Li, and α is the system size scaling coefficient. Of these, only L is observed. We do, however, have local estimates of DMDOHi and DMDH2Oi, which we model as

DMDOHî=DMDOHi+EOH,
(10)
DMDH2Oî=DMDH2Oi+EH2O,
(11)

where EOHN(0,SMDOĤi+τOHi2) and EH2ON(0,SMDH2Ôi+τH2Oi2) represent deviations from the idealized local diffusion coefficients. The error variance is modeled via two components: the posterior variance from the local model of Sec. II C 1 (SMDOĤi,SMDH2Ôi) and the excess variances τOHi2 and τH2Oi2 representing trajectory-specific idiosyncratic deviations not reflected by the within-trajectory estimates. We take the square roots of the excess variances to be generated by 0-truncated normal distributions, i.e., τOHTN0(μOH,γOH2),τH2OTN0(μH2O,γH2O2), with weakly informative standard half-Cauchy priors on the μ* and γ* parameters. We observe that this prior structure can be seen as flexibly generalizing several standard regression-like models: in the limit as E* → 0, we recover a model akin to weighted least squares, with weights based on the locally estimated variances; when E*S*̂i but γ* → 0, recover a model akin to a standard homoskedastic regression; and when γ* ≫ 0, we obtain a robust regression with a heavy-tailed error distribution. Finally, we take α to be a priori uniform on (0, 0.75) (as the 0.75 is expected to be strictly larger than the value of α for the conditions studied here).

Given the above, we perform posterior simulation using the No-U-Turn Hamiltonian Monte Carlo algorithm54 from the Stan library.55,56 Four chains were employed for each condition, with 105 burn-in iterations per chain followed by 105 additional iterations from which 1000 were retained (i.e., a thinning interval of 100) for a final sample size of 4000 draws per condition. Convergence was assessed with R̂.57 Posterior means and 95% posterior intervals were obtained for D·OH and DH2O for each condition for subsequent analysis, as discussed below.

We also examine the radial distribution functions [RDF, g(r)] of ·OH in its interactions with H2O under both normal conditions and the low temperature/high pressure regime. We simulate trajectories of ·OH in TIP4P/2005 water as described in Sec. II B, with the following modifications. Simulations were performed using an ∼50 Å cubic water box, at conditions of either 298 K and 1 atm pressure or 273 K and 10  000 atm pressure. Following initial adjustment, simulations were run for a total of 50 ns, with one frame sampled for analysis every 50 ps (for a total of 1000 frames). When estimating g, we consider distributions of selected atoms at specified distances from a focal atom, collecting the counts of all atoms of the requisite type in bins of width 0.2 Å out to a maximum distance of 10 Å.

To examine ·OH/H2O interactions, we estimate gOH,H2O for H–H, H–O, O–H, and O–O pairs (with distances denoted, respectively, rH,H, rH,O, rO,H, or rO,O); for purposes of comparison, we likewise estimate gH2O,H2O under equivalent conditions. For each trajectory and atom type pair, we obtain the counts of atoms of the selected type within each distance interval of the focal atom. Let C be the vector of such counts, with Ci being the total count in the ith distance bin and Cij being the count in the ith bin arising from frame j. It is natural to approximate the counts within each bin at each frame by

CijPois(λi),

where λi is the per-frame expected count for bin i. Here, we employ a Jeffreys prior for λi [i.e., p(λi) ∝ 1/λ2], taking λia priori iid and approximating counts across bins as conditionally independent given λ. Having observed total count ci from n frames, our posterior is then

λi|ciGamma(ci+1/2,1/n),

with posterior mean Eλi|ci = (ci + 0.5)/n and variance Var(λi|ci) = (ci + 0.5)/n2. (Posterior intervals can likewise be obtained directly from the quantiles of the Gamma distribution.) Finally, to obtain our estimate of g within the interval of bin i, we normalize the expected atom count by the corresponding expected count for bin i under a uniform distribution of atom locations. This process is repeated for each bin and each pair type to obtain each respective estimated g function. Finally, to facilitate visualization, we apply spline interpolation to the local posterior mean estimates for each bin, centering the estimate at the bin midpoint.

Figure 3 shows the estimated bulk diffusion coefficients for ·OH and H2O at the simulated temperatures and pressures. Although D·OH and DH2O approximately coincide under ambient conditions, we observe some differences in their response to temperature and pressure. While TIP4P/2005 reproduces the experimentally observed increase in H2O diffusion rate under moderate pressure, ·OH diffusion slows with pressure at moderate temperatures (eventually gaining very modest acceleration in the low temperature regime). At pressures approaching 10000 atm, both water and the hydroxyl radical slow considerably (as expected from the well-known increase in water viscosity in this regime), with D·OH diffusing more slowly than water despite its smaller size.

FIG. 3.

Posterior mean estimates and 95% posterior intervals for DH2O and D·OH. (In some cases, the interval width is narrower than point sizes.) Both behave similarly at low and very pressures, but D·OH lacks the accelerated diffusion seen in H2O at moderately high pressure.

FIG. 3.

Posterior mean estimates and 95% posterior intervals for DH2O and D·OH. (In some cases, the interval width is narrower than point sizes.) Both behave similarly at low and very pressures, but D·OH lacks the accelerated diffusion seen in H2O at moderately high pressure.

Close modal

To better summarize the posterior surface of D·OH in response to temperature and pressure and to provide a simplified approximation to the diffusion coefficient that can be used in other applications, we fit a least-squares log–log polynomial approximation to the posterior means. The resulting summary model [R2 = 0.995, standard error 0.032 on 13 residual degrees of freedom (log scale), root-mean-square error (RMSE) 0.0035 Å2 ps−1 on the phenomenological scale] is given by

DOĤ=exp4.375997×102+1.464854×102T+8.101061×101P+1.227679×101T2+2.817980×102P2+1.560849×101TP+2.622117×103P3Å2ps1
(12)

where T′ and P′ are the dimensionless quantities T=logTK1 and P=logPatm1. The resulting surface is shown in Fig. 4; the equivalent surface for DH2O is shown in Fig. S5. As expected from the known coincidence of D·OH and DH2O under ambient conditions, D·OH has broadly similar behavior to DH2O. However, ·OH shows less accelerated diffusion at moderate pressures and slows considerably at very high pressures (even in excess of H2O). We show this quantitatively in Fig. 5, which displays the ratio of estimated D·OH to be DH2O for the TIP4P/2005 system. D·OH and DH2O remain within a few percent of each other for pressures below ∼100 atm, with deviations being governed primarily by temperature. As pressure increases, however, temperature has less relative impact on the ratio of diffusion coefficients, with pressure having the dominant role.

FIG. 4.

DOĤ as a function of temperature and pressure, based on posterior inference from simulations in TIP4P/2005 water. ·OH diffusion is slightly faster at moderate pressures when temperatures are low (although less so than H2O) but slows considerably at high pressures.

FIG. 4.

DOĤ as a function of temperature and pressure, based on posterior inference from simulations in TIP4P/2005 water. ·OH diffusion is slightly faster at moderate pressures when temperatures are low (although less so than H2O) but slows considerably at high pressures.

Close modal
FIG. 5.

DOĤ/DH2Ô as a function of temperature and pressure based on posterior inference from simulations in TIP4P/2005 water. Diffusion rates are nearly equal at low pressure and either high or low temperature, with D·OH declining in relative terms as pressure increases.

FIG. 5.

DOĤ/DH2Ô as a function of temperature and pressure based on posterior inference from simulations in TIP4P/2005 water. Diffusion rates are nearly equal at low pressure and either high or low temperature, with D·OH declining in relative terms as pressure increases.

Close modal

It may be conjectured that the difference in the high-pressure behavior of ·OH vs its behavior under ambient conditions may be related to the observation of an MD study in BJH water at relatively high temperature (310 K) and 1 atm that ·OH tends to occupy cavities in the H2O hydrogen bonding network.36 These transient regions of empty space are not necessarily simple spheres, but can take a variety of shapes,58 and are balanced by regions of higher density liquid with more water molecules in tetrahedral arrangements.59 At high pressure, occupancy of such cavities would be expected to be highly favorable, and to retard diffusion. To investigate this possibility, we consider the radial distribution functions for ·OH/water interactions, as shown in Fig. 6. Immediately striking is the fact that there is very little change in the RDFs for interactions of either ·OH atom with solvent (top four panels), even across four orders of magnitude increase in pressure (and a temperature reduction of 25 K). This is in contrast to H2O, which shows marked shifts in the high-pressure regime (particularly a flattening of the O–O distribution that is compatible with compressive disruption of the hydrogen bond network). (We note in passing that the H2O O–O RDF is compatible with a recent study of TIP4P/2005 water by Camisasca et al.;60 they do not see an equivalent collapse of the 4.5 Å peak at low temperature but consider only pressures of 1 atm. This suggests that the effect is due to pressure rather than temperature.) The lack of change in local ·OH interactions suggests that enhanced occupancy of cavities at high pressure is not likely to explain the observed dynamics. However, a simple alternative explanation is that the flattening of the water O–O RDF indicates a sharp reduction in density of available cavities to occupy, leading ·OH to experience longer time between transitions. This would slow ·OH diffusion, without implying substantial changes in its local environment. We hypothesize that the interplay between temperature and pressure effects may impact the reactivity of ·OH in unexpected ways. At low temperature, the increased viscosity lengthens the lifetime of cavities formed by defects in the hydrogen bond network, which may impact the photochemistry of the initial reaction by trapping otherwise unstable excited states. A similar phenomenon is observed in the anomalous fluorescence signal measured for concentrated solutions of NaOH and HCl, which has been attributed to the stabilization of hydrated ion species within transient solvent cavities.61 High pressure is expected to work against this mechanism by collapsing the cavities; the nature of the excited states that may be stabilized and the pressure/temperature regime over which this occurs is likely to be fruitful topic for further investigation.

FIG. 6.

Posterior mean estimates of RDFs for interactions between ·OH and H2O (top panels) and self-interactions within H2O (bottom panels) by atom pair type. (Raw estimates have been smoothed with interpolating splines; 95% posterior intervals shown as shaded regions but are generally thinner than line widths.) Red lines indicate ambient (298 K, 1 atm) conditions, while blue lines indicate low temperature/high pressure scenario (273 K, 10000 atm). ·OH interactions change little with temperature and pressure, while water interactions vary markedly.

FIG. 6.

Posterior mean estimates of RDFs for interactions between ·OH and H2O (top panels) and self-interactions within H2O (bottom panels) by atom pair type. (Raw estimates have been smoothed with interpolating splines; 95% posterior intervals shown as shaded regions but are generally thinner than line widths.) Red lines indicate ambient (298 K, 1 atm) conditions, while blue lines indicate low temperature/high pressure scenario (273 K, 10000 atm). ·OH interactions change little with temperature and pressure, while water interactions vary markedly.

Close modal

Diffusion rates for the hydroxyl radical differ substantially at low temperatures and high pressures. For the most part, the pattern of variation in D·OH qualitatively follows that of DH2O (as would be expected for a small molecule in aqueous solution), but we do not observe the same extent of enhanced diffusion at moderate pressures (particularly at higher temperatures), and diffusion slows more substantially as pressures approach 10000 atm. This may have an impact on the efficiency of radical-scavenging mechanisms in high-pressure terrestrial environments and hypothetically in the even higher pressure conditions predicted to obtain within some outer solar system oceans.

To estimate the ·OH diffusion coefficient, we combined atomistic molecular dynamics with hierarchical Bayesian inference that allows us to easily pool information between solvent and radical and across multiple trajectories of varying size. The modular structure of Bayesian models lends itself naturally to this application, as does the ease with which one can, e.g., leverage a priori information on model parameters or account for variance across trajectories. This strategy (which builds on other recent work on statistical inference for the diffusion coefficient46) is readily adapted to the estimation of diffusion coefficients for other systems or to other related physical quantities.

Additional figures relating to the reported work are available as the supplementary material.

This work was supported by NASA Award No. 80NSSC20K0620 to R.W.M. and C.T.B. We thank George Miller and Vladimir Mandelshtam for helpful discussions about the chemistry of water radiolysis and nuclear quantum effects in water, respectively.

The authors have no conflicts to disclose.

Data sharing is not applicable to this article as no new data were created or analyzed in this study. CHARMM-compatible parameters for the ·OH model are included in the main text.

1.
Y.
Blanco
,
G.
de Diego-Castilla
,
D.
Viúdez-Moreiras
,
E.
Cavalcante-Silva
,
J. A.
Rodríguez-Manfredi
,
A. F.
Davila
,
C. P.
McKay
, and
V.
Parro
, “
Effects of gamma and electron radiation on the structural integrity of organic molecules and macromolecular biomarkers measured by microarray immunoassays and their astrobiological implications
,”
Astrobiology
18
,
1497
(
2018
).
2.
T. C.
Onstott
,
D. P.
Moser
,
S. M.
Pfiffner
,
J. K.
Fredrickson
,
F. J.
Brockman
,
T. J.
Phelps
,
D. C.
White
,
A.
Peacock
,
D.
Balkwill
,
R.
Hoover
,
L. R.
Krumholz
,
M.
Borscik
,
T. L.
Kieft
, and
R.
Wilson
, “
Indigenous and contaminant microbes in ultradeep mines
,”
Environ. Microbiol.
5
,
1168
1191
(
2003
).
3.
J. D.
Tarnas
,
J. F.
Mustard
,
B.
Sherwood Lollar
,
M. S.
Bramble
,
K. M.
Cannon
,
A. M.
Palumbo
, and
A.-C.
Plesa
, “
Radiolytic H2 production on Noachian Mars: Implications for habitability and atmospheric warming
,”
Earth Planet. Sci. Lett.
502
,
133
145
(
2018
).
4.
L.
Ojha
,
S.
Karunatillake
,
S.
Karimi
, and
J.
Buffo
, “
Amagmatic hydrothermal systems on Mars from radiogenic heat
,”
Nat. Commun.
12
,
1754
(
2021
).
5.
J. D.
Tarnas
,
J. F.
Mustard
,
B.
Sherwood Lollar
,
V.
Stamenković
,
K. M.
Cannon
,
J.-P.
Lorand
,
T. C.
Onstott
,
J. R.
Michalski
,
O.
Warr
,
A. M.
Palumbo
, and
A.-C.
Plesa
, “
Earth-like habitable environments in the subsurface of Mars
,”
Astrobiology
21
,
741
(
2021
).
6.
O.
White
,
J. A.
Eisen
,
J. F.
Heidelberg
,
E. K.
Hickey
,
J. D.
Peterson
,
R. J.
Dodson
,
D. H.
Haft
,
M. L.
Gwinn
,
W. C.
Nelson
,
D. L.
Richardson
,
K. S.
Moffat
,
H.
Qin
,
L.
Jiang
,
W.
Pamphile
,
M.
Crosby
,
M.
Shen
,
J. J.
Vamathevan
,
P.
Lam
,
L.
McDonald
,
T.
Utterback
,
C.
Zalewski
,
K. S.
Makarova
,
L.
Aravind
,
M. J.
Daly
,
K. W.
Minton
,
R. D.
Fleischmann
,
K. A.
Ketchum
,
K. E.
Nelson
,
S.
Salzberg
,
H. O.
Smith
,
J. C.
Venter
, and
C. M.
Fraser
, “
Genome sequence of the radioresistant bacterium Deinococcus radiodurans R1
,”
Science
286
,
1571
1577
(
1999
).
7.
R.
Cavicchioli
, “
Extremophiles and the search for extraterrestrial life
,”
Astrobiology
2
,
281
292
(
2002
).
8.
A. C.
Munteanu
,
V.
Uivarosi
, and
A.
Andries
, “
Recent progress in understanding the molecular mechanisms of radioresistance in Deinococcus bacteria
,”
Extremophiles
19
,
707
719
(
2015
).
9.
J. A.
LaVerne
, “
OH radicals and oxidizing products in the gamma radiolysis of water
,”
Radiat. Res.
153
,
196
200
(
2000
).
10.
M. S.
Matheson
, “
The formation and detection of intermediates in water radiolysis
,”
Radiat. Res., Suppl.
4
,
1
23
(
1964
).
11.
D.
Ghosal
,
M. V.
Omelchenko
,
E. K.
Gaidamakova
,
V. Y.
Matrosova
,
A.
Vasilenko
,
A.
Venkateswaran
,
M.
Zhai
,
H. M.
Kostandarithes
,
H.
Brim
,
K. S.
Makarova
,
L. P.
Wackett
,
J. K.
Fredrickson
, and
M. J.
Daly
, “
How radiation kills cells: Survival of Deinococcus radiodurans and Shewanella oneidensis under oxidative stress
,”
FEMS Microbiol. Rev.
29
,
361
375
(
2005
).
12.
R.
Roots
and
S.
Okada
, “
Estimation of life times and diffusion distances of radicals involved in x-ray-induced DNA strand breaks or killing of mammalian cells
,”
Radiat. Res.
64
,
306
320
(
1975
).
13.
J.
Du
and
J. M.
Gebicki
, “
Proteins are major initial cell targets of hydroxyl free radicals
,”
Int. J. Biochem. Cell Biol.
36
,
2334
2343
(
2004
).
14.
M. J.
Davies
,
S.
Fu
, and
R. T.
Dean
, “
Protein hydroperoxides can give rise to reactive free radicals
,”
Biochem. J.
305
,
643
649
(
1995
).
15.
K. J.
Barnham
,
C. L.
Masters
, and
A. I.
Bush
, “
Neurodegenerative diseases and oxidative stress
,”
Nat. Rev. Drug Discovery
3
,
205
214
(
2004
).
16.
H. A.
Schwarz
, “
Applications of the spur diffusion model to the radiation chemistry of aqueous solutions
,”
J. Phys. Chem.
73
,
1928
1937
(
1969
).
17.
M.
Huerta Parajon
,
P.
Rajesh
,
T.
Mu
,
S. M.
Pimblott
, and
J. A.
LaVerne
, “
H atom yields in the radiolysis of water
,”
Radiat. Phys. Chem.
77
,
1203
1207
(
2008
).
18.
B. G.
Ershov
and
A. V.
Gordeev
, “
A model for radiolysis of water and aqueous solutions of H2, H2O2 and O2
,”
Radiat. Phys. Chem.
77
,
928
935
(
2008
).
19.
K. A.
Omar
,
K.
Hasnaoui
, and
A.
de la Lande
, “
First-principles simulations of biological molecules subjected to ionizing radiation
,”
Annu. Rev. Phys. Chem.
72
,
445
465
(
2021
).
20.
Z.-H.
Loh
,
G.
Doumy
,
C.
Arnold
,
L.
Kjellsson
,
S. H.
Southworth
,
A.
Al Haddad
,
Y.
Kumagai
,
M.-F.
Tu
,
P. J.
Ho
,
A. M.
March
,
R. D.
Schaller
,
M. S.
Bin Mohd Yusof
,
T.
Debnath
,
M.
Simon
,
R.
Welsch
,
L.
Inhester
,
K.
Khalili
,
K.
Nanda
,
A. I.
Krylov
,
S.
Moeller
,
G.
Coslovich
,
J.
Koralek
,
M. P.
Minitti
,
W. F.
Schlotter
,
J.-E.
Rubensson
,
R.
Santra
, and
L.
Young
, “
Observation of the fastest chemical processes in the radiolysis of water
,”
Science
367
,
179
182
(
2020
).
21.
I.
Jordan
,
M.
Huppert
,
D.
Rattenbacher
,
M.
Peper
,
D.
Jelovina
,
C.
Perry
,
A.
von Conta
,
A.
Schild
, and
H. J.
Wörner
, “
Attosecond spectroscopy of liquid water
,”
Science
369
,
974
979
(
2020
).
22.
S.
Thürmer
,
M.
Ončák
,
N.
Ottosson
,
R.
Seidel
,
U.
Hergenhahn
,
S. E.
Bradforth
,
P.
Slavíček
, and
B.
Winter
, “
On the nature and origin of dicationic, charge-separated species formed in liquid water on x-ray irradiation
,”
Nat. Chem.
5
,
590
596
(
2013
).
23.
X.
Ren
,
E.
Wang
,
A. D.
Skitnevskaya
,
A. B.
Trofimov
,
K.
Gokhberg
, and
A.
Dorn
, “
Experimental evidence for ultrafast intermolecular relaxation processes in hydrated biomolecules
,”
Nat. Phys.
14
,
1062
1066
(
2018
).
24.
B.
Kirchner
,
J.
Stubbs
, and
D.
Marx
, “
Fast anomalous diffusion of small hydrophobic species in water
,”
Phys. Rev. Lett.
89
,
215901
(
2002
).
25.
S. T.
Roberts
,
P. B.
Petersen
,
K.
Ramasesha
,
A.
Tokmakoff
,
I. S.
Ufimtsev
, and
T. J.
Martinez
, “
Observation of a Zundel-like transition state during proton transfer in aqueous hydroxide solutions
,”
Proc. Natl. Acad. Sci. U. S. A.
106
,
15154
15159
(
2009
).
26.
D.
Marx
,
A.
Chandra
, and
M. E.
Tuckerman
, “
Aqueous basic solutions: Hydroxide solvation, structural diffusion, and comparison to the hydrated proton
,”
Chem. Rev.
110
,
2174
2216
(
2010
).
27.
I. M.
Svishchev
and
A. Y.
Plugatyr
, “
Hydroxyl radical in aqueous solution: Computer simulation
,”
J. Phys. Chem. B
109
,
4123
4128
(
2005
).
28.
C.
Knight
and
G. A.
Voth
, “
The curious case of the hydrated proton
,”
Acc. Chem. Res.
45
,
101
109
(
2012
).
29.
M. E.
Tuckerman
,
D.
Marx
, and
M.
Parrinello
, “
The nature and transport mechanism of hydrated hydroxide ions in aqueous solution
,”
Nature
417
,
925
9299
(
2002
).
30.
M.
Chen
,
L.
Zheng
,
B.
Santra
,
H.-Y.
Ko
,
R. A.
DiStasio
, Jr.
,
M. L.
Klein
,
R.
Car
, and
X.
Wu
, “
Hydroxide diffuses slower than hydronium in water because its solvated structure inhibits correlated proton transfer
,”
Nat. Chem.
10
,
413
419
(
2018
).
31.
M.
Ceriotti
,
W.
Fang
,
P. G.
Kusalik
,
R. H.
McKenzie
,
A.
Michaelides
,
M. A.
Morales
, and
T. E.
Markland
, “
Nuclear quantum effects in water and aqueous systems: Experiment, theory, and current challenges
,”
Chem. Rev.
116
,
7529
7550
(
2016
).
32.
L. M.
Dorfman
and
G. E.
Adams
,
Reactivity of the Hydroxyl Radical in Aqueous Solutions
(
U.S. Department of Commerce, National Bureau of Standards
,
Washington, DC
,
1973
).
33.
R. B.
Best
,
X.
Zhu
,
J.
Shim
,
P. E. M.
Lopes
,
J.
Mittal
,
M.
Feig
, and
J. A. D.
Mackerell
, “
Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone ϕ, ψ and side-chain χ(1) and χ(2) dihedral angles
,”
J. Chem. Theory Comput.
8
,
3257
3273
(
2012
).
34.
J. L. F.
Absacal
and
C.
Vega
, “
A general purpose model for the condensed phases of water: TIP4P/2005
,”
J. Chem. Phys.
123
,
234505
(
2005
).
35.
I. N.
Tsimpanogiannis
,
O. A.
Moultos
,
L. F.
Franco
,
M. B. M.
Spera
,
M.
Erdös
, and
I. G.
Economou
, “
Self-diffusion coefficient of bulk and confined water: A critical review of classical molecular simulation studies
,”
Mol. Simul.
45
,
425
453
(
2019
).
36.
A.
Pabis
,
J.
Szala-Bilnik
, and
D.
Swiatla-Wojcik
, “
Molecular dynamics study of the hydration of the hydroxyl radical at body temperature
,”
Phys. Chem. Chem. Phys.
13
,
9458
9468
(
2011
).
37.
P.
Bopp
,
G.
Jancsó
, and
K.
Heinzinger
, “
An improved potential for non-rigid water molecules in the liquid phase
,”
Chem. Phys. Lett.
98
,
129
133
(
1983
).
38.
M. H.
Kalos
and
P. A.
Whitlock
,
Monte Carlo Methods, Volume I: Basics
(
John Wiley and Sons
,
New York, NY
,
1986
).
39.
G. J.
Martyna
,
D. J.
Tobias
, and
M. L.
Klein
, “
Constant pressure molecular dynamics algorithms
,”
J. Chem. Phys.
101
,
4177
4189
(
1994
).
40.
S. E.
Feller
,
Y.
Zhang
,
R. W.
Pastor
, and
B. R.
Brooks
, “
Constant pressure molecular dynamics simulation: The Langevin piston method
,”
J. Chem. Phys.
103
,
4613
4621
(
1995
).
41.
J. C.
Phillips
,
R.
Braun
,
W.
Wang
,
J.
Gumbart
,
E.
Tajkhorshid
,
E.
Villa
,
C.
Chipot
,
R. D.
Skeel
,
L.
Kalé
, and
K.
Schulten
, “
Scalable molecular dynamics with NAMD
,”
J. Comput. Chem.
26
,
1781
1802
(
2005
).
42.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
, “
VMD: Visual molecular dynamics
,”
J. Mol. Graphics
14
(
33–38
),
27
28
(
1996
).
43.
J. V.
Ribeiro
,
B.
Radak
,
J.
Stone
,
J.
Gullingsrud
,
J.
Saam
, and
J.
Phillips
, “psfgen Plugin for VMD,” Software File (2020).
44.
L.
Martínez
,
R.
Andrade
,
E.
Birgin
, and
J. M.
Martínez
, “
PACKMOL: A package for building initial configurations for molecular dynamics simulations
,”
J. Comput. Chem.
30
,
2157
2164
(
2009
).
45.
S.
von Bülow
,
J. T.
Bullerjahn
, and
G.
Hummel
, “
Systematic errors in diffusion coefficients from long-time molecular dynamics simulations at constant pressure
,”
J. Chem. Phys.
153
,
021101
(
2020
).
46.
J. T.
Bullerjahn
,
S.
von Bülow
, and
G.
Hummel
, “
Optimal estimates of self-diffusion coefficients from molecular dynamics simulations
,”
J. Chem. Phys.
153
,
024116
(
2020
).
47.
I.-C.
Yeh
and
G.
Hummer
, “
System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions
,”
J. Phys. Chem. B
108
,
15873
(
2004
).
48.
M. A.
González
and
J. L. F.
Abascal
, “
The shear viscosity of rigid water models
,”
J. Chem. Phys.
132
,
096101
(
2010
).
49.
S.
Du
and
J. S.
Francisco
, “
The OH radical-H2O molecular interaction potential
,”
J. Chem. Phys.
124
,
224318
(
2006
).
50.
A.
Gelman
,
J. B.
Carlin
,
H. S.
Stern
, and
D. B.
Rubin
,
Bayesian Data Analysis
, 2nd ed. (
Chapman and Hall
,
London
,
2003
).
51.
R Core Team
,
R: A Language and Environment for Statistical Computing
(
R Foundation for Statistical Computing
,
Vienna, Austria
,
2021
).
52.
J. C.
Nash
,
Compact Numerical Methods for Computers. Linear Algebra and Function Minimization
(
Adam Hilger
,
1990
).
53.
L.
Scrucca
,
M.
Fop
,
T. B.
Murphy
, and
A. E.
Raftery
, “
mclust 5: Clustering, classification and density estimation using Gaussian finite mixture models
,”
R J.
8
,
289
317
(
2016
).
54.
M. D.
Homan
and
A.
Gelman
, “
The No-U-Turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo
,”
J. Mach. Learn. Res.
15
,
1593
1623
(
2014
); available at http://jmlr.org/papers/v15/hoffman14a.html.
55.
Stan Development Team
, “RStan: The R Interface to Stan” (2020), R package version 2.21.2.
56.
Stan Development Team
, “Stan: A C++ Library for Probability and Sampling” (2020), Software Library.
57.
A.
Gelman
and
D. B.
Rubin
, “
Inference from iterative simulation using multiple sequences (with discussion)
,”
Stat. Sci.
7
,
457
511
(
1992
).
58.
G. C.
Sosso
,
S.
Caravati
,
G.
Rotskoff
,
S.
Vaikuntanathan
, and
A.
Hassanali
, “
On the role of nonspherical cavities in short length-scale density fluctuations in water
,”
J. Phys. Chem. A
121
,
370
380
(
2017
).
59.
N.
Ansari
,
R.
Dandekar
,
S.
Caravati
,
G. C.
Sosso
, and
A.
Hassanali
, “
High and low density patches in simulated liquid water
,”
J. Chem. Phys.
149
,
204507
(
2018
).
60.
G.
Camisasca
,
H.
Pathak
,
K. T.
Wikfeldt
, and
L. G. M.
Pettersson
, “
Radial distribution functions of water: Models vs experiments
,”
J. Chem. Phys.
151
,
044502
(
2019
).
61.
A. M.
Villa
,
S. M.
Doglia
,
L. D.
Gioia
,
L.
Bertini
, and
A.
Natalello
, “
Anomalous intrinsic fluorescence of HCl and NaOH aqueous solutions
,”
J. Phys. Chem. Lett.
10
,
7230
7236
(
2019
).

Supplementary Material