Reactivity of arsenoplatin complex versus water and thiocyanate: a DFT benchmark study

Seven different density functionals, including GGAs, meta-GGAs, hybrids and range-separated hybrids, and considering Grimme’s empirical dispersion correction (M06-L, M06-2X, PBE0, B3LYP, B3LYP-D3, CAM-B3LYP, ωB97X) have been tested for their performance in the prediction of molecular structures, energies and energy barriers for a class of newly developed antitumor platinum complexes involving main group heavy elements such as arsenic. The calculated structural parameters, energies and energy barriers have been compared to the available experimental data. The results show that range-separated hybrid functionals CAM-B3LYP and ωB97X give good results in predicting both geometrical parameters and isomerization energies and barrier heights and are promising new tools for the theoretical study of novel platinum(II) arsenic compounds.


Introduction
Since the discovery of the antitumor activity of cisplatin [1], transition metal complexes, an essential group of compounds in diverse fields of chemistry, have been gaining an increased importance in medicinal chemistry [2][3][4][5]. Most of the studies in bioinorganic chemistry have been devoted to the structure-based optimization of cisplatin or its analogs to gain more potent and less toxic anticancer drugs [6][7][8]. The well-known toxicity of cisplatin is seemingly ascribed to the facile in vivo substitution of the chloride ligands by endogenous targets. Indeed, their replacement with less reactive chelate ligands in carboplatin and oxaliplatin has successfully yielded compounds with a higher therapeutic index [9][10][11][12].
Another major hurdle for the use of platinum compounds in cancer treatment is due to the intrinsic and acquired resistance to this class of drugs. Several innovative strategies have been developed to overcome this issue, such as the design of new types of platinum complexes like trans Pt(II) compounds, polynuclear Pt(II) compounds, Pt(IV) prodrugs or new nanomolecular or nanoparticle drug delivery approaches. Inspired by the identification of a new Pt-As bond in a nanoparticulate delivery system for loading arsenic trioxide (As 2 O 3 )-an inorganic compound with a wellknown anticancer activity-and cisplatin, a new complex Cl-Pt(NH = C(CH 3 )-O) 2 As(OH) 2 , arsenoplatin (AP), has been recently developed showing a promising anticancer activity in drug-resistant cell lines.
Indeed, arsenoplatin, an adduct bearing a square-planar Pt(II) center resembling the coordinative structure of cisplatin coordinated by an arsenic trioxide moiety [13], shows a synergistic effect of its two components, displaying an antitumor activity superior to both parent drugs As 2 O 3 and cisplatin in the majority of cell lines [13,14].
From a structural point of view, the AP complex and its derivatives disclose notable features: (i) The As(III) ligand acts in AP concomitantly as both a Lewis base through the sigma donation of a As lone pair to Pt and as a Lewis acid through the acceptance of two oxygen lone pairs from the imido ligands; (ii) As(III) displays a trigonal-bipyramidal geometry with two OH groups out of the Pt coordination plane (Scheme 1). This behavior is quite different from that shown by arsenous acid in other heterometallic complexes with transition metals, such as those with Ni(II) and Pd(II), where the As(III) center operates as a Lewis base with, respectively, a trigonal-pyramidal geometry [15], and a more tetrahedral character [16,17].
Arsenoplatin is quite stable in solution, and several studies have shown that it readily undergoes chloride exchange reactions at the Pt(II) center with retention of the Pt-As core and stereochemistry. Due to the strong trans effect of the arsenic moiety, the substitution of chloride is generally faster than in cisplatin. For instance, AP undergoes hydrolysis with a half-life t 1/2 = 2 h at 37 °C, significantly faster than cisplatin showing a half-time of ca. 2.5-3 h, and the substitution by small ions such as thiocyanate occurs immediately at room temperature [14].
The mechanism of action of AP is still unclear, and different hypotheses have been proposed on the basis of its reactivity with cell DNA and model proteins [13]. The presence of an easily removable chloride on Pt(II), as in most Pt-based anticancer compounds, suggests that AP could easily undergo hydrolysis followed by the facile binding of the resulting aquo complex to a DNA nucleobase to give monofunctional adducts. However, the lack of a second easily removable ligand in cis position prevents the formation of the intrastrand crosslinks between two adjacent nucleobases, assumed to be the critical lesions responsible for cisplatin cytotoxicity. Moreover, AP has been shown to bind to histidine residues of model proteins, maintaining the Pt-As core, suggesting that protein binding could be involved in its mechanism of action. Finally, it has been found that in the biological milieu, several reactions can occur and trigger the destabilization and, eventually, the dissociation of the Pt-As bond that eventuates in the release of cytotoxic and therapeutically relevant As(III) species.
A deeper insight at the molecular level into the reactivity of AP toward potential biomolecular targets would be of particular benefit in elucidating the mechanism of action of this new class of compound.
In this context, the structure-based control of substitution reaction occurring on either cisplatin derivatives or, more generally, on analogous Pt(II) square-planar complexes is pivotal in the design of new active compounds. Density functional theory (DFT) is a proper tool to investigate reaction mechanisms of transition metal complexes due to their usually ample size, thus resulting to be particularly suitable in the study of potential metallodrugs, and has been of great benefit to elucidate the initial steps of the mechanism of action of cisplatin [18][19][20][21][22][23][24][25][26] and Pt-based drugs [27][28][29][30][31][32]. Although most of these DFT approaches have been based on a few widespread density functionals developed in the 1990s, in particular B3LYP which has been shown to give reasonably good results for this class of compounds [18-20, 22, 25-30, 32], in the last decade these methods have been recognized to show some pitfalls-especially when the complexity of the considered systems increases and electron delocalization or non-covalent interactions play an important role-and new density functionals have been developed to improve the accuracy of DFT calculations.
It is nowadays commonly accepted that evaluation of the performance of density functionals through assessment on different realistic chemical models is a crucial step prior to the investigation of new systems. Though there is a plethora of studies on the performance of DFT functionals [33][34][35] on main group elements, only more recently density functionals have been benchmarked in many studies involving big datasets to assess their predictive capability in describing the most common transition metal reactions [36][37][38]. However, the accuracy of the wide range of newly developed density functionals has been found to depend both on the kind of considered transition metals (early, late, first or second and third row, etc.) and on the kind of properties taken into account (geometrical structure, reaction energy, energy barriers, excitation energies, etc.), so that no universal rule can be given for the choice of the "best" functionals. While many and extensive benchmarking studies have been performed on compounds containing either transition metals or main group heavy elements, assessments on the  [39][40][41][42][43]. A quick search in the literature for DFT functionals more frequently used in theoretical investigations on Pt-based anticancer complexes confirmed that B3LYP is still the most commonly employed functional both for geometrical optimizations and for the subsequent electronic energy calculations, despite in the last years more and more theoretical investigations employed more recent functionals, such as the Minnesota ones, though often without a systematic assessment study. A DFT study has recently been carried out on chloride substitution by nucleobases in arsenoplatin employing the B3LYP functional, although a preliminary benchmark with M06, M06L and MPW1PW functional was carried out to assess only the geometry structure [27].
In the present study, we assessed seven among the routinely employed and recently developed density functionals for their capability to predict the structure and reactivity of Cl-Pt(NH = C(CH 3 )-O) 2 As(OH) 2 , arsenoplatin (AP), in order to point out the adequate level of theory to use in future investigations on its mechanism of action. In particular, we considered the widely used hybrid GGA functional B3LYP [44,45] and its variant including Grimme's empirical dispersion correction B3LYP-D3 [46], the meta-GGA functional originally designed to study transition metal complexes M06-L [47], the Minnesota hybrid meta-GGA functional M06-2X [48], the one-parameter hybrid GGA functional PBE0 [49] and two more recent range-separated hybrid functionals which have been shown to provide good results for transition metal systems CAM-B3LYP [50] and ωB97X [51], which are known to give a good description of geometries and reaction profiles for transition or heavy metal-containing compounds [52][53][54][55][56], including structure and reactivity of transition metal complexes [57,58], and particularly antitumor Pt-based compounds [32,[59][60][61]. We also considered the influence of the use of three basis sets differing in their size and for the number of polarization and diffuse functions.
The performance of the selected functionals was checked on one side against the geometrical parameters from the X-ray structures of arsenoplatin and its Pt-substituted derivatives X-Pt(NH = C(CH 3 )-O) 2 As(OH) 2 (X = Cl − , SCN -) and on the other side against the experimentally determined semi-quantitative data for the energy barrier for the chlorine substitution in arsenoplatin by water and by thiocyanate SCNion. We also examined the accuracy of the considered functionals for predicting the quantitative enthalpy and free energy data for the thiocyanate to isothiocyanate isomerization determined by NMR spectroscopy [14]. The thiocyanate ion is omnipresent in mammalian biology, reaching up to mM concentrations in extracellular fluids [62], and also exhibits strong antimicrobial and anti-inflammatory effects [63]. This anion can coordinate to Pt(II) center of arsenoplatin with either sulfur or nitrogen, thus presenting an easy model for investigation of arsenoplatin interactions with various S-(cysteine, methionine) and N-containing (histidine) biomolecular targets, see Scheme 1.
The outcomes of this validation study would be of great value in future calculations on the prediction of the AP reactivity with potential biological targets and contributing to shed light on the molecular events involved in either its activation or target binding.
The C-PCM continuum solvent method was used to describe the solvation [65]. It has recently been shown to give considerably smaller errors than other continuum models for aqueous free energies of solvation for neutrals, cations and anions and to be particularly effective for the computation of solution properties requiring a high accuracy of solution free energies [73]. Solvation free energies, taken as the difference between the solution energies and the gas phase energies, were added to the gas phase enthalpies and free energies values to obtain the corresponding values in the aqueous solution. Experimental values have been used for small molecules and ions, − 6.3 kcal/mol and − 77.0 kcal/ mol for the solvation energies of water and chlorine anion, respectively [74].
For the geometry structure benchmark, geometrical optimizations have been carried out in the gas phase with all the three considered basis sets. For the energy benchmark, geometry optimizations were carried out in gas and in solution with basis sets BS1 and BS2 followed by single-point electronic and solvation energy calculations carried out with basis set BS3.
Frequency calculations were performed to verify the correct nature of the stationary points and to estimate zero-point energy (ZPE) and thermal corrections to thermodynamic properties. Intrinsic reaction coordinate (IRC) calculations were employed to locate reagents and products minima connected with the transition states for each considered reaction step.
The barrier for AP aquation was roughly estimated from the experimental data on reaction half-time by means of In order to evaluate the errors, we calculate the mean unsigned error (MUE), which is the average of two standard deviation-normalized MUE values for distances and angles:

Results and discussion
In this study, we investigated thermodynamics and kinetics of the substitution reactions in which the chloride ligand on AP is replaced by either water or thiocyanate ion, to identify the best computational approach among the density functional theory schemes most often applied in bioinorganic chemistry. A preliminary assessment of the considered density functionals was performed by comparing optimized and experimental geometries of AP and its chlorine-substituted derivatives X-Pt(NH = C(CH 3 )-O) 2 As(OH) 2 (X = Cl − , SCN − ), whose molecular geometries have been experimentally characterized by X-ray analysis [14] (Fig. 1). A reliable description of the potential energy surface, leading to an accurate prediction of minimum and transition state structures, is indeed a prerequisite for the benchmark of thermodynamics and kinetics parameters. Significant bond distances and bond angles calculated with all considered functionals and basis sets are compared with the corresponding experimental values in Table 1, where the accuracy of each density functional and basis set in the geometry prediction was evaluated through the mean unsigned errors for bond lengths only or with the inclusion of bond angles. structure of this class of compound. The slightly but significantly higher errors found for the BS1 set, especially for the metal ligands bond lengths, suggest that even when the large size of the system requires a small basis sets-as for instance in the AP interaction with proteins-at least the atoms directly bound to platinum should be described with a better basis.
Thermodynamics and kinetics of the chloride substitution in the AP complex with water, i.e., the aquation reaction, were investigated using all considered functionals and basis sets BS1 and BS2, in gas phase and in water solution, to locate the geometries of reactants, products, possible intermediates and transition states and carrying out singlepoint calculations of the electronic energy with the larger basis set B3. Solvation energy was calculated both with the basis set used for the optimizations, BS1 or BS2, and with the larger basis set BS3 employed in the single points. To disentangle the effect of the considered functionals on the electronic energy from that on the geometry, we also carried out single-point calculations of the electronic energy with the larger basis set B3 for all the considered density functionals on the geometries optimized with BS1 basis set and the best performing ωB97X functional. The well-established associative interchange mechanism was assumed for the chloride substitution that implies a transition state with a distorted trigonal-bipyramidal coordination of the Pt center (Scheme 2). Thus, our assessment of AP aquation (Table 2) was carried out through optimization and subsequent singlepoint energy calculation on separated reactants (R), reactantadduct (RA), transition state (TS), product-adduct (PA) and separated products (P). The reaction enthalpies and free energies were calculated as the difference between R and P, while the activation enthalpies and free energies were calculated as the difference between corresponding values for TS and the lowest between R and RA.
The reaction and activation free energies for the AP aquation, calculated for all considered levels of theory, are reported in Table 2, while the corresponding enthalpies can be found in Table S1. The only available kinetic experimental data for the AP aquation are its half-time, t 1/2 = 2 h at 37 °C and 4 mM Cl − [14,75], from which an approximate estimate of the experimental activation free energy of 23.9 kcal/mol can be obtained (see Method section).
Reasonable values of the activation free energy have been obtained for the majority of the considered density functional/basis set/optimization phase combinations. We first notice that, in spite of its relative inaccuracy in reproducing some bond length, the results obtained using geometries optimized with basis set BS1 are comparable with those obtained using geometries optimized with basis set BS2. Analogously, using the more accurate geometry evaluated with the best performing ωB97X functional in the single-point calculations with the different functionals does not significantly affect the calculated activation energies, confirming that the accuracy of reaction and activation energies are not very sensible to small geometry differences. When we compare the results for the activation free energy calculated through energy optimization with the more reliable basis set BS2 in the gas phase and in solution, we notice a small decrease of the free energy barrier for all the functionals with significantly better results for the optimization in solution for almost all other functionals. Focusing on the results obtained through energy optimization in solution, we first notice that the results are little affected from the use of BS2 or BS3 in the calculation of the solvation energy, with differences within 0.3 kcal/mol, in agreement with the literature, showing that solvation energy is not much sensitive to basis set [76], and we will then consider only the results obtained using BS3. We then observe that M062X gives the worst results with an error of − 5.3 kcal/mol, while the other functionals give all results within 3.0 kcal/mol within the experimental value. Among them, ωB97X, B3LYP-D3 and, to a lesser extent, CAM-B3LYP yield the most accurate results with errors of − 0.2, − 0.1 and + 1.4 kcal/mol, respectively. When we consider the free energy barriers from optimization in the gas phase, the best results are obtained with the same ωB97X, B3LYP-D3 and CAM-B3LYP functionals, with slightly larger errors of 1.0, 1.5 and 2.7 kcal/mol, respectively. The substitution of chloride on AP with thiocyanate ion leading to either AP-NCS or AP-SCN isomers was then investigated, assuming the same associative interchange mechanism considered for AP aquation (Scheme 2). We employed only the CAM-B3LYP/BS3//CAM-B3LYP/BS2 (C-PCM) level of theory, which gave the best results for the chloride substitution by water as discussed above.
Although no experimental quantitative kinetic data are available, by taking into account that in water this process occurs immediately at room temperature [14], we can roughly estimate that the activation free energy for this reaction is presumably lower than 19 kcal/mol, corresponding to a half-time of ca. 3 s at 37 °C.
The results of our calculations for both the Pt-SCN and Pt-NCS isomers, reported in Fig. 2, indicate the formation of the thiocyanate-substituted AP complex to be exothermic and exergonic, presumably, because of the more pronounced soft character of this nucleophile compared to the leaving chloride. As shown, the reported activation free energy (RA-> TS) of 17.8 kcal/mol for the most favorable path leading to the Pt-SCN isomer is fully consistent with a fast chloride substitution, in semi-quantitative agreement with the available experimental evidences [14].
T h e r e l a t i v e s t a b i l i t y o f P t ( S C N ) (NH = C(CH 3 )-O) 2 As(OH) 2 , AP-SCN, with respect to its coordinative isomer Pt(NCS)(NH = C(CH 3 )-O) 2 As(OH) 2 , AP-NCS, for which accurate experimental data are available, was then assessed against all the seven density functionals, and the results are reported in Table 3. The van′t Hoff plot of the temperature-dependent NMR spectra shows for the AP-SCN ⇆ AP-NCS equilibrium the thermodynamic parameters ∆H° = −3.8 kcal mol −1 , ∆S° = −13.7 cal mol −1 K −1 and ∆G° = 0.3 kcal/mol, indicating that the AP-NCS isomer in solution is enthalpically slightly favored, while the AP-SCN isomer is barely favored by free energy [14]. We focused on the enthalpy difference between the two isomers since the free energy values include the entropy difference for which the statistical thermodynamics approach employed in our calculation, based on ideal gas model with fully uncoupled translation, rotation, vibration and electronic motion, is intrinsically not accurate enough to accurately reproduce the experimental ∆S° value of few cal mol −1 K −1 . First of all we notice that values of ∆H° in reasonable agreement with the experimental value (within slightly more than 3 kcal/mol) were obtained for almost all functionals. As expected, the results of ∆G° are more irregular and less accurate presumably because of the inaccuracy in the calculation of entropy, which appears to be systematically underestimated.
Focusing on the more reliable results obtained through energy optimization with basis set BS2 in solution and using the largest basis set BS3 in the calculation of solvation energy, the most accurate estimates of AP-SCN ⇆ AP-NCS interconversion enthalpy were obtained with the CAM-B3LYP and B3LYP functionals, with − 1.2 and − 1.5 kcal/ mol deviations, respectively, the other functionals showing errors in the range 2.8-3.5 kcal/mol.

Conclusions
In this work we assessed the performance of seven different density functionals, including GGAs, meta-GGAs, hybrids and range-separated hybrids, and considering Grimme's empirical dispersion correction (M06-L, M06-2X, PBE0, B3LYP, B3LYP-D3, CAM-B3LYP, ωB97X) in predicting molecular structures, energies and energy barriers for a class of newly developed antitumor platinum complexes involving main group heavy elements such as arsenic.
In particular we considered structure and reactivity of Cl-Pt(NH = C(CH 3 )-O) 2 As(OH) 2 , arsenoplatin (AP), a novel highly successful agent for the treatment of cancer, in order to point out the most adequate level of theory to use in future investigations on its mechanism of action since this molecule has been demonstrated to display an antitumor activity superior to both its parent drugs.
The accuracy of the considered functionals was checked against the prediction geometrical parameters of some structurally characterized AP derivatives and the reaction and activation free energies for its reaction with water and the thiocyanate SCN − ion, for which semi-quantitative thermodynamic and kinetic data are available. We also examined the accuracy of the considered functionals for predicting the precise enthalpy and free energy data for the thiocyanate to isothiocyanate isomerization determined by NMR spectroscopy. The results show that range-separated hybrid CAM-B3LYP and ωB97X and the dispersion-corrected B3LYP-D3 functionals, when used with double zeta basis sets including polarization and diffuse functions on the main group atoms in the optimization and solvation energy calculations and triple-zeta basis sets including polarization and diffuse functions on hydrogen and metal atoms, give the best results in predicting both geometrical parameters and isomerization energies and barrier heights and are promising new tools for the theoretical study of novel platinum(II) arsenic compounds. Indeed, bond distances have been calculated within 0.01-0.02 Å from the experimental values and barriers for chloride substitution within 1 kcal/mol. However, the traditional B3LYP or its dispersion-corrected B3LYP-D3 variant and PBE0 hybrid