Predicting Stability Constants for Terbium(III) Complexes with Dipicolinic Acid and 4-Substituted Dipicolinic Acid Analogues using Density Functional Theory.

The relative stability constants of Tb(III) complexes exhibiting binding to a series of 4-substituted analogues of dipicolinic acid (2,6-pyridinedicarboxylic acid) (DPA) were calculated using density functional theory (DFT) with the standard thermodynamic cycle. DFT calculations showed that the strengths of the stability constants were modified by the substituents in the following (decreasing) order: -NH2 > -OH ∼ -CH2OH > -imidazole ∼ -Cl ∼ -Br ∼ -H > -F > -I, with the differences among them falling within one to two log units except for -NH2. Through population and structural analysis, we observed that the -NH2, -OH, -CH2OH, and halide substituents can donate electrons via resonance effect to the pyridine ring of DPA while inductively withdrawing electrons with different strengths, thus resulting in the different binding strengths of the 4-substituted DPAs to the Tb(III) ions. We believe that these observations possess utility not only in the ongoing development of luminescent probes for bioanalytical studies but also for more recent cross-industrial efforts to enhance reservoir surveillance capabilities using chemical tracers within the oil and gas sector.

Introduction

Since the demonstration of energy transfer from aromatic ligands to bound Tb(III) or Eu(III) results in dramatically enhanced luminescence from the bound lanthanide ions,[1−4] extensive studies have been carried out on such systems. Because of the spectrally narrow and unique emission obtained from lanthanide ions and the millisecond-range excited-state lifetimes, organic donor–lanthanide acceptor complexes have been used in a wide variety of applications including biosensors,[5] fluorescence immunoassays,[6−8] luminescent probes for time-resolved microscopy and bioimaging,[9−14] displays,[15−17] and telecommunications.[18,19] In recent years, the use of Tb(III)-sensitized luminescence in detecting ultralow levels of dipicolinic acid (DPA)[20−28] has attracted much attention due to its relevance in microbial diagnostics[29,30] and biodefense.[31−34] Dipicolinic acid is a major constituent of the bacterial endospore structure,[35] released from the spore upon germination; the exploitation of this naturally occurring system phenomenon lends itself toward utility as biological indicators for routine environmental monitoring, markers of sterilization efficiency,[36,37] and, in the most inauspicious application, as first responder assays for Bacillus anthracis spores that have been employed as delivery vehicles in anthrax attacks.[38−41] In addition, trace detection of DPA has been used in studying panspermia, defined by aeronautics and space-related life detection, and life in extreme environments.[42−45]

The specific motivation of this work stems from recent efforts to deploy DPA-based molecules as subsurface oil reservoir tracers.[46−48] Defined as injecting tracers contained within injection fluids and monitoring for their appearance in neighboring production wells, tracer tests provide a plethora of information about reservoir heterogeneity, well-to-well communication, and fluid dynamics.[49−52] This information can then be used for the optimization of improved oil recovery (waterflooding or chemical EOR).[53] Our laboratory has been actively synthesizing and characterizing DPA derivatives, which altogether serve as “barcodes” that enable large-scale oil reservoir tracer campaigns including a multitude of injection and production wells[47,48] (here, barcodes are referred to as uniquely identifiable molecules, which can be traced back to their specific original injection wells). To this end, it would be beneficial to develop computational methods for the prediction and prescreening of new ligand designs. Tb(III)-sensitized luminescence has been selected as the primary detection modality for DPA-based ligand tracers due to its high sensitivity. The first step of detection comprises the complexation of the ligands with the terbium (Tb3+) ions, which have conventionally been quantified by the stability constants (formation constants).[54] In spite of the developments in computational methods and the vast increase in computational power over the past decades, reliable calculations of aqueous stability constants have been difficult. Absolute errors of |log K1exp – log K1calc| > 101 (with K1exp and K1calc are the experimental and calculated stability constants) have been commonplace, especially for multivalent anion ligands (due to the difficult solvation calculations of anion ligands’ diffusive electronic structures).[55−57] Nevertheless, the relative errors comparing the stability constants of ligand–ion complexes with similar structures are more constrained. Using the cleverly designed thermodynamic cycles, the systematic errors within the computational protocols may be largely canceled.[58−68] In this work, we attempted to calculate the relative stability constants for the complexation of terbium ions with a series of 4-substituted DPA analogues to preliminarily predict their efficiency in the Tb(III)-sensitized luminescence assays. In addition, population and structural analyses were carried out to elucidate the detailed electronic structures and the relative binding strengths of the complexes. We believe that the findings presented in this work not only inform better designs of the oil field ligand tracers but also offer meaningful insight toward the ongoing development of luminescent bioprobe applications.

Results and Discussion

Absolute Stability Constants for Nonsubstituted DPA and Tb(III) Complexes

We first benchmark the different density functional theory (DFT) methods by calculating the absolute stability constants for the nonsubstituted DPA and Tb(III) complexes before transitioning to the (more robust) relative stability constant calculations for the 4-substituted DPA and Tb(III) complexes. Figure shows the representative equilibrium geometries of the terbium ion with eight coordinated water molecules, [Tb(H2O)8]3+, DPA2– ligand, and their complex, [Tb(DPA)(H2O)8]+. We choose these water coordination numbers in the complexes because they have been predicted to have the most stable structures.[58,69,70]Table summarizes the computed stability constants for the nonsubstituted DPA ligand and Tb(III) complexes (Figure ) using DFT with different density functionals, basis sets, and solvation models. In Table , we also included the available experiment value.[26] As shown, the computed stability constants deviated from the experiment value significantly. However, it has been identified that the prediction of “absolute” stability constants is still out of reach with current computationally efficient methods; instead, the “relative” values comparing the trends between the related structures using the same computation protocols are more trustworthy (next section)[55−57,59] (It is noted that more computationally intensive methods such as explicit solvent quantum calculations may improve the predictions for absolute stability constants further;[71] nevertheless, they are beyond the scope of this work). Finally, in Table , we observed that the selections of basis sets (6-31+G* or 6-311++G**) have the least effect on the calculated stability constants (with the calculated log K1 varied by ∼ ±1), while the selections of solvation models (COSMO or SMD) and density functionals (B3LYP or M06) have much larger effects (with the calculated log K1 varied by ∼ ±10 and ±30, respectively). Given the limitations of the current computationally efficient methods, it would be very difficult to extract more information from the absolute stability constant calculations in Table .

Figure 1

Representative equilibrium geometries of the terbium ion with eight coordinated water molecules, [Tb(H2O)8]3+, DPA2– ligand, and their complex, [Tb(DPA)(H2O)8]+.

Table 1

Calculated and Experimental Stability Constants (log K1) for the DPA Ligand and Tb(III) Complexes

	B3LYP		M06
	COSMO	SMD	COSMO	SMD	exp[26]
6-31+G*	75.87	68.36	45.77	37.96	7.01
6-311++G**	74.82	66.94	45.34	37.34

Relative Complexation Free-Energy Changes and Relative Stability Constants for 4-Substituted DPA and Tb(III) Complexes

In this work, we studied eight 4-substituted analogues of DPAs (4-X-2,6-pyridinedicarboxylic acid; XDPA2– in short) with X = −F, −Cl, −Br, −I, −OH, −NH2, −CH2OH, and −imidazole. The equilibrium geometries for selected XDPA2– and [Tb(XDPA)(H2O)8]+ complexes are shown in Figure S1. These substitutes have been considered as alternative bioprobes[2,4] as well as different barcodes in subsurface oil reservoir tracers. We define the relative free-energy changes in the gas phase and in aqueous solution (Figure A) aswith the brackets behind ΔG represent the free-energy changes of the specified complexations. There were multiple interesting observations in Figure A. First, the relative free-energy changes calculated in the gas phase were much higher than in solutions, with ΔΔG°g ∼ ±10 kcal/mol and ΔΔGaq ∼ ±2 kcal/mol (except for the −NH2 substitute). The reduced energy variations in solution compared to in the gas phase have been observed in prior studies (where the solvated ions included transition metals as well as rare earth elements).[55,57,72] Second, the trends of ΔΔG°g and ΔΔGaq generally agree with each other (e.g., positive values for both ΔΔG°g and ΔΔGaq for the halide and imidazole substitutes (using COSMO solvation method) and negative values for both ΔΔG°g and ΔΔGaq for the NH2 substitute); however, contradictory results also exist (e.g., for OH and CH2OH substitutes where we see positive ΔΔG°g but negative ΔΔGaq). As a result, caution must be taken when correlating solution stability constants with only the gas-phase calculations.[73−75] Finally, the different solvation models may give qualitatively different results (Cl, Br, and imidazole substitutes), which indicated that a consistent description of ion solvation still remains a challenge in this area.

Figure 2

(A) Relative free-energy changes in the gas phase, ΔΔG°g, and in aqueous solution, ΔΔGaq and (B) relative stability constants, Δlog K1, for the 4-substituted DPA ligands complexing to Tb(III).

Figure B shows the relative stability constants of the different 4-substituted DPAs complexing to Tb(III), which were defined bywith the brackets again representing the specified complexations. Our calculations suggested that the strengths of the stability constants were modified by the substituents in the following (decreasing) order: −NH2 > −OH ∼ −CH2OH > −imidazole ∼ −Cl ∼ −Br ∼ −H > −F > −I, with the differences among them falling within one to two log units except for −NH2. We note that the calculated order was consistent with available experimental data from the literature comparing quantum yields of 4-substituted DPA and Tb(III) complexes (with order: −NH2 > −OH > −NHAc > −Cl > −H ∼ −Br),[4] which partially validated our DFT calculations since the strong complexations of the ligands with the ions are the first step for strong sensitization and luminescence.

Population and Structural Analysis

To further understand the different binding strengths of the Tb(III) and XDPA2– complexes, we performed population and structural analysis for the [Tb(DPA)(H2O)8]+ and [Tb(XDPA)(H2O)8]+ complexes. Figure shows the identifying labels of the key atoms within the complexes.

Figure 3

Identifying labels of the key atoms within the [Tb(XDPA)(H2O)8]+ complexes.

It has been understood that the substitutes on rings can donate electrons via resonance effect to the rings as well as inductively withdrawing electrons.[76] The order of electron-withdrawing strengths on the X position (Figure ): halide > −OH > −NH2 substitutes on the pyridine ring of XDPA2– result in the following order of partial charges on the N1 atom (Figure ): halide > −OH > −NH2 substitutes, which was clearly observed from the partial charge analysis (Figures B and 5B). In addition, the −CH2OH and −OH substitutes were observed to have similar partial charges on N1 (Figures B and 5B). The order of electrostatic potential energy between Tb and N1 atoms (approximated by the product of their Mulliken charges): halide > −OH ∼ −CH2OH > −NH2 substitutes, thus, result in the order of stability constant strengths: −NH2 > −OH ∼ −CH2OH > halide substitutes (c.f. Figure ). On the other hand, the partial charges of Tb, O1, and O2 were seem unaffected with the 4-substitution (Figures A,C,D and 5A,C,D).

Figure 4

Partial charges predicted by the Mulliken population analysis in the gas phase for (A) Tb, (B) N1, (C) O1, and (D) O2 atoms in [Tb(XDPA)(H2O)8]+ complexes.

Figure 5

Partial charges predicted by the Mulliken population analysis in solution phase for (A) Tb, (B) N1, (C) O1, and (D) O2 atoms in [Tb(XDPA)(H2O)8]+ complexes.

The natural electron configuration of Tb in the complexes was analyzed by the Mulliken population analysis with 4f shell occupancy between 8.13 and 8.15 in [Tb(H2O)8]3+ complexes and between 8.15 and 8.18 in [Tb(DPA)(H2O)8]+ and [Tb(XDPA)(H2O)8]+ complexes for both the gas and solution phases. Indeed, the 0.13–0.18 electron excess on Tb inert 4f orbitals may be due to the primitive Mulliken analysis since electron donation to lanthanides occurs primarily to 5d and 6s orbitals. The natural 6s, 4f, and 5d electron configurations for Tb in different complexes are summarized in Table S4. It was observed that the 4-substitutions on DPA ligands have minimum effects on the relative electron occupations for Tb in the various complexes.

Figure shows the selected interatomic distances in the [Tb(XDPA)(H2O)8]+ complexes. It was observed that Tb–N1 distances for the strongest binding X = −NH2 complex were also the shortest (about 0.03–0.04 Å shorter then X = −H or halide substitutes; Figure A), and the next strongest binding X = −OH and −CH2OH complexes were the next shortest (about 0.01–0.02 Å shorter then X = −H or halide substitutes; Figure A). On the other hand, Tb–O1 and Tb–O2 distances were seem unaffected with the 4-substitution (Figure B,C).

Figure 6

Interatomic distances in [Tb(XDPA)(H2O)8]+ complexes for (A) Tb–N1, (B) Tb–O1, and (C) Tb–O2 pairs.

Conclusions

In the present study, we performed DFT calculations to study the absolute and relative stability constants for the Tb(III) complexes with DPA as well as a series of 4-substituted DPA analogues. Using different basis sets, density functionals, and solvation methods, we observed that the predictions of the absolute stability constants have remained challenging, with the selections of density functionals have the largest effect. We note that the same conclusion has been pointed out by prior studies,[55−57,59] and more research in this area is imminent. On the other hand, our calculations on the relative stability constants of the 4-substituted DPA ligands suggested that their binding strengths to Tb(III) were modified in the following order: −NH2 > −OH ∼ −CH2OH > −imidazole ∼ −Cl ∼ −Br ∼ −H > −F > −I, with most of the differences among them falling within one to two log units. In the context of our original motivation for selecting the ligands as oil field tracers, the results suggest that these substitutes should all be suitable for the application since the ligand–ion complexation efficiencies (as the first step for the Tb(III)-sensitization assays) are only moderately affected. Nevertheless, the predicted order of the modified stability constants remains to be experimentally verified. Last but not the least, we performed population and structural analyses on the [Tb(DPA)(H2O)8]+ and [Tb(XDPA)(H2O)8]+ complexes. We observed that the −NH2, −OH, −CH2OH, and halide substituents can donate electrons via resonance effect to the pyridine ring of DPA while inductively withdrawing electrons with different strengths resulting in different partial charges on N1 as well as different distances for Tb–N1, thus resulting in the different binding strengths.

Methods

Thermodynamic Cycle for Calculating Stability Constants

In this work, we focused on the binding of DPA ligands to Tb(III) for the formation of 1:1 complexes. For the 1:1 metal/ligand complexes with equilibrium M + L ⇌ ML, the experimental stability constant K1 is related to the Gibbs free-energy change for the reaction in solution, ΔGaq, bywhere aML, aM, and aL are the activities of the species involved in the equilibrium. First-principles methods have been developed to compute ΔGaq with the thermodynamic cycle shown in Figure .[57,59,60] In this process, ΔG°g represent the free-energy changes of the metal and ligand bindings in the gas phase, and ΔG*solv represent the free energies of solvation for the transfer of 1 mole of solutes from the gas phase to the aqueous phase. Note that because ΔG°g were usually calculated with the standard state of an ideal gas at 1 atm whereas ΔG*solv were calculated with the standard state of 1 M, a conversion of an ideal gas from 1 atm (24.46 mol/L) to 1 M (1 mol/L) was applied withFinally, with water, an additional correction, Gaq = Gaq* + RT ln([H2O]), is needed if the pure solvent H2O(l) is adopted as the reference state for the solvent. Here, RT ln([H2O]) = 2.38 kcal/mol is a free-energy change associated with moving a solvent from a standard-state solution-phase concentration of 1 M to a standard state of pure liquid, 55.34 M. Note that common logarithm (log) was used to express stability constants following experimental conventions, while natural logarithm (ln) was used in the rest of the thermodynamic equations.

Figure 7

Thermodynamic cycle used to calculate ΔGaq.

Computational Details

In this work, we studied [Tb(H2O)8]3+, DPA2–, XDPA2–, [Tb(DPA)(H2O)8]+, and [Tb(XDPA)(H2O)8]+ complexes with X = −F, −Cl, −Br, −I, −OH, −NH2, −CH2OH, and −imidazole. The proposed structures are shown in Figures and S1. Density functional theory (DFT) calculations reported in this work were performed using the NWChem 6.8.1 program package.[77] We used both B3LYP and M06 density functionals to examine the variations in the selection of different functionals.[78−80] The Stuttgart relativistic small core (RSC) segmented/ECP set was used for terbium.[81] The ECP on the terbium atom accounts for scalar relativistic effects by replacing 28 electrons with a relativistic pseudopotential. The 6-31+G* or 6-311++G** basis sets with diffuse functions were employed for the light atoms.[82−85] All geometry optimizations and frequency calculations were performed in the gas phase. Most of the geometry optimizations and frequency calculations were performed at the B3LYP/6-311++G** and M06/6-311++G** level except for nonsubstituted DPA complexes where 6-31+G* basis set was also considered. Frequency calculations were used to verify that geometries were minima on the potential energy surface and to compute thermal corrections to enthalpies and total entropies (at T = 298.15 K), which were then combined with the total DFT energies to calculate the gas-phase free energies for individual species and their differences, ΔG°g (Figure ).

Using the gas-phase geometries, we calculated single-point aqueous solvation free energies, ΔG*solv, with both the conductor-like screening model (COSMO)[86,87] and the SMD model.[88] The free-energy changes in the aqueous phase, ΔGaq, and the stability constants, log K1, can then be calculated by completing the thermodynamic cycle (Figure ). We note that similar approaches have been used to predict the stability constants for other metal/ligand complexes including transition metals and rare earth elements with good results.[56−58]

The Mulliken population analysis included in the NWChem software was performed to acquire the partial charges and the natural electron configuration of the studied complexes.