Insights into the coordination chemistry of antineoplastic doxorubicin with 3d-transition metal ions Zn²⁺, Cu²⁺, and VO²⁺: a study using well-calibrated thermodynamic cycles and chemical interaction quantum chemistry models

Abstract

We present a computational strategy based on thermodynamic cycles to predict and describe the chemical equilibrium between the 3d-transition metal ions Zn²⁺, Cu²⁺, and VO²⁺ and the widely used antineoplastic drug doxorubicin. Our method involves benchmarking a theoretical protocol to compute gas-phase quantities using DLPNO Coupled-Cluster calculations as reference, followed by estimating solvation contributions to the reaction Gibbs free energies using both explicit partial (micro)solvation steps for charged solutes and neutral coordination complexes, as well as a continuum solvation procedure for all solutes involved in the complexation process. We rationalized the stability of these doxorubicin-metal complexes by inspecting quantities obtained from the topology of their electron densities, particularly the bond critical points and non-covalent interaction index. Our approach allowed us to identify representative species in solution phase, infer the most likely complexation process for each case, and identify key intramolecular interactions involved in the stability of these compounds. To the best of our knowledge, this is the first study reporting thermodynamic constants for the complexation of doxorubicin with transition metal ions. Unlike other methods, our procedure is computationally affordable for medium-sized systems and provides valuable insights even with limited experimental data. Furthermore, it can be extended to describe the complexation process between 3d-transition metal ions and other bioactive ligands.

Introduction

It is well-documented that complexing bioactive compounds with transition metal ions often results in a reduction of their undesired secondary effects, while keeping or even enhancing their beneficial therapeutic properties [1]. Consequently, determining the thermodynamic complexation constants in the solution phase of such species is an active research field in pharmaceutical science and bioinorganic chemistry. However, this is not necessarily a straightforward task, as experimentalists may encounter several difficulties, including the solubility of the drug and/or metal complex(es), the kinetic and quantitative profiles of the reaction under study, the formation of poly-nuclear species (some of which may be undesired)[2], among others [3]. It should be noted that the experimental techniques used for this purpose often provide limited information about the structural properties of the formed complexes, for which it is common to appeal to additional spectroscopic studies.

Nowadays, theoretical protocols based on the quantum chemistry methods, are considered powerful auxiliary tools for studying several properties of organic and inorganic metal complexes, for instance, to estimate their complexation constants in gas and solution phases, to reveal their representative (most stable) coordination geometries in a given medium [4,5,6,7] and to discern molecular properties including some of their chemical reactivity features [8], among others. Particularly, the determination of the Gibbs free energies of complexation in solution phase by theoretical strategies, faces several drawbacks that can be classified in two main branches:

• 1. The quality of the electronic wave function used to find the equilibrium geometries and compute molecular properties. In computational chemistry, it is crucial to calibrate a suitable theoretical protocol (a method to determine the electronic structure of a chemical species) to ensure meaningful predictions, characterizations, and rationalizations of chemical properties or processes of interest. Among other considerations, a computationally affordable method is sought, which provides comparable results to experimental findings. If the problem under study lacks experimental information, it is common to appeal to results obtained from more sophisticated theoretical methods, such as Coupled Cluster (CC, the gold standard in quantum chemistry) or Configurations Interaction (CI), which, unfortunately, are computationally demanding and in general intractable for molecules constituted by more than a hundred of electrons. We emphasize that modern quantum chemistry offers a new and computationally affordable alternative, the Domain-Based Local Pair Natural Orbital (DLPNO) method which constitutes a variant of the standard CC formulation (the theoretical basis of this approach can be found in [9]); although, it is not suitable for large size molecules (with several hundreds of electrons), it can serve as a calibration standard to find an adequate theoretical protocol for the study of systems constituted by a large amount of electrons, for instance, transition metal complexes of bioactive compounds.

• 2. The models used for computing thermochemical properties in gas and solution phase. In most cases, gas phase thermochemistry is done under the ideal gas approximation, which often (not always) yields reasonable results. On the contrary, a more elaborated treatment is required for solution phase contributions, and it is common to use combined strategies, comprising at least by one of the following approaches: models based on a polarizable continuum [10, 11], molecular dynamics simulations [4, 12, 13], the Quantum-Mechanics Molecular-Mechanics (QM-MM) [14, 15] and, of particular interest in this work, Thermodynamic Cycles (TC) [5, 6, 16]. Recent studies have shown that TCs are promising strategies for the computation of solution phase quantities, especially for the reaction Gibbs free energy in solution phase. TCs are less time-consuming than Molecular Dynamics or QM-MM approaches, and, more importantly they can be designed to facilitate a systematic cancellation of various sources of errors, such as solvation patterns of solutes, proper charge screening by solvent molecules, and even those arising from the ideal gas approximation. This cancellation often translates into a substantial increase in the accuracy of the predicted results.

As an additional benefit, a proper theoretical study about the complexation equilibrium constant in solution phase, provides a set of outcomes (equilibrium geometries, wave and partition functions, electron densities, etc.) which, after further treatment, can be used to infer new and valuable information about the electronic and physico-chemical properties of the formed complex(es), for instance (and of interest in this study), key inter and intra-molecular interactions involved in the stability of the species.

This work explores the advantages of using these new affordable theoretical approaches and strategies to accurately describe and predict the chemical equilibrium between transition metal ions and bioactive molecules. Specifically, we investigate the complexation reaction between the antineoplastic doxorubicin and three transition metal ions, Cu²⁺, Zn²⁺ and VO²⁺, forming coordination species with 1:1 and 1:2 (metal-ligand) stoichiometries. To the best of our knowledge, there is no experimental data available for the corresponding equilibrium constants in aqueous solution. First, we generated all possible coordination geometries for each metal-doxorubicin complex and filled any coordination vacancies with water molecules. Next, we analyzed the stability of each of these geometries by computing the Gibbs free energy of complexation in solution phase using theoretical protocols based on thermodynamic cycles (which are explained in the following section). From these results, we identified plausible representative coordinated species for each complex in aqueous solution. Finally, we identified metal-ligand, ligand-ligand and solvent-ligand interactions that contribute to the stability of the most significant coordination geometries, using the Bond Critical Points (BCP) and the Non-Covalent Interaction Index (NCII) derived from the topological analysis of the electron density. Through this work, we demonstrate that even if experimental information about the complexation process between a bioactive molecule and a transition metal ion is scarce (or absent), a well-calibrated theoretical protocol can provide valuable chemical and thermodynamic information about the formation of the involved complexes in solution phase.

Background

Doxorubicin and doxorubicin transition metal complexes

Doxorubicin (HDOX) or Adriamycin (shown in Fig. 1) is a potent antineoplastic drug belonging to the anthracyclines family, known for its effectiveness against various solid tumors (as breast and lymphomas) and hematological neoplasms (as leukemia) [17, 18]. Unfortunately, as an adverse effect, its clinical use increases the risk of heart failure by a factor of eight [19,20,21,22,23]. Reports indicate that such undesired behavior can been attenuated by encapsulating HDOX inside liposomes [24] or nanoparticles [25,26,27], by using cardioprotective agents such as dexrazoxane [28,29,30,31,32] and by binding the free drug to transition metal ions [33,34,35,36,37]. Studies have reported several M-HDOX complexes (M = Co(II), Ni(II), Fe(III), Mn(II), Cu(II), Zn(II) and Mg(II)) able to induce significant apoptosis, implying that these metal complexes may exhibit enhanced therapeutic properties compared to free HDOX [37].

Fig. 1

Chemical structure of doxorubicin (HDOX). α and β correspond to the orientations where metal ions can bind the molecule

In this work we consider pertinent to analyze the chemical stability of the M-HDOX complexes particularly constituted by the Cu²⁺, Zn²⁺, and VO²⁺ ions due to the following reasons. Firstly, these three transition metals (copper, zinc, and vanadium) are highly bioavailable in the human body [38]. Secondly, the complexation of HDOX with Cu(II) has been shown to reduce the incidence of heart failure in both murine and human cell lines, while maintaining the same antitumor properties of the free drug [33]. Thirdly, Zn(II)-based drugs have demonstrated a significant effect in killing malignant cells in cancer therapy, with reports suggesting that complexes based on this metal show promise for the development of new bioactive species with favorable therapeutic effects [39,40,41,42]. Finally, several V(IV) and V(V) based compounds have shown promising antitumoral effects in in vitro studies [43,44,45,46,47,48,49,50]. Interestingly, even though HDOX and other vanadium complexes are effective against a similar spectrum of cell lines [43, 45, 47, 48, 50, 51], there has been no previous exploration of the structural or therapeutic properties of the V(IV)-HDOX species.

Experimental observations indicate that in aqueous solutions with pH values close to neutrality, where HDOX is predominantly in its neutral form (pK_a values of 7.84 and 10.04) [52, 53], transition metal ions bind to HDOX to form complexes with stoichiometric ratios of 1:1 and 1:2, depending on the oxidation state of the metal center and the reactant concentrations [34, 37, 54,55,56,57,58,59]. The complexation process proceeds with the release of a proton from each HDOX molecule, acidifying the medium [53]. Thus, the representative reaction between a specific metal ion M^q+ (q is the net charge on the metal) and n HDOX ligands in aqueous solution at a pH value beyond 7.84 and up to 10.04, forming the $\left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n } \right]_{\left( {aq} \right)}^{q - n}$ species, can be written according to the following chemical equation,

$${\text{M}}_{\left( {aq} \right)}^{q + } + n{\text{HDOX}}_{\left( {aq} \right)} \to \left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n } \right]_{\left( {aq} \right)}^{q - n} + n{\text{H}}_{\left( {aq} \right)}^+$$

(1)

It has been established that the two adjacent oxygen atoms (O) belonging to the quinone and phenol moieties (rings B and C in Fig. 1) of HDOX are the reactive centers involved in the coordination process with transition metal ions. Note that a HDOX molecule may bind a metal ion using O atoms from the alpha (α) or from the beta (β) sides (see Fig. 1).

Thermodynamic models

For the sake of clarity, at this stage we consider pertinent to specify the notation and nomenclature used through the remainder of the document. In $\left[ {{\text{M}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{q + }$, M is the metallic center, m corresponds to the number of water molecules coordinated to M, and q indicates the charge. In $\left[ {{\text{M}}^{\alpha ,\beta } \left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{q - n}$, n is the number of DOX⁻ species in the complex and l indicates the number of coordinated water molecules, while α and β indicates the orientation from where the doxorubicin molecule is attached to M.

Solvation of solutes by the cluster solvation cycle

To study the chemical equilibrium in solution phase from a theoretical perspective, it is essential to properly describe the solvation process of all solutes involved in the reaction. To address this issue, TCs commonly use a hybrid approach that combines a continuum model (where solvent is mainly represented by its static dielectric constant $\varepsilon$) to account for bulky contributions to the solvation Gibbs free energy, with an explicit "micro-solvation" step where solvent molecules are properly arranged around each solute to consider specific solute–solvent (and even solvent–solvent) interactions at the molecular level [60]. Continuum models provide a reasonable description of the solvation Gibbs free energy of neutral solutes and they often display a good performance when used as part of a thermodynamic cycle. Thus, a micro-solvation step offers an additional benefit. If charged solutes are involved in the chemical process being studied, explicit solvent molecules can screen their electric charges, contributing to a better description of the solvation process through a continuum approach.

In this work, we employed the “Cluster Solvation Cycle” shown in Fig. 2. This model considers that a single cluster composed of x solvent molecules spontaneously solvates a particular solute. Strategies based on the cluster cycle are adequate to maximize errors cancellation and have been recently utilized for accurately describing the thermodynamics of formation of some Cu(II) and Zn(II) complexes in water and ethanol solutions, respectively [7, 60].

Fig. 2

Thermodynamic cluster solvation cycle used for the computation of the solvation Gibbs free energy of charged solutes. The water molecules are arranged in a cluster fashion

In the cluster cycle, the solvation Gibbs free energy of a solute ${\text{A}}^{{\text{q}} \pm }$ is computed according to the following scheme,

$$\Delta_{\_c} G_{solv}^* ({\text{A}}^{{\text{q}} \pm } ) = \Delta G_g^o + \Delta G_{solv }^* \left( {\left[ {{\text{A}}\left( {{\text{H}}_{2} {\text{O}}} \right)_x } \right]^{{\text{q}} \pm } } \right) - \Delta G_{solv }^* \left( {{\text{H}}_{2} {\text{O}}} \right)_x - \Delta G_{corr}^o$$

(2)

where,

$$\Delta G_g^o = G\left( {\left[ {{\text{A}}\left( {{\text{H}}_{2} {\text{O}}} \right)_x } \right]^{{\text{q}} \pm } } \right) - G\left( {{\text{A}}^{{\text{q}} \pm } } \right) - G\left( {{\text{H}}_{2} {\text{O}}} \right)_x$$

(3)

(the subscript "g" will be used over the text to indicate a gas phase thermodynamic quantity) is the gas phase free energy contribution coming from the explicit solvation process, while $\Delta G_{solv }^* \left( {\left[ {{\text{A}}\left( {{\text{H}}_{2} {\text{O}}} \right)_x } \right]^{{\text{q}} \pm } } \right)$ and $\Delta G_{solv}^* \left( {{\text{H}}_{2} {\text{O}}} \right)_x$ are the continuum contribution to the solvation Gibbs free energies of the (screened) solute and of the solvent cluster constituted by x water molecules, respectively. The term,

$$\Delta G_{corr}^o = \Delta G^{o \to *} + RT\ln \left( {\left[ {{\text{H}}_{2} {\text{O}}} \right]/x} \right)$$

(4)

in Eq. (2) introduces proper corrections for the solute and solvent thermodynamic states; here, the $\Delta G^{ \circ \to \ast }$ quantity is the free energy correction due to the standard state conditions of solute, going from the ideal gas standard state (1 atm) to the solution standard state (1 M) at 298.15 K and it is estimated to be 1.89 kcal·mol⁻¹ under the ideal gas approximation. The thermodynamic correction corresponding to the solvent thermodynamic state is computed from the $RT\ln \left( {\left[ {{\text{H}}_{2} {\text{O}}} \right]/x} \right)$ term, which takes one 1 mol of the $\left( {{\text{H}}_{2} {\text{O}}} \right)_x$ species $\left( {55.34{\text{ M}}/x} \right)$ to the 1 M standard state.

Gibbs free energy of reaction in solution phase.

The formation of complexes between any of the metal centers under consideration (Zn²⁺, Cu²⁺ and VO²⁺) with n HDOX ligands in aqueous solution is represented by the generalized thermodynamic cycle shown in Fig. 3, and is described by the following equation,

$$\Delta_r G_{s,n}^* = \Delta_r G_g^o + \Delta \Delta G_{solv}^*$$

(5)

where,

$$\Delta_r G_g^o = G_g^o \left( {\left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n } \right]^{{\text{q}} - {\text{n}}} } \right) + nG_g^o \left( {{\text{H}}^+ } \right) - G_g^o \left( {{\text{M}}^{q + } } \right) - nG_g^o \left( {{\text{HDOX}}} \right)$$

(6)

is the Gibbs free energy of reaction in gas phase; $G_g^o \left( {{\text{H}}^+ } \right)$ is the gas phase free energy of the proton, and its value has been recently recalculated using quantum statistical mechanics under various thermodynamic conditions. At 1 atm and 298.15 K, it is recommended to be −6.28 kcal·mol⁻¹ [61]. The term,

$$\begin{gathered} \Delta \Delta G_{solv}^* = \Delta_{\_c} G_{solv}^* \left( {\left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n } \right]^{2 - n} } \right) + n\Delta G_{solv}^* \left( {{\text{H}}^+ } \right) \\ - \Delta_{\_c} G_{solv}^* \left( {{\text{M}}^{q + } } \right) - n\Delta G_{solv}^* \left( {{\text{HDOX}}} \right) \\ \end{gathered}$$

(7)

corresponds to the solvation Gibbs free energy of overall the involved solutes. $\left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n } \right]^{2 - n}$ and ${\text{M}}^{q + }$ are (charged) complexes which require solvent molecules to fill up their coordination vacancies and thus, their solvation Gibbs free energy will be computed using the thermodynamic solvation cycle; ${\text{HDOX}}$ is in its neutral form and a simple continuum solvation model will be applied for this species. The solvation Gibbs free energy of the proton $\Delta G_{solv}^* \left( {{\text{H}}^+ } \right)$ was set equal to − 265.9 kcal·mol⁻¹, a numerical value obtained from high level first-principles calculations [62]. Thus, using Eq. (2) in Eqs. (7) and (3) in Eq. (6), the following expression is obtained for the Gibbs free energy of reaction in aqueous solution ($\Delta_r G_{s,n}^*$, Eq. (5)),

$$\Delta_r G_{s,n}^* = \Delta_{r\_c} G_g^o + \Delta \Delta_{\_\varepsilon } G_{solv}^* + RT\ln \left[ \frac{l}{m} \right],$$

(8)

(recall that the "m" and "l" integers are used to indicate the number of solvent molecules considered for the explicit solvation of the metal centers and the $\left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n } \right]^{2 - n}$ species, respectively). The term $\Delta_{r\_c} G_g^o$ in Eq. (8) contains all the gas phase free energy contributions to the reaction, including those related to the explicit solvation, i.e., due to the formation of the solvated species $\left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{q - n}$ and $\left[ {{\text{M}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{q + }$, thus,

$$\begin{gathered} \Delta_{r\_c} G_g^o = G_g^o \left( {\left[ {M\left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_2 {\text{O}}} \right)_l } \right]^{q - n} } \right) + G_g^o \left( {{\text{H}}_2 {\text{O}}} \right)_m - G_g^o \left( {\left[ {{\text{M}}\left( {{\text{H}}_2 {\text{O}}} \right)_m } \right]^{q + } } \right) \\ - nG_g^o \left( {{\text{HDOX}}} \right) - G_g^o \left( {{\text{H}}_2 {\text{O}}} \right)_l + \left( { - 6.28{\text{ kcal}} \cdot {\text{mol}}^{ - 1} } \right) \\ \end{gathered}$$

(9)

while the term $\Delta \Delta_{\_\varepsilon } G_{solv}^*$ is given by,

$$\begin{gathered} \Delta \Delta_{\_\varepsilon } G_{{\text{solv}}}^* = \Delta G_{{\text{solv}}}^* \left( {\left[ {M\left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_2 {\text{O}}} \right)_l } \right]^{q - n} } \right) + \Delta G_{{\text{solv}}}^* \left( {{\text{H}}_2 {\text{O}}} \right)_m - \Delta G_{{\text{solv}}}^* \left( {\left[ {M\left( {{\text{H}}_2 {\text{O}}} \right)_m } \right]^{q + } } \right) \\ - n\Delta G_{{\text{solv}}}^* \left( {{\text{HDOX}}} \right) - \Delta G_{solv}^* \left( {{\text{H}}_2 {\text{O}}} \right)_l + \left( { - 265.9{\text{ kcal}} \cdot {\text{mol}}^{ - 1} } \right), \\ \end{gathered}$$

(10)

and contains all the solvation free energy contributions to the Gibbs free energy of reaction in solution, under the continuum media approximation. In Appendix A we provide a detailed guide on how to utilize this approach to compute the solvation Gibbs free energy of complexation in solution phase.

Fig. 3

Generalized thermodynamic cycle for the computation of the reaction Gibbs free energy $\Delta_r G_{s,n}^*$ in aqueous solution

It should be mentioned that the use of continuum solvation models, such as the ones employed in this study, can potentially lead to an overestimation of the loss of the entropy of reaction, particularly in non-equimolar processes due to solvent cage effects. To address this potential error, the Okuno model,[63] based on the free volume approximation [64], is commonly employed for reactions occurring in non-polar solvents [65, 66]. For instance, this model provides a correction of approximately 2.6 1 kcal·mol⁻¹ for 2:1 stoichiometric reactions, under these circumstances. While this approach has been used for the computation of Gibbs free energy in reactions involving aqueous solutions and other polar solvents [67, 68], we have chosen not to incorporate it into our calculations. It is important to acknowledge that this decision means that our predicted exergonicities may be slightly lower than the actual values, and this effect may have a greater impact on the formation of the 1:2 doxorubicin-metal complexes compared to the 1:1 complexes. Nevertheless, it is still relevant to assess how our model, which involves a cancellation of various sources of error, may contribute to the overestimation of the loss of the entropy of reaction. Thus, we computed the enthalpy contribution to the Gibbs free energies of all the studied reactions and are reported in the Supplementary Material.

Computational details

According to the thermodynamic strategy under consideration, the Gibbs free energy of reaction in solution phase between free HDOX molecules and a transition metal ion must be obtained by summing up contributions from gas and continuum phase calculations (Eq. (8)). Gas-phase contributions, those required for the complexation reaction and for the explicit solutes screening (Eq. (9)) were obtained from a standard gas-phase thermochemistry theoretical protocol (under the ideal gas approximation); we used the PBE [69], ωB97X-D [70], M06 [71] and PBE0 [72, 73] density functional approximations (DFAs) conjointly with the DGDZVP [74, 75] basis set, always ensuring that each of the optimized geometries corresponded to a minimum on the potential energy surface (PES) by an analysis of vibrational frequencies. These DFAs were selected as they have shown a good performance in the description of diverse chemical properties of transition metal complexes [76,77,78,79]. The $\left( {{\text{H}}_{2} {\text{O}}} \right)_x$ clusters involved in the explicit solvation process of solutes (Eqs. (3) and (10)), were taken from reference [80] and re-optimized according to our theoretical protocol. Continuum contributions to the solvation Gibbs free energy of the $\left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{q - n}$, $\left[ {{\text{M}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{q + }$, $\left( {{\text{H}}_{2} {\text{O}}} \right)_x$ and HDOX species (Eq. (10)) were computed as the energy difference between two single point electronic structure calculations, one at the M05-2X/6–31+G(d,p) and the other at M05-2X/6–31+G(d,p)/SMD [11] levels, both calculations were done on the previous gas phase optimized geometries (with each of the DFAs under consideration). All our calculations were performed using the Gaussian 09 quantum chemistry computational package [81] (optimized structures of the most stables complexes found in this study are reported in the Supplementary Material).

As a calibration step, we computed the electronic energy variation of the following reaction (gas phase formation of 1:1 stoichiometry complex for each metal):

$$\left( {{\text{H}}_{2} {\text{O}}} \right)_l + \left[ {{\text{M}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{q + } + {\text{ HDOX}} \to { }\left[ {{\text{M}}\left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{q - 1} + \left( {{\text{H}}_{2} {\text{O}}} \right)_m + {\text{H}}^+$$

(11)

at the DLPNO-CCSD/aug-cc-pVDZ [82] level of theory. Then the most suitable DFA for the study of the chemical equilibrium of each of the metal-doxorubicin complexes under investigation was selected and used for further analysis. The multireference wave function T1 diagnostic was performed at the DLPNO-CCSD/aug-cc-PVDZ level of theory for the most stable 1:1 stoichiometric metal-doxorubicin complexes found in this study. The computed norm of the single substitution amplitudes vector t₁, (divided by the square root of the number of correlated electrons) is reported in Table S1 of the Supplementary Material. As can be corroborated, these values were below the accepted threshold equal to 0.05 for compounds constituted by transition metals belonging to the 3d series [83]. We used the ORCA computational package [84, 85] for our Coupled-Cluster calculations.

The inter and intramolecular interactions among the HDOX species, solvent molecules and/or the metallic center were analyzed through the Quantum Theory of Atoms in Molecules (QTAIM), proposed by Bader [86]. We used two of the most outstanding QTAIM-based quantities, the Bond Critical Points (BCPs) and the Non-Covalent Interactions Index (NCII) [87]. The corresponding calculations were done using the GPUAM computational package [88, 89], building a grid of 10 points per atomic unit and using as input the wave function obtained with the ωB97X-D functional and DGDZVP basis. We selected this functional since it has shown significant improvements over several empirical dispersion-corrected DFAs when it comes to describing non-covalent interactions, [70] and thus, as indicated in previous work [90, 91], is an adequate theoretical method for the study of weak interactions using the mentioned QTAIM based tools. The NCII index is based on the density (ρ) and its derivatives. It enables the identification of non-covalent interactions from the reduced density gradient (s). To calculate “s” only the electron density is required. We computed both the 3D and 2D representations of the NCII indexes. The former was obtained from the corresponding molecular density, using an isosurface value of 0.5. The 2D-NCII maps were computed with the NCIPLOT software, using the pro-molecular approximation [92,93,94]. The sign of the second density Hessian eigenvalue ($\lambda_2$) was analyzed by evaluating sign $\left( {\lambda_2 } \right)\rho$. This analysis allows us to identify attractive interactions as regions where $\lambda_2 < 0$ and repulsive interactions as regions where $\lambda_2 > 0$. Additionally, weak van der Waals (vdW) interactions can be characterized by $\lambda_2 \approx 0$. To visualize these regions, it is common to appeal to the following color code also based on sign $\left( {\lambda_2 } \right)\rho$: blue for strong attractive interactions ($\lambda_2 < 0$, as strong hydrogen bonding), green for weak vdW interactions ($\lambda_2 \approx 0$, π–π stacking or van der Waals contacts) and red for strong repulsive interactions ($\lambda_2 > 0$, nonbonding, steric clashes). BCPs (Bond Critical Points) were employed to determine the coordination number of the metal complexes under investigation, as well as to identify hydrogen bonding interactions. On the other hand, Non-Covalent Interaction Index (NCIs) were utilized to locate specific molecular regions where significant non-covalent interactions occur.

Results and discussion

Before delving into the analysis of each family of doxorubicin-metal complexes investigated, it is important to highlight that the computed entropy contribution to the solution phase Gibbs free energy of reaction displayed values close to zero for all the studied chemical processes (computed values were lower than 0.1 kcal·mol⁻¹·K⁻¹). This can be corroborated by comparing the results in Tables 2, 4, and 6 with those in Tables S2, S3, and S4 of the Supplementary Material, respectively. Nonetheless, it should be noted that the contribution of the total entropic term $\left( {T\Delta S} \right)$ to the Gibbs free energy of reaction in the solution phase ranges between −6 and −18 kcal·mol⁻¹, and it certainly contributes to the exergonicity of the studied reactions.

Cu(II)-doxorubicin complexes

Among the three families of metal-complexes investigated in this work, those comprised by Cu(II) ions have been extensively studied. Both, the 1:1 and 1:2 stoichiometries are known to be stable in aqueous solution, but there is some controversy about which one is predominant at near-neutral pH conditions. Dabrowiak and Greenaway observed the 1:2 Cu(II):DOX complex over a range of pH values from 7 to 8.2 [34]; conversely, Garnier-Suillerot et al. observed that the 1:1 species is the one predominant at a neutral pH condition, although none of these studies provided information about the corresponding thermodynamic equilibrium constants [95, 96]. In any case, the square planar coordination geometry was considered as the most prone to be observed in solution phase. It is known that, in general, Cu(II) complexes may display different coordination numbers, from four to six [97,98,99]. For instance, by UV–vis titration experiments, Salmon et al.[100] found that the hexahydrate $\left[ {{\text{Cu}}\left( {{\text{H}}_{2} {\text{O}}} \right)_6 } \right]^{2 + }$ complex is the most abundant in aqueous solution, while Tabouli et al.[101] reported the $\left[ {{\text{Cu}}\left( {{\text{NH}}_{3} } \right)_5 } \right]^{2 + }$ (square pyramid) and $\left[ {{\text{Cu}}\left( {{\text{NH}}_{3} } \right)_6 } \right]^{2 + }$ (octahedral) species, as the most stable coordination geometries for cooper-ammonia complexes, using a quantum chemistry theoretical protocol. We thus consider necessary to explore all the coordination numbers in our computations, for both the free Cu(II) (from tetrahydrate to hexahydrate complexes, m = 4, 5, 6) and the $\left[ {{\text{Cu}}\left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{q - n}$ species, to find the chemical process which provides the lowest Gibbs free energy of complexation. An equivalent exploration will be performed for the remainder of the metal complexes under investigation.

Table 1 presents the results of the calibration study (Eq. (1)) for copper-doxorubicin metal-complexes. In addition to the electronic energy variations, we also report the relative reaction energies (in parentheses) referenced to the lowest energy complexation process, to facilitate comparison across different density functional approximations (DFAs). Remarkably, all the functionals here tested predicted the same trend as our DLPNO-CC standard; octahedral Cu(II)-DOX complexes (l = 4) are the most stable in gas phase, followed by the pentahedral (l = 3) and the tetrahedral (l = 2). For the octahedral and tetrahedral complexes, the HDOX molecule is preferable attached to the metallic center from the beta side, while for the pentahedral, the alpha orientation is preferred. Formation of complexes was generally enhanced when cooper tetrahydrate (square planar, (m = 4)) was used as starting material. It is worth noting that in some cases, our octahedral geometries evolved towards the penta-coordinated species after the geometry optimization process, due to a water molecule shifting towards the first solvation sphere. This was the case of the $\left[ {{\text{Cu}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_4 } \right]^+$ and the $\left[ {{\text{Cu}}\left( {^\alpha {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_4 } \right]^+$ species, particularly when their geometries were optimized with the PBE and PBE0 DFAs, respectively. The PBE0 functional displayed the closest reaction energies to our DLPNO-CC standard, with an average deviation (AD) of 14.21 kcal·mol⁻¹ followed by ωB97X-D (AD = 21.23 kcal·mol⁻¹), M06 (AD = 24.41 kcal·mol⁻¹) and PBE (AD = 43.93 kcal·mol⁻¹), while the M06 DFA displayed the best (closest to our standard) relative energy differences (quantities in parentheses in Table 1) with an average deviation of 9.71 kcal·mol⁻¹, followed by ωB97X-D (12.92 kcal·mol⁻¹), PBE0 (18.83 kcal·mol⁻¹) and PBE (20.40 kcal·mol⁻¹). Therefore, results from the PBE0 and M06 functionals will be prioritized in the analysis of the complexation process of copper-doxorubicin species in solution phase. Nonetheless, it is important to acknowledge that even with functionals that exhibit low averaged deviations, there may still be cases where the predicted relative stability trend deviates from our CC standard. For instance, the M06 functional predicts a higher stability for the alpha pentahedral complex compared to its beta counterpart. Additionally, in our calculations using the PBE0 DFA, we encountered difficulties in obtaining a minimum energy structure for the $\left[ {{\text{Cu}}\left( {^{\alpha ,\beta } {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^+$ complex, which is the second most stable species predicted by our CC protocol. Thus, for the analysis of chemical processes in the solution phase, we always prioritized the use of at least two different density functional approximations (DFAs), which enriched the discussion and provided more robust insights. Figure 4a, b, c–f show the BCPs and the NCII plots for the $\left[ {{\text{Cu}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_4 } \right]^+$ (octahedral) and $\left[ {{\text{Cu}}\left( {^\alpha {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]^+$ (square planar) coordination geometries, respectively. The most favored species displays two OH-O interactions (green dots, Fig. 4a) against the only one observed for the square planar complex (Fig. 4c). The 3D-NCII surfaces calculated for both conformers (Fig. 4c, d), conjointly with the 2D-NCII maps (Fig. 4e, f) indicates the presence of more attractive strong non-covalent interactions (blue surfaces) for the $\left[ {{\text{Cu}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_4 } \right]^+$ species and are mainly located at the coordination sphere for both compounds. This suggest that the metal–ligands attractive interactions may be stronger for the octahedral than for the tetrahedral complex.

Table 1 Computed energy variation corresponding to the reaction between $\left[ {{\text{Cu}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{2 + }$ and HDOX ligand, forming $\left[ {{\text{Cu}}\left( {^{\alpha ,\beta } {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^+$ complexes, in gas phase, using Eq. (11) relative values respect the most stable complexation process are reported in parentheses

Fig. 4

a, b BCPs, c, d 3D-NCII plots and e, f 2D-NCII maps obtained for the $\left[ {{\text{Cu}}\left( {^\alpha {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]^+$ and $\left[ {{\text{Cu}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_4 } \right]^+$ species, respectively, at the ωB97X-D/DGDZVP level of theory. Only BCPs between H₂O molecules and DOX⁻ are plotted, while BCPs between the metallic center and DOX⁻ or H₂O molecules were replaced by turquoise solid lines. “s” is the reduced gradient of the electron density

In Table 2 we report the Gibbs free energy of reaction in solution phase for the Cu(II)-DOX complexes. One may observe that both 1:1 and 1:2 stoichiometric species are prone to (co)exist in aqueous solution at the specified pH conditions (all the tested reactions, analyzed through the set of DFAs under consideration, were exergonic). Furthermore, these results indicate that the representative species at both stoichiometries may have a square planar geometry (l = 2 for n = 1 and l = 0 for n = 2). The formation of the former (1:1 stoichiometry) is particularly favored when the ${\text{DOX}}^-$ species is attached to the metal from the beta side, $\left[ {{\text{Cu}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]^+$, and also, if the free Cu(II) ion is found in its hexahydrate (m = 6) coordination geometry prior to the complexation process. It is worthy to recall that precisely the Cu(II) hexahydrate is the most abundant in aqueous solution. For the 1:2 complexes, results obtained with the PBE and PBE0 functionals indicate that the beta orientation is the one preferred by the HDOX molecule, $\left[ {{\text{Cu}}\left( {^\beta {\text{DOX}}} \right)_2 } \right]$, while results from the M06 and the ωB97X-D DFAs rather suggest that is the alpha orientation the one attached to the metallic center, $\left[ {{\text{Cu}}\left( {^\alpha {\text{DOX}}} \right)_2 } \right]$. In any case, the Gibbs free energy difference between both orientations is less than 4 kcal·mol⁻¹ and conjointly, these results suggest that for the 1:2 stoichiometry situation, both are prone to coexist (the reaction energy difference between the process with m = 5 and m = 6, using the M06 functional, is less than 0.3 kcal·mol⁻¹). Results obtained with the PBE0 DFA, suggest that the 1:1 stoichiometry is slightly more stable than the 1:2 (by barely 4 kcal·mol⁻¹). However, it is important to note that experimental validation is necessary to confirm this prediction, for instance, it is crucial to consider the potential for a higher overestimation of the loss of reaction entropy in the formation of the 1:2 complex. Figure 5a shows the geometry and BCP of the $\left[ {{\text{Cu}}\left( {^\alpha {\text{DOX}}} \right)_2 } \right]$ species. The complex is stabilized by five intermolecular BCPs, with one located on the left and four on the right-hand side. The 3D-NCII analysis shown in Fig. 5b supports this observation and shows that, both weak (green surface) and/or strong (blue surfaces) attractive non-covalent interactions are abundant in these molecular moieties. Interestingly, despite the size of the doxorubicin ligands, the 2D-NCII plot (Fig. 5c) indicate that repulsive interactions are the less abundant for this complex.

Table 2 Solution phase Gibbs free energy of complexation between $\left[ {{\text{Cu}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{2 + }$ and n HDOX ligands, forming $\left[ {{\text{Cu}}\left( {^{\alpha ,\beta } {\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{2 - n}$ complexes

Fig. 5

a BCPs, b 3D-NCII and c 2D-NCII plots obtained for the $\left[ {{\text{Cu}}\left( {^\alpha {\text{DOX}}} \right)_2 } \right]$ species at the ωB97X-D/ DGDZVP level of theory. Only BCPs between H₂O molecules and DOX⁻ are plotted, while BCPs between the metallic center and DOX⁻ or H₂O molecules were replaced by turquoise solid lines. “s” is the reduced gradient of the electron density

Zn(II)-doxorubicin complexes

Zn(II) species are commonly found in their tetra, penta and/or hexa coordination modes in aqueous solution, for instance, the most abundant coordination geometries for the zinc-aquo complexes in solution phase are the tetrahedron and the octahedron, while the pentahydrate is practically absent [102,103,104,105,106,107,108,109,110,111,112,113]. Jabłońska-Trypuć et al. [37] reported that under near-neutral pH conditions, the Zn(II)-DOX complexes are predominantly found in a 1:2 stoichiometry, based on their partial molar ratio measurements. No further reports were found about the formation of these species. Therefore, like the analysis of Cu(II)-DOX complexes, we will explore all coordination patterns for both hydrated zinc and zinc-doxorubicin species, to propose plausible complexation processes in aqueous solution.

Our theoretical gas-phase benchmarking for zinc complexes (Table 3) reveals that the PBE0 functional provides the closest results to our DLPNO-CCSD standard (AD = 6.24 kcal·mol⁻¹), followed by ωB97X-D (AD = 13.09 kcal·mol⁻¹), M06 (AD = 14.65 kcal·mol⁻¹) and PBE (AD = 21.61 kcal·mol⁻¹); while, the M06 functional displays the best relative energy differences with an average deviation of 3.90 kcal·mol⁻¹, followed by ωB97X-D (AD = 8.64 kcal·mol⁻¹), PBE0 (AD = 12.59 kcal·mol⁻¹) and PBE (AD = 16.76 kcal·mol⁻¹). Thus, as previously considered for Cu(II)-DOX species, for zinc-doxorubicin metal complexes results obtained from PBE0 and M06 functionals will be prioritized in our analysis of results computed in solution phase. Octahedral Zn(II)-DOX geometries were predicted to be the most stable 1:1 stoichiometric complex under gas phase conditions, particularly if the ${\text{DOX}}^-$ species is attached to the metal from the beta side. Also note that the spontaneity of these reactions is enhanced when the tetracoordinated Zn(II)-aquo species is considered as reactant (m = 4). In Fig. 6a, b we show the BCPs of the $\left[ {{\text{Zn}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_4 } \right]^+$ and the $\left[ {{\text{Zn}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]^+$ complexes, respectively, for comparison purposes. The complex with the higher stability displays two OH–O bonds between the H₂O molecules that coordinates the metal center and the O atoms of the DOX⁻ species. Conversely, the less stable tetrahedral complex only exhibited one of these OH-O bonds. Likewise, both the 3D-NCII plots (Fig. 6c, d) and the 2D-NCII maps (Fig. 6e, f) reveal that the attractive non-covalent interactions (blue and green surfaces) are substantially less abundant in the less stable complex, with particular emphasis on the region nearby the coordination sphere.

Table 3 Computed energy variation corresponding to the reaction between $\left[ {{\text{Zn}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{2 + }$ and HDOX ligand, forming the $\left[ {{\text{Zn}}\left( {^{\alpha ,\beta } {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^+$ complex, in gas phase, using Eq. (11)

Fig. 6

a, b BCPs, c, d 3D-NCII plots and e, f 2D-NCII maps obtained for the $\left[ {{\text{Zn}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_4 } \right]^+$ and $\left[ {{\text{Zn}}\left( {^\beta {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]^+$ species, respectively, at the ωB97X-D/ DGDZVP level of theory. Only BCPs between H₂O molecules and DOX⁻ are plotted, while BCPs between the metallic center and DOX⁻ or H₂O molecules were replaced by turquoise solid lines. “s” is the reduced gradient of the electron density

The Gibbs free energy of complexation of the Zn(II)-DOX species in solution phase (Table 4) generally exhibit positive or slightly negative values for both stoichiometries, indicating that the formation of zinc complexes in aqueous solution is less favorable than that of copper complexes. The tested functionals consistently suggest that the 1:1 stoichiometric species may preferentially adopt a pentahedral coordination geometry (l = 3), and, according to results obtained with the ωB97XD-D, M06 and PBE0 functionals, the corresponding formation process is almost equally favored if the Zn(II)-aquo complex is in its penta- or hexahydrate forms (m = 5, 6; note that the latter is reported as the most stable in solution phase). All tested functionals indicate that this process is almost insensitive to whether the DOX⁻ molecule is attached to the metal center from the alpha or beta side, and thus, both binding patterns are prone to be observed. Based on results obtained from the M06 functional, the formation of the octahedral complex is more likely, especially when the zinc pentahydrate (m = 5) is considered as reactant; nonetheless, one must recall that this Zn(II)-aquo complex is practically absent in aqueous solution, and thus, we rather preferred to consider the pentahedral zinc-doxorubicin complex as the representative 1:1 stoichiometric species, which can be formed from the hexahydrate species (which is abundant in aqueous solution), and has a complexation energy difference less than 1.0 kcal·mol⁻¹ respect the octahedron.

Table 4 Solution phase Gibbs free energy of complexation between $\left[ {{\text{Zn}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{2 + }$ and n HDOX ligands, forming $\left[ {{\text{Zn}}\left( {^{\alpha ,\beta } {\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{2 - n}$ complexes

Our findings suggest that the formation of the 1:2 stoichiometric Zn(II)-DOX complexes are not too favored; the PBE0 and ωB97X-D functionals predict that none of the complexation processes analyzed for these species will occur (all are endergonic), while M06 and PBE indicate that the octahedral (l = 2) or the tetrahedral (l = 0) complexes may exist in solution, respectively, with the drug attached to the metal from the alpha side, in any case. Based on our benchmarking study (Table 3), where we concluded that results obtained with the M06 functionals are of high priority, we selected the $\left[ {{\text{Zn}}\left( {^\alpha {\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]$ complex as the representative Zn(II)-DOX 1:2 species. This complex can be formed from three different zinc hydrates (m = 4, m = 5, m = 6). However, it should be noted that the pentahydrate is less abundant in aqueous solution, and more importantly, that the formation from the tetrahydrate is the most exergonic process. Now that we have identified the representative species for both stoichiometries, we can analyze their relative stability using the results obtained from the M06 functional. Both, the $\left[ {{\text{Zn}}\left( {^\alpha {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_3 } \right]^+$ and $\left[ {{\text{Zn}}\left( {^\alpha {\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]$ complexes are expected to be formed from the zinc tetrahydrate, which is abundant in solution phase. The relative stability of these species is thus directly determined by the difference in their molecular Gibbs free energies. Our calculations show that such energy difference is about 3 kcal·mol⁻¹, favoring the 1:1 stoichiometric species and is too small to draw definitive conclusions about the predominance of one species over the other, for instance, recall that a higher overestimation of the loss of entropy is expected for the formation of the 1:2 complex. Figure 7a, b depict the geometry of the $\left[ {{\text{Zn}}\left( {^\alpha {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_3 } \right]^+$ and the $\left[ {{\text{Zn}}\left( {^\alpha {\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]$ species conjointly with the BCPs surrounding their corresponding coordination centers. The former display two BCP flanking the coordination sphere, while the second displays seven BCPs. These findings suggest that intramolecular interactions are more important for the stability of the 1:2 stoichiometry species, than in their 1:1 counterpart.

Fig. 7

BCPs for the a $\left[ {{\text{Zn}}\left( {^\alpha {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_3 } \right]^+$ and b $\left[ {{\text{Zn}}\left( {^\alpha {\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]$ species, obtained at the ωB97X-D/ DGDZVP level of theory. Only BCPs between H₂O molecules and DOX⁻ are plotted, while BCPs between the metallic center and DOX⁻ or H₂O molecules were replaced by turquoise solid lines

VO ²⁺ doxorubicin complexes

It is documented that the VO²⁺ ion is found rather in acidic aqueous solutions at pH values up to 6.0 [114]. Therefore, the formation of the oxovanadium-doxorubicin complexes must occur at pH values below 6.0, however, the protonated (cationic) form of doxorubicin predominates under such conditions. Consequently, both doxorubicin and VO²⁺ species will carry positive electric charges, and their complexation may be hindered by electrostatic repulsion forces unless doxorubicin gets first deprotonated to its neutral form (for instance, due to the interaction with the metallic center). The whole complexation process can be thus described by the following Hess law,

$$\begin{gathered} \left[ {{\text{VO}}\left( {{\text{H}}_{2} {\text{O}}} \right)_5 } \right]^{2 + } + n{\text{HDOX}} \to \left[ {{\text{VO}}\left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{2 - n} + n{\text{H}}^+ \\ n{\text{H}}_2 {\text{DOX}}^+ \to n{\text{HDOX}} + n{\text{H}}^+ \\ \_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\\ \left[ {{\text{VO}}\left( {{\text{H}}_{2} {\text{O}}} \right)_5 } \right]^{2 + } + n{\text{H}}_2 {\text{DOX}}^+ \to \left[ {{\text{VO}}\left( {{\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{2 - n} + 2n{\text{H}}^+ \\ \end{gathered}$$

(12)

The Gibbs free variation of such process is given by,

$$\begin{gathered} \Delta_r G_{aq,tot}^* = \Delta_r G_{aq,n}^* + nRT \cdot {\text{p}}K_a \\ = \Delta_r G_{aq,n}^* + n\left[ {4.65{\text{ kcal}} \cdot {\text{mol}}^{ - 1} } \right] \\ \end{gathered}$$

(13)

where $\Delta_r G_{aq,n}^*$ is the solution phase Gibbs free energy of reaction computed through our strategy based on TCs. In Eq. (13) we used the experimental pK_a value of 7.84 corresponding to the deprotonation of the amine moiety of ${\text{H}}_2 {\text{DOX}}^+$, while the value of the thermodynamic temperature T was set equal to 298.15 K. It should be noted that the addition of HDOX species to the VO(IV) center is more restricted than for Cu(II) and Zn(II), due to the presence of the oxo ligand. The coordination vacancies must be accurately classified to unambiguously identify the generated coordination geometries, which were established to be equatorial or axial, as depicted in Fig. 8. A HDOX molecule can be linked to the VO(IV) center through two equatorial vacancies or one equatorial and one axial. To indicate the former case, a superscript ^E will be included in our notation, whereas the superscript ^A will be used for the latter.

Fig. 8

Molecular structure of the $\left[ {{\text{VO}}\left( {{\text{H}}_{2} {\text{O}}} \right)_5 } \right]^{2 + }$ species displaying the axial (ax) and equatorial (eq) coordination sites

Table 5 presents the results of our benchmarking study for the VO(IV)-doxorubicin complexes. All the theoretical methods here tested predicted that the octahedral complex $\left[ {{\text{VO}}\left( {^{\beta^E } {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_3 } \right]^+$ is the most stable 1:1 stoichiometric species in gas phase. DOX⁻ was observed to prefer the beta orientation, occupying two of the equatorial coordination vacancies of the vanadium center. The PBE0 functional displayed the closest reaction energies to our DLPNO-CC standard, with an average deviation (AD) of 6.88 kcal·mol⁻¹ followed by ωB97X-D (AD = 14.81 kcal·mol⁻¹), M06 (AD = 15.62 kcal·mol⁻¹) and PBE (AD = 81.40 kcal·mol⁻¹). The ωB97X-D DFA in turn displayed the best (closest to our standard) relative energy differences (quantities in parentheses in Table 1) with an average deviation of 2.49 kcal·mol⁻¹, followed by PBE0 (4.53 kcal·mol⁻¹), M06 (5.28 kcal·mol⁻¹) and PBE (5.98 kcal·mol⁻¹). Figure 9a–c show the BCP, 3D-NCII and 2D-NCII plots for the $\left[ {{\text{VO}}\left( {^{\beta^E } {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_3 } \right]^+$ species, respectively displaying two attractive hydrogen bonding interactions (indicated by blue surfaces) that stabilize the complex. The most pronounced attractive non-covalent interactions are observed in the vicinity of the coordination sphere, primarily resulting from the interaction between solvent molecules and the metal center. This finding highlights the importance of adequately filling coordination vacancies to determine the most stable oxovanadium complex in solution phase. By ensuring the appropriate coordination environment, the complex can engage in favorable interactions with the surrounding solvent molecules, enhancing its stability and overall behavior.

Table 5 Computed energy variation corresponding to the reaction between $\left[ {{\text{VO}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{2 + }$ and HDOX species, forming $\left[ {{\text{VO}}\left( {^{\alpha ,\beta } {\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{2 - n}$ complexes, in gas phase, using Eq. (11)

Fig. 9

a BCPs, b 3D-NCII and c 2D-NCII plots obtained for the $\left[ {{\text{VO}}\left( {^{\beta^E } {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_3 } \right]^+$ species at the ωB97X-D/ DGDZVP level of theory. Only BCPs between H₂O molecules and DOX⁻ are plotted, while BCPs between the metallic center and DOX⁻ or H₂O molecules were replaced by turquoise solid lines. “s” is the reduced gradient of the electron density.

The Gibbs free energies of complexation of VO(IV) with doxorubicin in acidic aqueous solution were computed using Eq. (12) and are reported in Table 6. Except for the PBE functional, all the methods predicted that the $\left[ {{\text{VO}}\left( {^\alpha {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]$ complex, which has a square pyramid coordination geometry, is the only predominant 1:1 stoichiometric species in solution. As in the gas phase, the oxygen atoms of DOX occupy equatorial vacancies of the oxovanadium center, but in solution phase, the alpha orientation was preferred by the drug. For the 1:2 stoichiometry case, we observed that the formation of the square pyramidal $\left[ {{\text{VO}}\left( {^\beta {\text{DOX}}} \right)_2 } \right]$ species was highly exergonic only when the drug binds to the oxovanadium center through the beta side. Nonetheless, our results strongly suggest that the 1:1 complex predominates over the 1:2 complex, as indicated by a significant Gibbs free energy difference of 53.78 kcal·mol⁻¹ (computed using the PBE0 functional). This substantial energy difference is unlikely to be counteracted by errors arising from the overestimation of the loss of entropy. Figure 10a, b show the coordination geometries of the $\left[ {{\text{VO}}\left( {^\alpha {\text{DOX}}} \right)\left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]$ and $\left[ {{\text{VO}}\left( {^\beta {\text{DOX}}} \right)_2 } \right]$ species. Both structures show only one BCP, which may be indicating that the intramolecular interactions may not be too relevant for the stability of VO(IV)-doxorubicin complexes in solution phase.

Table 6 Solution phase Gibbs free energy of complexation between $\left[ {{\text{VO}}\left( {{\text{H}}_{2} {\text{O}}} \right)_m } \right]^{2 + }$ and n HDOX ligands, for $\left[ {{\text{VO}}\left( {^{\alpha ,\beta } {\text{DOX}}} \right)_n \left( {{\text{H}}_{2} {\text{O}}} \right)_l } \right]^{2 - n}$ complexes

Fig. 10

BCPs obtained at the ωB97X-D/aug-cc-pVDZ level. (a) $\left[ {{\text{VO}}\left( {^\beta {\text{DOX}}} \right)({\text{H}}_2 {\text{O}})_2 } \right]$ and (b) $\left[ {{\text{VO}}\left( {^\beta {\text{DOX}}} \right)_2 } \right]$. Only BCPs between H₂O molecules and DOX⁻ are plotted, while BCPs between the metallic center and DOX⁻ or H₂O molecules were replaced by turquoise solid lines

Conclusions

In this study, we used a hybrid thermodynamic cycle approach to analyze the complexation equilibrium between the antineoplastic doxorubicin and three transition metal ions, Cu²⁺, Zn²⁺, and VO²⁺, in aqueous solution. To improve accuracy, we performed a preliminary benchmarking study to minimize errors resulting from the approximated electronic structure of the involved species. Additionally, we used explicit solvation with solvent clusters to enhance errors cancellation between the solvation patterns of the participating solutes. Our analysis revealed the representative 1:1 and 1:2 stoichiometric species, as well as the most probable complexation process for each metal–complex family. For instance, we found that copper-doxorubicin complexes have a squared planar coordination geometry at both stoichiometries and are formed from Cu(II) in its hexahydrate form; doxorubicin is linked to the metal center from oxygens located at the beta side. The 1:1 complex was found to be slightly more stable than the 1:2 complex but the validity of this result strongly depends on the overestimation of the loss of entropy going from gas to solution phase. For zinc, we found that the most favored complexes were the penta- and tetrahedral species for the 1:1 and 1:2 stoichiometries, respectively. These complexes are formed from the tetra and hexahydrates of free zinc and chelated by oxygens from the alpha side of doxorubicin. We did not find evidence of one stoichiometry predominating over another. In the case of oxovanadium-doxorubicin complexes, both 1:1 and 1:2 stoichiometric species exhibited a squared pyramidal geometry and were formed from the free octahedral VO²⁺ species in acidic aqueous solution, with doxorubicin exhibiting the alpha and beta orientations, respectively. Our results suggest that the 1:1 stoichiometric oxovanadium complex is predominant, while the 1:2 species may be scarce in aqueous solution.

Interestingly, our thermodynamic approach predicted entropy contributions close to zero to the Gibbs free energy of reaction in all the studied processes. This consistent pattern across the different complexes would be attributed to the thermodynamic model employed in estimating solution phase quantities. Our computational strategy is carefully designed to minimize errors and facilitate the cancellation of potential inaccuracies. However, it is important to acknowledge that the cancellation of minority Gibbs free energy components can also impact the computed reaction entropies. The accurate determination of the reaction entropy in the solution phase remains a challenging task and can be one of the major sources of error in strategies that rely on continuum solvation models and thermodynamic cycles for calculating (contributions to) solvation Gibbs free energies of reactions in polar solvents. Unfortunately, the only way to quantify this error is by comparing the predicted values with experimental measurements.

The coordination numbers inferred in our thermodynamic study were corroborated by the analysis of corresponding BCPs, which also served as indicators of relevant inter- and intramolecular interactions involved in the stability of the studied complexes. NCI indexes allowed us to build plausible discussions regarding the stability of one complex over another, in terms of the non-covalent interactions exhibited by each compound. Together, they constituted valuable tools to rationalize the thermodynamic stability of the set of metal-doxorubicin complexes investigated here.

Among the four DFAs tested in our study, PBE0 consistently produced Gibbs free energies of reactions that were the closest to our DLPNO/aug-cc-pVDZ standard for all three families of metal complexes investigated, followed by ωB97X-D. Relative energies were reasonably well predicted with the M06 functional, while ωB97X-D was the second-best performing functional. Likewise, in solution phase, these three functionals consistently provided the same trends among the sets of metal complexes studied, which also can be attributed to the cancellation of errors propitiated by our hybrid thermodynamic approach. Based on these results, we conclude that PBE0, ωB97X-D and M06 are suitable for studying the complexation process between the three metal ions considered in this study and biomolecules containing the anthracycline moiety, in both, gas and solution phase, particularly using the hybrid thermodynamic approach described herein. Conversely, our findings suggest that the PBE functional may not provide reliable results in further studies and must be avoided.

Finally, our study demonstrates that with the aid of a proper thermodynamic theoretical protocol and adequate theoretical tools, valuable information about the complexation between metal ions and bioactive molecules can be inferred, even in the absence of experimental data. Our results also suggest that reasonable predictions can be made, which can subsequently be compared with experimental determinations to gain a deeper understanding of the coordination chemistry of drugs bound to transition metal ions. Overall, our study highlights the potential of computational methods based on thermodynamic cycles in the field of bioinorganic chemistry and provides a useful framework for future studies in this area.

Appendix A

In this section, we will provide a detailed guide on how to utilize the cluster solvation cycle (Fig. 2) for calculating the Gibbs free energy of reaction in solution phase ($\Delta_r G_{{\text{aq}},n}^*$) for the interaction between transition metal ions and ligands of interest. As an illustrative example, we will examine the complexation between Zn(II) and HDOX resulting in the formation of the $\left[ {{\text{Zn}}\left( {{\text{DOX}}} \right)_2 } \right]$ species. The chemical equilibrium for this process can be described as follows,

$${\text{Zn}}\left( {{\text{II}}} \right) + 2{\text{HDOX}} \to \left[ {{\text{Zn}}\left( {{\text{DOX}}} \right)_2 } \right] + 2{\text{H}}^+$$

(14)

Upon filling up coordination vacancies for free Zn(II) and ${\text{Zn}}\left( {{\text{DOX}}} \right)_2$, two compounds are formed: the hexahydrate species and the tetrahedral complex, respectively. The representative equilibrium is,

$$\left( {{\text{H}}_{2} {\text{O}}} \right)_2 + \left[ {{\text{Zn}}\left( {{\text{H}}_{2} {\text{O}}} \right)_6 } \right]^{2 + } + 2{\text{HDOX}} \to \left[ {{\text{Zn}}\left( {{\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right] + \left( {{\text{H}}_{2} {\text{O}}} \right)_6 + 2{\text{H}}^+$$

(15)

The Gibbs free energy of complexation for this process can be calculated using Eq. (8), with a temperature of 298.15 K. It should be noted that for this case, with ligand coordination number (l) of 2 and metal coordination number (m) of 6, the equation takes the following form,

$$\Delta_r G_{{\text{aq}},2}^* = \Delta_{r\_c} G_g^o + \Delta \Delta_{\_c} G_{solv}^* + \left( {0.592{\text{kcal}} \cdot {\text{mol}}^{ - 1} } \right)\ln \left[ \frac{2}{6} \right]$$

(16)

The term $\Delta_{r\_c} G_g^o$ contains all the Gibbs free energy contributions in the gas phase, including those related to the explicit solvation of solutes, and is computed from Eq. (9)

$$\begin{gathered} \Delta_{r\_c} G_g^o = G_g^o \left( {\left[ {{\text{Zn}}\left( {{\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]} \right) + G_g^o \left( {{\text{H}}_{2} {\text{O}}} \right)_6 - G_g^o \left( {\left[ {{\text{Zn}}\left( {{\text{H}}_{2} {\text{O}}} \right)_6 } \right]^{2 + } } \right) \\ - 2G_g^o \left( {{\text{HDOX}}} \right) - G_g^o \left( {{\text{H}}_{2} {\text{O}}} \right)_2 + 2 \cdot \left( { - 6.28{\text{ kcal}} \cdot {\text{mol}}^{ - {1}} } \right) \\ \end{gathered}$$

(17)

All the remainder quantities in Eq. (17) can be obtained from standard thermochemistry theoretical calculations. On the other hand, $\Delta \Delta_{\_c} G_{solv}^*$ contains all the bulky contributions to the solvation Gibbs free energy of solutes,

$$\begin{gathered} \Delta \Delta_{\_\varepsilon } G_{solv}^* = \Delta G_{solv}^* \left( {\left[ {{\text{Zn}}\left( {{\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]} \right) + \Delta G_{solv}^* \left( {{\text{H}}_{2} {\text{O}}} \right)_6 - \Delta G_{solv}^* \left( {\left[ {{\text{Zn}}\left( {{\text{H}}_{2} {\text{O}}} \right)_6 } \right]^{2 + } } \right) \\ - 2\Delta G_{solv}^* \left( {{\text{HDOX}}} \right) - \Delta G_{solv}^* \left( {{\text{H}}_{2} {\text{O}}} \right)_2 + 2 \cdot \left( { - 265.9{\text{ kcal}} \cdot {\text{mol}}^{ - 1} } \right) \\ \end{gathered}$$

(18)

Each of the solvation free energies at the right-hand side of Eq. (18) were obtained from a continuum approach, i.e., as the difference between the electronic energies computed in both, gas ($EE_g$) and continuum ($EE_c$) phases,

$$\Delta G_{solv}^* \left[ {{\text{Zn}}\left( {{\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right] = EE_c \left[ {{\text{Zn}}\left( {{\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right] - EE_g \left[ {{\text{Zn}}\left( {{\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2 } \right]$$

(19)

$EE_c$ can be obtained from any continuum solvation model; the one selected in this work was SMD. Finally, substituting each free energy contribution with the corresponding numerical values one gets $\Delta_r G_{{\text{aq}},2}^* = 0.17{\text{ kcal}} \cdot {\text{mol}}^{ - 1}$ for the formation of the ${\text{Zn}}\left( {{\text{DOX}}} \right)_2 \left( {{\text{H}}_{2} {\text{O}}} \right)_2$ complex in aqueous solution. This procedure can be easily implemented in a simple python (or other language) script or in a standard spreadsheet.