Computational insight into manganese(II) complexes comprising macrocyclic ligands for magnetic resonance imaging

Abstract

The increasing interest in advances in the field of magnetic resonance imaging (MRI) and the promise of extensive development of Mn(II)-based MRI contrast agents have motivated us to investigate thoroughly manganese(II) complexes comprising macrocyclic ligands. Herein, 23 mononuclear Mn(II) complexes were selected and classified into five families based on the structural motif of the parent macrocycle. DFT methods were applied together with two implicit solvation models, the conductor-like polarizable continuum model (CPCM) and the solvation model based on density (SMD). The stability constants of Mn(II) complexes (log K_MnL), A(¹⁷O) hyperfine coupling values of aqua ligand(s), the zero-field splitting (ZFS) of the sextet ground state, and the thermodynamics of water dissociation were addressed and compared to the parameters applied in the analysis of the experimental data. Moreover, ZFS parameters were calculated using the multireference CASSCF/NEVPT2 method.

---

Introduction

Magnetic resonance imaging (MRI) is one of the most important non-invasive imaging modalities used for accurate visualization of soft tissues or body functions. Unlike computed tomography (CT), positron emission tomography (PET), or single photon emission computed tomography (SPECT), MRI does not require the use of X-rays or any other ionizing radiation; instead, it utilizes a strong external magnetic field combined with the application of radio frequency pulses (RF pulses) to generate detailed images.^1,2 However, to enhance MR image resolution and improve diagnostic accuracy, the use of contrast agents (CAs) is often required. The CAs currently used in clinics are mainly T₁-CAs – namely coordination compounds containing Gd(III).^3,4 T₁-CAs affect primarily the T₁ relaxation time of the ¹H water protons, which consequently results in positive contrast on T₁-weighted MRI images.⁵ The efficacy of Gd(III) complexes as T₁-relaxation agents is attributed mainly to the magnetic dipole–dipole interactions between nuclear and electron spins, which are particularly efficient due to the presence of the seven unpaired electrons and a relatively long electronic relaxation time.^6,7 The use of Gd(III)-based CAs is considered to be safe; however, several cases of adverse effects have been reported,⁸ and some of the utilized Gd(III) compounds have been associated with nephrogenic systemic fibrosis (NSF).⁹ Moreover, it has been reported that in patients who received multiple doses of Gd-based CAs, long term-accumulation of Gd(III) occurred.¹⁰ Lately, concerns have also been raised about the ecological impact of Gd-based MRI agents, which have become detectable in natural waters around the world.¹¹ With respect to all these factors, there has been an active search for appropriate MRI CA replacements.

As a result, complexes with paramagnetic transition metal ions, particularly Mn(II) or Fe(III) (S = 5/2),^12–14 have been extensively investigated. The 3d⁵ high-spin electron configuration of Mn(II), along with its slow electronic relaxation and fast water exchange, makes Mn(II) the most promising alternative to Gd(III).^15–17 Additionally, the shorter Mn–H_water bond distance compared to Gd–H_water helps compensate for the effect of the lower electron spin number on the relaxivity of Mn(II)-based CAs.¹⁸

Manganese is also an essential biological element acting, e.g., as a cofactor in several enzymes; thus, living organisms can handle a slight excess of free Mn(II). However, in large concentrations, it may result in a neurodegenerative disorder called “manganism,” which manifests with symptoms similar to Parkinson's disease.¹⁹

Therefore, Mn(II) must be strongly chelated to secure safe in vivo use. In the case of Mn(II), penta- or hexadentate ligands, either cyclic or non-cyclic, have been examined for the complexation of Mn(II), providing coordination numbers (C.N.) of 6 or 7 with at least one free site for the coordination of water molecules.^15,20

The overall efficiency of CAs is described as relaxivity, r₁, defined as the paramagnetic enhancement of the longitudinal water proton relaxation (T₁) in the presence of 1 mM concentration of a paramagnetic CA (eqn (1)).

The relaxivity is determined by the properties of the paramagnetic complex and originates from inner-, second- and outer-sphere contributions (eqn (2)).⁵

r_i = r^IS_i+ r^SS_i+ r^OS_i (2)

The inner-sphere relaxivity contribution is governed by microscopic parameters such as the number of water molecules coordinated to the metal center (described as the hydration number, q), the water exchange rate (k_ex), the rotational motion of the chelate described by the rotational correlation time (τ_r) and the relaxation rates of the electron spin of the metal ion, T_1,2e. The values of these microscopic parameters governing relaxivity are obtained after simultaneous analysis of ¹H NMRD and variable-temperature ¹⁷O NMR data using a set of analytical equations (see the SI) based on the Solomon–Bloembergen–Morgan theory of paramagnetic relaxation.²¹ The equations important for the consequently calculated parameters are given below (eqn (3)–(5)).

Appropriate ligand design allows fine-tuning of the microscopic parameters, thus significantly affecting the resulting relaxivity.

Moreover, the suitability of paramagnetic complexes as efficient CAs depends on several other factors,²² such as thermodynamic stability and kinetic inertness. The thermodynamic stability must be sufficient enough to prevent the dissociation and release of the metal ions, so it can be considered the key factor regarding the long-term persistence of the complexes in a living organism. The importance of kinetic inertness lies in the initial biological exposure to a contrast agent and thus is related to the safe in vivo use of CAs.²³

Recently, several attempts have been made to find suitable computational tools for predicting thermodynamic stability, to reduce synthetic efforts and help with ligand design. In one of the first attempts to correlate the structure and stability of Gd(III) complexes, topological and quantum chemical descriptors were used to characterize the structures of a small data set of about 20 compounds.²⁴ Moreover, in 2014, a quantitative structure–property relationship method based on machine learning was developed to predict the stability constants of Gd(III) complexes (log K_GdL), relying on parametrized functions encoding the 2D structure of the ligand.²⁵

More recently, other empirical correlations have been proposed to predict and rationalize the thermodynamic stabilities of Gd(III) and Mn(II) complexes with respect to their structures.^26,27 These correlations implied that the thermodynamic stability constants and pGd/pMn (pM = −log([M]ⁿ⁺_{free metal})) values can be approximated using structural descriptors (ligand motifs) that are easily identifiable, e.g. the number and type of donor groups, the size of the macrocyclic platform, the type of chelate ring, etc. The contribution of each structural descriptor was determined through a least-squares fitting procedure applied to stability data available in the literature. The overall stability constants and pGd/pMn values were then obtained by summing the individual contributions of these structural descriptors. This methodology, developed for both Gd(III) and Mn(II) complexes, is remarkably accurate. However, there are still small deviations when compared to the experimental stability constants, arising from the use of different media during measurements and their ionic strength. In the case of Mn(II) complexes, there are several limitations to the applicability of this approach, such as the need for stability constant datasets for specific donor atoms or ligand archetypes; thus, it does not anticipate novel chelators like bispidine.^28,29 It is also possible to observe huge differences between theoretical vs. experimental stability constants of Mn(II) complexes (log K_MnL) for 15-pyN₃O₂ derivatives.^30,31

The prediction of other important parameters, in the context of MRI CA development, namely the water exchange rate (k_ex), has been investigated using DFT (density functional theory) and wave function analysis for selected Gd(III) complexes.³² The examination of Gd–O_water bond distances, electron density (ρ_BCP), and electron localization function (ELF) at the bond critical points then allowed us to find a correlation between the experimentally observed water exchange rates and the calculated ρ_BCP and ELF values.

Additionally, more recent studies on Gd(III) and Mn(II) complexes also suggested the possible impact of transient and static ZFS on electron relaxation.³³ The effect of the electron spin relaxation on relaxivity was investigated in Mn(II) complexes with triazacyclononane derivatives that lack a water molecule coordinated in the inner sphere (q = 0); hence, the observed relaxivity involved only outer-sphere contribution.³⁴ The study compared the ZFS parameter values obtained from ¹H NMRD profiles with those calculated using CASSCF wave functions, but despite the effort, the rationalization of electron relaxation still remains difficult.

The aim of this work was to design a new computational model based on the structures of Mn(II) complexes with selected macrocyclic ligands (Scheme 1). Generally, macrocyclic ligands tend to be more suitable than acyclic chelators in terms of more sufficient thermodynamic stability, higher kinetic inertness and reasonable relaxivity values.^22,35 Therefore, for the purposes of this work, five different ligand groups containing various structural cyclic motifs were employed (Groups I–V), and their reported experimental values were compared to the calculated (theoretical) ones. This computational model aims to predict complex stability constants (log K_MnL) using several approaches (methods A–E), taking into account the geometry of the selected compounds as well as estimating other important parameters such as the water exchange rate (k_ex) and the metal–water distance (d(Mn–O_H2O)). Such predictions of complex stability can serve as a powerful tool to aid ligand design and to rationalize various factors influencing complex stability in solution. Furthermore, theoretical insights into several microscopic parameters can provide preliminary information regarding potential CA efficiency.

Scheme 1 Structural formulas of the selected macrocyclic ligands of the respective Mn(II) complexes employed in the computational study, distinguished by colors: Group I (red), Group II (blue), Group III (green), Group IV (purple) and Group V (orange).

Results and discussion

The primary challenge in theoretical calculations, regarding metal complexes, is the selection of an appropriate theoretical (computational) method. In this study, the calculations are based on cost-effective density functional theory (DFT).³⁶ Consequently, the second important issue is the selection of the DFT functional and the basis set. As we wanted to study Mn(II) complexes in water solutions, we tested two hybrid functionals, B3LYP and PBE0, and two (hybrid) meta-GGA functionals, r²SCAN and TPSSh. Concerning the basis set, we tested the Ahlrichs basis set def2-TZVP and its relativistic alternative ZORA-def2-TZVP.

First, the geometry of the {[Mn(H₂O)₆]·12H₂O} complex, including the second coordination sphere in analogy with previous work,³⁷ was optimized, and the average Mn–O distance was compared to the value of 2.165 Å derived from EXAFS measurements of Mn(II) water solutions.³⁸ The best agreement was achieved for the r²SCAN/ZORA-def2-TZVP combination (Table S1); therefore, the reported calculations were carried out in the same manner.

Next, we focused on the DFT calculation of the stability constants of Mn(II) complexes (log K_MnL) – see Table S2a. Generally, several different approaches can be adopted for this purpose.^39,40 First, the molecular geometries of reactants and products were optimized in vacuo. The molecular vibration calculations confirmed that the final geometries corresponded to the local minima due to the absence of imaginary frequencies. Then, solvation energies were calculated in four different ways using implicit solvation models, CPCM (Conductor-like Polarizable Continuum Model) and SMD (Universal Solvation Model based on Solute Electron Density) – see Scheme 2. In the case of methods A and B, the solvation energies were calculated with CPCM and SMD, respectively, on fixed geometries resulting from the vacuum geometry optimization. Considering two other methods C and D, the solvation energies were calculated with CPCM and SMD, respectively, on relaxed geometries. Finally, for method E, the geometry optimization, frequency calculations, and thus Gibbs energy calculations were carried out directly using the CPCM implicit solvation model for water – Scheme 2.

Scheme 2 Computational workflow for log K_MnL values according to methods A–E.

The molecular structures of the Mn(II) complexes obtained by method E are depicted in Fig. S1 and the XYZ coordinates are provided in the SI. The comparison of experimental and calculated log K_MnL values is presented in Fig. 1 for methods A, C and E related to CPCM, and in Fig. S2 for methods B and D related to SMD. Both figures also include linear fits to log K_MnL. Overall, these methods provide similar outputs/trends. It is evident that, across all methods, Group I displays the greatest deviation from the linear relationship.

Fig. 1 The comparison of the calculated and experimental log K_MnL values of the complexes depicted in Scheme 1 using methods A, C and E.

Interestingly, it was found that geometry optimization in vacuo (methods A–D) resulted in significantly longer Mn–O_H2O distance(s) compared to optimization with CPCM (Fig. 2) or SMD (Fig. S3) solvation models. This trend was observed for the majority of the complexes studied, with the exception of complexes 7, 10 and 22 in the case of CPCM, and complexes 10 and 22 when using SMD. In contrast, Mn–O_H2O distances obtained with CPCM and SMD solvation models are very similar (Fig. 2), indicating that both methods produce very similar, hence comparable, molecular geometries. To support this, we calculated the RMSD parameter⁴¹ between the geometries, of the selected Mn(II) complexes optimized in vacuo, with CPCM and SMD (excluding hydrogen atoms), which confirmed significant changes in the molecular geometries when optimized in vacuo and with solvation models – Table S3a.

Fig. 2 The comparison of the distances between the Mn^II ion and the oxygen atoms of coordinated aqua ligands as deduced from the DFT-optimized molecular structures in vacuo and with CPCM (top), and CPCM and SMD (bottom) – methods A–E.

Furthermore, methods A–E were modified by adding two extra molecules of water per aqua ligand for Mn(II) complexes 1–23, as it was already pointed out by other authors that it can improve the description of their properties.^37,42 Herein, these methods are labelled as A2–E2, and the results of the log K_MnL calculations are shown in Fig. S4 and Table S2b. Again, we observe that the match between the experimental and the calculated values is not perfect, and the overall trends are very similar to those obtained with methods A–E. The molecular structures of the complexes obtained by method E2 are depicted in Fig. S1 and the XYZ coordinates are provided in the SI. The comparison of RMSD parameters and Mn–O_H2O distances is available in Table S3b. Notably, methods A2–E2 provided more consistent and comparable Mn–O_H2O distance(s) upon geometry optimization, both in vacuo and when using CPCM or SMD models – Fig. S5, relative to the unmodified methods A–D described above.

The results discussed above suggest that applying the solvation model during geometry optimization is preferable. Therefore, we will further focus on methods E and E2. It is evident that calculations of log K_MnL values do not provide perfect correlation with experiments. This can be expected due to the large structural diversity of macrocyclic ligands, and consequently, the resulting Mn(II) complexes, as well as the imperfection of the solvation models.

Nevertheless, considering the future application of computational chemistry in this area of research, it would be useful to provide, at the very least, reasonable predictions of thermodynamic stability – specifically, log K_MnL values within groups/families sharing a similar molecular scaffold. This could assist synthetic chemists in rationalizing their synthetic efforts toward the design and development of more efficient (Mn-based) MRI CAs. Therefore, the results for Group II and Group IV, each having the largest number of members – six, are discussed in detail – Fig. 3.

Fig. 3 The detailed comparison of the calculated and experimental log K_MnL values of the complexes from Groups II and IV using method E (top) and method E2 (bottom).

In the case of Group IV, based on the 15-pyN₃O₂ macrocycle, the results are promising, as methods E and E2 correctly predict an increase of log K_MnL compared to the experimental values. Within Group II, there is also a linear trend except for complex 6. Complex 6 is a bit exceptional among all studied complexes, because it has three carboxylate groups and no coordinated water molecule. These groups could possibly interact with solvent water molecules through strong hydrogen bonds, a feature that is not captured by the approaches used herein and is not likely perfectly reproduced with solvation models.

Among other important parameters used in the analysis of potential Mn(II)-based MRI agents is also the hyperfine splitting (coupling constant) of ¹⁷O in aqua ligands. In Fig. 4, a comparison of the results from methods E and E2 with respect to the applied values for the analysis of the experimental data is presented.

Fig. 4 The comparison of DFT-calculated hyperfine splitting values of ¹⁷O of aqua ligand(s) using method E (top) and method E2 (bottom) vs. the experimentally applied values.

It seems that the experimentally applied values of |A(¹⁷O)| are usually fixed to the range of ca. 5–7 MHz, whereas the theoretically calculated values span a larger interval of 5–12 MHz. Concerning the comparison of methods E and E2, it was found that method E2 can provide both smaller or larger values of |A(¹⁷O)| than method E – Fig. S6. Moreover, it was found that the calculated |A(¹⁷O)| data reasonably correlate with the donor–acceptor Mn–O_H2O distance for method E; more specifically, for longer Mn–O_H2O distances, the smaller |A(¹⁷O)| values were calculated – Fig. 5. In contrast, method E2 does not provide such correlation – Fig. S7.

Fig. 5 The comparison of DFT-calculated hyperfine splitting values of ¹⁷O of aqua ligand(s) using method E vs. respective Mn–O_H2O distances.

Furthermore, another parameter of importance is Δ², which is the square of the trace of the zero-field splitting (ZFS) tensor. The ZFS describes the splitting of the ground spin state S = 5/2 into three Kramers doublets caused by the ligand field and the spin–orbit coupling according to spin Hamiltonian (eqn (6)):

Ĥ = D(Ŝ_z² − Ŝ/3) + E(Ŝ_x² − Ŝ_y²) + μ_BBgŜ (6)

where D and E parameters represent axial and rhombic zero-field splitting parameters, respectively. Hence, Δ² is then calculated as (2/3D² + 2E²). Herein, DFT was also applied to calculate these ZFS parameters, but it is well known that more sophisticated methods, like the post-Hartree method CASSCF/NEVPT2, are more suitable for metal complexes. The comparison is depicted in Fig. S8, showing the experimentally applied values of Δ² vs. DFT-calculated values (methods E and E2). Evidently, DFT values are too low and do not cover the whole range of experimental values. In contrast, CASSCF/NEVPT2 provides much larger values of the Δ² parameter; however, the comparison to the experimental data is rather scattered (Fig. 6).

Fig. 6 The comparisons of Δ² values related to the ZFS of Mn^II complexes applied in the analysis of the experimental data and calculated using CASSCF/NEVPT2 or DFT (method E).

The rather unclear or scattered correlation between theoretical and experimental data may arise from the use of “fixed” parameters such as A_O/ħ, or, in special cases, Δ².⁴³ However, Δ² strongly influences the shape of the low-field region of the ¹H NMRD profile and is typically extracted from experimental data; therefore, it is usually not treated as a fixed parameter. A more implicit explanation for the observed deviations, considering experimental vs. theoretical values of Δ², is the intrinsic limitations of the theoretical treatment of ZFS. In particular, the computed values can reflect, with rather reasonable accuracy, the contribution of the static ZFS, but the transient ZFS contribution may differ significantly, since rapid and short-lived distortions of the coordination polyhedron may lead to a broad distribution of ZFS values over time.^33a,44,45 A more rigorous evaluation of these transient effects requires molecular dynamics simulations that capture the ultrafast fluctuations of the coordination environment through e.g. ab initio MD studies.⁴⁶

As already mentioned in the Introduction above, the correlation between the experimentally observed water exchange rates (k_ex) and the calculated electron density (ρ_BCP) or electron localization function (ELF) of Gd–O_aq bond distances was reported.³² Thus, we also employed Quantum Theory of Atoms in Molecules (QT-AIM), and analyzed the properties of Mn–O_aq bonds. More specifically, the bond critical points (BCPs) were identified, and ρ_BCP values were evaluated for methods E and E2 and are displayed in Fig. S9. In contrast to Gd(III) complexes, the studied Mn(II) complexes did not provide any clear correlations.

These differences may originate from the fact that the vast majority of Gd(III) complexes undergo dissociatively activated water exchange reactions (with only a few exceptions).^21,47 Therefore, the correlation discussed above for Gd(III) compounds is rather not unexpected. In contrast, for Mn(II) complexes, the situation is different, since both associatively and dissociatively activated mechanisms of reaction have been reported.^15,21,48 Moreover, in the case of an associatively activated reaction mechanism, the departure of the coordinated water molecule—and thus the strength of the Mn–O_H2O bond—is no longer the rate-limiting step for the exchange process.

Finally, we also investigated the thermodynamic parameters of the dissociation of the aqua ligand from the respective Mn(II) complexes according to eqn (7):

[MnL(H₂O)_n]^2−x(aq) → [MnL(H₂O)_n−1]^2−x(aq) + H₂O(l) (7)

following the protocol for method E. Herein, we attempted to correlate the Gibbs free energy of this chemical reaction (ΔG_diss) with the (largest) Mn–O_H2O bond distance. There appears to be a correlation between the Mn–O_H2O bond distance and the dissociation free energy (Fig. 7). Generally, longer Mn–O_H2O distances correspond to lower (more negative) ΔG_diss, suggesting easier dissociation. However, there are some outliers, namely complexes 1–3 (Group I) and 22–23 (Group V). A possible explanation could be that these macrocyclic ligands are rather small and relatively rigid, thus able to easily accommodate the loss of water without significantly modifying the shape of the complexes.

Fig. 7 The comparison of the Gibbs free energy of the water dissociation according to eqn (7) with the d(Mn–O_H2O) distance.

Having information about ΔG_diss at our disposal, we tried to evaluate whether there is any correlation with experimental k_ex values – Fig. 8.

Fig. 8 The comparison of the Gibbs free energy of the water dissociation according to eqn (7) with experimental k_ex values. The coordination numbers (C.N.) of [MnL(H₂O)_n]^2−x are depicted using squares and circles for C.N. equal to 6 and 7, respectively.

Generally, various mechanisms of water exchange are proposed for these Mn(II) complexes, varying from associative interchange (I_a) through interchange (I) up to dissociative (D). If the Mn(II) complex has a C.N. of seven, it can be safely presumed that the D mechanism should be operative. In such a case, the k_ex should correlate with ΔG_diss. The dataset is significantly scattered (Fig. 8 – circles), but with the exception of 18 and 19, it holds that more exergonic dissociation reactions have larger k_ex experimental values.

Experimental

The ORCA 5.0.4 software was used for DFT calculations.⁴⁹ The starting molecular structures were either extracted from the Cambridge Crystallographic Database (CSD)⁵⁰ or built using program Avogadro 1.2.⁵¹ The following DFT functionals were used: B3LYP,⁵² PBE0,⁵³ r²SCAN⁵⁴ and TPSSh.⁵⁵ The atomic basis set used was either def-TZVP or its relativistic version ZORA-def2-TZVP.⁵⁶ The calculations were sped up using the RIJCOSX approximation.⁵⁷ The largest integration grid (DefGrid3) and tightSCF convergence criteria were used in all calculations. The implicit solvation models CPCM^58,59 or SMD⁶⁰ were used during geometry optimization based on methods A–E or A2–E2. The calculated Gibbs free energies were corrected for the change of 1 mol of an ideal gas from 1 atm to 1 M using the term RT ln(24.46) = 1.89 kcal mol⁻¹ at T = 298.15 K.⁶¹ Within the DFT framework, the zero-field splitting parameters were calculated using a coupled-perturbed method,⁶² including both spin–spin and spin–orbit contributions (DTensor keyword set to “ssandso”), and the former part was dealt with UNO orbitals.⁶³ To obtain ZFS parameters using the post-Hartree–Fock method CASSCF⁶⁴ using CAS(5e,5o) with NEVPT2 correction,^65,66 ORCA 6.0 was used⁶⁷ by taking into account 1 sextet, 24 quartets, and 75 doublets, because the newly implemented AVAS procedure was found very useful for the preparation of the active space.^68,69 The QT-AIM analysis was carried out with AIMAll software.⁷⁰

Conclusions

Thorough DFT calculations were performed for 23 Mn(II) macrocyclic-based complexes. The use of CPCM solvation during geometry optimization (method E) proved to be the most reliable, as it provided a reasonable description of the aqua and ligand(s) coordination. Our main motivation was to provide a reasonable estimation of stability constants log K_MnL for the predictive purpose of chemical design of novel MRI contrast agents. While absolute agreement with experimental data varied, especially for complexes with multiple carboxylate groups, the method captured the main trends within at least some ligand families.

Additionally, the parameters applied during the evaluation of the relaxivity were also investigated, namely, hyperfine splitting |A(¹⁷O)| and zero-field splitting Δ². It was shown that these parameters usually span a larger interval than is assumed during the analysis of experimental data. Attempts to correlate the water exchange rate (k_ex) with the electron density at the bond critical point of the Mn–O_H2O bond were unsuccessful. However, it seems that there is some correlation with the dissociative Gibbs free energy for heptacoordinate Mn(II) complexes. These outcomes should motivate further research in this area, as a comprehensive computational framework for predicting MRI-relevant properties, namely, stability constants (log K_MnL, thermodynamics) and water exchange rates (k_ex, kinetics), would be valuable and of great interest.

Author contributions

Radovan Herchel: conceptualization, formal analysis, funding acquisition, investigation, methodology, visualization, writing – original draft, and writing – review & editing. Marie Pražáková: investigation and writing – original draft. Bohuslav Drahoš: conceptualization and writing – original draft.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the SI. See DOI: https://doi.org/10.1039/d5dt01381j.

Acknowledgements

This research was funded by Palacký University in Olomouc through projects IGA_PrF_2023_007, IGA_PrF_2024_009, and IGA_PrF_2025_012. Computational resources were provided by the e-INFRA CZ project (ID:90254), supported by the Ministry of Education, Youth and Sports of the Czech Republic.

We would like to express our gratitude to all the reviewers for their insightful comments and suggestions, which have contributed to enhancing the quality of this theoretical study.

References

1. G. Yan and R. Zhuo, Chin. Sci. Bull., 2001, 46, 1233.
2. J. Wahsner, E. M. Gale, A. Rodriguez-Rodriguez and P. Caravan, Chem. Rev., 2019, 119, 957.
3. S. Aime, M. Botta and E. Terreno, Gd(III)-based Contrast Agents for MRI. in Advances in Inorganic chemistry, Elsevier, 2005, 57, pp 173–237.
4. G. Tircso, E. Molnar, T. Czupacz, Z. Garda, R. Botar, F. K. Kalman, Z. Kovacs, R. Brucher and I. Toth, Gadolinium(III)-Based Contrast Agents for Magnetic Resonance Imaging. A Re-Appraisal. in Metal-Ions in Bio-Imaging Techniques, ed. A. Siegel, E. Freisinger and R. K. O. Sigel, De Gruyter, 2021, pp 39–70.
5. E. Terreno, D. D. Castelli, A. Viale and S. Aime, Chem. Rev., 2010, 110, 3019.
6. L. Helm, Prog. Nucl. Magn. Reson. Spectrosc., 2006, 49, 45.
7. C. Luchinat, Magn. Reson. Chem., 1993, 31, S145.
8. M. Rogosnitzky and S. Branch, BioMetals, 2016, 29, 365.
9. T. Grobner, Nephrol., Dial., Transplant., 2006, 21(4), 1104.
10. M. Le Fur and P. Caravan, Metallomics, 2019, 11, 240.
11. Z. Zhang, W. Jiang, T. Gu, N. Guo, R. Sun, Y. Zeng, Y. Han and K. Yu, Environ. Pollut., 2024, 359, 124740.
12. N. Kuznik and M. Wyskocka, Eur. J. Inorg. Chem., 2016, 445.
13. E. A. Kras, E. M. Snyder, G. E. Sokolow and J. R. Morrow, Acc. Chem. Res., 2022, 55, 1435.
14. S. Chen, L. An and S. Yang, Molecules, 2022, 27, 4573.
15. B. Drahos, I. Lukes and E. Toth, Eur. J. Inorg. Chem., 2012, 1975.
16. S. Daksh, A. Kaul, S. Deep and A. Datta, J. Inorg. Biochem., 2022, 237, 112018.
17. A. Gupta, P. Caravan, W. S. Price, C. Platas-Iglesias and E. M. Gale, Inorg. Chem., 2020, 59, 6648.
18. J. S. Troughton, M. T. Greenfield, J. M. Greenwood, S. Dumas, A. J. Wiethoff, J. Wang, M. Spiller, T. J. McMurry and P. Caravan, Inorg. Chem., 2004, 43(20), 6313.
19. S. Rivera-Mancia, C. Rios and S. Montes, BioMetals, 2011, 24, 811–825.
20. C. Henoumont, M. Devreux and S. Laurent, Molecules, 2023, 28, 7275.
21. A. Merbach and E. Toth, The chemistry of contrast agents in Medical Magnetic Resonance Imaging, Wiley, 2013.
22. P. Caravan, J. Ellison, T. McMurry and R. Lauffer, Chem. Rev., 1999, 99, 2293.
23. G. A. Rolla, C. Platas-Iglesias, M. Botta, L. Tei and L. Helm, Inorg. Chem., 2013, 52, 3268.
24. Y. H. Qi, Q. Y. Zhang and L. Xu, Chem. Inf. Comput. Sci., 2002, 42, 1471.
25. F. Dioury, A. Dupart, G. Dreyfus, C. Ferroud and J. Cossy, J. Chem. Inf. Model., 2014, 54, 2718.
26. R. Uzal-Varela, A. Rodriguez-Rodriguez, H. Wang, D. Esteban-Gomez, I. Brandariz, E. M. Gale, P. Caravan and C. Platas-Iglesias, Coord. Chem. Rev., 2022, 467, No.214606.
27. R. Uzal-Varela, F. Perez-Fernandez, L. Valencia, A. Rodriguez-Rodriguez, C. Platas-Iglesias, P. Caravan and D. Esteban-Gomez, Inorg. Chem., 2022, 61, 14173.
28. (a) D. Ndiaye, M. Sy, A. Pallier, S. Meme, I. Da Silva, S. Lacerda, A. M. Nonat, L. J. Charbonniere and E. Toth, Angew. Chem., Int. Ed., 2020, 59, 11958; (b) D. Ndiaye, P. Cieslik, H. Wadepohl, A. Pallier, S. Meme, P. Comba and E. Toth, J. Am. Chem. Soc., 2022, 144, 22212; (c) M. Sy, D. Ndiaye, I. Da Silva, S. Lacerda, L. J. Charbonniere, E. Toth and A. M. Nonat, Inorg. Chem., 2022, 61, 13421.
29. P. Cieslik, P. Comba, B. Dittmar, D. Ndiaye, E. Toth, G. Velmurugan and H. Wadepohl, Angew. Chem., Int. Ed., 2022, 61, e202115580.
30. M. Prazakova, D. Ndiaye, E. Toth and B. Drahos, Dalton Trans., 2023, 52, 7936.
31. M. Prazakova, D. Ndiaye, E. Toth and B. Drahos, Inorg. Chem., 2025, 64, 8205–8221.
32. M. Regueiro-Figueroa and C. Platas-Iglesias, J. Phys. Chem. A, 2015, 119, 6436.
33. (a) S. Rast, A. Borel, L. Helm, E. Belorizky, P. H. Fries and A. E. Merbach, J. Am. Chem. Soc., 2001, 123(11), 2637–2644; (b) C. Platas-Iglesias, D. Esteban-Gomez, L. Helm and M. Regueiro-Figueroa, J. Phys. Chem. A, 2016, 120(32), 6467–6476.
34. R. Uzal-Varela, L. Valencia, D. Lalli, M. Maneiro, D. Esteban-Gomez, C. Platas-Iglesias, M. Botta and A. Rodriguez-Rodriguez, Inorg. Chem., 2021, 60, 15055–15068.
35. R. E. Mewis and S. Archibald, Coord. Chem. Rev., 2010, 254, 1686–1712.
36. F. Neese, Coord. Chem. Rev., 2009, 253, 526–563.
37. D. Esteban-Gomez, C. Cassino, M. Botta and C. Platas-Iglesias, RSC Adv., 2014, 4, 7094–7103.
38. H. Konieczna, D. Lundberg and I. Persson, Polyhedron, 2021, 195, 114961.
39. F. Jensen, Introduction to Computational Chemistry, John Wiley & Sons, Incorporated, Newark, UNITED STATES, 3rd edn, 2017.
40. C. J. Cramer, Essentials of Computational Chemistry: Theories and Models, John Wiley & Sons, Incorporated, Chichester, Great Britain, 2nd edn, 2004.
41. Version rmsd 1.3.2, see https://github.com/charnley/rmsd for citation information.
42. V. Patinec, G. A. Rolla, M. Botta, R. Tripier, D. Esteban-Gómez and C. Platas-Iglesias, Inorg. Chem., 2013, 52, 11173–11184.
43. A. de Sa, C. S. Bonnet, C. F. G. C. Geraldes, E. Toth, P. M. T. Ferreira and J. P. Andre, Dalton Trans., 2013, 42, 4522–4532.
44. E. Belorizky and P. H. Fries, Phys. Chem. Chem. Phys., 2004, 6, 2341–2351.
45. J. C. Miller and R. R. Sharp, J. Phys. Chem. A, 2000, 104, 4889–4895.
46. A. Lasoroski, R. Vuilleumier and R. Pollet, J. Chem. Phys., 2014, 141, 014201.
47. L. Helm and A. E. Merbach, Coord. Chem. Rev., 1999, 187, 151–181.
48. M. Botta, F. Carniato, D. Esteban-Gómez, C. Platas-Iglesias and L. Tei, Future Med. Chem., 2019, 11, 1461–1483.
49. F. Neese, Comput. Mol. Sci., 2022, 12, e1606.
50. I. J. Bruno, J. C. Cole, P. R. Edgington, M. Kessler, C. F. Macrae, P. McCabe, J. Pearson and R. Taylor, Acta Crystallogr., Sect. B: Struct. Sci., 2002, 58, 389–397.
51. M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek and G. R. Hutchison, J. Cheminf., 2012, 4, 17.
52. (a) C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter, 1988, 37, 785; (b) A. D. Becke, J. Chem. Phys., 1993, 98, 1372; (c) A. D. Becke, J. Chem. Phys., 1993, 98, 5648; (d) P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623.
53. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
54. J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew and J. Sun, J. Phys. Chem. Lett., 2020, 11, 8208–8215.
55. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
56. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297.
57. F. Neese, F. Wennmohs, A. Hansen and U. Becker, Chem. Phys., 2009, 356, 98.
58. M. Garcia-Rates and F. Neese, J. Comput. Chem., 2020, 41, 922.
59. V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995.
60. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378.
61. V. S. Bryantsev, M. S. Diallo and W. A. Goddard III, J. Phys. Chem. B, 2008, 112, 9709–9719.
62. F. Neese, J. Chem. Phys., 2007, 127, 164112.
63. C. Riplinger, J. P. Y. Kao, G. M. Rosen, V. Kathirvelu, G. R. Eaton, S. S. Eaton, A. Kutateladze and F. Neese, J. Am. Chem. Soc., 2009, 131, 10092–10106.
64. P. A. Malmqvist and B. O. Roos, Chem. Phys. Lett., 1989, 155(2), 189–194.
65. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger and J. P. Malrieu, J. Chem. Phys., 2001, 114(23), 10252–10264.
66. C. Angeli, R. Cimiraglia and J. P. Malrieu, Chem. Phys. Lett., 2001, 350(3–4), 297–305.
67. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2025, 15, e70019.
68. E. R. Sayfutyarova and S. Hammes-Schiffer, J. Chem. Theory Comput., 2019, 15, 1679–1689.
69. E. R. Sayfutyarova, Q. Sun, G. K.-L. Chan and G. Knizia, J. Chem. Theory Comput., 2017, 13, 4063–4078.
70. T. A. Keith, AIMAll (Version 19.10.12), TK Gristmill Software, Overland Park KS, USA, 2019, (https://www.aim.tkgristmill.com).

---