Molecular Dynamics Simulation of the Interaction between Common Metal Ions and Humic Acids

Abstract

Humic acids (HAs) have important environmental and geochemical effects on soil, water environments and sediment. HAs strongly complex some metal ions, which affects the migration of metal ions and the colloidal aggregation of HA. Here, the complexation of Ca²⁺ and Mg²⁺ with HA in aqueous solution under neutral conditions has been systematically studied by molecular dynamics (MD) simulation. The results show that the aggregation of HA is caused by the complexation of HA and metal ions, mainly due to the intermolecular bridging between Ca²⁺, Mg²⁺ and COO⁻ groups. Monodentate and bidentate coordinations have been found between Ca²⁺ and COO⁻ groups of different HA molecules in the same simulation system. Mg²⁺ only has a monodentate coordination with COO⁻ group.

1. Introduction

Humic acids (HAs) are a common chemical component in soil, groundwater and sediment. The acids are ubiquitous in the environment and are mainly formed by the weathering and decay of plants [1,2,3]. Complexes of HAs with inorganic and organic compounds have different chemical and biological stability and play an important role in the purification and advanced treatment of petrochemical wastewater. HA is considered to be the main pollution source of ultrafiltration (UF), nanofiltration (NF) and reverse osmosis (RO) membranes, which directly leads to the formation of the surface regulation layer (“biofouling”) to which microorganisms are attached [4,5,6]. HA can also combine with metal ions and form metal complexes in natural water and petrochemical wastewater, thus affecting the solubility, fluidity, toxicity and bioavailability of metals [7]. Some research has reported the effects of pH [8], metal species [9] and temperature [10] on the aggregation of HA. However, the effect of metal ion types on the aggregation behaviors of HA has not yet been clearly discussed. The strong complexation between HA and metal ions plays a fundamental role in the treatment of metal pollution and the elimination of HA in solution [11,12].

Most studies on aggregation have only focused on macroscopic experimental phenomena due to different experimental variables, such as the different aggregate sizes and electrokinetic properties of HA colloids. The relationship between different critical coagulation concentration (CCC) values and the adhesion efficiency of HA aggregates due to the presence of typical metal ions has rarely been studied [13]. Moreover, different types of metal ions may have different effects on the aggregation of HA, even if they have the same valence state. Recent analysis revealed that the aggregation behaviors and binding interactions of colloids rely on metal cation valence as well as cation types. For example, Stirpe et al. [14] found that Cu²⁺ demonstrates greater influence than Zn²⁺ in modifying the kinetics of protein aggregation. Kalinichev et al. [15] studied the molecular dynamics (MD) simulation of the interaction between natural organic compounds (NOM) and metal ions (Na⁺, Cs⁺, Mg²⁺, Ca²⁺). The results showed that monovalent ions form very weak outer spherical complexes with NOM, and Mg²⁺ has a weak interaction with NOM because of its strong hydration shell; the binding of Ca²⁺ with NOM was the strongest of all studied ions. The above results were obtained only in terms of the radial distribution functions (RDFs) and coordination number, However, the interaction energy was not calculated. Ai et al. [16] studied the aggregation mechanism of HA with different heavy metal ions (Ag⁺, Cd²⁺, Cr³⁺) and common metal ions (Na⁺, Ca²⁺, Al³⁺) by MD simulation. However, the aggregation of HA with the same valence state has been less widely discussed. In addition, the internal mechanism of the aggregation process between HA and metal aggregates, such as an intramolecular or intermolecular binding interaction, has not been studied thoroughly. A comprehensive study of the binding behavior of metal ions with HA, especially different aggregations of HA, will increase our understanding of their impacts in petrochemical wastewater treatment. In the present study, Ca²⁺ and Mg²⁺ with the same valence state in concentrated petrochemical wastewater were chosen to evaluate the interaction of HA under neutral conditions [17]. To characterize the aggregation process of HA metal complexes, the interaction energy, RDFs, solvent accessible surface areas (SASAs) and hydrogen bonds of different metal ion complexes were analyzed by using MD simulation. These analyses enhance our understanding the structure, kinetics and aggregation characteristics of HA, which is crucial to evaluate the application of HA in the removal process of metal ions from petrochemical wastewater and to reduce membrane pollution.

2. Models and Computational Methods

The Temple-Northeastern-Birmingham (TNB) [18,19] model is usually used as a typical representative of HA in experiments and is successfully used in simulation because TNB is similar to HA and its MD simulation is simple [15,16,20]. The molecular weight of the TNB model is 753 Daltons, which is in the range of the most common HA molecules [21]. As shown in Figure 1a, the TNB model contains three carboxyl groups, three carbonyl groups, two phenolic groups, two amino groups and four other R-OH alcohol groups [16,21,22]. Figure 1a shows the size, composition and molecular structure of the TNB model (C₃₆H₄₀O₁₆N₂), and Figure 1b shows the morphology of the TNB model in simulation software (Materials Studio 2019, Accelrys, San Diego, USA). As shown in

Table 1. Experimental data of several humic acid (HA) molecular models.

	C (wt.%)	H (wt.%)	O (wt.%)	N (wt.%)	S (wt.%)	P (wt.%)	Reference
TNB	57.4	4.9	34.0	3.7	-	-	[18]
Steelink	59.5	5.2	33.5	1.8	-	-	[24]
CSHA#9	53.3	4.2	37.9	1.4	3.0	-	[25]
Schulten	51.5	4.0	41.8	2.0	0.6	-	[26]
Exp. (SRHA)	52.5	4.2	42.7	1.1	0.6	0.02	[27]

, the experimental characterization of TNB model is extremely consistent with that of other HA models [23].

In the simulation, functional groups of HA were shown to have different protonation states under different pH conditions. The experimental results show that the three carboxyl groups of TNB (pKa = 4~5) are completely deprotonated under neutral conditions, while the two phenolic groups are completely protonated under neutral conditions (pKa ≈ 9) [27]. Because the amino group of the TNB model is located on the benzene ring, which is similar to aniline, according to the pKa value of aniline (pKa = 4.6), the amino group of TNB model can always remain electrically neutral (NH₂) [16]. These deprotonated carboxyl groups are the main source of negative charge in the simulated system at neutral pH and are also important binding sites of metal ions in the aggregation state because of the charge attraction to deprotonated sites. In this study, the force field used was COMPASS [26], which is commonly used in organic compounds. The force field can be regarded as the empirical expression of the potential energy surface, which is the basis of molecular mechanics, MD and the Monte Carlo simulation. In MD simulation, the force field method is used to calculate the energy of the system according to the position of atoms. Compared with the previous quantum mechanical method, it can save a large amount of time and can be used to perform calculations for systems containing tens of thousands of particles. The COMPASS force field can accurately simulate various organic and inorganic substances in a gaseous state, liquid state or solid state [28,29]. Firstly, energy minimization was performed for TNB and water molecules, respectively [26]. The optimized O-H bond length of water molecules was 0.957 Å, and the H-O-H angle was 104.526°. A simple point charge (SPC) was used by building a water model to describe the water molecules in the simulation system, and then four TNB molecules were randomly inserted into a water box (density: 997 kg/m³ [30]) by using an AC module. The simulation system had 12 negative charges. Then, six divalent cations (Ca²⁺/Mg²⁺) were added into the above frame to compensate for excess negative charges. The Smart method was used to minimize the energy of the system. The Smart method integrates the steepest descent method, conjugate gradient method and Newton method. It can shorten the time required for configuration optimization and improve the optimization accuracy. An NPT ensemble was used in the simulation system. The temperature was controlled at 298K by Andersen [26]. The initial velocity of each molecule was randomly generated by the distribution. The Ewald method was used to calculate the long-range electrostatic interaction, and an atom-based method was used to calculate the van der Waals interaction. The time step was 1 fs, and the trajectory information was recorded every 1 ps. The total simulation time was 1 ns. All MD simulation systems were completed by using the Materials Studio 2019 package. In the simulation system, the nitrogen and oxygen in the second period of the periodic table of elements, directly connected to the hydrogen atom, were hydrogen bonded, including -OH, -COO⁻, -NH₂ and water molecules. The interaction energy and the number of hydrogen bonds were calculated using a self-made program.

In order to study the interaction between metal ions and TNB to reveal the difference in the translational mobility of ions around TNB molecules, different ion situations were chosen to analyze the diffusion coefficients of ions and TNB in the first 200 ps. The mean square displacement (MSD) [31,32] and diffusion coefficient (D) [33] were obtained by fitting the curves. MSD was determined by

$$6Dt={<| \mathrm{r}}_j(t{)-\mathrm{r}}_j{(0) |}^2>$$

(1)

and (D) was defined as

$$D=\underset{t\to \infty }{\lim }\frac{1}{6t}{<| \mathrm{r}}_j(t{)-\mathrm{r}}_j{(0) |}^2>$$

(2)

where r_j(t) is the atomic position of the target molecule at time t and < > indicates that all ion position parameters in the simulated system are averaged.

The interaction between TNB and metal ions was investigated because RDFs enable a general understanding of the distance and strength of the interaction between atoms. Therefore, in this study, RDFs analysis was carried out. The coordination number (n) [34,35] refers to the number of coordination atoms around the central atom in a compound, and the n of atomic interaction was obtained by the following equation:

$$n={\displaystyle {\int }_{R_1}^{R_2}4\mathsf{\pi }{r}^2}\rho \mathrm{g}(r)\mathrm{d}r$$

(3)

where n is the number of specific atoms in a peak, r is the distance from the atom to the target atom, ρ is the number density of the atom in the main phase and g(r) is the RDFs from the atom to the target atom.

3. Results and Discussion

3.1. System Equilibrium and Interaction Energy

The completed simulation system is shown in Figure 2a,b. The distance between the metals and TNB is very small, and they become aggregated. Important criteria determining whether the NPT simulation system was in equilibrium were as follows: (1) temperature balance—the standard deviation of temperature change was required to be less than 5%; (2) energy balance—the energy was required to be constant or to fluctuate along a constant value; and (3) density balance—the density change was required to be less than 5%. Through calculation, it was found that the temperature, energy and density deviation of the last 200 ps were less than 3%, and the system reached equilibrium. Therefore, the subsequent analysis could be carried out, and the data were valid.

Figure 3 shows the MSD of metal ions and TNB in different systems. It can be seen that TNB in water molecules has a low migration ability, which is mainly caused by the strong hydrogen bonds and electrostatic interaction between TNB and water. The corresponding D of molecules and ions are calculated and listed in

Table 2. Diffusion coefficients of TNB and ions in different systems.

	Ca2+	Mg2+	TNB (Ca2+-TNB)	TNB (Mg2+-TNB)
D/10−9 m2 s−1	1.179	0.709	0.399	0.230
R2	0.998	0.999	0.996	0.983

. A comparison between Ca²⁺-TNB and Mg²⁺-TNB shows a larger D value of Ca²⁺ and a stronger diffusion ability of Ca²⁺, which is mainly because Mg²⁺ is more electronegative than Ca²⁺. The stronger electronegativity of Mg²⁺ leads to a smaller ion radius. Therefore, the charge of Mg²⁺ nucleus is higher. Mg²⁺ has a stronger electric field strength and a greater interaction with water molecules [13] (Figure 4).

3.2. RDFs and Coordination Number

RDFs and coordination numbers are calculated and summarized in Figure 5 and

Table 3. Radii of the first and second coordination shells around Ca²⁺/Mg²⁺ and coordination numbers of association with COO⁻ and OH of TNB and water molecules.

System	Hydrophilic Group	Peak1/Å	Peak2/Å	n1	n2	Total
Ca2+-TNB	COO−	2.1	3.7	1.59	0.10	1.69
	OH	4.1	/	0.27	/	0.27
	H2O	2.3	4.5	0.57	0.30	0.87
Mg2+-TNB	COO−	1.7	3.5	0.92	0.31	1.23
	OH	3.7	/	0.25	/	0.25
	H2O	1.9	4.3	0.74	0.27	1.01

. As shown in Figure 5a, the first sharp peaks of “Ca²⁺-carboxyl”, “Ca²⁺-water” and “Ca²⁺-hydroxyl” are located in d_Ca,carboxyl ^I = 2.1 Å, d_Ca,water ^I = 2.3 Å and d_Ca,OH = 4.1 Å, respectively (d_A,B ^I: distance of first sharp peaks between A and B). Each Ca²⁺ has 1.59 O atoms of the COO⁻ groups, 0.27 O atoms of hydroxyl and 0.57 H₂O molecules in the first coordination shell on average. The second peaks of “Ca²⁺-carboxyl” and “Ca²⁺-water” are located in d_Ca,carboxyl ^II = 3.7 Å and d_Ca,water ^II = 4.5 Å. Each Ca²⁺ has 0.10 O atoms of the COO⁻ groups and 0.30 H₂O molecules in the second coordination shell. The RDFs of Mg²⁺ are shown in Figure 5b; the first sharp peaks of “Mg²⁺-carboxyl” and “Mg²⁺-water” are located in d_Mg,carboxyl ^I = 1.7 Å and d_Mg,water ^I = 1.9 Å. Each Mg²⁺ has 0.92 O atoms (COO⁻) and 0.74 H₂O molecules in the first coordination shell on average. There are 0.31 O atoms (COO⁻) and 0.27 H₂O molecules in the second coordination shell of Mg²⁺ on average (d_Mg,carboxyl ^II = 3.5 Å, d_Mg,water ^II = 4.3 Å). Because the electronegativity of Mg²⁺ is greater than that of Ca²⁺, the shell of Mg²⁺ hydration (d_Mg,water ^I = 1.9 Å, d_Mg,water ^II = 4.3 Å) is thinner than that of Ca²⁺ (d_Ca,water ^I = 2.3 Å, d_Ca,water ^II = 4.5 Å) [13]. The number of water molecules in the first hydration shell is an effective parameter for evaluating the interaction strength between polar groups and water [34]. The total coordination number of Ca²⁺ and H₂O is 0.87, which is smaller than that of Mg²⁺ and H₂O (1.01). The binding of water molecules to Ca²⁺ is not as close as other combinations because the charge radius ratio of Ca²⁺ is smaller than that of Mg²⁺, which is consistent with the above analysis. It can also be seen that the distance between different shell layers of the hydration shell is about 2.2~2.4 Å, which is very close to the length of the hydrogen bond. Therefore, it is thought that the water molecules in the first water shell can form hydrogen bonds with the water molecules in the second water shell. In other words, the first and second water shells are connected by hydrogen bonds [31]. The coordination number of Mg²⁺ and the hydroxyl group is 0.25 (d_Mg,OH = 3.7 Å); this is smaller than that of Ca²⁺ and the hydroxyl group (0.27), and the binding ability of these two ions to the hydroxyl group is low.

All the above results can be further verified by the geometric analysis of a snapshot from the MD simulation, which is shown in Figure 6. The snapshot of simulated trajectories in Figure 6a shows that both monodentate (2.133 Å) and bidentate (2.029 Å and 2.154 Å) coordination were present between Ca²⁺ and COO⁻ groups of different TNB molecules. The geometry of the Ca²⁺-TNB complex reveals the strong bidentate coordination of COO⁻ groups of two different TNB molecules to Ca²⁺. This shows that Ca²⁺ is complexed with two or more carboxylic acids, and this structure is consistent with previous simulation results [26,36]. In Figure 6b, it is seen that a monodentate coordination exists in the Mg²⁺-TNB complex, and the coordination distance (1.663~1.719 Å) is shorter than that of the Ca²⁺-TNB system (2.029~2.154 Å).

3.3. Solvent Accessible Surface Areas and Hydrogen Bonds

As shown in Figure 7, the SASAs of TNB molecules in the studied systems under a neutral condition are calculated, and the solvent radius is set to 1.4 Å [26]. Hydrophilic functional groups of TNB, such as carboxyl, hydroxyl and amino groups, play an important role in the system, and the hydrophilic region accounts for the majority of total SASAs. The interaction energy between TNB and Mg²⁺ decreased to −668.28 kcal/mol at the end of the simulation. However, in the simulated Ca²⁺-TNB system, the interaction energy decreased slowly and finally decreased to only −282.72 kcal/mol, indicating that Ca²⁺ could not interact with TNB effectively (Figure 8). The SASA of Ca²⁺-TNB was 2300.52 Ų, and that of Mg²⁺-TNB was 2145.10 Ų. A smaller interaction between TNB and Ca²⁺ leads to an increase of the hydrophilic area of the Ca²⁺-TNB system because dispersive TNB will expose more hydrophilic groups. The above analysis shows that the interaction between TNB and Ca²⁺ is much weaker than that of Mg²⁺-TNB.

The surface of TNB molecules easily forms intermolecular or intramolecular hydrogen bonds, which accelerates the aggregation of TNB molecules. The number of hydrogen bonds between functional groups of TNB was obtained using a self-made program, and the process under a neutral pH condition was analyzed. The aggregation process of the TNB molecule in the metal ion solution with hydrogen bonds between functional groups is shown in Figure 9. These results show that the hydrogen bonds between different functional groups of TNB molecules may be an important factor for aggregation under a neutral condition. The hydrogen bonds generated by the hydroxyl group and carboxyl group amounted to 25.3 in Ca²⁺-TNB, and the hydrogen bonds generated by the hydroxyl and carboxyl groups amounted to 14.8 in Mg²⁺-TNB; these values are greater than those of the other three groups. Hydroxyl and carboxyl groups are the weakest of all groups in the system. Hydroxyl and carboxyl groups are dominant in the aggregation process of HA in all simulated metal ion solutions under a neutral condition.

4. Conclusions

The aggregation mechanism of HA in aqueous solutions containing Ca²⁺ and Mg²⁺ was studied by computational molecular dynamics under neutral conditions. The TNB model was used to represent HA. Under neutral conditions, the carboxyl group of TNB is deprotonated to COO⁻, and the negative charge of the system is neutralized by the addition of metal cations. Compared with Ca²⁺-TNB and Mg²⁺-TNB, the interaction between Mg²⁺ and water is stronger, and so the diffusion coefficient of Mg²⁺ is lower. The addition of metal ions reduces the electrostatic repulsion of TNB and coordinates with the COO⁻ and OH of TNB. Different coordination numbers have different effects on the aggregation behavior of TNB under neutral conditions. Metal ions can form intramolecular or intermolecular bridges with different functional groups through electrostatic synergy. Ca²⁺ and COO⁻ of TNB can form bidentate (2.029 Å and 2.154 Å) complexes, while Mg²⁺ can only form monodentate (2.133 Å) complexes with the COO⁻ of TNB. The interaction between TNB and Ca²⁺ is weaker, so the dispersed TNB in Ca²⁺-TNB will expose more hydrophilic groups. The SASAs of TNB in Ca²⁺-TNB (2300.52 Ų) were found to be larger than those of TNB in Mg²⁺-TNB (2145.10 Ų). The hydrogen bonds between hydroxyl groups and carboxyl groups are dominant in the aggregation process of HA in all simulated metal ion solutions under neutral conditions.