Zinc- and copper-porphyrins in aqueous solution – two similar complexes with strongly contrasting hydration

Abstract

We present detailed analysis of the hydration behavior of zinc and copper bound porphyrins treated via ab initio quantum mechanical charge field molecular dynamics which agrees well with available experimental data. The computed metal–water coordination in the case of zinc bound porphyrin demonstrates a strong association of water with zinc compared to the copper–water interaction which correlates well with the calculated free energy of binding. The H-bond dynamics in these hydrated systems yield weaker H-bond interactions as compared to that observed in the case of metal-free porphyrin; nevertheless, the effect of metal association with porphyrin resulted in shifts in the vibrational frequencies. These characteristic data suggest a contrasting behavior between these metalloporphyrins in solution which could also serve to correlate with the properties of biological systems.

---

1 Introduction

The liquid state is an integral component of a large number of chemical and virtually all biochemical systems, since most of the processes taking place in these systems require the presence of a liquid. Investigations of the underlying dynamics in solution, which may consist of hydrated ions and composite solutes as well as aqueous bioorganic/bioinorganic complexes and macromolecules in water, require a rigorous treatment that can be achieved by the applications of sophisticated experimental techniques. On the other hand, ab initio based theoretical methods have been shown to provide valuable insights not only into structural but also dynamical properties of the solvated systems.¹ In particular the dynamics of chemical systems is known to also play a major role in the understanding of structure–function relationship of various chemical and biochemical systems. A number of experimental techniques such as nuclear magnetic resonance (NMR) and X-ray crystallography provide a means for the profiling of overall structural attributes of the system but very little information on ultrafast dynamical properties at the sub-nano second time scale can be retrieved.² The investigation of dynamics is challenging due to insurmountable technical issues associated with the experimental setups. On the other hand, the evaluation of dynamics of the solute–solvent interactions in aqueous ionic solutions and investigation of the mean lifetimes of bonds and coordination in pure water were made possible by the application of incoherent quasielastic neutron scattering (IQENS) experiments and femtosecond laser pulse spectroscopy, respectively.^3,4 The dynamical properties of (bio)chemical coordination complexes in an aqueous environment are known as more intricate targets, since experimental methods were found to be sensitive to both the size and number of molecules as well as the timescale of dynamical processes.

The limitation of experimental methods to investigate the dynamical properties of complex systems led to the development of theoretical methods, and ab initio based simulation studies are amongst the most advanced tools capable of reproducing time-dependent chemical and biochemical processes at the femtosecond time-scale.^4,5 Thus, detailed investigations of the structure and dynamics of complex systems in aqueous solution were made possible by employing modern theoretical methods coupled with statistical simulation techniques such as molecular dynamics (MD).⁶ These methods are generally based on potential functions derived either from empirical potential data or from increasingly accurate ab initio calculations of potential energy surfaces. Potential functions are approximations, and thus, in the case of stronger interactions between particles, sophisticated three- or even four-body interaction potentials are required, which are in many cases too time consuming and troublesome to construct. Further important features related to interactions such as polarization and charge transfer effects are in many cases not taken into account by classical approaches, for instance in the case of complex biological systems in which chemical/biochemical processes take place at the electronic level.

To take n-body, polarization and charge transfer effects into account, the application of quantum mechanics (QM) in conjunction with the established molecular dynamics (MD) framework proved helpful but the treatment of a representative complex system consisting of several hundred particles via a QM method is not within the scope of contemporary computational equipment. In this context two solutions were put forward to investigate solvated systems and solutions: the first approach known as Car–Parinello MD (CPMD) is based on the utilization of simplified density functional theory (DFT) methods at the generalized gradient approximation (GGA) level that was applicable only for systems having a limited number of particles.⁷ On the other hand, hybrid quantum mechanical/molecular mechanical (QM/MM) techniques partition the system into the most relevant subsystem containing a solute and its immediate solvation layers treated by QM whereas the interactions in the remaining part of the system are accounted for via adequate force field techniques.⁸ A variety of QM/MM protocols have been developed and the quantum mechanical charge field molecular dynamics (QMCF-MD) approach was explicitly designed for the treatment of solvated systems.^9–11 The QMCF-MD formalism was successfully applied to numerous solvated systems ranging from small aqueous ions and composite anions to organic and complex biochemical systems in aqueous solution.¹² Recently, a study was reported on the QMCF-MD simulation results of structural and dynamical properties of the hydrated apo- and magnesium-porphyrin.¹³

Considering the vigorous need to obtain insights into the structure and dynamics of solvent molecules in complex biochemical systems at the atomic level, an extended effort was made by applying the QMCF-MD technique to zinc- and copper-bound metalloporphyrins in aqueous solution. The association of porphyrins and metals leads to different forms of metalloporphyrins with specific reactivity in many proteins and enzymes. This has stimulated considerable interest to probe biological and chemical properties of these molecules via experimental methods. Despite the widespread utilization of experimental techniques, atomic level descriptions concerning the structure and dynamics of these molecules could not be achieved. This was complemented by employing different theoretical approaches based on various levels of theories ranging from classical to ab initio methods.¹⁴ A number of force field based simulation studies were reported on metalloporphyrins for the investigation of their structures and dynamics related to biochemical processes.¹⁵ Metalloporphyrin–water interactions are also a much debated area of theoretical research that are investigated to obtain an accurate description of solvation properties, since the majority of the biological porphyrins exist in an aqueous environment.

One of the less investigated metalloporphyrin systems is the zinc protoporphyrin (ZnPP), a metabolite formed in trace amounts during heme production in the case of lead poisoning or iron deficiency.¹⁶ The analysis of the fluorescence of the ZnPP molecule was mimicked using zinc porphyrin (ZnP) which signifies its importance as a prototype in many cases in understanding the photochemical properties of conjugated systems. A number of ZnP species were also shown to play an essential role in various chemical processes due to their characteristic structural properties and excitation energies of ZnP systems in aqueous solution determined by using various theoretical chemistry methods.^17,18 The meso-substitution of ZnP species imparts solubility to these molecules thus leading to their aggregation in aqueous solution.¹⁹ The coordination or binding properties are also of considerable interest since energy decomposition analysis showed that the ZnPP molecule has greater affinity for nitrogen based ligands compared to ligands binding via oxygen atoms.^20,21 A further metalloporphyrin of great significance is copper porphyrin (CuP) which was identified as a naturally occurring red color pigment and an isolate from deep-sea sediments as terrestrially derived oxidized organic matter.^22,23 The biological role of CuP was implicated in bronze baby syndrome (BBS), an uncommon side-effect of phototherapy.²⁴ CuP and its various derivatives are involved in photoprocesses that are of interest because of their DNA-binding interactions which have been probed by both experimental and theoretical methods.²⁵ Like ZnP, meso-substituted CuP molecules also show a tendency to form aggregates in aqueous solution.¹⁹ As far as axial binding of ligands to these porphyrins is concerned, it was reported almost without exception that CuP is four-coordinated in-plane lacking the axial ligands whereas ZnP tends to be penta-coordinated.²⁶

All four nitrogen atoms of porphyrin systems were found to coordinate with metal ions. However, the latter are displaced from the plane due to the high conformational freedom of porphyrin molecules.¹³

2 Methods

Considering the importance of ZnP and CuP in aqueous solution, an ab initio QMCF-MD simulation study was carried out to elucidate their structural and dynamical properties. The QMCF-MD formalism^9–11 basically follows a similar partitioning scheme on which the conventional QM/MM approaches are based, but a multi-level description was required to explicitly treat the solution systems, since the solvents in the QM and MM regions diffuse across the two regions. The QM region in the case of the QMCF ansatz is further divided into the inner core and the solvation layer zone (both treated at the quantum mechanical level) whereas interactions between atoms in the remaining MM region are treated via force field potentials. For the quantum mechanical calculation, the Hartree–Fock level of theory along with double-ξ + polarization basis sets for non-metals with some additional basis sets in particular for metals was found to be a good compromise between accuracy and computational effort,^9–11 also yielding data in good agreement with experiment in a previous study of porphyrin and MgP.¹³

Typically, the accuracy of theoretical levels is benchmarked using small gas-phase model systems (e.g. metal–water clusters with an increasing number of ligands) comparing the predicted energy and structural data.^27–29 However, due to the reduced number of ligands in this case (i.e. the porphyrin system plus one or two water molecules) and the complex electronic structure of transition metal ions (in some cases preferring five- or even four-fold coordination^30–36), cluster computations do not permit an unambiguous benchmarking of the theoretical level via gas-phase computations. The application of polarizable continuum models may improve this situation; however, due to the different numbers of ligands, the volume and surface area of the associated cavity are not constant, which introduces further uncertainties in the evaluation. In such a case a treatment in the bulk environment (i.e. a water box under periodic boundary conditions) at finite temperature as done in QM/MM MD simulations provides a more reliable means to describe the system and to assess the quality of the theoretical treatment. This was for instance also the case in the previous study of porphyrin and MgP¹³ discussed above.

The description of interactions of QM atoms with MM particles is the most challenging task when formulating a QM/MM approach. In the case of QMCF-MD, the MM partial charges were included in the core Hamiltonian for the QM region via the perturbation term V′:

where M is the number of partial charges q_J, and r_iJ corresponds to the distance vector between electron i and partial charge J. This treatment often referred to as electrostatic embedding³⁷ simultaneously accounts for the interaction between atoms in the QM zone and their interactions with the field of the surrounding charge distribution of the MM atoms. Thus, quantum effects such as charge transfer, polarization and many-body interactions occurring in the QM zone are visible in the QM/MM coupling. Since the electrostatic embedding only accounts for the Coulombic contributions, non-Coulombic interaction potentials have to be provided for the QM/MM coupling as well (i.e. potential embedding). If the QM region is sufficiently large (diameter of approx. 15 Å or beyond), the separation of atoms located in the central core QM region^9–11 such as the Cu(II) and Zn(II) ions from the MM atoms is beyond the typical non-Coulombic cutoff distance in the range of 6–8 Å. Since atoms in the QM layer region^9–11 (i.e. atoms of the porphyrin ligand and QM solvent molecules) are much closer to the QM/MM interface, non-Coulombic potentials are explicitly required to adequately account for the respective interactions. Since a full layer of hydration is explicitly considered in the QM region (making the computation rather demanding), the main contribution from the non-Coulombic QM/MM coupling results from QM solvent–MM solvent interactions. With water being the most investigated substance, a large number of suitable interaction potentials are available in the literature. A further advantage of including a full layer of water molecules is a large distance between atoms of the porphyrin ligand and MM particles. Due to these large separations, the contributions of the respective non-Coulombic potentials are small compared to the Coulombic contributions resulting from the embedding.

The QM and the MM region is interfaced through a special treatment for the smooth exchange of particles.^9–11 A small layer of a typical thickness of 0.2 Å between the two major regions is defined at the QM/MM interface region. This typical QM/MM interfacial distance was found to be optimum for the smooth transitions of molecules as well as to avoid the occurrence of significant transition artifacts. Further details of the simulation technique are discussed in previous reports.^9–11

Simulation protocol

The ab initio QMCF-MD simulations of Zn(II)- and Cu(II)-porphyrin have been carried out in a simulation cube with a side length of 39.28 Å (1000 water molecules) subject to periodic boundary conditions. Both systems were modelled by replacing the ion of the previously studied aqueous magnesium porphyrin simulation¹³ with Zn(II) and Cu(II), respectively. The Hartree–Fock (HF) level of theory and 6-31G** basis sets for water and porphyrin atoms, and Stuttgart RSC 1997 and LANL2DZ basis sets coupled with effective core potentials for zinc and copper atoms were employed, respectively. As mentioned earlier, the HF method used in this study has already been successfully applied to numerous aqueous systems ranging from single ions to composite anions to complexes in aqueous solution,^38–40 while DFT methods have been shown to yield less accurate results for aqueous systems⁴¹ due to an inadequate description of hydrogen bonding. The SPC/E water model (extended simple point charge model) was employed to describe water in the MM region in both systems.^42,43 The equations of motions were integrated using a velocity-Verlet algorithm⁴⁴ with a time step of 2 fs. The simulations were performed in the NPT ensemble, and temperature and pressure were maintained at 298.15 K and 1 atm employing the Berendsen thermostat and manostat⁶ with a relaxation time of 0.1 ps and 0.5 ps, respectively. The reaction field method was employed using a cutoff distance of 18.0 Å to treat long-range electrostatics and all bond lengths involving hydrogen atoms of the systems were constrained employing the M-RATTLE algorithm,⁴⁵ similar to the M-SHAKE⁴⁶ technique. The entire QM zone (QM core plus QM layer region) with a radius of 8.5 Å was treated at the respective quantum mechanical level. To realise the QM/MM coupling according to the QMCF protocol^9–11 outlined above the QM region was further subdivided into the core and the solvation layer zone. Only the metal (i.e. zinc and copper) was included in the QM core region (radius of 0.5 Å). The QM layer region comprises the remaining part of the QM zone, i.e. all atoms of the porphyrin ligand plus solvent molecules. The smoothing region was defined to interface QM and MM regions via a separation of 0.2 Å to ensure smooth diffusion of solvent particles between the QM and MM regions.^9,10 Both systems were equilibrated at 298.15 K followed by a sampling of simulation trajectories for 10 ps using a time step of 2 fs.

3 Results and discussion

The radial distribution function g(r) enables the analysis of solute–solvent relationships in the systems by providing information on the distribution of particles as a function of the respective pair distance. The evaluation of water molecules around the central atom enables the detection of the overall hydration shell surrounding the entire solute. The distribution of water oxygen (O_w) and hydrogen (H_w) with respect to the central metal ion in both ZnP and CuP has been evaluated (Fig. 1). The Zn–O_wg(r) profile depicts a pronounced peak between 1.91 and 2.93 Å with a split pattern observed at the tip of the peak that indicates variations in the binding of water molecules to the zinc ion.

Fig. 1 Radial distribution functions of water molecules with respect to (a) zinc and (b) copper atoms of the ZnP and CuP molecules, respectively.

The splitting of the first peak is in part due to statistical noise. However, the presence of other ligands in the vicinity of the ions (as seen also in the distance plots depicted in Fig. 4) has a direct influence on the binding characteristics leading to small shifts in the respective metal–water distances. From this perspective the splitting of the first peak reveals significant information. In order to address this point in more detail, much longer simulations have to be carried out, which is at present beyond the capacities of available computational equipment. Despite these limitations, the presented simulation data are adequate to highlight the differences between the two presented systems.

The Zn–O_w peak was separated from the second peak (between 3.48 and 5.01 Å) highlighting a strong coordination of water at the zinc coordination site. In the case of the Cu–O_w profile the intensity of this peak is strongly diminished with a low water population near the Cu(II) ion. The distinction between the first and second peaks is barely visible, corresponding to a very labile character of the first shell water molecules, most likely displaying ultrafast ligand exchanges with the second shell. The second peaks in both g(r) profiles were continued to broad peaks ranging from 5.31 to 8.72 and from 6.62 to 8.65 Å for the ZnP and CuP case, respectively. These peaks correspond to the hydration layer covering the entire solutes including all water molecules involved in the coordination to the metal and hydrogen bond formation.

The peaks correspond to a slight ordering of water molecules surrounding the entire porphyrin moiety. From the classical pre-equilibration of the MgP system it has been observed that beyond a distance of 8.5 to 9.0 Å no distinct ordering effect is observed in the pair distribution.¹³ Based on this finding the size of the QM region was chosen, which was also applied in the current work. The slight peak at 8.5 Å in the Zn–O pair distribution can be regarded as noise, since no pronounced peak is observed in the Zn–H RDF or the CuP case.

It should be mentioned that a main benefit of the embedded treatment inherent to the QMCF procedure is that QM/MM transition artifacts are greatly reduced over a conventional QM/MM treatment employing a subtractive approach (i.e. treating the full system via MM, then subtracting the MM contributions of the isolated high level region and finally adding the QM contributions of the high level zone). In the latter case some systems showed pronounced QM/MM artifacts which was one of the reasons to develop the more advanced QMCF simulation technique.^9–11

Based on the probability functions of the water oxygens in the case of ZnP, the mean distance between the oxygen atom of the axially coordinated water molecule and the zinc ion was found to be 2.13 Å whereas in the case of hydrated CuP, the weak coordination of water molecules to the copper ion prevents the analysis of a particular ion–oxygen distance. The running integration associated with the Zn–O_wg(r) functions points toward the coordination of approximately one water molecule to the zinc ion whereas the probability of a water molecule in direct coordination to the copper ion was negligible, thus indicating penta- and tetra-coordination for Zn(II) and Cu(II), respectively.⁴⁷ The respective five- and four-fold coordination in the hydrated ZnP and CuP molecules is in good agreement with experimental estimations.²⁶ The difference in the coordination behavior of water molecules towards zinc and copper ions also depends on the degree of coordinative unsaturation of ZnP and CuP molecules, which results due to the electronic properties of the host molecule to accommodate a nucleophilic molecule at its coordination site.⁴⁸

RDFs poorly specify pairwise contributions in non-spherical systems, and thus, the relative positions and orientations of molecules must be characterized in more detailed distribution functions. Having this objective in mind, angular-radial distribution (ARD) functions were plotted to obtain detailed insight into the structural properties, since these types of distributions include the angular components along with distances of the two interacting atoms.¹¹Fig. 2a–d illustrates the ARD plots of O_w atoms surrounding the whole ZnP and CuP molecules complementing the structural data deduced from the g(r)-functions. Contour plots of the hydrated solutes provide information on different regions of space occupied by water ligands. A single peak near the central metal was observed in both cases, thus indicating the presence of axially bound water molecules (cf.Fig. 2a and b). In the case of ZnP, the much higher intensity clearly indicates the presence of a strongly bound water molecule in the axial position of the Zn(II) ion, whereas in the case of hydrated CuP, a much lower intensity is visible highlighting the strongly reduced possibility of finding axial water at the copper ion (cf.Fig. 2b). The coordination of water towards the metal in both solutes correlates well with an earlier simulation study on the hydration of magnesium bound porphyrin reporting preferential localization of one water molecule bound to the magnesium ion.¹³ On the other hand, further information on water molecules involved in hydrogen bonding with the nitrogen atoms of ZnP and CuP molecules was based on separate ARD functions generated for water molecules, thereby excluding the axial one, since the high density of the latter suppresses the density distributions of all surrounding water molecules. Symmetric density distributions were observed for water ligands involved in H-bonding with these solutes above and below the plane of the ring as assessed from the contour plots (cf.Fig. 2c and d). In the case of Zn–O_w density, off-axis peaks corresponding to second shell water molecules that are H-bonded to first shell water were observed. In addition structural organization above and below the porphyrin-plane is clearly visible in the ARD, while the rather uniform distribution visible in the CuP case implies that no pronounced structural motif occurs throughout the simulation.

Fig. 2 Angular radial distribution contour plots of water oxygens with respect to (a) zinc of ZnP, (b) copper of CuP showing the probability of axial water, and in the absence of axial water in the case of (c) ZnP and (d) CuP.

The structural features discussed so far describe the potential interaction between the solute and water in different hydration shells, and thus an evaluation of the free energy of binding associated with metal–water interactions in both solutes enables quantification of the attractive/repulsive character of water binding. For this purpose, the respective potentials of mean force (PMF), W(r), were calculated for the water molecules as a function of the radial separation from the metallic center in the ZnP and CuP molecules using the following expression;⁴⁹

W(r) = −RT ln(g(r)) (2)

with R being the universal gas constant and T the temperature. Fig. 3a and b display the PMF plots associated with the Zn–O_w and Cu–O_wg(r) functions, thus yielding a measure of the free energy of binding of the water–metal interaction. The respective binding energies of −0.75 and +2.39 kJ mol⁻¹ for the Zn–O_w and Cu–O_w coordination clearly indicate the strong contrasting character of water binding towards these metals. The negative value of the binding energy for the zinc–water interaction points towards a strong affinity of O_w atoms to zinc(II) which is comparable to the metal-free porphyrins having a similar value of the free energy of binding. In that case, the H_w atoms orient themselves towards the nitrogen atoms to form H-bonds which on the other hand are hindered by the presence of metals bound to porphyrins.¹³ In the case of CuP on the other hand, a dominant repulsive interaction was deduced from the positive magnitude of the free energy of binding that correlates well with the structural information on the Cu–O_w interaction. The water repulsion by CuP is assumed to enable the formation of a H-bond network between hydration layers that could also be assessed by plotting O_w–O_wg(r) functions. The difference in the free energies of binding for these two solutes is attributed only to the electronic structure of the two transition metals interacting with porphyrin and water molecules.^50,51

Fig. 3 Potential of mean force of the water molecules derived as a function of Zn–O_w and Cu–O_wg(r) functions in (a) ZnP and (b) CuP, respectively.

The preferred orientation of O_w atoms towards the metals has a competing effect from the H_w atoms, thus lowering the H-bond expectancy between the water hydrogens and the nitrogen atoms of these solutes. This effect along with the direct coordination of water molecules to Zn(II) and Cu(II) was quantified in the form of dynamical properties by calculating the mean residence time (MRT,τ) for potential H-bonds between solute and solvent based on the minimum cut-off distance of about 2.5 Å.⁵²

The MRT based on the direct method provides values for the time t* (a minimum time span taken by the exchanging ligand between hydration shells) of 0.0 and 0.5 picosecond that may be denoted as τ_0.0 and τ_0.5.^3,4,52Table 1 lists the estimates of MRTs for the hydrated ZnP and CuP in comparison to the dynamical data for pure water determined employing a similar technique (1.3 ps).³⁸ The low MRT values for H-bonds in the hydrated metal-bound porphyrins is a direct result of the presence of the metal ions that preferentially attract water oxygens. This behavior is in contrast to that observed in the hydrated metal-free porphyrin demonstrating much stronger N–H_w bonds. The dynamical data for the Zn–O_w coordination indicate no ligand exchange between hydration shells which is in fair agreement with the strong coordination observed in the corresponding g(r) profile, thus confirming the penta-coordinated structure of ZnP in aqueous solution. In contrast to the ZnP molecule, water exchanges take place at the copper site as illustrated in the distance plot, in agreement with the native four-coordinate structure of CuP in aqueous solution (cf.Fig. 4). The ligand exchange dynamics in the case of CuP also correspond to that of an earlier simulation study on Cu(II) in aqueous solution as well as to experimental data.^50,53,54 The MRT values for the H-bond dynamics were computed as 0.5 and 0.4 ps for the ZnP and CuP molecules which are close to that in the case of hydrated magnesium bound porphyrin but the ligand exchange data between 0 and 1 demonstrate weak H-bonding. Variations in the structural and dynamical properties of aqueous metalloporphyrins could also be due to a combined effect of the inherent flexibility of the porphyrin ring, a competing behavior of the water oxygen for metal coordination and water exchanges between hydration shells.^55,56

Table 1 Comparative H-bond dynamics for the water molecules interacting with nitrogen atoms in ZnP and CuP molecules and in pure water

H-bond (system)	CN	N ^0.5ex	N ^0.0ex	τ 0.5 (ps)
N–Hw (metal-free porphyrin)^13	∼1	∼1	264	4.7
N–Hw (MgP)^13	0.1	—	39	0.6
N–Hw (ZnP)	0.1	—	36	0.5
N–Hw (CuP)	∼0.1	—	34	0.4
H2O^38	4.1	20	131	1.3

---

Fig. 4 Metal–O_w distance plots for (a) ZnP and (b) CuP; dashed lines show boundaries of hydration layers.

The substantially decreased affinity observed in water coordination towards ZnP and CuP is in line with the results or previous QM/MM MD simulation studies of copper(II) and zinc(II) ammine complexes. When compared to the respective aquo-complexes^51,57 a substantial weakening of the ion–water interaction along with a decrease in the coordination number upon an increase of the number of ammonia ligands has been observed.^{30,31,33,34,36}

The reduced conformational flexibility of the porphyrin systems arising due to metal coordination is further investigated by a vibrational analysis. The prediction of the respective power spectra is based on the Fourier transform of the velocity autocorrelation function.^58,59 The resulting vibrational modes associated with the hydrated ZnP and CuP systems excluding all bonds involving hydrogen atoms are presented in Fig. 5. It was reported earlier that the incorporation of divalent metals into the porphyrin ring leads to a spectral shift from 940 cm⁻¹ representing the vibrational frequency for metal-free porphyrin to a higher frequency. The ZnP and CuP molecules also demonstrate the property of spectral shift which was deduced from the peaks between 119 and 1866 cm⁻¹. This characteristic property of metalloporphyrins is attributed to the strong metal–ligand interaction, and was previously observed in a similar QMCF-MD simulation study of magnesium-bound porphyrin as well.^13,60

Fig. 5 Power spectra in cm⁻¹ of the ZnP and CuP molecules in aqueous solution.

The porphyrins also display residual hydrophilicity/hydrophobicity which seems to vary upon binding to different metal ions. The evaluation of the hydrophilic/hydrophobic character of hydrated ZnP and CuP enabled the determination of the relative H-bonding potential between the porphyrin nitrogens and water hydrogens. This was accomplished by the estimation of the time-averaged solvent accessible surface area (SASA) using VMD⁶¹ (cf.Table 2). Comparison of the reported SASA for the ligand-free and magnesium-bound porphyrin with those calculated for the ZnP and CuP provides a complete picture of the influence of different metal ions on H-bond formation in the hydrated metalloporphyrins. The relative magnitudes of the SASA for these solutes show that the CuP molecule displays the highest hydrophobic character among all porphyrin systems listed in the table. This further confirms that the CuP has low capacity for H-bond formation with water molecules which also correlates with the structural and dynamical properties evaluated for these metalloporphyrins, thus establishing that the complexation of porphyrins to metal ions alters the H-bonding potential.¹³

Table 2 Time-averaged solvent accessible surface area in nm². Hydrophobic and hydrophilic surfaces are decomposed to show their individual contribution

Porphyrins	Hydrophobic	Hydrophilic	Total
a Values are computed with standard deviations.
Metal-free porphyrin^13	4.58	0.29	4.87
Magnesium porphyrin^13	4.62	0.21	4.83
Zinc porphyrin^a	4.68 ± 0.11	0.18 ± 0.04	4.87 ± 0.11
Copper porphyrin^a	4.73 ± 0.20	0.20 ± 0.11	4.90 ± 0.19

---

4 Conclusion

In this study, detailed insight into the hydration behavior of ZnP and CuP complexes has been provided, stressing the capability of the advanced ab initio QMCF-MD formalism for the investigation of bioinorganic complexes in aqueous solution. The inclusion of quantum mechanics in the molecular dynamics framework to describe all inter- and intramolecular interactions between solute and solvent makes the technique highly competent for the evaluation of the structural and dynamical properties of complex chemical systems in aqueous solution. The results show that metalloporphyrins in aqueous solutions demonstrate very characteristic hydration patterns which were analyzed via a number of structural and dynamical parameters. Besides the evaluation of the structure and dynamics of water molecules, the computations of free-energies of water binding to the metal ions and relative hydrophilicity/hydrophobicity of the aqueous metalloporphyrins are also of considerable significance which would indeed be helpful for understanding the increased role of metal-bound porphyrins in aqueous biological systems. The simulation results also encourage us to unravel the dynamical features of several other more complex metallo-organic systems by using both theoretical and experimental methods.

Acknowledgements

Financial support for this work by the ATC-AIC-APC project from BMWF ASEA-UNINET Austria for Syed Tarique Moin is gratefully acknowledged. This work was supported by the Austrian Ministry of Science, Research and Economy BMWFW UniInfrastrukturprogramm as part of the Research Focal Point Scientific Computing at the University of Innsbruck. This work was supported by the ASEAN-European Academic University Network (ASEA Uninet) funded by the Austrian Federal Ministry of Science, Research and Economy (BMWFW).

References

1. B. M. Rode and T. S. Hofer, Pure Appl. Chem., 2006, 78, 525–539.
2. H. Ohtaki and T. Radnai, Chem. Rev., 1993, 93, 1157–1204.
3. A. J. Lock, S. Woutersen and H. J. Bakker, J. Phys. Chem. A, 2001, 105(8), 1238–1243.
4. A. J. Lock, S. Woutersen and H. J. Bakker, Femtochemistry and Femtobiology, Word Scientific, Singapore, 2001.
5. V. Sundström, Annu. Rev. Phys. Chem., 2008, 59, 53–77.
6. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
7. J. Hutter, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 604–612.
8. A. Warshel and M. Levitt, J. Mol. Biol., 1976, 103, 227–249.
9. B. M. Rode, T. S. Hofer, B. R. Randolf, C. F. Schwenk, D. Xenides and V. Vchirawongkwin, Theor. Chem. Acc., 2006, 115, 77–85.
10. T. S. Hofer, A. B. Pribil, B. R. Randolf and B. M. Rode, Adv. Quantum Chem., 2010, 59, 213–246.
11. A. K. H. Weiss and T. S. Hofer, RSC Adv., 2014, 3, 1606–1635.
12. M. Hitzenberger, T. S. Hofer and A. K. H. Weiss, J. Chem. Phys., 2013, 139, 114306.
13. S. T. Moin and T. S. Hofer, Mol. BioSyst., 2014, 10, 117–127.
14. S. Shaik and S. P. D. Visser, Cytochrome P450, Springer, 2005, pp. 45–85.
15. J. A. Shelnutt, J. Porphyrins Phthalocyanines, 2000, 4, 386–389.
16. C. J. Martin, C. L. Werntz and A. M. Ducatman, Occup. Med., 2004, 54, 587–591.
17. P. Fan, M. Valiev and K. Kowalski, Chem. Phys. Lett., 2008, 458, 205–209.
18. N. Govind, M. Valiev, L. Jensen and K. Kowalski, J. Phys. Chem. A, 2009, 113, 6041–6043.
19. R. F. Pasternack, L. Francesconi, D. Raff and E. Spiro, Inorg. Chem., 1973, 12, 2606–2611.
20. O. Schalk, Y. Liang and A. Unterreiner, Z. Phys. Chem., 2013, 227, 35–48.
21. H. Imai, K. Misawa, H. Munakata and Y. Uemori, Chem. Pharm. Bull., 2008, 56, 1470–1472.
22. A. H. Church, Philos. Trans. R. Soc. London, 1869, 159, 627–636.
23. S. E. Palmer and E. W. Baker, Science, 1978, 201, 49–51.
24. A. F. McDonagh, J. Pediatr., 2011, 158, 160–164.
25. J. Li, Y. Wei, L. Guo, C. Zhang, Y. Jiao, S. Shuang and C. Dong, Talanta, 2008, 76, 34–39.
26. T. Kaufmann, B. Shamsai, R. Lu, R. Bau and G. M. Miskelly, Inorg. Chem., 1995, 34, 5073–5079.
27. T. S. Hofer, H. Scharnagl, B. R. Randolf and B. M. Rode, Chem. Phys., 2006, 327, 31–42.
28. A. D'Incal, T. S. Hofer, B. R. Randolf and B. M. Rode, Phys. Chem. Chem. Phys., 2006, 8, 2841–2847.
29. S. T. Moin, T. S. Hofer, A. B. Pribil, B. R. Randolf and B. M. Rode, Inorg. Chem., 2010, 49, 5101–5106.
30. C. F. Schwenk and B. M. Rode, Phys. Chem. Chem. Phys., 2003, 5, 3418–3427.
31. M. Q. Fatmi, T. S. Hofer, B. R. Randolf and B. M. Rode, Phys. Chem. Chem. Phys., 2006, 8, 1675–1681.
32. T. S. Hofer, B. R. Randolf, A. A. Shah, B. M. Rode and I. Persson, Chem. Phys. Lett., 2007, 445, 193–197.
33. M. Q. Fatmi, T. S. Hofer, B. R. Randolf and B. M. Rode, J. Phys. Chem. B, 2007, 111, 151–158.
34. M. Q. Fatmi, T. S. Hofer, B. R. Randolf and B. M. Rode, Phys. Chem. Chem. Phys., 2008, 112, 5788–5794.
35. T. S. Hofer, B. R. Randolf, B. M. Rode and I. Persson, Dalton Trans., 2009, 1512–1515.
36. M. Q. Fatmi, T. S. Hofer and B. M. Rode, Phys. Chem. Chem. Phys., 2010, 12(33), 9713–9718.
37. H. Lin and D. G. Truhlar, Theor. Chem. Acc., 2007, 117, 185–199.
38. T. S. Hofer, B. M. Rode, A. B. Pribil and B. R. Randolf, Adv. Inorg. Chem., 2010, 62, 143–175.
39. T. Sakwarathorn, S. Pongstabodee, V. Vchirawongkwin, L. R. Canaval, A. O. Tirler and T. S. Hofer, Chem. Phys. Lett., 2014, 595, 226–229.
40. A. O. Tirler and T. S. Hofer, J. Phys. Chem. B, 2014, 118(45), 12938–12951.
41. D. Xenides, B. R. Randolf and B. M. Rode, J. Chem. Phys., 2005, 122, 4506–4515.
42. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91(24), 6269–6271.
43. S. Chatterjee, P. G. Debenedetti, F. H. Stillinger and R. M. Lynden-Bell, J. Chem. Phys., 2008, 128(12), 124511.
44. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Clarendon press, Oxford, 1987.
45. H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura and T. Head-Gordon, J. Chem. Phys., 2004, 20(20), 9665–9678.
46. V. Kräutler, W. F. van Gunsteren and P. H. Hünenberger, J. Comput. Chem., 2001, 22(5), 501–508.
47. M. Q. Fatmi, T. S. Hofer and B. M. Rode, Phys. Chem. Chem. Phys., 2010, 12, 9713–9718.
48. P. Muller, Pure Appl. Chem., 1994, 66, 1077–1184.
49. D. Chandler, Introduction to modern statistical mechanics, Oxford University Press, 1987.
50. S. T. Moin, T. S. Hofer, A. K. H. Weiss and B. M. Rode, J. Chem. Phys., 2013, 139, 014503.
51. M. Q. Fatmi, T. S. Hofer, B. R. Randolf and B. M. Rode, J. Chem. Phys., 2005, 123, 054514.
52. T. S. Hofer, H. T. Tran, C. F. Schwenk and B. M. Rode, J. Comput. Chem., 2004, 25(2), 211–217.
53. D. H. Powell, L. Helm and A. E. Merbach, J. Chem. Phys., 1991, 95, 9258–9265.
54. L. Helm and A. E. Merbach, Coord. Chem. Rev., 1999, 187, 151–181.
55. W. R. Scheidt and Y. J. Lee, Metal Complexes with Tetrapyrrole Ligands I, Springer, 1987, pp. 1–70.
56. M. O. Senge, Chem. Commun., 2006, 243–256.
57. C. F. Schwenk and B. M. Rode, J. Chem. Phys., 2003, 119, 9523–9531.
58. U. Balucani, J. P. Brodholt and R. Vallauri, J. Phys.: Condens. Matter, 1996, 8, 6139–6144.
59. R. N. Bracewell, Sci. Am., 1989, 260, 86–89.
60. L. J. Boucher and J. J. Katz, J. Am. Chem. Soc., 1967, 89, 1340–1345.
61. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.

---