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.

## I. INTRODUCTION

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 chemotrophy^{2} 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,

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

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,

These reactions serve as the starting point for the formation of a complex mixture of reactive oxygen species^{18} and/or reactions of these primary products with biomolecules.^{19} The recent advent of experimental techniques enabling direct detection of the long-hypothesized H_{2}O^{+} intermediate^{20} and observation of attosecond dynamics in liquid water^{21} as well as the discovery of new reaction pathways^{22,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.

## II. METHODS

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.

### A. Parameterization of the ·OH model

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 H_{2}O, an extremely unfavorable reaction under normal conditions. Although more subtle quantum effects (like those hypothesized to lead to enhanced diffusion of OH^{−}^{29–31}) are possible, the experimental observation that the diffusion coefficient of ·OH under ambient conditions is close to that of water^{32} 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-based^{33} 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

where *R*_{AB} is the distance between atoms *A* and *B*, *ɛ*_{A} and *ɛ*_{B} are species-specific well-depth parameters, and *r*_{min A}/2 and *r*_{min 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}, *r*_{min O}/2, and *r*_{min 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, *r*_{min O}/2 ∈ (0, 3) Å, and *r*_{min H}/2 ∈ (0, 2) Å using a four-dimensional Halton sequence (bases 2, 3, 5, and 7). (Halton sequences are low-discrepancy sequences^{38} 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 $\u224820$ Å 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 Hummel^{45} 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 Hummel^{46} 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 H_{2}O under these simulation conditions. These raw estimates were then corrected for finite sample sizes using the analytical correction factor of Yeh and Hummer,^{47}

where *D*_{∞} is the diffusion constant at infinite size, *D*_{MD} 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/m^{3} from González and Abascal^{48} 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.23\u2212D\u22c5OH\u0302|/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 H_{2}O 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, *r*_{min O}/2 ≈ 0.802, and *r*_{min 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, *r*_{min O}/2 ∈ (0.6, 1.1) Å, and *r*_{min 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–H_{2}O orientation.)

The final parameters for the ·OH model (including both predetermined and calibrated parameters) are *r*_{OH} = 0.9751 Å, *E*_{OH} = 475.6 kcal/mol/Å^{2}, *m*_{O} = 15.994 Da, *m*_{O} = 1.008 Da, *ɛ*_{O} = −0.235 kcal/mol, *ɛ*_{H} = −0.306 8 kcal/mol, *r*_{min O}/2 = 0.781 85 Å, and *r*_{min H}/2 = 0.746 939 Å. The relatively small value of *r*_{min 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 Francisco^{49} based on DFT-based potential energy scanning calculations (with both differing from^{36}), 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 $D\u22c5OH\u0302\u22480.222$ Å^{2}/ps for an ARE of ∼4% vs the target value of 0.23 Å^{2}/ps.

### B. Diffusion simulation

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.

### C. Inference for the diffusion coefficient

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 *D*_{MD} using a modified Brownian motion process model and then integrating these local estimates via a hierarchical model that combines estimates of $DMD\u22c5OH$ and $DMDH2O$ across systems of varying size to obtain final estimates of *D*_{·OH}.

#### 1. Local estimation of *D*_{MD}

Although the covariance estimator of Bullerjahn, von Bülow, and Hummel^{46} 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 Hummel^{46}). Here, we thus use a strategy of Bayesian estimation for the local (size uncorrected) diffusion constant *D*_{MD}, 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* = (*Y*_{1}, *Y*_{2}, …) at times 1, 2, …, the process may be defined in one dimension by

where *Z*_{t} is iid *N*(0, *σ*^{2}) and *W*_{t} is iid *N*(0, *a*^{2}). 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*(*Y*_{i}, *Y*_{j}) = min(*i*, *j*)*σ*^{2} + *a*^{2}*I*(*i* = *j*), and hence, the likelihood is given by (conditioning on and centering the first observation)

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; likewise^{46} it is reasonable to expect *a*^{2} to be comparable to (or possibly smaller than) *σ*^{2}. We thus use independent half-Gaussian priors for *σ*^{2} and *a*^{2}, 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 *L*_{2} 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 approximation^{50} to obtain posterior standard deviations. We perform estimation using a custom R^{51} script, with direct optimization of the log posterior surface using Broyden-Fletcher-Goldfarb-Shanno (BFGS);^{52} the mclust package^{53} was used for efficient calculation of the multivariate Gaussian log-likelihood. The posterior standard deviation of $D\u0302$ 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 H_{2}O 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 *D*_{MD} for both ·OH and H_{2}O 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 H_{2}O 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:

We begin with the observation that if *D*_{·OH} and $DH2O$ are the respective bulk diffusion constants for ·OH and H_{2}O, then

where $DMD\u22c5OHi$ and $DMDH2Oi$ are the diffusion coefficients for a PBC system with length scale *L*_{i}, and *α* is the system size scaling coefficient. Of these, only *L* is observed. We do, however, have local estimates of $DMD\u22c5OHi$ and $DMDH2Oi$, which we model as

where $E\u22c5OH\u223cN(0,SMD\u22c5OH\u0302i+\tau \u22c5OHi2)$ and $EH2O\u223cN(0,SMDH2O\u0302i+\tau 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 ($SMD\u22c5OH\u0302i,SMDH2O\u0302i$) and the excess variances $\tau \u22c5OHi2$ and $\tau 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., $\tau \u22c5OH\u223cTN0(\mu \u22c5OH,\gamma \u22c5OH2),\tau H2O\u223cTN0(\mu H2O,\gamma 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*\u226bS*\u0302i$ 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 algorithm^{54} from the Stan library.^{55,56} Four chains were employed for each condition, with 10^{5} burn-in iterations per chain followed by 10^{5} 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\u0302$.^{57} Posterior means and 95% posterior intervals were obtained for *D*_{·OH} and $DH2O$ for each condition for subsequent analysis, as discussed below.

### D. Simulation and inference for radial distribution functions

We also examine the radial distribution functions [RDF, *g*(*r*)] of ·OH in its interactions with H_{2}O 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/H_{2}O interactions, we estimate $g\u22c5OH,H2O$ for H–H, H–O, O–H, and O–O pairs (with distances denoted, respectively, *r*_{H,H}, *r*_{H,O}, *r*_{O,H}, or *r*_{O,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 *C*_{i} being the total count in the *i*th distance bin and *C*_{ij} being the count in the *i*th bin arising from frame *j*. It is natural to approximate the counts within each bin at each frame by

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 *λ*_{i}*a priori* iid and approximating counts across bins as conditionally independent given *λ*. Having observed total count *c*_{i} from *n* frames, our posterior is then

with posterior mean **E***λ*_{i}|*c*_{i} = (*c*_{i} + 0.5)/*n* and variance Var(*λ*_{i}|*c*_{i}) = (*c*_{i} + 0.5)/*n*^{2}. (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.

## III. RESULTS AND DISCUSSION

Figure 3 shows the estimated bulk diffusion coefficients for ·OH and H_{2}O 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 H_{2}O 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.

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 [*R*^{2} = 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

where *T*′ and *P*′ are the dimensionless quantities $T\u2032=logTK\u22121$ and $P\u2032=logPatm\u22121$. 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 H_{2}O). 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.

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 H_{2}O 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 H_{2}O, 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 H_{2}O 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 $\u22484.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.

## IV. CONCLUSION

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 coefficient^{46}) is readily adapted to the estimation of diffusion coefficients for other systems or to other related physical quantities.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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.

## REFERENCES

_{2}production on Noachian Mars: Implications for habitability and atmospheric warming

_{2}, H

_{2}O

_{2}and O

_{2}

_{2}O molecular interaction potential