Molecular dynamics simulation of humic substances

Abstract

Humic substances (HS) are complex mixtures of natural organic material which are found almost everywhere in the environment, and particularly in soils, sediments, and natural water. HS play key roles in many processes of paramount importance, such as plant growth, carbon storage, and the fate of contaminants in the environment. While most of the research on HS has been traditionally carried out by conventional experimental approaches, over the past 20 years complementary investigations have emerged from the application of computer modeling and simulation techniques. This paper reviews the literature regarding computational studies of HS, with a specific focus on molecular dynamics simulations. Significant achievements, outstanding issues, and future prospects are summarized and discussed.

Introduction

Humic substances (HS) consist of a large variety of natural organic molecules that originate from the decomposition, and related microbial activity, of dead biological material, especially plant tissues [1]. HS are ubiquitous in the natural environment where they contribute to the regulation of many crucial ecological and environmental processes. For example, HS sustain plant growth and terrestrial life in general, and control the fate of environmental contaminants by acting as sorbents for toxic metal ions, radionuclides, and organic pollutants [2]–[4]. Furthermore, HS account for most of the planet’s organic material, and represent the most abundant reservoir of carbon [1],[5]. In fact, HS are receiving growing attention in recent years because of their potential role in land management strategies aimed at promoting carbon sequestration, to ultimately reduce atmospheric CO₂ and hence help tackle climate change [5].

Despite much research carried out over many decades now, the detailed nature of HS is still not fully understood. While main molecular building blocks have long been identified as hydrocarbon, quinone, phenol, catechol, and sugar groups (Figure 1), the crucial issue of how exactly these chemical moieties are organized at the molecular and supramolecular levels is still debated. A complicating factor is represented by the high variability of HS, as their composition depends on the specific ecosystem where they originate, in terms for instance of vegetation, climate, and topography. Following a polymer analogy, a popular early hypothesis described HS as collections of organic macromolecules, with molecular weights of up to several tens of thousands or even hundreds of thousands of grams per mole [1],[6],[7]. However, the current consensus, supported by the most recent experimental evidence, describes HS as heterogeneous supramolecular mixtures of relatively small molecules, with molecular weight of a few thousands of grams per mole, which associate dynamically through weak (non-covalent) interactions, especially hydrogen bonds and hydrophobic forces [8]–[14]. According to this hypothesis, HS are also capable of self-assembling into micellar structures, whereby an inner hydrophobic core is shielded from outer water through interfacial hydrophilic regions [15]–[17].

Figure 1

Common HS chemical groups. Chemical structures of the main molecular building blocks forming HS.

Over the past 20 years, traditional experimental investigation of HS has been compounded by various computer modeling and simulation approaches. This review focuses mainly on computational studies of HS conducted using the molecular dynamics method. After an introduction to the main methodological aspects, the available literature is categorized, summarized, and critically discussed.

Review

The molecular dynamics simulation method

Molecular dynamics (MD) is a computer simulation technique which is widely used in science and engineering, and is employed to obtain equilibrium and transport properties for collections of discrete particles. MD is a powerful method to simulate matter at the molecular scale; applications can be found for a wide range of systems, from simple gases and liquids [18]–[22] to various complex materials including proteins [23]–[26], lipid membranes [27]–[35], polymers [36]–[39], and carbon nanostructures [40]–[42]. Popular computer programs that implement MD include LAMMPS [43],[44], GROMACS [45], AMBER [46], GROMOS [47], DL_POLY [48], and CHARMM [49]. In this section, the main aspects of the MD method are summarized; more details can be found in dedicated books [50]–[54] and review articles [55]–[58].

The key components of the main MD algorithm are reported in Algorithm 1.

The first stage typically involves initializing the calculation by supplying the computer program with the coordinates of all atoms in the system (x), together with the models (V) which determine how the atoms interact. Such models are typically called potentials, or force fields. It should be noted that the focus of this summary is on fixed-charge biomolecular/organic force fields, as these are used predominantly in the simulation of humic substances. However, several other types of force fields exist, as documented extensively in the literature; in particular, significant progress has been recently made on polarizable models [59].

The standard form of a force field is

$$V=V_{\text{bonded}}+V_{\text{nonbonded}}$$

(1)

where V_bonded defines the (intramolecular) interactions between atoms covalently bonded to each other, and V_nonbonded defines the intermolecular interactions. In particular, V_bonded typically contains simple harmonic terms, whereby for instance a covalent bond is modeled with a potential representing a mechanical spring:

$$V\left(l\right)=\frac{1}{2}k{(l-l_{\text{eq}})}^2$$

(2)

with k the rigidity constant, l the actual bond length, and l_eq the equilibrium bond length. Additional terms are used to control the angles within groups of three consecutively bonded atoms, as well as torsions involving groups of four bonded atoms. Rigidity constants and equilibrium values are typically optimized to reproduce vibrational and conformational properties from experiments or ab initio quantum mechanics calculations. The nonbonded components of a force field normally describe van der Waals and electrostatic interactions. The van der Waals interaction between each pair of atoms a distance r apart is modeled with the Lennard-Jones potential:

$$V\left(r\right)=4\epsilon \left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]$$

(3)

with σ defining the collision distance and ε the attractive energy. Atoms interacting through the Lennard-Jones potential can be thought of as spheres which repel each other at short range (a feature that mimics overlap between electron clouds) and attract each other at long range (corresponding to attractive dispersion forces). The Lennard-Jones parameters σ and ε are normally optimized to reproduce thermodynamic data from experiment, including liquid densities and enthalpies of vaporization. Electrostatic forces between each pair of atoms i,j located a distance r apart are modeled with the Coulomb potential:

$$V\left(r\right)=\frac{Q_i{Q}_j}{4\pi {\epsilon }_{0}r}$$

(4)

with Q _i and Q _j the corresponding charges and ε₀ the permittivity of free space. In general, charges are assigned empirically to reproduce experimental observables such as known multipole moments or thermodynamic properties.

The second stage in Box 1 represents the main computational loop of a molecular dynamics simulation. The first part of the loop involves computing the force on each atom, which is obtained from the gradient of the potential V. In the second part of the loop, the forces are used to move each atom forward in time. This is done by solving numerically Newton’s equations of motion. For example, considering one of the most widely used algorithms [60], given the force f(t) and velocity v(t) at the current time t, each atom is moved one timestep, from position x(t) at time t to position x(t+Δ t) at time t+Δ t, according to

$$\mathbf{x}(t+\mathrm{\Delta t})=\mathbf{x}\left(t\right)+\mathrm{\Delta t}\left[\mathbf{v}\right(t)+\mathrm{\Delta t}\mathbf{f}(t)/2m]$$

(5)

where m is the atom’s mass. Each iteration of this second stage advances the system in time by a typically small timestep (Δ t=10⁻¹⁵ s), and thus complete simulations normally require up to 10⁶−10⁹ iterations.

The third stage in Box 1 refers to the output data generated by the simulation. In particular, a trajectory is obtained consisting of consecutive snapshots of the system taken at regular time intervals during the simulation. The output trajectory is typically analyzed using statistical mechanics to obtain various thermodynamical and dynamical properties of interest, such as energy terms, average and local densities, diffusion and viscosity coefficients, mechanical parameters, and electrical potentials.

Simulations of humic substances

In this section, a number of representative molecular dynamics investigations of HS reported in the literature are reviewed. The studies considered are organized into different subsections corresponding to different types of systems investigated. It should be noted that in the literature HS are sometimes referred to as ‘NOM,’ ‘SOM,’ or ‘DOM’ [61],[62]. These acronyms stand for natural organic matter (NOM), soil organic matter (SOM), and dissolved organic matter (DOM). Specifically, NOM refers to a complex mixture of organic material that is found in water, soils, and sediments [63],[64]. SOM refers to all carbon-containing substances in soils [65]. DOM is defined as the portion of NOM which passes through a filter of 0.45- μm pore size [66]. For all of these three categories, HS represent a major constituent [7],[65].

Modeling HS and their fundamental properties

Computational molecular models for HS began to appear in the 1990s. One of the first and most significant contributions involved the development of the TNB model (from Temple-Northeastern-Birmingham) [67]–[70]. The TNB model was aimed at representing a typical, ‘average’ HS molecule, with a chemical composition determined on the basis of analytical measurements. Specifically, the TNB model comprises three carboxylic groups, three carbonyl groups, two phenolic groups, two amine groups, and four other R-OH alcohol groups, for a total chemical formula of C ₃₈H₃₉O₁₆N₂ and a molecular weight of 753 g mol ⁻¹. The molecular structure of the TNB molecule is reported in Figure 2. Earlier simulations of the TNB molecule were performed in vacuo, meaning that no other substance (such as water) was included in the system. In particular, these calculations aimed at finding the most realistic (optimized) geometrical arrangements by minimizing the potential energy of the system [67],[70].

Figure 2

The TNB model. Molecular structure of the Temple-Northeastern-Birmingham (TNB) model for HS. Reprinted with permission from Sein et al. [70]. Copyright 1999 American Chemical Society.

Alvarez-Puebla and Garrido [71] studied the effect of pH on the aggregation of the TNB humic model [68]–[70]. By simulating the aggregation process, it was found that the molecular size increased with increased pH values due to intramolecular electrostatic repulsion, while the size of the aggregates decreased with increased pH because of increased repulsive intermolecular interactions [71]. Alvarez-Puebla et al. [72] subsequently developed a modified version of the TNB model [68]–[70] aimed at better representing a set of experimental data on HS composition. A series of simulations were conducted to investigate HS aggregation as a function of the model’s ionic states, both in vacuo and in aqueous solution [72].

Leenheer et al. [73] developed a model of HS and used it for the interpretation of experimental data on metal-HS association. Kubicki and Apitz [69] later used the Leenheer model to predict equilibrium structures through classical molecular mechanics and quantum calculations and to test the effect of the specific computational methodology on the structures obtained. The Leenheer model [73] was also adopted by Porquet et al. [74] to investigate hydrogen bonding and clustering of neutral HS molecules in water.

While HS molecular models are typically constructed by assembling atoms manually into the desired compositions and geometries, an interesting alternative approach was developed by Diallo et al. [75], who proposed a series of structural models for soil HS by processing an extensive set of experimental data through an automated algorithm, implemented into specifically designed computational software. The model molecules obtained were relatively small, with an average molecular weight of ≈1,000 g mol ⁻¹. As opposed to the more traditional approaches, the method of Diallo et al. [75] has the advantage that only the appropriate isomers are selected when multiple structures can be deduced from the same set of analytical data.

The specific role of water in its interactions with HS was investigated by Aquino et al. [76],[77]. In this work, HS were represented by simple hydrocarbon chains containing hydrophilic (carboxyl) groups. The MD simulations showed that distant hydrophilic groups can be cross-linked by water molecular bridges [76],[77].

In general, the HS models reviewed so far include molecules characterized by relatively low numbers of atoms, on the order of 100, yielding molecular weights of ≈1,000 g mol ⁻¹. However, the development of molecular models comprising substantially larger numbers of atoms has also been reported. In particular, significant work in this context has been carried out by Schulten and coworkers [61],[65],[78], who proposed model HS molecules comprising over 1,000 atoms and corresponding molecular weights of up to and over 10,000 g mol ⁻¹.

Specifically, Schulten and Schnitzer [78] designed a SOM molecule by hydrogen bonding a humic structure to a hexapeptide and a trisaccharide, obtaining a compound with molecular formula C ₃₄₂H₃₈₈O₁₂₄N₁₂ and corresponding molecular weight of 6,651 g mol ⁻¹. Based on the SOM model, Schulten [61] subsequently proposed a model for DOM and investigated complexes with xenobiotic substances. In general, xenobiotics are substances such as pollutants or pesticides, which are found in the environment yet are not naturally expected to be present. Schulten [61] performed molecular simulations of systems including DOM, water, and the xenobiotic pentachlorophenol (a pesticide), atrazine (a herbicide), and DDT (an insecticide). Geometry optimization calculations were performed to analyze energetics and hydrogen bonds. It was found that van der Waals forces and hydrogen bonds were the main contributors to the temporary retention of xenobiotic substances in DOM [61].

By considering the compounds observed in pyrolysis and other experimental studies of organic matter, Schulten further refined a molecular model for DOM [65]. In particular, a prototypical DOM molecule was obtained taking into account the most frequently occurring molecular building blocks, as well as an averaged elemental composition. Various organic functional groups found in HS were employed, including aromatic, alkyl, carboxyl, ketone, quinone, phenol, alcohol, ether, amine, amide, and heterocyclic N and S functional groups. The resulting molecule contained 1,262 atoms, with molecular formula C ₄₈₇H₄₉₂O₃₀₆N₁₅S₂ and corresponding molecular weight of 11,515 g mol ⁻¹; a 3D representation of the Schulten DOM molecule [65] is reported in Figure 3. In the simulations, 35 water molecules were also added to the surface of the DOM molecule. Energy minimization calculations were carried out to study the contributions from van der Waals and electrostatic terms, as well as to characterize hydrogen bonds [65].

Figure 3

The Schulten DOM molecule. Snapshot of a 3D representation of the Schulten DOM model [65]. Color codes for atom types are as follows: carbon (cyan), hydrogen (white), oxygen (pink), nitrogen (blue), and sulfur (yellow). Reprinted with permission from Schulten [65]. Copyright 1999 Elsevier.

Sutton et al. [62] refined the Schulten DOM molecule, obtaining a compound with molecular formula C ₄₄₇H₄₂₁O₂₇₂N₁₅S₂ and corresponding molecular weight of 10,419 g mol ⁻¹. This model was then simulated under conditions of increased hydration, typical of natural soil and water environments. In particular, the DOM molecule was surrounded by water molecules corresponding to a hydration layer of ≈5 Å thickness. The systems were shown to reproduce experimental physical and chemical properties of HS for several characteristic environmental conditions of soil. Specifically, results were obtained for density, hydrogen bonds, radius of gyration, and the Hildebrand solubility parameter [62].

For the HS molecular models considered in this review, Table 1 reports the corresponding chemical formula and molecular weight.

Table 1 Properties of model HS molecules

HS in complex with soil minerals

In general, HS are stabilized by their association with soil minerals, which prevent microbial attack and resulting rapid decomposition of HS [79],[80]. As a consequence, HS adsorption to minerals regulates the presence of carbon in soils [79]. An improved understanding of organo-mineral interactions is thus highly desirable, as it could lead to new strategies for soil carbon retention and sequestration through stabilization of HS.

A pioneering MD simulation study in this area was carried out by Teppen et al. [81], who investigated trichloroethylene (C ₂HCl₃), taken as a basic compound representative of organic material, adsorbed on clay mineral surfaces (kaolinite and pyrophyllite) in the presence of water. By considering different levels of hydration, it was found that water can outcompete C ₂HCl₃ for adsorption at the clay surface [81].

Shevchenko et al. [82] simulated organo-mineral aggregates in water using a NOM model based on an oxidized lignin-carbohydrate complex. MD simulations were conducted using the simulated annealing approach, whereby structural optimization is obtained by cooling-heating cycles which allow energy barriers to be overcome, eventually leading to optimized geometries.

To investigate further the nature of HS-mineral interactions, Sutton and Sposito [79] simulated complexes comprising the Schulten DOM molecule [65] and Ca-montmorillonite, a clay mineral of the smectite group. In particular, two DOM-Ca-montmorillonite systems were constructed. The first system included a C ₄₄₇H₄₂₁O₂₇₂N₁₅S₂ protonated DOM molecule, a 32-unit Ca-montmorillonite clay layer, 12 Ca ²⁺ ions, and 543 interlayer water molecules. The second system comprised a C${}_{447}{\mathit{\text{H}}}_{345}{\mathit{\text{O}}}_{272}{\mathit{\text{N}}}_{15}{\mathit{\text{S}}}_2^{76-}$ DOM polyanion, a 32-unit Ca-montmorillonite clay layer, 50 Ca ²⁺ ions, and 852 interlayer water molecules. In both systems, the DOM molecule was inserted into one of the clay interlayers. From their simulations, Sutton and Sposito [79] were able to ascribe the stabilization of organo-mineral systems to significant direct hydrophobic and hydrogen bonding interactions between organic and mineral groups. A simulation snapshot from this work is reported in Figure 4.

Figure 4

DOM-montmorillonite system. Snapshot from a simulation of protonated DOM-montmorillonite system. Water molecules are represented with cylinders, whereas DOM and clay are represented as balls and sticks. Color codes for atom types are as follows: carbon (gray), hydrogen (white), oxygen (red), nitrogen (blue), sulfur (yellow), and calcium ions (brown). Reprinted with permission from Sutton and Sposito [79]. Copyright 2006 Elsevier.

Petridis et al. [83] modeled an Al ₂O₃ mineral surface in contact with the organic compounds stearic acid and glucose. The aim of this work was to study the mechanism by which glucose accumulates in a layer between Al ₂O₃ and stearic acid, as observed experimentally. The simulations conducted revealed that glucose deposits onto Al ₂O₃ driven by a lower entropic penalty with respect to stearic acid [83].

HS and metal ions

Understanding the association between HS and metal ions is a crucially important issue, as this process controls the speciation, solubility, and toxicity of trace metals [84],[85].

Sutton et al. [62] simulated systems including the Schulten DOM molecule [65], water, and the Na ⁺ and Ca ²⁺ ions. It was found that Ca ²⁺ ions associate more strongly than Na ⁺ with the carboxylate groups of the humic molecule. Moreover, Ca ²⁺ was shown to promote better hydration of the humic molecule [62].

Alvarez-Puebla et al. [86] studied the interaction between brown humic acid (BHA) with Cu ²⁺, Ni ²⁺, and Co ²⁺ ions. BHAs are the most polar and soluble components of HS, because of their high content in carboxylic and phenolic acidic groups. The BHA structure was developed based on the TNB model [68]–[70]. Specifically, a BHA polymer was obtained by concatenating 13 TNB monomer units. Due to the high computational cost of simulating such a large molecule, Alvarez-Puebla et al. [86] did not include solvating water, although its effect was approximated by introducing frictional forces through a Langevin scheme [87]. The BHA was observed to display higher affinity for Cu ²⁺ (most reactive), followed by Co ²⁺, and then by Ni ²⁺ (most inert). This behavior was attributed to electrostatic retention, a mechanism consistent with both experimental and simulation results [86].

Xu et al. [88] performed molecular dynamics simulations of complexes involving NOM and metal ions. In particular, they investigated the interactions between Cs ⁺ and Cl ⁻ with NOM in water. For NOM, they adopted the TNB model [68]–[70]. Several simulations were performed for a range of metal ion concentrations. A representative simulation snapshot is reported in Figure 5. The data obtained showed that Cs ⁺ associates with NOM through rapid exchange with the bulk solution, whereas Cl ⁻ does not significantly associate with NOM; these results were found to be consistent with nuclear magnetic resonance experiments [88].

Figure 5

Hydrated TNB molecule interacting with ions. Simulation snapshot from a system comprising the TNB molecule [70] in 0.3 M Cs ⁺ aqueous solution. Color codes for atom types are as follows: carbon (brown), hydrogen (gray), oxygen (red), nitrogen (blue), Cs ⁺ (green). Water molecules are shown in transparent representation. Reprinted with permission from Xu et al. [88]. Copyright 2006 Elsevier.

The study by Xu et al. [88] was extended by Kalinichev and Kirkpatrick [89] and by Iskrenova-Tchoukova et al. [90], who considered the Na ⁺, Mg ²⁺, and Ca ²⁺ ions. It was found that metal-NOM binding is primarily driven by electrostatic attraction between the positive ions and the negatively charged carboxylate groups of the NOM molecule (whereas phenolic groups were not significant binding sites). Moreover, the propensity for metal-NOM aggregate formation was found to be correlated with the charge to radius ratio and the size of the ions [89].

A rather original methodological study was performed by Kalinichev et al. [91], who considered the effects of different models and system sizes on the simulation results for a NOM-Ca ²⁺ association process. In particular, they tested combinations of the force fields CVFF [92], CHARMM [93], and AMBER [94], with the water models SPC [95] and TIP3P [96]. The properties considered, which included radial distribution functions and potentials of mean force, were found to be fairly robust with respect to the different model parameters used [91].

HS and contaminants

Antimicrobials make up a large proportion of the contaminants detected in the environment [97]–[99]. The occurrence of antimicrobials in soil and water is caused by their widespread use in agriculture and medicine [100]–[103], as well as their presence in a wide range of healthcare and household goods [104]–[106]. The detrimental effects of antimicrobials include the disruption of key microbial processes in soil, toxicity to organisms, and the development of microbial resistance [107]–[110]. These problems are significantly mitigated when antimicrobials are adsorbed in organic matter, such as HS. To gain insights into the adsorption process, Aristilde and Sposito [111] carried out molecular dynamics simulations of the binding of the antimicrobial ciprofloxacin by HS. Ciprofloxacin is a frequently prescribed antibiotic commonly found in hospital wastewaters [112]. Regarding the HS component, Aristilde and Sposito [111] used the Schulten DOM model [62],[65]. The simulations showed that the ciprofloxacin-HS association involved the disruption of original hydrogen bonds within the DOM molecule and their replacement with intermolecular hydrogen bonds with ciprofloxacin [111].

Another class of ubiquitous contaminants is represented by polycyclic aromatic hydrocarbons (PAHs), which are highly toxic compounds that form as a result of the combustion of organic fuels such as coal, oil, and natural gas. It has been shown in a number of studies that organic matter can regulate the transport, fate, degradation, and bioavailability of PAHs [113]–[118]. Saparpakorn et al. [119] investigated by simulation the binding of PAHs to different HS models; in particular, they simulated Schulten’s SOM molecule [78] and implemented models for earlier molecules proposed by Buffle et al. [120] and by Stevenson [1]. The simulations performed aimed at quantifying the role of intermolecular interactions, as well as docking energies and binding modes [119].

Schulten et al. [121] modeled complexes of HS and the xenobiotic diethyl phtalate (DEP), with the objective of investigating the sorption process. Interactions were studied between a single HS molecule and an increasing number of DEP molecules, from 1 to 30. From their simulations, Schulten et al. [121] were able to quantify the sorption process in terms of the different contributions from electrostatic, van der Waals, and hydrogen bonding interactions. In particular, sorption inside free-volume pockets of HS was observed to take place between a single HS molecule and up to seven DEP molecules, whereas additional DEP molecules were adsorbed at the HS surface [121].

Another category of contaminants of increasing relevance is represented by carbon nanoparticles. The general use of carbon nanomaterials in industry is rapidly growing, raising health and environmental concerns which demand quantitative assessment. Wang et al. [122] used molecular simulations to investigate the interactions between DOM and fullerene (C ₆₀), a typical carbon nanoparticle. Fullerene plays a role in a wide range of industrial applications and is known to display some degree of toxicity [123]–[126]. Wang et al. [122] selected seven small organic molecules, representative of main DOM building blocks, and characterized their interaction with C ₆₀ in terms of adsorption energy and water solubility; it was found that the presence of DOM can stabilize C ₆₀[122]. Further insights into DOM-C ₆₀ systems were obtained by Sun et al. [127], who considered a range of low molecular weight organic acids as key components of DOM. By estimating adsorption energies, it was observed that aromatic acids interact more strongly with C ₆₀ than aliphatic acids [127]. Wu et al. [128] simulated an aggregate comprising ten C ₆₀ molecules associated with a small HS molecule. In particular, the HS model was constructed by connecting a benzoic group to a hydrocarbon tail. It was found that hydrophobic and π- π interactions were the two main mechanisms of association [128]; a simulation snapshot from this work is reported in Figure 6.

Figure 6

Interaction between a small HS molecule and fullerene. Snapshot from a MD simulation by Wu et al. [128]. Color codes for atom types are as follows: carbon (gray), hydrogen (white), oxygen (red). Reprinted with permission from Wu et al. [128]. Copyright 2006 Elsevier.

To obtain insights into the sorption of volatile organic compounds into HS, Shih et al. [129] studied the interaction between the TNB humic acid model [68]–[70] and toluene (representative volatile organic compound) in vacuo. Specifically, the diffusion coefficient of toluene was characterized as a function of temperature from 300 to 400 K. The results obtained are in qualitative agreement with experiment, in that diffusivities were observed to increase with temperature. However, the experimental data were slightly overestimated [129].

HS and water filtration

Satisfying the world’s population need for clean and drinking water is one of the greatest challenges of our time. To address this challenge, it is paramount to develop and optimize industrial processes aimed at filtering and desalinating sea water and municipal waste water. The currently most promising filtration technology relies on membranes operating in reverse osmosis plants. In these processes, the presence of HS is a fundamental aspect to consider. In fact, a key problem that greatly limits the efficiency of current filtration membranes is fouling, a phenomenon whereby particles deposit and accumulate on the membrane surface ultimately causing a reduction in the filtering performance. A major category of fouling agents is represented by organic substances, particularly HS [130]–[135].

A number of MD studies have been devoted to different aspects of the fouling process. Ahn et al. [136] investigated the effects of metal ions on the adsorption of a NOM model [68] onto the surface of polyethersulfone membranes. It was found that divalent ions (Mg ²⁺ and Ca ²⁺) induce fouling by promoting aggregation of NOM molecules [136]. However, the interactions between NOM and the filtration membrane were not explicitly investigated.

The fouling of a polyamide membrane was investigated by Hughes and Gale [137],[138]. Specifically, they considered glucose and phenol molecules as representative HS fouling agents, as both glucose and phenol are common building blocks of HS. Membrane-foulant interactions were quantified in terms of free energies and hydrogen bonding. It was found that both foulants bind strongly to the membrane surface, with phenol sometimes diffusing through the membrane pores [137],[138]. A simulation snapshot from this study is reported in Figure 7; a phenol molecule can be seen penetrating the polymeric membrane.

Figure 7

Interaction between phenol and polymeric membrane. Snapshot from a MD simulation by Hughes and Gale [138]. A phenol molecule (colored yellow) permeates into a polyamide membrane (colored purple). Reprinted with permission from Hughes and Gale [138]. Copyright 2012 Royal Society of Chemistry.

Myat et al. [139] investigated possible specific mechanisms of interaction between representative organic foulants. Specifically, they focused on the biopolymer bovine serum albumin (BSA) [140] and the polysaccharide sodium alginate, taken to be representative of high molecular weight compounds typically found in surface and waste waters. Moreover, they considered the TNB humic acid model [68]–[70] as representative of HS. No water was explicitly included. Simulations of a BSA-HS complex revealed the presence of various electrostatic and hydrophobic interactions, as well as hydrogen bonding. On the other hand, analysis of an alginate-HS complex highlighted the presence of exclusively ion-mediated interactions. The simulation results were found to be consistent with corresponding experimental data [139].

Conclusions

Achievements, issues, and future prospects

Over the past 20 years, a growing number of computer models have been developed and applied to study many important structures and processes involving humic substances (HS), including their basic molecular properties [62],[65],[67]–[70],[75], their aggregation behavior [71],[72],[74],[76],[77], their interaction with various substances including minerals [79],[81]–[83],[86], ions [62],[85],[86],[88]–[91], and contaminants [111],[119],[121],[122],[128],[129],[141],[142], and their fouling capability in relation to membrane-based water filtration technologies [136],[138],[139].

These investigations yielded considerable molecular-level insights into the structure and function of HS, as summarized in the previous sections of this review. However, a few issues should be considered. In particular, it is important to bear in mind that none of the HS models developed so far correspond to real humic molecules. Rather, the models represent putative compounds obtained by assembling molecular building blocks which are known experimentally to be most prevalent in HS. Furthermore, several investigations, especially among the earliest simulations reported, focused on energy minimization calculations, with the aim of finding the most energetically favorable (optimized) conformations for a molecule or molecular aggregate [61],[65],[67],[70],[121],[142],[143]. However, it should be noted that energy optimization methods yield properties corresponding to a temperature of 0 K, as only the potential energy is considered, while there is no kinetic energy in the system. When temperature and thermal motion are important, as is typically the case for systems of organic and biological molecules, full MD simulations, while computationally more demanding than optimizations, are to be preferred. A final issue to highlight involves the fact that many simulation studies of HS did not include hydrating water (in vacuo assumption) [67],[69]–[71],[73],[78],[82],[121],[128],[129],[139],[144]–[150]. As already pointed out elsewhere [72],[143],[149],[151], HS are hydrated in reality, and water interactions with HS are likely to influence important properties. For example, the large molecular dipole of water is expected to interact strongly with HS polar groups, and hydrogen bonds between water and HS are expected to be prevalent. The presence of appropriate amounts of water in MD simulations of HS is therefore recommended.

In terms of future prospects, there is an expectation that specific HS structures will be accurately identified from experiment, opening up opportunities for MD simulations of realistic HS compounds. As a result, simulated systems will likely become larger and more complex, and hence also more computationally expensive. While this could represent an obstacle, there are reasons to be optimistic. From a hardware perspective, the continuous increase in computational power will keep extending the attainable simulation times and sizes. Moreover, ongoing research in multiscale methods [152]–[157] promises to substantially improve simulation efficiency in the near future. Self-assembly simulations of large numbers of different HS molecules might soon become a reality, opening up the opportunity to study and quantify atomic-level properties within realistic HS supramolecular structures.

More generally, the study of HS in the foreseeable future will have great relevance for several areas of key global importance. Owing to the role of HS in controlling CO ₂ in the ecosystem, advances in HS research could lead to new solutions for carbon capture and storage, thus contributing to address the urgent global challenge of increasingly rapid climate change [5]. Moreover, a better understanding of HS can be instrumental in increasing food production to satisfy the needs of a growing population [158], as well as in optimizing filtration technologies to obtain clean and drinking water [159]. While experimental research will always be essential, in the years to come, molecular simulations of HS are expected to become increasingly useful, particularly for providing a more detailed understanding of experimental observations, for guiding the design of new experiments, and for predicting properties and phenomena at the molecular scale.