Voronoi polyhedra probing of hydrated OH radical †

We present 3D-characterization of the solvated c OH radical in water at 37 (cid:1) C using the Voronoi polyhedron technique. Assuming that the c OH solvation cage is represented by the minimal-volume polyhedron we have calculated statistical distributions of the metric and topological properties of the solvation cage. The statistical means 10.1, 3.26 ˚ A, 27.2 ˚ A3, and 1.50, have been obtained for the solvation number, the face-weighted cage-radius, the free volume of c OH aq , and the asphericity factor, respectively. The mean value of the face-weighted radius coincides with the maximum of the radical-water radial distribution function. Similar properties of the polyhedra constructed for the hydrogen bonded radical and for the unbound one con ﬁ rm the mechanistic view of the localization of c OH in cavities in the hydrogen-bond network. The distribution of the asphericity factor reveals the existence of ice-like regions. A short review of the graphical software used for 3D-visualization includes MS Excel, Maple, and POV-Ray programs. A practical user guide is provided in the ESI. †


Introduction
Insight into a solvation cage (shell) of a solute provided by classical molecular dynamics (MD) and Monte Carlo (MC) simulations or by DFT-based Car-Parrinello Molecular Dynamics (CPMD) requires computational tools for visualization and statistical description of the simulated structures. 1,2he simplest and commonly used tool is a partial radial distribution function (RDF).RDF describes site-site solutesolvent correlation in space. 1 A serious drawback resulting from the one-dimensional nature of a RDF is that 3D correlations are averaged out.Although more thorough insight can be provided by using angular-dependent pair distribution or spatial distribution functions, 3,4 construction of the Voronoi polyhedron (VP) offers a unique qualitative tool for 3D-characteristics of the nearest neighbourhood of a solute. 5The Voronoi polyhedron associated with a given point (atom) is dened as the convex region of space closer to this point than to any other points (atoms).Mathematically VP follows a concept embodied in a Wigner-Seitz primitive cell in crystallography. 6Metric and topological properties of VP quantify the space owned by a solute in a solution.These properties may be used to describe size and shape of a solvation cage (shell) providing estimates for the free volume and the number of neighbouring solvent molecules (the solvation number).
In the present paper we employ a method of VP construction to provide the rst comprehensive statistical description of the metric and topological properties of the hydrated hydroxyl radical (cOH aq ).The problem of cOH solvation in aqueous media is important because of its relevance to atmospheric chemistry, biology, medicine, radiobiology, and radiation chemistry.Hydroxyl radicals are highly reactive and consequently short-lived.In biological systems, cOH radicals produced from the decomposition of hydro-peroxides, cause much cell damage. 7,8Hydroxyl radicals generated in the troposphere remove volatile organic compounds (VOCs) and methane from the air. 9In aqueous systems exposed to ionizing radiation cOH is the main oxidant formed. 10Taking into account a key role of the hydroxyl radical in biology and medicine we have recently carried out classical MD simulation for a diluted aqueous solution of cOH at the body temperature (37 C). 11 Using partial RDFs calculated for the solute-solvent and solvent-solvent sites we showed that cOH radical is coordinated by 13-14 solvent molecules and tends to occupy cavities in the hydrogen-bond (HB) network of water.Liquid water at 37 C resembles a gel-like HB network, where rotations cause individual hydrogen bonds to break and quickly re-form in new congurations. 12The localization in cavities is thus consistent with a small number of hydrogen bonds established by cOH.According to the geometrical denition of hydrogen bond, cOH radical was not hydrogenbonded in ca.30% of cavities. 11In other cases, cOH formed mostly one hydrogen-bond to the surrounding water molecules, usually acting as a proton-donating partner (H-donor).Compared to neat water the continuous lifetime of H-donor bond (0.033 ps) was an order of magnitude smaller, but the intermittent lifetime of a few picoseconds was similar. 13Our later study showed that although the mean number of waterwater hydrogen bonds is the same in solution and in neat water, the HB connectivity pattern of solvent molecules is different. 14Namely, in the presence of cOH ice-like patches, i.e. supramolecular structures of continuously connected four-bonded molecules, are noticeably smaller.
To develop our mechanistic view on localization of cOH in aqueous media we have elaborated a VP-based approach and processed a molecular dynamics simulation trajectory to calculate metric and topological properties of the molecular neighbourhood of cOH aq .At the same time properties of the constructed VPs characterize cavities in the HB network occupied by the hydroxyl radical.To provide the 3D-visualization we use graphical facilities of commercial MS Excel spreadsheet application, Maple computer algebra system, and freeware POV-Ray (Persistence of Vision Raytracer) ray-tracing program.Advantages and drawbacks of these visualization methods are shortly discussed.Technical and mathematical details are given in ESI.†

Methodology and results
A problem of VP construction has been formulated with respect to the oxygen atoms, nearly allocating mass centres of the cOH and H 2 O molecules.The oxygen atom of cOH is taken as the central point (centre) about which the Voronoi polyhedron is constructed.By denition, VP is the minimal polyhedron whose planar faces are perpendicular bisector planes of lines joining the centre to the oxygen atoms of the neighbouring H 2 O molecules.A region delimited by so determined polyhedron unambiguously denes the space owned by the radical in solution.Since there is one face for each of the nearest neighbours, the number of faces species the coordination (solvation) number.
To calculate the topological and metric properties of cOH aq we have processed the MD simulation trajectory obtained previously. 11The simulation was carried out for the system containing 400 water molecules and one radical molecule at the density 0.994 g cm À3 , corresponding to the density of neat water at 37 C. Radical and water molecules were described by the exible models two-site 11 and three-site, 15 respectively.The employed model potentials include short-range pair interactions of hydrogen atoms to better describe spatial hindrance resulting from the presence of HB network.The simulation time step of 0.1 fs was assumed.The stability of the total energy was 10 À6 < DE/E < 10 À5 .Aer the equilibration stage of about 40 ps, the simulation was extended up to 30 ps.The average temperature of the run was 313 AE 6 K.The equilibrium congurations were recorded every tenth simulation step.We have analyzed 2.5 Â 10 4 equilibrium congurations.The algorithm designed to construct VP is described below.Mathematical details can be found in ESI.†

Construction of Voronoi polyhedron
For a given conguration we have selected coordinates of water oxygen atoms (O i ¼ 1,2,.,N ) within an initially chosen spherical neighbourhood of the radical oxygen atom (central point C) and determined N perpendicular bisector planes to lines joining C with each of the O i atoms.For each plane we marked the region shared with point C.An intersection point of three planes located on the same side of the individual plane as the point C has been classied as VP vertex.To determine edges and faces vertices belonging to each bisector plain have been sorted and correspondingly numbered.Details of sorting method are given in ESI.† For any convex polyhedron the numbers of faces (N F ), edges (N E ), and vertices (N V ) are related by Euler relationship: Since each vertex is the intersection of three bisector planes, N E and N V are related by N E ¼ 3/2 N V , and the topological condition, eqn (1), can be written as:

Calculation of metric properties
Using the coordinates of the numbered vertices we calculate a surface area of each face S i (i ¼ 1, 2,., N F ): In eqn (3) index 1 i corresponds to the vertex #1 belonging to the i-th face and M i denotes the intersection point of CO i line and the i-th bisector plane.
The volume of VP can be calculated as: Taking the tabulation step of 0.01 Å we have varied the neighbourhood radius from 2.5 to 6 Å.The topological condition eqn (2) has not been satised for the radius smaller than 3.26 Å. Fig. 1 shows how the VP volume (Vol) and the number of faces (N F ) change with the increasing neighbourhood radius.Whilst N F increases Vol rapidly decreases approaching the lowest level at a certain radius R marked in the gure.Polyhedra constructed for higher radius are less sharp because of newly-built small faces corresponding to more distant solvent molecules.Appearance of these faces does not change the volume.Polyhedron that corresponds to the radius R has been identied with the solvation cage of cOH.Consequently, the number of faces N F (R) has been taken as the solvation (coordination) number of cOH aq .Metric and topological properties have been calculated for the minimalvolume VP with the smallest number of faces.

Calculation of topological properties
Apart from the cage-radius (R) we dene the face-weighted radius (R w ) as: In eqn ( 5) the contribution of molecules providing larger faces is enhanced by the weighting factor.Taking eqn (4) R w can be alternatively expressed as: So dened face-weighted radius is related with asphericity factor (a) commonly expressed as: 16,17 From eqn ( 6) and ( 7) we obtain: In order to improve a description of the topological properties we have also dened volume quotient Q vol , and distance quotient Q dist .Namely, for a given VP we have set vertex O max which is the most distant from the centre C and plane P perpendicular to the vector CO !max .Plane P divides VP into part A of volume V A , comprising O max , and part B of volume V B ¼ Vol À V A .Quotients Q vol and Q dist have been expressed as: In eqn (10) O max,B denotes the most distant vertex in part B.
Visualization 3D-presentation of VP is essential for both verication of the implemented computational procedure and visualization of a solvation cage.We have tested three visualization methods using computational tools provided by Microso Excel spreadsheet application, Maple computer algebra system, and Persistence of Vision Raytracer (POV-Ray) program.The two former are commercial soware, whereas POV-Ray is a freeware program.Advantages and drawbacks of the tested methods are shortly discussed below.Sample imagines are depicted in Fig. 2.
Input data for all the methods tested are the coordinates of sorted vertices belonging to each of the VP faces.Threedimensional presentation of a solvation cage by using the spreadsheet application is not straightforward.The user has to implement orthogonal projections applying the Gram-Schmidt process. 18The Maple soware offers a command-line utility and ready-to-use macros accepting basic graphical options (colour, transparency, line-style, etc.).It makes 3D-visualization of a solvation cage rather intuitive and easy.POV-Ray program creates photo-realistic images using an advanced rendering technique, called ray-tracing.The POV-Ray code is written in object-oriented C++, hence handling may be difficult for users less familiar with programming languages.Technical and mathematical details on three visualization methods are provided in ESI.†

Statistical handling of VP properties
Calculated probability distributions (normalized probability density functions) of the metric and topological properties of cOH solvation cage are displayed in Fig. 3-5.Statistical description of these distributions is presented in Table 1.We have calculated mean, standard deviation (s), median, mode,  and four dimensionless parameters: Skew 1 (eqn ( 11)), Skew 2 (eqn ( 12)), asymmetry coefficient g (eqn ( 13)), and Kurtosis (eqn ( 14)).Skew 1 and Skew 2 , describe to what extent a given distribution function leans of its mean.
For a perfectly symmetric unimodal distribution, mean ¼ median ¼ mode, and the skewness parameters are equal to zero.
The asymmetry coefficient g is dened as the third standardized moment of a given distribution function: In eqn (13) E and m 3 denote the expectation operator and the 3rd central moment, respectively.
Kurtosis is dened as the fourth standardized moment.It is a measure of the "peakedness" of a distribution and the "heaviness" of its shoulders.
For the normal distribution Kurtosis ¼ 3.

RDF versus VP method
Solute-solvent partial RDFs are commonly used to assess the size of a solvation cage and the coordination (solvation) number.Regarding the size, both the maximum and the minimum of the rst peak are considered, whereas the solvation number is dened as a running integration number spanning the rst peak: In eqn (15) r j is the number density of the j-th site, g ij is the partial RDF calculated for sites i and j, r ij(0) and r ij(m) delimit the position of the rst peak.The nearest neighbourhood of cOH is described by a set of four partial RDFs, O r O w , O r H w , H r O w , H r H w , where the subscripts r and w refer to the radical and water molecules.Features of the calculated RDFs are listed in Table 2. Taking the values of r ij(m) one may assess a size of the solvation cage as about 4.5 Å.This estimate is well above the mean cageradius R (see Table 1).Although the probability density function f(R), displayed in Fig. 3, shows slightly positive skew, the righthand tail extends up to 4.2 Å, i.e. noticeably below r ij(m) .It indicates that the position of the rst minimum substantially overestimates the size of cOH solvation cage.Compared to f(R) the probability distribution of the face-weighted radius R w , also depicted in Fig. 3, is more compact, noticeably shied to the le, and shows the negative Skew 1 .These differences suggest Fig. 3 The normalized probability distribution function of the cageradius R and the face-weighted radius R w .Inset shows the O r O w RDF (solid) 11 and the O w O w RDF (dash) obtained from the MD simulation of neat water at 37 C. 12 Table 1 Statistical description of the probability distributions of: the cage-radius R, the face-weighted radius R w , the volume Vol, the number of faces N F (the solvation number), the total area S tot , the asphericity factor a and the quotients   Å) is substituted for the upper integration limit in eqn (15), the running integration number gives 4.5.This value coincides with the n O w O w coordination number extracted from the simulation of neat water at 37 C. 12 Thus the spherical neighbourhood delimited by R w comprises closer located solvent molecules in the number characteristic for the structural properties of liquid water.It may suggest that cOH occupying deformations in the hydrogen-bond network takes place of H 2 O molecules.
To summarise this section, the conventional methods based on RDFs overestimate both the solvation number and the size of the solvation cage.Regarding the latter, the position of the RDF maximum seems to be more reliable estimate.

Metric and topological properties
Metric and topological properties are not available in the conventional analysis based on RDFs.These properties can be easily captured using the developed methodology.Fig. 5 presents the probability distribution of the metric properties of solvation cage: the volume (Vol) and the total number of faces (S tot ).Both f(Vol) and f(S tot ) are highly leptokurtic (Kurtosis > 3) and show noticeable positive skewness.The mean ¼ 27.24 Å3 provides an estimate for the free volume of cOH in aqueous solution.This cage capacity is reproduced by a sphere of radius 1.87 Å that is much smaller compared to mean R or R w .It may suggest asphericity or shape irregularity of solvation cage.To verify our supposition we have calculated the probability distributions of spatial anisotropy parameters (a, Q vol , Q dist ).The distribution of the asphericity factor a ranges from 1.2 to 2.3 and shows mean ¼ 1.50, median ¼ 1.46, mode ¼ 1.41, and highly positive asymmetry coefficient g.The calculated mode is close to the asphericity factor of a rhombic dodecahedron ð3 ffiffiffi 2 p =p z 1:35Þ, whereas mean is close to a ¼ 1.58 obtained for   the local environment of H 2 O molecules in ambient water. 17It was previously believed that the radical predominantly occupies distortions in the HB network. 11,14However, we have found that ca. 5% of solvation cages shows the asphericity comparable with the Wigner-Seitz cell of the I h ice lattice (a Ih ¼ 2.25).This result indicates that cOH is also located in ice-like regions (patches).At the same time it conrms the existence of patches in neat water 12 and in the diluted solution 14 at 37 C. Since noticeably smaller patches were found in solution 14 one may conclude that localization of cOH in patches perturbs the connectivity of fourbonded water molecules.Shape irregularity of the solvation cage is characterized by Q vol and Q dist .Some of the constructed VPs show distinct disproportion of faces making an impression of shape sharpening (cf.Fig. 2).Such irregularity indicates close proximity of some neighbours and compaction of the others at the opposite side of the radical.As dened by eqn (9), Q vol presents a volumetric ratio of the sharp and the more regular parts of the solvation cage, whereas Q dist in eqn (10) describes an axial asymmetry of these two parts.The probability distributions calculated for Q vol and Q dist show mean and median close to unity, but both are highly leptokurtic.For about one-h of cages Q vol and Q dist exceeded 1.5.In Fig. 6 the anisotropy parameters (a, Q vol , Q dist ) are presented as a function of the solvation number N F .As it can be seen, higher coordination numbers correspond to less aspherical and less anisotropic cages.We have found that for N F > 12 the volume quotient Q vol is less than unity.It means that the compact part of cOH cage takes a larger volume than the sharper one.

Hydrogen bonding issue
Water is a highly associated liquid.Extensive assembling of molecules via hydrogen bonds controls most of the solvent properties of water.Solvation may be considered as a compromise between water-water and solute-water interactions, minimizing Gibbs free energy of the system.Although the Voronoi polyhedra are constructed on the basis of the oxygen atoms, the hydrogen bonding problem can be addressed.We have compared the properties of VP constructed for the hydrogen bonded cOH radical and for the unbound one.To distinguish these two species we have followed the extended energetic denition of H-bond, used previously. 11,14According to this denition a pair of the molecules is hydrogen bonded if a distance between the hydrogen atom of the H-donating partner and the oxygen atom of the H-acceptor is less than 2.5 Å, an angle between the O-H bond of the H-donor and the line connecting the oxygen atoms is not larger than 30 , and the pair interaction energy is at least equal to À8 kJ mol À1 .
The shape of the VP constructed the H-bonded cOH have been analysed for 2 Â 10 3 congurations.Analysis of the molecular neighbourhood of the unbound radical has been performed for 2 Â 10 4 congurations.The statistical parameters of the distributions of metric and topological properties calculated for the H-bonded cOH and for the unbound radical differ by less than 2%.We have found that the asphericity factor a obtained for the solvation cage of the H-bonded cOH is slightly higher, whereas the mean and mode of distribution of other properties (R, R w , N F Vol, S tot , Q vol and Q dist ) are slightly lower.Such a small difference between the H-bonded and the unbound cOH can be expected assuming the cavity mechanism of localization. 11,14Given that, the constructed VPs reveal the shape of cavities of the HB network.

Summary and conclusion
In this paper we report an analysis of the hydrated cOH using computational techniques based on Voronoi polyhedra to process MD simulation trajectory.Construction of Voronoi polyhedra offers a qualitative tool for 3D-characterization of the molecular neighbourhood.Assuming that the solvation cage of cOH is represented by the minimal-volume Voronoi polyhedron (constructed around the radical oxygen atom) we have calculated the statistical distribution of the metric and topological properties of cOH aq at the biologically important temperature (37 C).Our calculations show that the conventional methods based on RDFs overestimate both the solvation number and the size of the solvation cage.The statistical mean of the number of faces of VP indicates that cOH is coordinated by 10 water molecules.This number is noticeably smaller compared to the solvation number (13-14) resulting from the analysis based on partial RDFs.Metric and topological properties of the constructed Voronoi polyhedra reveal features that cannot be captured by RDFs.The most important results that improve understanding of the solvation of cOH in water are listed below.unbound radical show very similar properties as it can be expected assuming localization of cOH in cavities in HB network.

Fig. 1
Fig. 1 Volume (left axis) and the number of faces (right axis) of VP constructed for the water oxygen atoms within a varying spherical neighbourhood.The constructed VP is visualized for the selected radii (3.29, 3.30, 3.32, and 3.43 Å).The solvation cage is represented by the minimal-volume VP with the smallest number of faces (the coordination number).

Fig. 2
Fig. 2 Three-dimensional imagines of cOH solvation cage created by Microsoft Excel, Maple, and POV-Ray.The position of cOH is marked by a central dot.
some asymmetry of cOH solvation cage discussed in the next section.As shown in the inset, the mean value of R w (3.26 Å) coincides with the maximum position of the O r O w RDF.At the same time it coincides with the rst minimum of the O w O w RDF calculated from the simulation of pure water at 37 C. 12 Accuracy of partial coordination number n ij listed in Table 2 depends on a quality of the respective RDF.If the rst minimum is broad and poorly dened n ij is subjected to a signicant uncertainty.The coordination numbers n O r O w n O r H w n H r O w n H r H w indicate that the solvation cage of cOH comprises 13-14 water molecules.The VP-based approach provides a statistical distribution of the number of faces N F , identied with the solvation number.The probability density function f(N F ) displayed in Fig. 4 is not Gaussian, although its kurtosis is near 3.A small negative skewness is indicated by Skew 1 , Skew 2 parameters, and asymmetry coefficient g.The mean value (10.10) and the standard deviation (1.08) of f(N F ) show that the typical solvation cage contains 9, 10, or 11 water molecules.As shown in the inset mean N F is reproduced by the O r O w running integration number at about 3.6 Å, corresponding to the mean value of cage-radius R. If, in turn, the mean value of face-weighted radius R w (3.26

Fig. 4
Fig. 4 Probability distribution of the cOH coordination number, i.e. the number of faces N F .Inset shows the O r O w running integration number versus the upper integration limit in eqn (15).

Fig. 5
Fig. 5 Probability distribution of the metric properties: (a) VP volume Vol; and (b) total area of faces S tot .

( 1 )
We have found that smaller solvation numbers correspond to more aspherical and more anisotropic solvation cages.(2) The face-weighted radius (3.26 Å) delimits the neighbourhood of closer located solvent molecules in the number characteristic for the structural properties of liquid water.(3) The mean volume of Voronoi polyhedra (27.2 Å3 ) provides an estimate for the free volume of cOH aq .(4) The distribution of the asphericity factor reveals the existence of ice-like regions.(5) The Voronoi polyhedra constructed for the H-bonded cOH and for the

Fig. 6
Fig. 6 Spatial anisotropy parameters versus the cOH solvation number (the number of faces N F ): the asphericity factor a (circles), the volume quotient Q vol (squares), and the distance quotient Q dist (diamonds).

Table 2
Features of the partial RDFs describing the nearest neighbourhood of the hydrated cOH radical at 37 C. Positions of the firstpeak maximum (r ij(max) ) and the first-peak minimum (r ij(min) ) are given in Å. Integration over the first peak provides site-site coordination numbers (n ij ) (see eqn(15))