Prediction of stability constants for cu2+ complexes with organic fluorescent ligands using thermodynamic cycle in combination with dft theory and smd solvent model

Hue University Journal of Science: Natural Science Vol. 129, No. 1D, 15–23, 2020 pISSN 1859-1388 eISSN 2615-9678 DOI: 10.26459/hueuni-jns.v129i1D.5947 15 PREDICTION OF STABILITY CONSTANTS FOR Cu2+ COMPLEXES WITH ORGANIC FLUORESCENT LIGANDS USING THERMODYNAMIC CYCLE IN COMBINATION WITH DFT THEORY AND SMD SOLVENT MODEL Mai Van Bay1,2, Nguyen Khoa Hien3, Hoang Kim Thanh4, Pham Cam Nam5, Duong Tuan Quang1* 1 University of Education, Hue University, 34 Le Loi St., Hue, Vietnam 2 The University of Danang−University of Science and Education, 459 Ton Duc Thang St., Danang, Vietnam 3 Mientrung Institute for Scientific Research, VAST, 321 Huynh Thuc Khang St. Hue City, Vietnam 4 Danang University of Medical Technology and Pharmacy, 99 Hung Vuong St., Da Nang, Vietnam 5 The University of Danang−University of Science and Technology, 54 Nguyen Luong Bang St., Danang, Vietnam * Correspondence to Duong Tuan Quang (Received: 02 August 2020; Accepted: 14 August 2020) Abstract. Accurately predicting the stability constant (β) of the Cu2+ complex with organic fluorescent ligands provides an important basis to design molecular fluorescent sensors for selective detection of Cu2+. With appropriate reference complexes, the calculated stability constants are in good agreement with experimental values. The logβ values of the predicted stability constants of Cu2+ complexes with Calcein blue (H3Cb) and FluoZin-1 (H2Fz) are 13.33 (exp. 14.27) and 6.59 (exp. 6.01), respectively. More importantly, the results could be applied to the investigation of complexes. Keywords: fluorescent, stability constant, complex, thermodynamic cycle, DFT 1 Introduction The complex interaction between metal ions and organic fluorescent ligands is one of the important approaches to design fluorescent sensors for the detection of metal ions as well as other analytical species, such as anions and biothiols through complex exchange reactions [1]. Besides a requirement of optical properties, a fluorescent sensor for detection of metal ions needs sufficiently strong interaction with the target metal ion, usually through complexation reactions [2]. Predicting this complex formation is very necessary for designing sensors [3]. Understanding the complexation characteristics of metal ions in aqueous solutions is the basis for predicting and controlling the behavior of metal ions in the environment, biological systems, and other industrial applications [4]. The key point of this problem is to estimate the affinity of metal ions with the ligands through the stability constant of complexes [5, 6]. By definition, the stability constant of complexes between metal ion M and ligand L ( βMLn ) in the solution is determined according to Eq. (1). M(aq) + nL(aq) = MLn(aq) βMLn = aMLn aMaL n (1) where "a" is the activity at equilibrium, and "aq" indicates the state in the aqueous solution. The stability constant can be approximately evaluated through concentration instead of activity; therefore, the stability constant can be expressed as in Eq. (2). Mai Van Bay et al. 16 βMLn ≈ [MLn] [M][L]n (2) and the relationship between βMLn and Gibbs free energy of complex formation reaction (∆𝐺aq 0 ) is expressed in Eq. (3). ∆𝐺aq 0 = −𝑅𝑇 ln βMLn (3) In principle, it is possible to determine the stability constant of a complex by calculating the theoretical Gibbs energy of the reaction in solutions [6]. However, there are still various major obstacles in accurately assessing the Gibbs free energy value of reaction solutions [7]. This may be related to determining the true form of metal ions in aqueous solutions. For example, Cu2+ ions can exist in solution in possible stable forms such as [Cu(H2O)4] 2+ , or [Cu(H2O)5] 2+ , or [Cu(H2O)6] 2+ [8]. In addition, another difficulty is to assess the solvation energy of the substance in water. Especially for ions, when using common solvent models such as PCM (Polarizable Continuum Model) and COSMO (Conductor-like Screening Model), the calculated results are quite different from the experimental value [9]. A large number of stability constants for metal complexes have been experimentally determined, forming a database to serve relevant applications [10]. However, such data are not sufficient because numerous complexes exist between metal ions and ligands. There have been various attempts to theoretically predict the stability constants of metal complexes [6]. On the basis of the correlation between the stability constant and the properties of metal ions (i.e., ion radius, charge, electronegativity, and ionization potential), the equations for determining the stability constant from the empirical database are formulated [11, 12]. This approach is useful but requires a large amount of suitable empirical data that are suitable for the structure of the studied complex. So far, numerous research groups have made significant efforts in finding quantum computational methods to accurately predict the stability constants of metal complexes [13]. These studies show that a good control of the calculation models may lead to the results that are close to experimental values [14, 15]. In this study, we report an approach of using a thermodynamic cycle in a combination with the DFT theory and SMD solvent model to predict the stability constants of Cu2+ complexes with two organic fluorescent ligands, namely, Calcein blue (H3Cb) and FluoZin-1 (H2Fz) solutions. Two reference ligands, adipic acid (H2A) and iminodiacetic acid (H2Ia, are used. The chemical structure of the ligands is shown in Fig. 1. Fig. 1. Structural formula of fluorescent ligands and reference ligands Hue University Journal of Science: Natural Science Vol. 129, No. 1D, 15–23, 2020 pISSN 1859-1388 eISSN 2615-9678 DOI: 10.26459/hueuni-jns.v129i1D.5947 17 2 Methods 2.1 Thermodynamic methods Because Cu2+ forms 1:1 stoichiometry complexes with the four selected ligands [16-19], the metal complexes with ligands in a 1:1 molar ratio with the presence of water molecules are considered in this study. Therefore, the complexation is modeled in solution according to Eq. (4). [M(H2O)m](aq) 𝑥 + L(aq) 𝑦 ∆Gaq,ML 0 → [ML(H2O)n](aq) 𝑥+𝑦 + (𝑚 − 𝑛)H2O(aq) (4) Similar complexation of reference ligands (Lref) follows Eq. (5). [M(H2O)m](aq) x + Lref(aq) 𝑧 ∆Gaq,MLref 0 → [MLref(H2O)k](aq) 𝑥+𝑧 + (𝑚 − 𝑘)H2O(aq) (5) Eq. (4) and (5) are combined to give Eq. (6). [MLref(H2O)k](aq) x+z + L(aq) 𝑦 ∆Gaq 0 → [ML(H2O)n](aq) 𝑥+𝑦 + Lref(aq) z + (𝑘 − 𝑛)H2O(aq) (6) Three Eqs. (4), (5), and (6) can be combined to obtain Eq. (7). ∆𝐺aq 0 = ∆𝐺aq,ML 0 − ∆𝐺aq,MLref 0 (7) According to thermodynamics, ∆𝐺aq,ML 0 is related to the equilibrium concentration of substances according to Eq. (8). ∆𝐺aq,ML 0 = −𝑅𝑇 ln [[ML(H2O)n]][H2O] 𝑚−𝑛 [M(H2O)m][L] (8) Compared with the definition of stability constant, Eq. (9) is obtained. ∆𝐺aq,ML 0 = −𝑅𝑇 ln(βCuL [H2O] 𝑚−𝑛) (9) Similarly, for the reference complex, Eq. (10) is derived. ∆𝐺aq,MLref 0 = −𝑅𝑇 ln(βCuLref [H2O] 𝑚−𝑘) (10) Eqs. (7), (9), and (10) are combined to obtain Eq. (11). ∆𝐺aq 0 = −𝑅𝑇 ln ( βCuL βCuLref [H2O] 𝑘−𝑛) (11) Eq. (11) can be rewritten in the form of Eq. (12). log βCuL = − ∆𝐺𝑎𝑞 0 𝑅𝑇 ln(10) + log βCuLref − (𝑘 − 𝑛) log[H2O] (12) The concentration of water under the standard state is 55.56 M [20, 21]. From Eq. (12), it can be seen that using a reference complex with known empirical stability constant value (βCuLref) can avoid estimating the existing form of metal ions in the water environment ( [M(H2O)m] 𝑛+ ). Moreover, the choice of a complex reference structure that is more similar to the studied complex can partially eliminate the systematic error due to the calculation method [22]. βCuL can be determined from Eq. (12) by using the Gibbs energy of reaction in Eq. (6), which can be found from the thermodynamic cycle in Fig. 2. Mai Van Bay et al. 18 Fig. 2. The relationship between ∆𝐺aq 0 , ∆𝐺g 0, ∆𝐺solv 0 , and ∆𝐺g→aq 0 where “g” indicates the gas phase; ∆𝐺g→aq 0 is the Gibbs free energy change when transferring a mole of a substance from the standard condition in the gas phase (1 atm or 24.46 L.mol−1) into solution (1 L.mol−1). This change is equivalent to the process of compressing one mole of an ideal gas with a volume of 24.46 L to 1 L at 298.15 K, and the change of Gibbs free energy is calculated according to Eq. (13). ∆𝐺g→aq 0 = ∫ 𝑉d𝑝 1 24.46 = −𝑅𝑇 ∫ d𝑉 𝑉 1 24.46 = −𝑅𝑇 ln 1 24.46 = 𝑅𝑇 ln 24.46 (13) Thus, for a substance transferred from the gas phase to the solution, the Gibbs free energy is calculated following Eq. (14). 𝐺aq 0 = 𝐺g 0 + ∆𝐺solv 0 + 𝑅𝑇 ln 24.46 (14) And ∆𝐺aq 0 of a complex exchange reaction is determined from the thermodynamic cycle according to Eqs. (15) and (16). ∆𝐺aq 0 = (𝑘 − 𝑛)𝐺aq,H2O 0 + 𝐺aq,[CuL(H2O)n] 0 + 𝐺aq,Lref 0 − 𝐺aq,[CuLref(H2O)k] 0 − 𝐺aq,L 0 (15) ∆𝐺aq 0 = ∆𝐺g 0 + ∆∆𝐺solv 0 + ∆∆𝐺g→aq 0 (16) where ∆𝐺𝑔 0 is the Gibb free energy of reaction in the gas phase; ∆∆𝐺solv 0 is the solvation free energy of reaction, and ∆∆𝐺g→aq 0 is the Gibbs free energy for standard state change of reaction. 2.2 Computational methods All calculations were performed by using Gaussian16 software [23]. The optimized geometry of each substance is calculated at the theory level of PBE0/6-31+G(d) in the gas phase. Gibbs free energy in the solution of substances was calculated according to Eq. (17). where ε0 is the electronic structure energy obtained from the single-point energy calculation at PBE0/6311++G(d, p) from the corresponding optimized geometry at PBE0/6-31+G(d,p) in the gas phase. 𝐺corr is the thermal correction to free energy, which includes zero-point energy correction, and is determined at the same level as the geometry at PBE0/6-31+G(d). ∆𝐺solv is the change of solvation free energy calculated according to the experimentally parameterized model by Truhlar at M052X/631G(d)/SMD [24]. 𝐺s = 𝐺g + ∆𝐺solv + 𝑅𝑇 ln 24.46 = ε0 + 𝐺corr + ∆𝐺𝑠𝑜𝑙𝑣 + 𝑅𝑇 ln 24.46 (17) Hue University Journal of Science: Natural Science Vol. 129, No. 1D, 15–23, 2020 pISSN 1859-1388 eISSN 2615-9678 DOI: 10.26459/hueuni-jns.v129i1D.5947 19 3 Results and discussion 3.1 Determination of stable geometry for reference complexes The value of log of the two reference complexes [Cu(Ia)] and [Cu(Adp)] is 10.54 and 3.35, respectively [16-18]. The geometric structure of these reference complexes is needed for the calculation models. It is possible to predict their stable geometric structures according to Gibbs free energy change of the conversion between the existing forms of complexes according to Eq. (18). [CuLref]aq 𝑥+𝑧 + xH2Oaq ∆Gaq 0 → [CuLref(H2O)x]aq 𝑥+𝑧 (18) The more thermodynamically favorable the conversion reaction, the more stable [CuLref(H2O)x] complex (the more negative value of ∆𝐺aq 0 ). Possible stable geometries of the reference complexes are shown in Fig. 3. The calculated results of Gibbs free energy change for the conversion reactions of different complexes (Table 1) show that [Cu(Ia)] and [Cu(Adp)(H2O)] are the most stable forms of the two reference complexes. Therefore, these two complexes are used in subsequent calculations. Fig. 3. Stable geometry of the reference complexes Table 1. Gibbs free energy change for the conversion reaction of the reference complex forms Reaction ∆𝑮𝐚𝐪 𝟎 (𝐤𝐜𝐚𝐥.𝐦𝐨𝐥−𝟏) [𝐂𝐮(𝐈𝐚)]𝐚𝐪 + 𝐇𝟐𝐎𝐚𝐪 → [𝐂𝐮(𝐈𝐚)(𝐇𝟐𝐎)]𝐚𝐪 1.35 [𝐂𝐮(𝐈𝐚)]𝐚𝐪 + 𝟐 𝐇𝟐𝐎𝐚𝐪 → [𝐂𝐮(𝐈𝐚)(𝐇𝟐𝐎)𝟐]𝐚𝐪 8.58 [𝐂𝐮(𝐀𝐝𝐩)]𝐚𝐪 + 𝐇𝟐𝐎𝐚𝐪 → [𝐂𝐮(𝐀𝐝𝐩)(𝐇𝟐𝐎)]𝐚𝐪 −11.07 [𝐂𝐮(𝐀𝐝𝐩)]𝐚𝐪 + 𝟐 𝐇𝟐𝐎𝐚𝐪 → [𝐂𝐮(𝐀𝐝𝐩)(𝐇𝟐𝐎)𝟐]𝐚𝐪 −6.22 Mai Van Bay et al. 20 3.1 Determination of stability constant of Cu2+ complex with HCb2− and Fz2− The complex forms that were chosen for the investigation include [Cu(HCb)], [Cu(HCb)(H2O)], [Cu(HCb)(H2O)2] , [Cu(Fz)] , [Cu(Fz)(H2O)] , and Cu(Fz)(H2O)2] . The stable geometries of ligands and complexes are shown in Fig. 4 and 5. Fig. 4. Stable geometries of fluorescent ligands and reference ligands Fig. 5. Stable geometries of complex forms between Cu2+ and HCb2− or Fz2− Hue University Journal of Science: Natural Science Vol. 129, No. 1D, 15–23, 2020 pISSN 1859-1388 eISSN 2615-9678 DOI: 10.26459/hueuni-jns.v129i1D.5947 21 The calculated results in Table 2 indicate that the stability constants of Cu2+ complexes with HCb2− and Fz2− when using the reference complexes [Cu(Ia)] and [Cu(Adp)(H2O)] decrease in the following order [Cu(HCb)(H2O)] > [Cu(HCb)(H2O)2] > [Cu(HCb)] and [Cu(Fz)(H2O)] > [Cu(Fz)] > [Cu(Fz)(H2O)2]. Therefore, the most stable complexes between Cu2+ and HCb2− or Fz2− are [Cu(HCb)(H2O)] and [Cu(Fz)(H2O)]. The calculated stability constants of these two complexes are then compared with the experimental values. For the [Cu(HCb)(H2O)] complex, the empirical log βexp is 14.27 [19] and the calculated log βcalc when using the reference complexes [Cu(Ia)]) and [Cu(Adp)(H2O)] are 13.33 and 11.39, respectively. The results show that when using the reference complex [Cu(Ia)], the predicted stability constant of the Cu2+ complex with HCb2− is closer to the experimental value. This might be due to the similar structure of the chelating center [25]. Both ligands have a chelating center with three atoms containing nonbonding pairs of electrons on the two oxygen atoms in the two carboxyl groups (COO) and the nitrogen atom. For Adp2− reference ligands, the similarity between the two chelating centers is lower due to the absence of nitrogen atom in Adp2−. Hence, the calculated constant is less consistent with its experimental value. For the [Cu(Fz)(H2O)] complex, the empirical log βexp is 6.01 [19], and the calculated values of log βcalc when using the reference complexes [Cu(Ia)]) and [Cu(Adp)(H2O)] are 4.25 and 6.59, respectively. As a result, Adp2− is a better reference ligand than Ia2− for calculating the stability constant of the Cu2+ complex with Fz2−. There is a similarity in the chelating center between the Adp2− and Fz2− ligands. The geometric structures of the complexes in Figures 3 and 5 show that the chelating center Adp2− has three oxygen atoms containing nonbonding pairs of electrons in two carboxyl groups. Fz2− also has three oxygen atoms containing nonbonding pairs of electrons in two carboxyl and methoxy groups (OCH3), and a nonbonding pair of electrons in the nitrogen atom. However, unlike the nitrogen atom in Ia2-, the nitrogen atom in Fz2− is conjugated with the aromatic ring, reducing the density of the electron pair in nitrogen; therefore, it may reduce the ability to form a coordination bond between the nitrogen atom and Cu2+. Table 2. ∆𝐺aq 0 (kcal.mol−1) and predicted log β of different complex forms Reference complex [𝐂𝐮(𝐈𝐚)] [𝐂𝐮(𝐀𝐝𝐩)(𝐇𝟐𝐎)] ∆𝐺aq 0 log βref log βcalc ∆𝐺aq 0 log βref log βcalc [Cu(Cb)] 5.17 10.54 6.75 −4.31 3.35 4.77 [Cu(Cb)(H2O)] −1.43 10.54 13.33 −11.70 3.35 11.93 [Cu(Cb)(H2O)2] 6.37 10.54 9.35 −3.91 3.35 7.96 [Cu(Fz)] 10.07 10.54 3.17 −5.42 3.35 5.59 [Cu(Fz)(H2O)] 10.95 10.54 4.25 −4.42 3.35 6.59 [Cu(Fz)(H2O)2] 20.47 10.54 −0.99 5.09 3.35 1.36 Mai Van Bay et al. 22 4 Conclusions In this study, we propose an approach that uses a thermodynamic cycle in combination with the density functional function theory and the SMD solvent model to predict the stability constant of Cu2+ complexes with organic fluorescent ligands. The calculated results show that there is a good agreement between the theoretical stability constant and the experimental value. The predicted stability constants ( log β ) of the Cu2+ complex with Calcein blue and FluoZin-1 are 13.33 (exp. 14.27) and 6.59 (exp. 6.01), respectively. The results also indicate that the selection of reference ligands is a very important task to calculate the stability constants of the target ligands. The more similar the reference and chelating center of a ligand are, the more accurate it is to predict the stability constant of complexes. The method for predicting stability constant presented in this work could be applied to many other complexes. Funding statement This research is funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 104.06-2017.51 (D.T.Q.). References 1. Lo KK. Molecular design of bioorthogonal probes and imaging reagents derived from photofunctional transition metal complexes. Accounts of Chemical Research. 2020;53(1):32-44. 2. Thomason JW, Susetyo W, Carreira LA. Fluorescence studies of metal-humic complexes with the use of lanthanide ion probe spectroscopy. Applied Spectroscopy. 1996;50(3):401-408. 3. Pan X, Jiang J, Li J, Wu W, Zhang J. Theoretical design of near-infrared Al3+ fluorescent probes based on salicylaldehyde acylhydrazone schiff base derivatives. Inorganic Chemistry. 2019;58(19):12618- 12627. 4. Bistri O, Reinaud O. Supramolecular control of transition metal complexes in water by a hydrophobic cavity: a bio-inspired strategy. Organic & Biomolecular Chemistry. 2015;13(10):2849-2865. 5. Roy LE, Martin LR. Theoretical prediction of coordination environments and stability constants of lanthanum lactate complexes in solution. Dalton Transactions. 2016;45(39):15517-15522. 6. Vukovic S, Hay BP, Bryantsev VS. Predicting stability constants for uranyl complexes using density functional theory. Inorganic Chemistry. 2015;54(8):3995-4001. 7. Kim M, Sim E, Burke K. Ions in solution: Density corrected density functional theory (DC-DFT). The Journal of Chemical Physics. 2014;140(18):18A528. 8. Galván-García EA, Agacino-Valdés E, Franco-Pérez M, Gómez-Balderas R. [Cu(H2O) n ]2+ (n = 1–6) complexes in solution phase: a DFT hierarchical study. Theoretical Chemistry Accounts. 2017;136(3). 9. Klamt A. The COSMO and COSMO‐RS solvation models. WIREs Computational Molecular Science. 2017;8(1). 10. The IUPAC stability constants database. Chemistry international - Newsmagazine for IUPAC. 2006;28(5). 11. Shiri F, Salahinejad M, Momeni-Mooguei N, Sanchooli M. Predicting stability constants of transition metals; Y3+, La3+, and UO2 2+ with organic ligands using the 3D-QSPR methodology. Journal of Receptors and Signal Transduction. 2020;41(1):59- 66. 12. Ghasemi JB, Salahinejad M, Rofouei MK. Review of the quantitative structure–activity relationship modelling methods on estimation of formation constants of macrocyclic compounds with different guest molecules. Supramolecular Chemistry. 2011;23(9):614-629. 13. Chen H, Shi R, Ow H. Predicting stability constants for terbium(III) complexes with dipicolinic acid and 4-substituted dipicolinic acid analogues using density functional theory. ACS Omega. 2019;4(24):20665-20671. 14. Mohammadnejad S, Provis JL, van Deventer JS. Computational modelling of gold complexes using density functional theory. Computational and Theoretical Chemistry. 2015;1073:45-54. 15. Devarajan D, Lian P, Brooks SC, Parks JM, Smith JC. Quantum chemical approach for calculating Hue University Journal of Science: Natural Science Vol. 129, No. 1D, 15–23, 2020 pISSN 1859-1388 eISSN 2615-9678 DOI: 10.26459/hueuni-jns.v129i1D.5947 23 stability constants of mercury complexes. ACS Earth and Space Chemistry. 2018;2(11):1168-1178. 16. Lukeš I, Šmídová I, Vlček A, Podlaha J. Study of bis (iminodiacetato) cuprates(II) and tetrakis (iminodiacetato) cuprates(II). A Chemical Papers. 1984;38(3):331-339. 17. Das AK. Stabilities of ternary complexes of cobalt(II), nickel(II), copper(II) and zinc(II) involving aminopolycarboxylic acids and heteroaromaticN-bases as primary ligands and benzohydroxamic acid as a secondary ligand. Transition Metal Chemistry. 1990;15(5):399-402. 18. Casasnovas R, Ortega-Castro J, Donoso J, Frau J, Muñoz F. Theoretical calculations of stability constants and pKa values of metal complexes in solution: application to pyridoxamine–copper(II) complexes and their biological implications in AGE inhibition. Physical Chemistry Chemical Physics. 2013;15(38):16303. 19. Pandey R, Kumar A, Xu Q, Pandey DS. Zinc(II), copper(II) and cadmium(II) complexes as fluorescent chemosensors for cations. Dalton Transactions. 2020;49(3):542-568. 20. Pliego JR. Reply to comment on: ‘Thermodynamic cycles and the calculation of pKa’ [Chem. Phys. Lett. 367 (2003) 145]. Chemical Physics Letters. 2003;381(1-2):246-247. 21. Bryantsev VS, Diallo MS, Goddard III WA. Calculation of solvation free energies of charged solutes using mixed cluster/continuum models. The Journal of Physical Chemistry B. 2008;112(32):9709- 9719. 22. Rebollar-Zepeda AM, Campos-Hernández T, Ramírez-Silva MT, Rojas-Hernández A, Galano A. Searching for computational strategies to accurately predict pKas of large phenolic derivatives. Journal of Chemical Theory and Computation. 2011;7(8):2528-2538. 23. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, et al. Gaussian 16 Rev. A.03. Wallingford, CT2016. 24. Marenich AV, Cramer CJ, Truhlar DG. Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. The Journal of Physical Chemistry B. 2009;113(18):6378-6396. 25. Alexander MD. Chelate ring conformations and substitution rates of cobalt(III) complexes. Inorganic Chemistry. 1966;5(11):2084-2084.