Theoretical study on the hydrophobic and hydrophilic hydration on large solutes: The case of phthalocyanines in water

THE JOURNAL OF CHEMICAL PHYSICS 143, 044502 (2015) Theoretical study on the hydrophobic and hydrophilic hydration on large solutes: The case of phthalocyanines in water Elisa I. Martín,1,2 José M. Martínez,1and Enrique Sánchez Marcos1,a) 1Departamento de Química Física, Universidad de Sevilla, E-41012 Sevilla, Spain 2Departamento de Ingeniería Química, Universidad de Sevilla, E-41012 Sevilla, Spain (Received 27 January 2015; accepted 6 July 2015; published online 23 July 2015) A theoretical study on the hydration phenomena of three representative Phthalocyanines (Pcs): the metal-free, H2Pc, and the metal-containing, Cu-phthalocyanine, CuPc, and its soluble sulphonated derivative, [CuPc(SO3)4]4−, is presented. Structural and dynamic properties of molecular dynamics trajectories of these Pcs in solution were evaluated. The hydration shells of the Pcs were defined by means of spheroids adapted to the solute shape. Structural analysis of the axial region compared to the peripheral region indicates that there are no significant changes among the different macrocycles, but that of [CuPc(SO3)4]4−, where the polyoxoanion presence induces a typically hydrophilic hydration structure. The analyzed water dynamic properties cover mean residence times, translational and orientational diffusion coefficients, and hydrogen bond network. These properties allow a thorough discussion about the simultaneous existence of hydrophobic and hydrophilic hydration in these macrocycles, and indicate the trend of water structure to well define shells in the environment of hydrophobic solutes. The comparison between the structural and dynamical analysis of the hydration of the amphipathic [CuPc(SO3)4]4−and the non-soluble Cu–Pc shows a very weak coupling among the hydrophilic and hydrophobic fragments of the macrocycle. Quantitative results are employed to revisit the iceberg model proposed by Frank and Evans, leading to conclude that structure and dynamics support a non-strict interpretation of the iceberg view, although the qualitative trends pointed out by the model are supported. C2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4927003] I. INTRODUCTION Hydrophobic effects are crucial in several chemical and biological processes.1–7They are generally known as hydrophobic hydration and hydrophobic interaction. The first term is referred to the process of nonpolar solutes dissolved in water (thermodynamic properties, structure, and dynamics), while the second one describes the tendency of nonpolar groups to mutually associate in aqueous solution, this aggregation leading to minimize the total hydrophobic surface exposed to water.8,9Both effects are associated with molecular properties (shape, size, charge, etc.) of the non-polar solute.10 In 1945, Frank and Evans11 proposed a model to understand these observations. The water molecules surrounding a hydrophobic solute must be arranged in the direction of greater crystallinity, forming icelike cages around the solute, which theycalledicebergs, wherethe water-waterinteractions should be reinforced. This simple model has been widely invoked in aqueous solutions containing hydrophobic solutes. On the iceberg model, Kauzmann8proposed the entropic origin of the hydrophobic attraction among caged hydrophobes in water, widely used to explain biochemical processes mediated by water.12–14 During the last two decades, many studies have attempted to confirm experimentally this model by means of different techniques such as neutron diffraction, dielectric relaxation, a)[email protected] and nuclear magnetic resonance,15–23 without reaching a clear evidence for the primitive model.12,24 The hydrophobicity has also been studied theoretically. Computer simulations including Molecular Dynamics (MD), Monte Carlo, and reference hypernetted-chain statisticalmechanic-based methods have addressed the question of characterizing the phenomenon, not only by examining structural properties of water molecules close to the non-polar solutes, but also determining their dynamic properties.25–28 The size and shape of non-polar solutes are related to the water network and its dynamic properties.24,29–38 For the time being, it is generally accepted that hydrophobic properties cannot be associated only with water structural properties, but also with the dynamic behaviour.12,24,26,28 All these studies are focused on the description of the structural and dynamic hydration around non-polar solutes and their comparison with the results found for the hydrophilic hydrationofpolarsolutes.35,39 But,what happenswhendealing withamphipathic solutes?How doesthe concurrenceofhydrophilic and hydrophobic groups on the same solute influence their respective hydration structures in the case of nano-size planar solutes? Phthalocyanines (Pcs) define a family of macrocycles which not only fulfill the previous structural requirements, but have also attracted much attention for their technological applications.40–43 They were recently identified as promising photosensitizers for photodynamic therapy (PDT) because their typical absorption band at ∼680 nm is stronger than that of the first-generation PDT photosensitizers: the porphyrin 0021-9606/2015/143(4)/044502/10/$30.00 143, 044502-1 ©2015 AIP Publishing LLC Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08 044502-2 Martín, Martínez, and Sánchez Marcos J. Chem. Phys. 143, 044502 (2015) family.44 However, Pcs are generally known to be insoluble in waterduetotheirstrongtendencytoaggregateinaqueoussolution.45–48 This fact decreases their potential use as photosensitizer in cancer treatment and therefore chemical substitutions in the macrocycle become an essential requirement to increase their water solubility.49 Thus, apart from metal cations inserted in the macrocycle center, several hydrophilic and amphiphilic groups have been introduced in the macrocycle periphery. These chemical changes have led to different substituted and soluble metal Phthalocyanines (MPcs).50–52 In a previous work,53 we presented a novel procedure to develop MPc-water intermolecular potentials, based on first-principles quantum mechanical information. The chosen Pcs, CuPc, H2Pc, and [CuPc(SO3)4]4−(Figure 1), allow us to monitor the changes from a typical example of hydrophobic hydration, as in free phthalocyanine (H2Pc), to a Pc derivative with solubilizing anionic groups, as [CuPc(SO3)4]4−, passing through an intermediate Pc derivative, where only two hydrogen atoms are replaced by the divalent metal cation Cu() in the macrocycle center (CuPc). Due to macrocycle dimensions, these chemical species cover the intermediate range between the typical small hydrophobic solutes and the large planar surface or cylinder-like biomolecules. Thus, [CuPc(SO3)4]4−is a representative molecule of amphiphilic compounds, bearing both a central planar topology provided by the macrocycle and solubilizing groups in the peripheral regions.26,28,52 Classical MD simulations of the three systems allowed us to study their interactions with water, and analyze the energetic and structural properties of their solutions. Interestingly, the analysis of the surrounding solvent structure for the three compounds does not show as dramatic changes as might be expected from their different solubility character. These structural results open new questions such as the relationships between the hydrophobic chemical nature of the planar complex FIG. 1. Structure and atomic types employed to develop intermolecular potentials between water and the three phthalocyanines studied: (a) CuPc, (b) H2Pc, and (c) [CuPc(SO3)4]4−(red circle denotes the non-equivalent sulphonate ligand). andthe dynamic propertiesof thewater moleculessurrounding the complex. The main goal of the present work is to get insight into the hydration phenomenon in amphiphilic solutes by turning the attention to the hydrogen bond network and dynamic properties of water in the closest macrocycle environment. To fulfill this aim, novel definition of hydration shells adapted to the molecular topology of the macrocycles must be introduced together with a detailed shell-by-shell analysis of the main structural and dynamic properties. This set of results is analyzed in the framework of previous results on hydrophobic hydrationof somenon-polarand polarsolutes of differentsizes and shapes.13,26–28,54 II. METHODOLOGY A. Molecular dynamics simulations MD simulations of CuPc, H2Pc, and [CuPc(SO3)4]4−in water carried out in a previous work53 were performed with the DLPOLY code55 in the canonical ensemble (NVT) using the Nosé-Hoover thermostat at 300 K and periodic boundary conditions. Pc complex-water interaction potentials were specifically developed in our previous study53 and combined with the SPC/E (extended single point charge) water model. Each simulation box contained one Pc complex and 1398, 1401, and 2085 water molecules for the CuPc, H2Pc, and [CuPc(SO3)4]4−systems, respectively. Box length was set to reproduce the water density 0.997 g/cm3at this temperature. B. Shell model definition Due to the oblate shape of the Pc macrocycles, the radial criterion to define the successive hydration shells is no longer valid. A more appropriate definition of a shell model based on spheroidal surfaces is proposed, Sspheroid =(x a)2 +(y b)2 +(z c)2 .(1) For a water molecule centered at (x, y, z), if Sspheroid ≤1 the water molecule is inside the volume defined by the previous surface, whereas if Sspheroid >1 the water molecule is outside the spheroid. If the molecular plane of the Pc macrocycle is defined as the xyplane, then zaxis is perpendicular to the molecular plane. Axes are limited by the geometry of the complex plus the cutoffdistances for the first and second hydration shells. In order to quantify these cutoffvalues we make use of the anglesolved radial distribution function (RDF) defined by an azimuthal angle θrelative to the zaxis in the range θ∈[0◦,30◦], for each atomic type of the complex with respect to the water molecule oxygen. Figure 2shows these angle-solved RDFs for the CuPc case together with the θangle definition. Quite similar values for the minima corresponding to the first and the second hydration shells, centered at 4.7 and 7.8 Å, respectively, are observed (dotted vertical lines in the figure). Anglesolved distribution functions have proven to be an excellent tool for the structural analysis, in particular, of regions which are not defined by a radial distribution, as those imposed by Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08 044502-3 Martín, Martínez, and Sánchez Marcos J. Chem. Phys. 143, 044502 (2015) FIG. 2. Angle-solved X–O distribution function (X =Cu, N1, N2, C1, C2, C3, HC) in the axial region, θ∈[0◦,30◦], for each atomic type of the CuPc complex. Dotted line indicates the cutoffdistance of the water molecules belonging to the first or second shell. (Inset depicts the θangle definition.) molecular planar geometry. Details about angle-solved RDFs can be found in previous papers dealing with square-planar complexes.53,56,57 Thus, the water molecules belonging to the first shell are those inside the spheroidal surface defined by the firstshell axis parameters (a1=b1=11.2 Å; c1=4.7 Å), i.e., they must satisfy the condition S1,spheroid ≤1 (orange region in Figure 3(b)). Second-shell water molecules are those enclosed between the spheroidal surface defined by the first-shell axis parameters, and that of the second one (a2=b2=14.3 Å; c2=7.8 Å). Therefore, they must satisfy two conditions. S1,spheroid >1 and S2,spheroid ≤1 (yellow region in Figure 3(b)). Figure 3(a) depicts a snapshot randomly taken from the simulation of the CuPc complex in water solution showing solvent molecules in the two first hydration shells. Due to the solute shape, the solvent dynamics in different zones inside the first shell can be examined, thus the axial and peripheral sub-regions were defined. From the decomposition of the Cu–O RDF for three θintervals (Figure 9 in Ref. 53), it may be concluded that the axial sub-region includes the first-shell water molecules located at x, y ≤3.0 Å (the axes origin defined at the center of the macrocycle). In the case of [CuPc(SO3)4]4−complex, an additional sub-region inside the first shell was defined to deal with the sulphonate hydration. Figure 4shows the water molecules belonging to the different first-shell sub-regions of the tetrasulphonated CuPc complex. FIG. 3. Snapshot randomly taken of the CuPc complex simulation. (a) Water molecules belonging to the first (red) and second (yellow) solvation shells. (b) Sketch of the spheroidal model showing the limits of the first (orange) and second (yellow) solvation shells. FIG. 4. First-shell water molecules belonging to the axial (blue), peripheral (green), and sulphonate (yellow) sub-regions for a snapshot randomly taken of the [CuPc(SO3)4]4−complex simulation. The sulphonate sub-regions centered on the sulphonate oxygen atoms have been defined from the well-defined peaks appearing in the O–OW (O being sulphonate oxygen) and O–HW RDFs (Figure S1 of the supplementary material).58 III. RESULTS In this section we present a shell-by-shell analysis of the MD simulations carried out for the CuPc, H2Pc, and [CuPc(SO3)4]4−complexes in aqueous solution. A. Orientational distribution To elucidate the way water molecules orientate nearby complex surface, orientational distribution functions of the water molecule inertia axes corresponding to the molecular dipole vector ( uµ), the H–H direction vector ( uHH), and the normal vector to the molecular plane ( u⊥) are considered. The following time-averaged orientational distribution functions with respect to the polar angles θ(azimuthal) defined in Figure S2 of the supplementary material58 can be calculated by means of the following equations: Pu(θ)dθ=⟨θ⟩φdθ. (2) The angle probability for a pure liquid with random orientation of the vectors is59,60 Pu(θ)=1 2sin(θ).(3) Left side of Figure 5shows the orientational distribution functions of the three inertia axes ( uµ, uHH, and  u⊥) for the water molecules belonging to the first (red line) and second shell (green line) relative to bulk water (black line) for the CuPc system. Second shell distributions almost match those of bulk water (green vs. black lines). However, the first-shell distributions are quite different from that of the bulk (red vs. black lines). PHH(θ)probability function presents a narrow peak centered at 90◦, this means that the most representative orientation of the H–H axis is parallel to the Pc Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08 044502-4 Martín, Martínez, and Sánchez Marcos J. Chem. Phys. 143, 044502 (2015) FIG. 5. Bulk-relative orientational distribution functions (θpolar component) of the CuPc system for the uµ,uHH, and u⊥vectors of the water molecules belonging to: left side first shell (red), second shell (green), and bulk (black); right side axial (yellow) and peripheral sub-region (blue). molecular plane, Pµ(θ)function shows a wider preferential distribution around 90◦, whereas the P⊥(θ)distribution presents two maxima around 0◦and 180◦. Altogether, it means that first-shell water molecule plane adopts an average orientation parallel to the macrocycle plane, as found for other authors when dealing with small non-polar solutes.10,61 The narrower peaks for the PHH(θ)distribution reveal that the parallel orientation of the H–H axis with respect to the macrocycle plane is the most restricted rotation for water molecules in this region. In order to check the heterogeneity degree of the orientational distribution along the first hydration shell of CuPc, its distribution was split into contributions due to the axial and peripheral sub-regions. Results are shown on the right side of Figure 5, only the axial water molecules (yellow lines) are preferentially orientated parallel to the molecular plane, whereas the distributions corresponding to the peripheral subregion (blue lines) follow a bulk behavior. Therefore, the axial water molecules are responsible for the peculiar shape of the first-shell distributions already found for the first shell on the leftpartofFig.5for PHH(θ)(redlines).Asimilarreasoningline may be derived when analyzing Pµ(θ), although the broader peak of the axial sub-region (yellow line in Pµ(θ)on the right side of Figure 5) causes a broader first-shell distribution of Pµ(θ)(red line on the left side of Figure 5). The shell-by-shell study of the solvent in the CuPc solution has shown that only the first-shell water molecules exhibit noticeable changes in their orientation against a random distribution. A structural comparison among the three Pc complexes of the orientational distributions, P(θ), of the different first-shell sub-regions has been made on Figure S3 of the supplementary material.58 For the distributions corresponding to the axial sub-region (Figure S3,58 left side), the P(θ)function behavior highlights that the parallel orientation of water moleculeabove andbelowthemacrocyclecenterregionadopts a slightly more parallel orientation for the tetrasulphonated Pc than for the CuPc or the H2Pc. In that concerning the peripheral sub-region (right side of Figure S3 in the supplementary material58), the behavior of H2Pc and CuPc indicates a marginal change of the orientational distribution with respect to the bulk (black and red lines). The significant change is only appreciated for the sulphonate sub-regions of [CuPc(SO3)4]4− (magenta lines) where the distributions correspond to the average orientation of water molecules hydrogen bonded to the sulphonate oxygen atoms. This well defined hydration structure around the solubilizing anionic ligands is also deduced from the O–OW and O–HW RDFs (Figure S1 in the supplementary material).58 From the previous results, it may be concluded that the Pc macrocycle does not impose a large orientational order in the peripheral sub-region of the complexes, except for the hydrophilic environment provided by the sulphonate anions. This sub-region of the soluble Pc macrocycle aside, the main orientational perturbation affects the axial sub-region, where a ratherparallelorientation ofwater molecules isobserved.Nevertheless, the small population of molecules in this sub-region is compatible with the fact obtained from neutron diffraction experiments for aqueous solutions containing small hydrophobic solutes.15,16 Neutron diffraction measurements do not show a particular enhancement in the disordered tetrahedral network of pure water. The small fraction of first-shell water molecules really affected by the non-polar solute has a marginal effect on the experimental diffractograms, which are dominated by bulk water molecules. Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08 044502-5 Martín, Martínez, and Sánchez Marcos J. Chem. Phys. 143, 044502 (2015) B. Dynamic properties Now we turn our attention to the dynamic properties of water in the different shells. The mean residence times (MRT) of water molecules in the first and second shells were obtained adapting the definition given by Impey et al.62 to describe solvated ions. Therefore, we implemented an algorithm to numerically evaluate a survival probability function, Pj(tn,t;t∗)(Eq.(4)).This binaryfunction takesthevalue ofone if the j-th water molecule has been in the solvation shell from a time tnto a time tn+t, without getting out in the interim of this interval more than t∗, and a value of zero otherwise. The average of this function over each time origin tnand the sum over all water molecules gives a survival time-correlation function, n(t), n(t)=1 Nt Nt  n=1 j Pj(tn,t;t∗).(4) Ntis the total number of configurations. Because the ascription of a value for t∗is not unique, two values, previously used by Garcia and Stiller63 and Impey et al.,62 were adopted in our analysis, 0 and 2 ps, in order to provide a reasonable range of MRT values. Thus, t∗=0 ps represents the strict criterion where the survival autocorrelation function is computed on the basis of the exclusive counting of water molecules which did not move from the considered shell. t∗=2 ps is the average lifetime of a water molecule in the hydration shell of another one. Therefore, the latter criteria are more relaxed and consequentlylead tolonger meanresidencetimes.Fromthe fittingof the function n(t)∼e−t/τ to a single-exponential, it is possible to extract a characteristic relaxation time, τ. This means the average time a water molecule resides in the shell before it is replaced by another water molecule coming from the bulk or from another solute hydration shell. MRT results for t∗=2 ps are collected in Table I. For comparative purposes regarding the impact of t∗on the MRT values, in Table S1 of the supplementary material,58 a table with the MRT values corresponding to t∗=0 ps and 2 ps was included. The three systems show similar trends: solvent molecules belonging to the second shell are more labile than water molecules in the first shell, even for the non-soluble H2Pc complex. The presence of copper ion and sulphonate groups slightly increases the residence time in the first shell. To quantify this qualitative observation, the residence time for the axial sub-region was also calculated. These MRT values TABLE I. Mean residence times (ps) of water molecules belonging to the first and second solvation shells and axial sub-region obtained from the MD simulations of the CuPc, H2Pc, and [CuPc(SO3)4]4−in aqueous solution for a t∗value of 2 ps. Values in parentheses indicate the standard deviation of the results. System 1st shell 2nd shell Axial sub-region CuPc 29.4 (0.4) 14.2 (0.1) 19.9 (0.6) H2Pc 28.0 (0.6) 14.0 (0.1) 13.3 (0.6) [CuPc(SO3)4]4−31.5 (0.1) 16.7 (0.1) 22.5 (0.6) indicate that the lack of the copper ion in the center of free phthalocyanine increases the lability of the water molecules in the central region of the macrocycle. In order to compare on the same foot the MRTs obtained for the axial sub-region with those of the peripheral one, the peripheral sub-region for the CuPc system was split in subregions having the same volume than that of the axial one. The MRT value obtained for the solvent molecules belonging to each of these split sub-regions was 12 ±1 ps. These MRTs are comparable with those obtained for several mono- and divalent cations.62,64,65 The translational self-diffusion coefficients of water moleculeswerecomputed accordingtoEinsteinrelation bycomputing the mean square displacement (MSD). Figure S4 of the supplementary material58 shows the global MSD curves for the first and second solvation shells of CuPc system, as well as their components in the x,y, and zdirections. The time windows taken into account for the computation of these functions were the MRTs of Table I. Similar MSD functions were obtained from the simulations of H2Pc and [CuPc(SO3)4]4−in aqueous solution. Table II collects the translational self-diffusion coefficients of water molecules by shells for the three Pcs, as well as the values obtained for the specific first-shell sub-region: axial, peripheral, and sulphonated regions. Diffusion coefficients for first and second shells agree with the behavior already observed for the MRTs. Hence, the water molecule mobility increases from first shell to bulk (the bulk SPC/E value is 2.7×10−5cm2s−1), the second shell having almost recovered the bulk water value. The only exception is found for the diffusion in the sulphonate sub-region, where the mobility is constrained by polyoxoanion interactions beyond the first shell, in part due to the charged character of the sulphonate ligands. TABLE II. Macrocycle-relative translational self-diffusion coefficients, Di(10−5cm2s−1), of the water molecules belonging to the first and second shells, axial, peripheral, and sulphonated sub-regions obtained from the simulations of CuPc, H2Pc, and [CuPc(SO3)4]4−in aqueous solution. The standard deviation was 0.1 (10−5cm2s−1) for all cases. 1st shell Axial sub-region Peripheral sub-region Sulphonated sub-region 2nd shell System Global xyzGlobal xyzGlobal xyzGlobal xyzGlobal CuPc 1.2 1.4 1.4 0.8 0.6 0.9 0.8 0.0 1.2 1.4 1.4 1.0 2.7 H2Pc 1.2 1.4 1.5 0.8 0.9 1.3 1.2 0.1 1.3 1.4 1.4 1.0 2.7 [CuPc(SO3)4]4−0.9 1.0 1.0 0.6 0.5 0.7 0.7 0.0 1.0 1.1 1.1 0.9 1.0a1.0a1.0a1.0a2.3 0.6b0.5b0.8b0.7b aValue obtained for the equivalent sulphonate groups. bValue obtained for the non-equivalent sulphonate group. Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08 044502-6 Martín, Martínez, and Sánchez Marcos J. Chem. Phys. 143, 044502 (2015) First-shell water molecules behave differently from those of the second shell. A homogeneous translation in the x,y, and zdirections is observed for the second shell, whereas for the first-shell water molecules a slowdown in the direction perpendicular to the molecular plane (i.e., z) compared to the xor ydirection is observed. This is the expected anisotropicdiffusion,reflectingthemoleculartopologyofthesolute, which was already found by different authors in the study of protein macromolecules and large hydrophobic surfaces.26,28 The Cu() presence in the center of the macrocycle slows down the mobility of the axial water molecules. The components of the translational self-diffusion coefficient indicate that the axial water molecules do not diffuse along the z axis in contrast with those of the peripheral sub-region. The water mobility in the peripheral sub-region also decreases with respect to pure water. When considering the sulphonate sub-regions, the values for the different components of the axes are similar and homogeneous, the values for the three equivalent sulphonate groups being larger than those of the non-equivalent one, because of the more restricted rotation in the latter case. The comparison among the peripheral regions of the three complexes shows how the sulphonate groups’ presence mostly affects the periphery of the Pc ring decreasing the translational dynamics of the solvent molecules in this sub-region. This was already observed in the protein surface where there are residues bearing quite different polarity. In general, it is observed that the water dynamics near the hydrophilic residues is slower than near the hydrophobic ones.12,28 The rotational properties of the water molecules close to the solute become an additional interesting feature of the solution dynamics. Figure 6shows the reorientational times, τl, µ (l=1 or 2) associated to the inertial axis of the dipole moment unit vector  uµof the water molecules belonging to the first and second shells and bulk for the aqueous solution of the three Pc complexes. Figure S5 in the supplementary material58 shows the corresponding reorientational times associated to  uHH and  u⊥. For the sake of comparison, the bulk SPC/E water FIG. 6. Reorientational times, τl, µ, (l=1 or 2) for the unit vector uµof the water molecules belonging to the first and second shells and bulk water obtained from the simulations of CuPc (red circle), H2Pc (green square), and [CuPc(SO3)4]4−(blue triangle) complexes in aqueous solution. The picture includes the values for 1-order (black dotted line) and 2-order (red dotted line) calculated from the bulk SPC/E water simulation. reorientationaltimeshave beenincluded in thefigure. The bulk values are similar to those published by other authors.66 The computed reorientational times confirm that the firstshell water molecules present a slower dynamics than that of the second shell, whose reorientational dynamics is similar to that of pure water.26,28,67 This behavior agrees with recent studies on the rotational motion of water molecules in the solvation shells of non-polar solutes using polarization-re- solved mid-infrared pump-probe spectroscopy.21,22 In these studies,reorientationaltimesof thewatermoleculesaroundthe solute are two or three times greater than those of pure water. Studies of nuclear magnetic resonance and dielectric relaxation indicate that the average mobility of water molecules in solutions containing hydrophobic solutes decreases.17–20 Since these methods measure the global response of water molecules, these techniques cannot distinguish between bulk water and the non-polar solute solvation shells. Recent classical MD simulations carried out by Laage et al.24 on the water reorientational dynamics around hydrophobic solutes with a medium-size and rather spherical-shape, like Me3NO or (C4H10)3COH, find the same type of partial immobilization detected here. Within the first shell, the values calculated for the water moleculesbelonging to theaxial sub-regionare twicethose obtained for the peripheral one, except for the free Pc (Table III). An interesting feature was found for the sulphonated subregions. The solvent molecules present a different tendency of the reorientational times for the unit vector  uµ(data of  uHH and  u⊥are given in Table S2 of the supplementary material58). In this case, the longest reorientational time is obtained for the H–H direction vector as an indicator of the coordination mode of these water molecules to the polyoxoanion oxygen atoms via donor hydrogen bonds. In addition, the restricted rotation of one of the sulphonated groups causes longer reorientational times of the water molecules solvating it. C. Hydrogen bonding (HB) The energetic, structural, and dynamic analysis of the solvent are completed by a shell-by-shell study of the HBs formed by water molecules. ThestructuralHB definition proposedby Luzar andChandler68,69 was adopted, on the basis of two geometrical parame- TABLE III. Reorientational times, τl, µ (l=1 or 2) for the unit vector uµ of the water molecules belonging to the axial, peripheral, and sulphonated sub-regions inside the first shell obtained from the simulations of CuPc, H2Pc, and [CuPc(SO3)4]4−complexes in aqueous solution. The standard deviation was less that 0.5 ps for all cases. Axial sub-region Peripheral sub-region Sulphonated sub-region System τ1, µ τ2, µ τ1, µ τ2, µ τ1, µ τ2, µ CuPc 15.1 5.0 7.6 2.6 H2Pc 9.9 3.2 7.6 2.7 [CuPc(SO3)4]4−16.8 5.7 8.6 2.9 9.9a3.2a 14.4b4.7b aValue obtained for the equivalent sulphonate groups. bValue obtained for the non-equivalent sulphonate group. Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08 044502-7 Martín, Martínez, and Sánchez Marcos J. Chem. Phys. 143, 044502 (2015) ters, rOO ≤3.5 Å and αH–O···O≤30◦(Figure S6 in the supplementary material58). The distribution functions of the two geometrical parameters, rOO and αH–O···O, did not show any differences for the first and second shells, compared to the bulk water for CuPc and H2Pc (Figure S7 in the supplementary material58). The only noticeable change in the geometrical parameter distribution appears in the case of the [CuPc(SO3)4]4−complex where the distribution for the sulphonate sub-regions differs from the rest of distributions (Figure S8 in the supplementary material58). The disruption is a consequence of the hydrophilic interactions established between the polyoxoanion groups and water molecules. These attracting anion-water interactions compete with the water-water ones, breaking them partially. Thus, the wider Pα(HO·O)distribution in the sulphonate sub-regions shows less orientated HBs as a consequence of the fact that firstshell water molecules accommodate to mainly interact with sulphonate oxygen atoms (see Figure 4, yellow molecules). An interesting feature of the HB network is the coordination number. Table IV shows the average number of waterwater HBs for a water molecule inside each shell for the three Pc solutions. For the bulk water, the coordination number obtained is ∼3.5. This is the SPC/E value for a simulation of pure water. The first shell shows a smaller number of HBs than the bulk as a consequence of the confined region explored, it is worth pointing out that the main disruption is observed for the tetra-sulphonated case. This is the most representative case where hydrophilic regions cause a decrease of the HB average number as a consequence of the attracting anion-water interactions which compete with water-water interactions forming anion-water HBs. The HB number evolution from the first shell to bulk is not regular. On the contrary, the value of the second shell is smaller than those of both the first shell and the bulk. This behavior may be understood on the basis of the particular first-shell structure which must reorganize due to the steric restrictions imposed by the solute and the less favorable average macrocycle-water interactions. The second shell becomes an intermediate zone which cannot follow a well-definedpattern.Interestingly,itisforthetetrasulphonated Pc solution where the number of HBs changes the least when passing from the first-shell to bulk, this is the case where the first shell is partially mediated by hydrophilic interactions. Figure 7complements Table IV giving the distribution of these HBs by shells, indicating for each shell the number of them formed exclusively between water molecules belonging to the same shell (solid bars) or involving water molecules of another shell (dashed bars). The HB fraction formed with outer-shell molecules is ∼25% for the first shell, ∼30% for the second shell, and less than 5% for the bulk (in order to simplify Figure 7, this small contribution has not been plotted). TABLE IV. Average number of HBs of water molecules in the different regions around the solute. System CuPc H2Pc [CuPc(SO3)4]4− 1st shell 3.1 3.1 2.8 2nd shell 2.6 2.6 2.7 bulk 3.4 3.4 3.5 FIG. 7. Average number of HBs formed by water molecules of a shell with molecules of their same shell (solid bars) and those established with water molecules of other shells (dashed bars) for the three Pc aqueous solutions. The hydration network can not only be characterized by the number of hydrogen bonds, but also by their bond energy. Figure 8shows the average HB energy computed for the water molecules in the different shells. Although the energy changes are small, one should bear in mind that we are dealing with average interaction energies, the general trend is extremely instructive. Xu and Berne find a similar behavior for the solvation shell of a 16-residue polypeptide.70 They estimate a relative increase in HB energy of 0.25 kcal/mol with respect to the bulk value. Looking at Figure 8, for the first shell this increase is in between 0.1 and 0.2 kcal/mol. It is worth pointing out that it is for the most hydrophilic solute, i.e., the [CuPc(SO3)4]4−, where this reinforcement is smaller. The second shell also reflects a higher HB energy, i.e., more stabilizing than the bulk value. On the contrary, the smaller values for the first-second intershell regions denote the difficulties associated to the HB arrangement of water molecules obeying different structural environments. The described differences of the binding energy among water molecules of the same or different hydration shells were already reported by Lazaridis13 in a study on hydrophobic hydration of small FIG. 8. Hydrogen bond energy per water molecule, EH B (kcal/mol), for the first and second shell and bulk water obtained from simulations of CuPc (black circle), H2Pc (red square), and [CuPc(SO3)4]4−(green triangle) complexes in aqueous solution. The black dotted line indicates the value calculated for pure water simulations. Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08 044502-8 Martín, Martínez, and Sánchez Marcos J. Chem. Phys. 143, 044502 (2015) non-polar solutes. Although less marked, the same behavior is observed in the second shell-bulk transition region. Therefore, the HB strength is one of the keys to support the old concept of hydrophobic hydration, although the small number of water molecules involved in the phenomenon with respect to the bulk may make their experimental evidence difficult. A femtosecond spectroscopic study of the solvation of hydrophobic solutes shows that the addition of the hydrophobic solutes does not shift the intramolecular OD-stretch vibration to the red.21 It is worth commenting that the use of a polarizable water model should magnify the observed behavior given that the additional degrees of freedom associated to the charge distribution of water molecules must optimize the intermolecular water-water interactions inside a given shell. In order to study the HB dynamics, we computed the HB lifetimes using the interrupted model.71 This analysis is based on a binary function h(t), which is 1 when a hydrogen bond is present and 0 otherwise. Luzar72 defines the time autocorrelation function C(t)based on h(t)function, C(t)=⟨h(0)·h(t)⟩ ⟨h(0)2⟩.(5) C(t)accounts for the probability that a HB persists at time tafter having been formed at time zero. The fitting of this function to a single-exponential form C(t)∼e−t/τ allows the computation of τ.26 Figure S9 in the supplementary material58 depicts the time autocorrelation functions C(t)calculated shell-by-shell for the CuPc simulation. The relaxation times, τ, for the three systems are collected in Table V. The bulk water value agrees with that of 4.4 ps obtained from a SPC/E pure water simulation. This value is similar to those obtained from SPC73 and TIP4P74 water simulations. The HB relaxation time within the first shell is about three times larger than the values for other shells. The axial sub-region relaxation times show (Table V) that the effect of introducing the copper ion in the center of the macrocycle increases the HB lifetime of water molecules in this sub-region. This correlates well with the diffusion coefficient and reorientational time previously discussed, given that it is reflecting a slow down of the water dynamics in this hydration zone. The average relaxation time computed for the HBs formed between the sulphonate oxygen atoms and the water hydrogen atoms belonging to its hydration shell, 26 ±5 ps, is similar to the value for the axial region. This indicates that residence times of water molecules in the hydrophilic region present similar values to HBs close to the central Cu() metal cation. TABLE V. Relaxation times, τ(ps), for HBs formed between two water molecules belonging to a particular region calculated from simulation of CuPc, H2Pc, and [CuPc(SO3)4]4−complex in aqueous solution. Values in parentheses indicate the standard deviation of the results. Region CuPc H2Pc [CuPc(SO3)4]4− 1st shell 14.0 (0.7) 12.4 (0.5) 16.9 (1.1) Axial (1st shell) 22 (4) 7 (4) 25 (5) Peripheral (1st shell) 11.3 (0.6) 12.8 (0.6) 11.7 (0.8) 2nd shell 3.9 (0.1) 3.7 (0.2) 3.9 (0.1) bulk 4.5 (0.1) 4.5 (0.1) 4.6 (0.1) Second-shell values are consistent with the average number of HBs given in Table IV and are smaller than those found for both the first shell and the bulk. IV. CONCLUDING REMARKS The set of phthalocyanines studied here can be considered as an archetype of medium-size solutes having well-defined non-polar fragments and, in one case, hydrophilic fragments. We are not aware of any previous theoretical studies on the hydration of nano-size planar molecules of amphipathic character. The MD simulations reported here were based on refined intermolecular potentials with site-site parameters fitted to quantum-mechanical information. The new strategy to explore hydration of large planar macrocycles makes use of a spheroidal shell model defined for solvent molecules that is based on solute shapes and MD structural information. Definition of sub-regions within the first hydration shell allowed us to monitor the hydrophobic and hydrophilic interactions in the regions close to the solute. We observed that the main geometrical parameters of the solvent structure are not significantly deviated from those of the pure water network. Only water molecules of the first shell present a parallel orientation of their H–H vector with respect to the macrocycle plane. As expected, the orientational distribution is different in the sulphonate regions where water molecules solvate the negatively charged groups. Structural changes were found in the HB network. The average number of HBs in the first shell is slightly smaller than that found in bulk, partly due to the presence of the macrocycle surface. However, this average number decreases slightly when the phthalocyanine macrocycle includes hydrophilic sulphonate groups. In that case, the structure of water is altered (disrupted) to a larger extent because of polyoxoanionwater HB formation. It is worth pointing out that the second shell exhibits the lowest average of HBs. The mean HB energy evolution tends to compensate this effect at least partially. The highest (most negative) value is found in the first shell and decreases up to the bulk, as reported in other structural cases, such as protein and micelle surfaces.28 The lowest interaction energy is observed in the first shell of the [CuPc(SO3)4]4− complex, which reflects the competition between water-water and sulphonate-water interactions. According to the structural andenergyanalysisof theHBnetwork,it canbeconcluded that the solvent adopts an average structure defined by shells where in-shell interactions prevail over the inter-shell interactions. The same behavior is found by Lazaridis13 in the study of methane. This suggests a similar energy pattern for the inter- and intra-shell solvent-solvent interaction energies, regardless of the solute size and shape. Translational and orientational dynamics of water molecules in the first shell clearly exhibit a general slowdown of their diffusional properties in an anisotropic way. The properties associated to the axis perpendicular to the molecular plane are more affected than those belonging to the axes defining the molecular plane. This behavior was previously observed in protein-solvent interface.26 HB lifetimes are larger for first-shell solvent molecules than for bulk, in agreement with the above-mentioned slower dynamics. Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08 044502-9 Martín, Martínez, and Sánchez Marcos J. Chem. Phys. 143, 044502 (2015) This global analysis leads us to conclude that the principal signature of hydrophobic hydration is found in the dynamics rather than any structural changes. Both phenomena would support the iceberg model proposed by Frank and Evans11 if their simple interpretation of the hydrophobic hydration was not literally assumed. As already pointed out by Galamba,38 from a structural point of view, rather than an ice-like structure of the hydration shell, the hydrophobic hydration acts like a local cooling of the sample. However, experimental confirmation would be extremely difficult considering that hydrophobic contactsconstituteonly asmallfraction ofthetotal. Asdemonstrated in other non-polar solutes,10,12 the only method capable of probing local dynamical effects is atomistic simulations. The amphipathic character of the solute has been examined in the sulphonated derivative. The impact of the coupling between the hydrophobic and hydrophilic fragments on the phthalocyanines has been determined by the comparative analysis of the water properties in the different first-shell subregions of CuPc and [CuPc(SO3)4]4−. The general outcome is the slight influence exerted by the sulphonate groups on the properties of the axially located solvent molecules, revealing a large degree of uncoupling of the axial region properties from those of the peripheral one. This conclusion seems to support the experimental finding of Nyokong et al.52 who observe that asymmetric substitution of solubilizing groups on Pcs tends to increase their water solubility. This experimental fact agrees with the weak coupling among regions, supporting that the main contribution to the Pc stabilization in aqueous solution comes from the increase of high-order multipole moments, provided by the localized contributions of the separated polar or charged solubilizing groups. Analysis of the dynamics is clearly able to characterize the solute hydrophobic effects, although literal adoption of the iceberg concept seems to be too strict. Hynes et al.24 concluded that in the case of a metallo-macrocycle there is some retardation of the diffusive properties rather than a total immobilization of water molecules. Our results support that conclusion, and indicate a formal cooling rather than a freezing.26,38 Pc derivative compounds hold the promise of important applications across chemistry, chemical engineering, and pharmacology. As Nyokong et al.51,52 have stressed, understanding their hydration behaviour is key to the synthesis of Pc derivatives in liquid media. We believe that this work gives a new perspective on Pc solution chemistry. ACKNOWLEDGMENTS We thank Dr. C. Domene (King’s College of London) and Professor K. Refson (Royal Holloway, University of London) for their useful comments and suggestions. Financial support from the Spanish Junta de Andalucia (No. P11-FQM7607) is acknowledged. E.I.M. acknowledges a predoctoral grant from the Spanish Minister of Sciences and Innovation (No. CTQ2008-05277). 1N. T. Southall, K. A. Dill, and A. D. J. Haymet, J. Phys. Chem. B 106, 521 (2002). 2D. Chandler, Nature 437, 640 (2005). 3K. P. Murphy, P. L. Privalov, and S. J. Gill, Science 247, 559 (1990). 4B. Widom, P. Bhimalapuram, and K. Koga, Phys. Chem. Chem. Phys. 5, 3085 (2003). 5B. Lee and G. Graziano, J. Am. Chem. Soc. 118, 5163 (1996). 6G. Hummer, S. Garde, A. E. Garcia, A. Pohorille, and L. R. Pratt, Proc. Natl. Acad. Sci. U. S. A. 93, 8951 (1996). 7S. K. Pal and A. H. Zewail, Chem. Rev. 104, 2099 (2004). 8W. Kauzmann, Adv. Protein Chem. 14, 1 (1959). 9P. L. Chau, Mol. Phys. 101, 3121 (2003). 10M. Kinoshita, J. Mol. Liq. 119, 47 (2005). 11H. S. Frank and M. W. Evans, J. Chem. Phys. 13, 507 (1945). 12P. Ball, Chem. Rev. 108, 74 (2008). 13T. Lazaridis, J. Phys. Chem. B 104, 4964 (2000). 14X. Huang, C. J. Margulis, and B. J. Berne, Proc. Natl. Acad. Sci. U. S. A. 100, 11953 (2003). 15P. Buchanan, N. Aldiwan, A. K. Soper, J. L. Creek, and C. A. Koh, Chem. Phys. Lett. 415, 89 (2005). 16J. Turner and A. K. Soper, J. Chem. Phys. 101, 6116 (1994). 17R. Haselmaier, M. Holz, W. Marbach, and H. Weingartner, J. Phys. Chem. 99, 2243 (1995). 18Y. Ishihara, S. Okouchi, and H. Uedaria, J. Chem. Soc., Faraday Trans. 93, 3337 (1997). 19U. Kaatze, H. Gerke, and R. Pottel, J. Phys. Chem. 90, 5464 (1986). 20W. Wachter, R. Buchner, and G. Hefter, J. Phys. Chem. B 110, 5147 (2006). 21Y. L. A. Rezus and H. J. Bakker, Phys. Rev. Lett. 99, 148301 (2007). 22A. A. Bakulin, C. Liang, T. La Cour Jansen, D. A. Wiersma, H. J. Bakker, and M. S. Pshenichnikov, Acc. Chem. Res. 42, 1229 (2009). 23C. Pangali, M. Rao, and B. J. Berne, J. Chem. Phys. 71, 2975 (1979). 24D. Laage, G. Stirnemann, and J. T. Hynes, J. Phys. Chem. B 113, 2428 (2009). 25G. M. Torrie and G. N. Patey, J. Phys. Chem. 97, 12909 (1993). 26A. R. Bizzarri and S. Cannistraro, J. Phys. Chem. B 106, 6617 (2002). 27M. Kinoshita, J. Solution Chem. 33, 661 (2004). 28B. Bagchi, Chem. Rev. 105, 3197 (2005). 29W. Blokzijl and J. B. F. N. Engberts, Angew. Chem., Int. Ed. Engl. 32, 1545 (1993). 30S. Okazaki, K. Nakanishi, H. Touhara, and Y. Adachi, J. Chem. Phys. 71, 2421 (1980). 31S. Okazaki, K. Nakanishi, H. Touhara, and Y. Adachi, J. Chem. Phys. 72, 4253 (1980). 32S.Okazaki,K.Nakanishi,H.Touhara, N. Watanabe,andY.Adachi,J.Chem. Phys. 74, 5863 (1981). 33P. L. Chau, T. R. Forester, and W. Smith, Mol. Phys. 89, 1033 (1996). 34C. Y. Lee, J. A. McCammon, and P. J. Rossky, J. Chem. Phys. 80, 4448 (1984). 35S. H. Lee and P. J. Rossky, J. Chem. Phys. 100, 3334 (1994). 36A. Wallqvist and B. J. Berne, J. Phys. Chem. 99, 2885 (1995). 37M. Mateus, N. Galamba, and B. Cabral, J. Chem. Phys. 136, 014507 (2012). 38N. Galamba, J. Phys. Chem. B 117, 2153 (2013). 39M. T. Stone, P. J. In’t Veld, Y. Lu, and I. C. Sánchez, Mol. Phys. 100, 2773 (2002). 40M. S. Xu, J. B. Xu, M. Wang, and L. Que, Appl. Phys. 91, 748 (2002). 41L. E. Norena-Franco and F. K. Vasnik, Analyst 121, 1115 (1996). 42A. Hagfeldt and M. Gratzel, Acc. Chem. Res. 33, 269 (2000). 43M. Calvete, G. Y. Yang, and M. Hanack, Synth. Met. 141, 231 (2004). 44R. Bonnett, Chemical Aspects of Photodynamic Therapy (Gordon and Breach Science Publishers, Amsterdam, 2000). 45C. C. Leznoffand A. B. P. Lever, Phthalocyanine-Properties and Applications (VCH Publishers, Germany, 1989). 46N. B. Mckeown, Phthalocyanine Materials-Synthesis Structure and Function (Cambridge University Press, UK, 1998). 47R. B. Hammond, K. J. Roberts, R. Docherty, M. Edmondson, and R. Gairns, J. Chem. Soc., Perkin Trans. 2 1996, 1527. 48S. Yim, S. Heutz, and T. S. Jones, Phys. Rev. B 67, 165308 (2003). 49M. C. DeRosa and R. J. Crutchley, Coord. Chem. Rev. 233, 351 (2002). 50N. L. Oleinick, A. R. Antunez, M. E. Clay, B. D. Rihter, and M. E. Kenney, Photochem. Photobiol. 57, 242 (1993). 51T. Nyokong, Coord. Chem. Rev. 251, 1707 (2007). 52F. Dumoulin, M. Durmus, V. Ahsen, and T. Nyokong, Coord. Chem. Rev. 254, 2792 (2010). 53E. I. Martín, J. M. Martínez, and E. Sánchez Marcos, J. Chem. Phys. 134, 024503 (2011). 54X. Huang, C. J. Margulis, and B. J. Berne, J. Phys. Chem. B 107, 11742 (2003). Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.182.116 On: Wed, 19 Oct 2016 14:02:08