Coarse-grained simulations of sulfanilamide and hexachlorobenzene mobility in soil organic matter

Abstract

Soil organic matter (SOM) is a complex and heterogeneous molecular system, crucial for soil health and ecosystem functioning. Therefore, the release of pollutant molecules into the environment poses a significant environmental threat. Mechanisms that are governed by the interactions of these pollutants with SOM at the molecular level remain largely unexplored. In this study, coarse-grained molecular dynamics simulations were employed to investigate the behavior of hexachlorobenzene (HCB) and sulfanilamide (SAA) in humic substance (HS) systems with varying compositions. Diffusion coefficients indicated a strong influence of water on SAA, with SAA displaying higher mobility, whereas HCB exhibited greater accumulation in the HS phase. Calculations of spatial distributions supported these observations, showing that SAA is predominantly situated in the water phase, while HCB's interaction was influenced by the hydrophobicity of the SOM system. Simulations in the microsecond range, which were possible by the coarse-grained representation, revealed temporary trapping of SAA in the SOM matrix. These were anti-correlated with water diffusion, while HCB's behavior was dominated by direct pollutant–SOM interactions. In general, the coarse-graining approach provides novel insights into the trapping processes of pollutants in SOM and offers a realistic representation of molecular interactions at larger spatial and temporal scales. The proposed method enhances the understanding of pollutant mobility in soil systems, thus enabling future studies on the ecological impact of pollutant–SOM interactions.

---

Environmental significance

Pollutants released into terrestrial ecosystems, including polar and non-polar organic pollutants, threaten soil functions and human health. Their chemical stability enables them to evade natural degradation; however, their behavior in soil is not well understood due to its structural and chemical complexity. Soil organic matter (SOM), which is rich in diverse functional groups, is a principal sorbent that governs pollutant mobility. In this study, we employ coarse-grained molecular dynamics to investigate the diffusion of hexachlorobenzene and sulfanilamide within SOM of different compositions. Simulations in the microsecond range provide mechanistic insight into interactions between pollutants and SOM and offer a transferable framework for investigating pollutant dynamics in SOM.

1 Introduction

Soil organic matter (SOM) is a significant soil component, comprising decomposed plant and animal residues, microbial biomass, and various organic compounds.¹ It contributes to improved soil structure stability, increasing its water infiltration and retention capabilities.^2,3 By enhancing the soil's nutrient-holding capacity, SOM plays a crucial role in maintaining soil fertility and supporting plant growth.⁴ Furthermore, SOM provides a habitat for soil microorganisms, which drive nutrient cycling and the functioning of the ecosystems.⁵

The soil's ecological functions are endangered by the release of chemical substances into the environment.^6,7 Persistent organic pollutants (POPs) are toxic chemical compounds with a long lifetime, resulting from their resistance to degradation through environmental processes.⁸ These compounds are widespread in the environment and can be transported over long distances through air and water currents.^9,10 Moreover, polar organic pollutants, such as pharmaceutical antibiotics and personal care products, are of particular concern due to their potential to cause harm to terrestrial ecosystems.^11–13 Numerous studies have been carried out to investigate the sorption of POPs and polar organic pollutants in soil using experimental and theoretical methods.^14–35

Understanding the interactions of pollutants within the soil matrix and the properties of SOM is crucial for optimizing remediation strategies, such as immobilization or bioremediation.³⁶ Traditionally, SOM has been closely studied by obtaining so-called humic substances (HSs) via alkaline extraction, which benefits from standardized protocols and reference samples provided by the International Humic Substances Society (IHSS).³⁷ Despite ongoing discussions about how accurately HSs can reflect native SOM,^38–41 they remain the subject of research aiming to systematically understand more about SOM formation, composition, and activities in diverse environments.^42,43

From a theoretical viewpoint, computational chemistry provides insight into processes at the molecular level and expands the understanding of the pollutant–SOM interaction. The modeling of SOM is challenging due to its heterogeneity and complexity, which has led to the ongoing development of various SOM models.^44–55 The Vienna Soil Organic Matter Modeler 2 (VSOMM2)^56,57 was recently released, which allows model systems of condensed phase HS to be generated, following a building block approach. Based on elemental and organic fractions as input parameters, heterogeneous systems of HS with different compositions can be investigated, thus approaching the properties of real SOM. These inputs (e.g., bulk C to N ratios and fractions of aromatic, aliphatic, or acidic functional groups) are commonly used to describe SOM heterogeneity and serve as constraints on the generated structures, consistent with established practices in supramolecular SOM modeling (see Gerzabek et al.⁵⁸).

While such systems consisting of larger, supramolecular structures can improve the accuracy of the models, the accessible time and length scales of simulations may be compromised. To reduce computational cost, coarse-graining techniques are often employed,^59–62 whereby multiple atoms are grouped into larger particles known as beads. This approach sacrifices atomistic detail for computational efficiency by reducing the number of interactions that must be evaluated. Coarse-graining strategies can follow either a top-down or bottom-up philosophy: in the top-down approach, interaction potentials are tuned to reproduce experimental observables or macroscopic properties, whereas in the bottom-up approach, they are derived systematically from atomistic simulations.⁶⁰ Furthermore, several coarse-grained force fields have emerged,^63–65 including the widely used Martini model,⁶⁶ which provides pre-parametrized beads to represent common chemical functionalities.

To apply the coarse-graining method to SOM systems, we have developed a tool^67,68 that converts systems generated by the VSOMM2 from an atomistic to a coarse-grained representation, using the Martini 3 framework.⁶⁹ On the one hand, this enables an abstraction of the SOM molecular structures, whose exact shape remains unknown. On the other hand, the simulation of diffusion processes in large SOM systems becomes possible, thus accounting for the heterogeneous nature of interaction sites of SOM. To the best of our knowledge, no molecular-level study has examined the interaction between pollutants and supramolecular SOM structures over extended time scales, with a focus on their diffusion and distinct interaction sites.

The present study fills this gap by investigating the mechanisms and key functional groups that drive the interaction and trapping of POPs and polar organic pollutants in SOM at a coarse-grained level of theory. Sulfanilamide (SAA) and hexachlorobenzene (HCB) were selected to represent distinct pollutant classes and interaction modes with SOM. SAA is a widely detected pharmaceutical with polar interaction sites, making it useful for probing site-specific interactions with carboxylic and phenolic moieties.¹⁵ In contrast, HCB is a legacy POP with hydrophobic character that primarily partitions into aromatic or black-carbon-like domains.⁵⁰ Together, this pair of pollutants represents complementary sorption mechanisms, reflects environmental relevance, and allows comparison to previous studies. Their interactions with SOM model systems were analyzed, the latter being represented by HSs from different samples provided by the IHSS. Coarse-grained molecular dynamics (CGMD) simulations were employed to examine the mobility of pollutants in the SOM systems and their interaction with SOM functional groups. The analysis of the diffusion and trapping of the pollutant molecules within the SOM matrix provided insights into the underlying processes that govern their behavior in soil environments.

2 Materials and methods

2.1 Coarse-grained modeling approach

The parametrization of the coarse-grained pollutants and SOM models was conducted within the framework of the Martini 3 force field, which offers a wide range of available topologies for, e.g., small molecules,⁷⁰ biomolecules,^71–74 or carbon materials.^75,76 In previous work,^67,68 coarse-grained SOM systems based on VSOMM2 structures were introduced, employing an automated parametrization approach for the generation of the models. Validation of the coarse-grained models involved matching the probability distributions of bonded interactions, using the Swarm-CG tool,⁷⁷ and reproducing thermodynamic observables, specifically transfer free energies.

The mapping of pollutant molecules was performed with careful consideration of chemical functionality and molecular symmetry. Martini 3 bead types were chosen to reflect the polarity and hydrogen-bonding capacity of functional groups, aiming to reproduce experimental octanol–water partition coefficients as a first-order validation of solute–solvent interactions. To this end, coarse-grained transfer free energies were calculated using the Bennett acceptance ratio method,⁷⁸ and compared against experimental data if available. While the Martini framework cannot fully capture directionality in hydrogen bonding, it allows to approximate the overall partitioning behavior relevant for environmental partitioning. Bonded interaction parameters were further refined using Swarm-CG, with atomistic MD trajectories of each pollutant in explicit water serving as a structural reference during iterative optimization. A detailed description of this protocol is included in Section S1 of the SI.

The primary focus of the investigation is the mobility and interactions of SAA and HCB in SOM systems, as captured by CGMD simulations. Using elemental and carbon distribution data on HS samples from the IHSS, five characteristically distinct HS systems were generated using VSOMM2 and subsequently converted into a coarse-grained representation. Guided by the principal component analysis of IHSS reference materials conducted by Escalona et al.,⁷⁹ we selected HS samples that span distinct regions along the acidity and aromaticity axes, thereby capturing the chemical diversity and inherent heterogeneity of SOM within a concise and manageable set of models. Please refer to Table 1, and Tables S1 and S2 (SI) for further details regarding the HS systems, which comprise information on the compositions and functional group fractions.

Table 1 HS samples from the IHSS, which are modeled and analyzed in this work. The elemental and organic fractions of the HS models were used as input parameters for the generation with the VSOMM2. These include carbon (C), nitrogen (N), carbonyl, carboxylic (carbox.), aromatic (arom.), acetal, hetero-aliphatic (hetero-aliph.) and aliphatic (aliph.) fractions

Identifier	HS sample	Characteristics	Elemental and organic fractions
			C	N	Carbonyl	Carbox.	Arom.	Acetal	Hetero-aliph.	Aliph.
SOM I	Elliot soil I (FA)	High carboxyl and carbonyl content	0.5057	0.0272	0.1212	0.2525	0.3031	0.0101	0.0909	0.2222
SOM II	Elliot soil IV (HA)	Low carboxyl and carbonyl content	0.5951	0.0390	0.0100	0.1100	0.4100	0.0600	0.1400	0.2700
SOM III	Leonardite (HA)	High aromatic content	0.6381	0.0123	0.0800	0.1500	0.5800	0.0400	0.0100	0.1400
SOM IV	Pahokee peat II (FA)	Intermediate carboxyl and aromatic content	0.5131	0.0234	0.0373	0.1936	0.4037	0.0621	0.1128	0.1905
SOM V	Suwannee river II (FA)	High hetero-aliphatic & aliphatic content	0.5234	0.0067	0.0495	0.1683	0.2178	0.0594	0.1584	0.3466

---

As initial configurations, equilibrated SOM structures derived from 20 μs long CGMD simulations from previous work were adopted, in which different water contents were considered.⁶⁸ To enhance the sampling of the interactions between pollutants and HSs, the SOM systems with the highest water content considered were utilized, corresponding to a heavy atom fraction of water of approximately 0.567. Here, one coarse-grained water bead represents four atomistic water molecules. Each of the five systems consists of 1600 VSOMM2 fragments, with five fragments per HS molecule and 5000 water beads.

The pollutants SAA and HCB were studied separately and positioned into the different SOM systems. Each simulation included a single pollutant molecule, representing a low environmental concentration of the pollutant. Additional simulations at higher pollutant concentrations are reported in Section S9 of the SI. To ensure sufficient sampling of the different interaction sites of the heterogeneous SOM structures, 20 simulations were conducted, in which the initial position of the pollutant was varied. These positions were further refined to ensure a minimum distance of 0.75 nm between the pollutant and HS molecules, thus preventing initial artificial trapping within a local minimum of the SOM matrix. Therefore, the pollutant molecules initially tended to interact with the HS molecules from the water phase.

2.2 Computational details

The CGMD simulations were performed using the GROMACS simulation package version 2019.4.⁸⁰ Each system underwent an energy minimization, followed by an NPT equilibration of 10 ns and the production run, which lasted 1.5 μs. The temperature and pressure in the systems were 298.15 K and 1 bar, respectively. Further computational details and information about the simulation parameters are provided in Section S2 (SI).

Python packages were utilized to carry out the analysis, including NumPy,⁸¹ Matplotlib,⁸² MDAnalysis^83,84 and IPython.⁸⁵ The molecular structures were visualized using VMD.⁸⁶

3 Results and discussion

3.1 Pollutant parametrization

The assigned mappings for SAA and HCB are depicted in Fig. 1. First, the available standard Martini bead types were used for the atomic groups, including the aromatic moieties (TC5) and the amino group (TN6d) of SAA. The bead types for the sulfonamide group of SAA (Fig. 1(a) in blue) and for the HCB groups (Fig. 1(b) in green) were then determined by optimizing the octanol–water transfer free energies (ΔG), calculated with the models, against the experimental energies. In addition, the solvent accessible surface area (SASA) and Connolly surfaces (Fig. 1(c) and (d)) were calculated from the atomistic pollutant models to account for the reproduction of the molecular shape. A comparison of the values is shown in Table 2.

Fig. 1 Applied mapping to SAA (a) and HCB (b). The different colors and labels respectively represent the corresponding mapped groups and Martini bead types. Connolly surfaces of SAA (c) and HCB (d) are depicted for the atomistic (cyan) and coarse-grained models (orange).

Table 2 Comparison of model predictions with experimental and reference values for SAA and HCB. Octanol–water transfer free energy values (ΔG) are compared to experimental values,⁸⁸ while the solvent-accessible surface area (SASA) values are compared to values from atomistic simulations. The energy values are in kJ mol⁻¹ and surface areas in nm². The experimental transfer free energies were derived from partition coefficients

	SAA		HCB
	Model	Reference	Model	Reference
ΔG	−5.27	−3.54	35.47	32.71
SASA	4.15	4.07	4.37	4.37

---

The octanol–water values of the models closely matched the experimental values in the range of k_bT ≈ 2.5 kJ mol⁻¹, which is considered to be a good agreement for a Martini molecule parametrization.⁸⁷ The resulting bead types for the sulfonamide group and the HCB groups were P5 (strong polarity) and X1 (halogen, lowest polarity), respectively, correctly reflecting the polarity of the corresponding groups. In accordance with the Martini 3 parametrization approach, no hydrogen bonding label was applied for the bead representing the sulfonamide group because it can act as both a donor and an acceptor.⁶⁹ In Martini, the hydrogen bonding labels can influence the interaction with other beads to mimic an overall reduced or increased attraction. However, it is important to note that the information about the directionality of such interactions cannot be reproduced by the coarse-grained model in any case.

The SASA values between the atomistic and coarse-grained models of the pollutants can be considered to be in sufficient agreement. The bond distance between the HCB groups obtained from the Swarm-CG optimization was increased by 12.7% to achieve a better agreement of the SASA values and to optimize the match between the Connolly surfaces of the atomistic and coarse-grained model. Further information can be found in Section S1 of the SI.

3.2 Diffusion coefficients

To assess the mobility of pollutants, diffusion coefficients were determined based on the mean squared displacements (MSDs) of SAA and HCB. Details regarding the calculation can be found in Section S3 (SI). Fig. 2 presents the diffusion coefficients of the water beads (blue), SAA (orange), and HCB (green) across different SOM types. Please note that the error bars in Fig. 2 and all following Figures containing bar plots are the standard errors arising from the averaging over the 20 individual trajectories for each case.

Fig. 2 Diffusion coefficients of the water beads (blue), SAA (orange), and HCB (green) for the different SOM types.

For each SOM system, the diffusion of SAA correlates with that of water, with the highest diffusion coefficient observed in the SOM II system. The strong influence of water diffusion on SAA suggests that SAA molecules either spend extended periods in the water phase or remain in continuous contact with water. The elevated water diffusion in the SOM II system can be attributed to the system's more hydrophobic nature (low acidity), which promotes local phase separation and results in a more interconnected water phase,⁶⁸ thereby facilitating higher water mobility.

The overall diffusion of HCB is slower, reflecting its high accumulation potential in SOM.^50,51 The lowest diffusion coefficient is observed in the SOM III system. This result aligns with a theoretical study on HCB, which reported relatively low diffusion in Leonardite humic acid (LHA),⁵⁴ corresponding to SOM III in this study. In contrast, in the SOM I system, HCB diffusion closely matches that of water, indicating that HCB is predominantly located in the water phase for this SOM type. Conversely, SOM II exhibits the greatest disparity between HCB and water diffusion, probably due to a strong accumulation of HCB in the HS phase.

It is important to note that Martini coarse-grained simulations are carried out without explicit frictional or stochastic forces, unlike methods based on the generalized Langevin equation.^89–91 This omission leads to artificially accelerated dynamics, as the loss of atomistic-level friction prevents proper dissipation of energy. While such acceleration can be advantageous for exploring long-timescale conformational changes, it compromises the accuracy of real-time dynamics compared to atomistic models. Furthermore, the degree of acceleration may vary depending on the coarse-graining resolution and mapping scheme. Several approaches have been developed to quantify or correct for this acceleration, including relative entropy methods,^92,93 excess entropy scaling relationships,^94–97 biasing potentials^98,99 and assessments based on molecular roughness.^100,101

Therefore, the absolute diffusion coefficients shown in Fig. 2 are not directly comparable to experimental values if such data were available. Because coarse-graining in Martini models typically accelerates dynamics by roughly a factor of 3–8,¹⁰² relative speed-ups of SAA and HCB are not expected to differ by orders of magnitude. Nevertheless, the precise magnitude of this artificial acceleration may vary between the two models. Moreover, the environmental heterogeneity in SOM (e.g., highly charged vs. non-polar domains) may alter the effective friction and thus the magnitude of the coarse-grained speed-up. Consequently, direct comparison of their diffusion constants should be limited to qualitative trends rather than precise quantitative differences.

3.3 Coordination numbers

To assess the overall partitioning of the pollutants within the SOM systems and its connection to diffusion behavior, the average coordination numbers of water and HS beads around SAA and HCB were calculated and illustrated in Fig. 3 (see Section S4 (SI) for details of the calculation). The coordination number is a measure of how many beads of water or HS the pollutant molecules are surrounded by. Note that one water bead is representative of four atomistic water molecules.

Fig. 3 Coordination numbers of water (W) and HS beads around SAA (orange) and HCB (green), averaged over the 20 trajectories for each pollutant and each SOM type.

SAA exhibits similar behavior across different SOM types, with a slightly higher value for the HS beads in the SOM III system. Compared to HS, the coordination number of water is approximately four to five times higher, indicating a relatively low but non-negligible SAA concentration within the HS phase. This explains the observed correlation between the calculated diffusion coefficients of SAA and water in Fig. 2. Furthermore, the results suggest that at high water contents, SAA diffusion is primarily influenced not by direct interactions with HS molecules but rather by the structural and spatial arrangement of HS molecules, which modulates the mobility and flow of the surrounding water beads.

In contrast, the coordination number of HS beads around HCB varies significantly depending on the SOM type. The highest value is observed for the SOM II (low acidity) and SOM III (high aromaticity) systems. A notable contribution is also evident in the highly aliphatic SOM (SOM V), while the lowest coordination number is found in SOM I, which exhibits high acidity. Overall, HCB has a higher concentration in the HS phase compared to SAA, suggesting that HCB has a stronger binding affinity to SOM. This distinct behavior of SAA and HCB in the SOM II system is illustrated in Fig. 4. Unlike SAA, HCB diffusion is strongly influenced by direct HCB–SOM interactions.

Fig. 4 Merged final frames from 20 independent replicas of SOM II. Colored spheres mark the pollutant bead positions in each replica's final frame; SAA: yellow in (a), HCB: green in (b). The brown (HS) and translucent white (water) density isosurfaces were obtained by averaging bead positions across the 20 final frames, then contouring at a fixed threshold to indicate regions frequently occupied by HS and water. The surfaces are replica-endpoint densities (not time averages); pollutant spheres are not averaged. HCB endpoints predominantly fall within HS-rich regions, whereas SAA endpoints are often in water.

Interestingly, due to the low coordination number of HS beads, SOM I constitutes an exception, displaying a relatively low probability of HCB residing in the HS phase. This explains why, in SOM I, HCB exhibits diffusion comparable to that of water (see Fig. 2). Across all SOM systems, HS and water molecules are not perfectly mixed; instead, larger interconnected water regions exist, as previously noted for SOM II. The size of these regions is likely to be determined by the specific composition and arrangement of HS molecules. Particularly in SOM I, HCB diffuses within this water phase rather than integrating into the matrix formed by HS molecules, as identified from simulation trajectories.

At first glance, this behavior appears contradictory given HCB's strong hydrophobicity. However, it is presumably caused by the high polarity of the abundant carboxylate groups in SOM I, which prevent HCB from being absorbed into the HS phase. Moreover, the variations in coordination number across the HS phase can be linked to the acidity of the SOM system. High-acidity SOM (e.g., SOM I) is associated with low hydrophobicity of HS molecules, while low-acidity SOM (e.g., SOM II) exhibits higher hydrophobicity. This emphasizes the crucial role of HCB's hydrophobicity in its accumulation within SOM.

3.4 Interaction energies with functional groups

In order to characterize the interaction of the pollutants with SOM in more detail, interaction energies between the pollutant molecules and HS functional groups were calculated and plotted for the groups with the highest contributions to the interaction (see Fig. 5), considering the whole 1.5 μs of each of the production runs.

Fig. 5 Interaction energies (time-averaged) between SAA (a) or HCB (b) and beads representing different functional groups (x-axis). Values are further averaged across the 20 independent trajectories for every pollutant and SOM type. Proportions of the sulfonamide group (blue), aromatic ring (gray), and amino group (orange) are depicted for SAA. Green represents the HCB interaction energies. If available for a group, “T” or “S” stands for the Martini size of the bead, namely tiny or small, respectively. In the case of aromatic beads, a distinction is also made between those representing electron-rich regions of a ring, which are denoted by the letter “e”.

The negatively charged carboxylate groups exhibit the strongest interaction with SAA in each system. Furthermore, ketone, hydroxyl, phenol, acid anhydride, and aromatic functional groups play an important role. Overall, the preference for these functional groups can be explained by the hydrogen bonding ability of SAA and the hydrophilic character of the sulfonamide group. Conversely, the aromatic nature of SAA is reflected in the attraction to aromatic beads of the HSs, which is the most hydrophobic functional group in Fig. 5 that appears for SAA. The overall order for these groups is mostly determined by their abundance within each SOM system.

Another notable feature of SAA is its site-specific adsorption to the HSs, which is most apparent in the SOM III case when comparing the contribution of the sulfonamide group (blue) across different functional groups. The site-specific adsorption is characterized by selective interactions of the aromatic ring or sulfonamide group with specific functional groups or sites within SOM. The ratio ξ = E_int,sulf./E_int,arom. of interaction energies of the sulfonamide group compared to the aromatic part with the amino group (gray + orange) in SOM III is as follows for the five depicted functional groups, quinone, acid anhydride, aromatic, phenol, and carboxylate: 1.89, 1.94, 0.90, 1.88, and 2.36. This indicates, for example, a stronger attraction of the sulfonamide group to the negatively charged carboxylate groups (ξ = 2.36), and a higher attraction of the aromatic ring to the aromatic parts of the HS molecules (ξ = 0.90). This feature of SAA in the SOM interaction was reported by theoretical and experimental investigations^51,52 and is reproduced by the coarse-grained model, thus further validating the present parameters.

In comparison, for HCB, significant differences in the distribution of functional groups are visible. The hydrophobic aromatic and alkane groups show high contributions in every SOM type. The SOM I case shows HCB exhibiting exceptional behavior, whereby higher proportions of hydrophilic groups, e.g., carboxylate and ketone groups, are also recognizable, but where the overall values are also low. The higher interaction energies occurring in SOM II and III underline the importance of aromatic and alkane groups for the interaction of HCB with the HS structures. Overall, the interaction energies reflect the strong hydrophobicity of HCB, which alters the distribution of preferred functional groups compared to SAA.

3.5 Pollutant trapping

To explore the potential of HS structures to fixate pollutant molecules, an analysis of the trapping of the pollutants was conducted. This could be realized using quantities such as mean-first-passage times¹⁰³ or residence time distributions.¹⁰⁴ However, these require systematically defined positions or volumes within the SOM systems, which is impractical in light of the large total number of trajectories. Instead, the time intervals along the CGMD trajectories wherein the pollutant molecules were fixated were determined using the pollutants' spatial displacements and a hidden Markov model¹⁰⁵ (please refer to Section S7 (SI) for details about the calculations containing values of the displacements and identified trapping times). This type of trapping interval, where SAA remains fixed in the SOM matrix for a certain time, is illustrated in the trajectory snapshots in Fig. 6. The ratio of trapped to the total simulation time (trapping probability) for the different SOM types and pollutants is depicted in Fig. 7.

Fig. 6 Visualization of trapping and release of an SAA molecule in the SOM I system. The SAA trajectory is color-coded from red to blue, representing the movement from start to end point in the respective time interval. The HS molecules are depicted in brown, whereas the water beads are omitted for clarity. (a) SAA migrates from the bulk aqueous phase (red) into the HS matrix (blue). (b) The molecule then remains virtually stationary, illustrating its trapped state within the matrix. (c) Finally, SAA escapes, moving from the trapped location (red) back into the water phase (blue).

Fig. 7 Trapping probabilities of SAA (orange) and HCB (green) for the different SOM types.

In the case of SAA, the trapping probabilities are relatively high for SOM I, III, and IV. These systems exhibit a relatively low water diffusion (please compare to Fig. 2), which correlates with a higher probability of the molecules being trapped. In contrast, a low trapping probability is visible for SOM II, which is most likely connected to the high diffusion of SAA in this system, caused by the higher water diffusion. It is therefore assumed that the strong influence of the diffusion of water on SAA, indicated by the diffusion coefficients, is directly reflected in the probability of SAA being trapped in the SOM matrix, where water diffusion and trapping probability are anti-correlated. This is further supported by a strong negative Pearson correlation (r = −0.90) between water diffusion coefficients and SAA trapping probabilities. The results further indicate that the molecules do not follow a continuous diffusion, but rather go through sequential intervals of trapping and free diffusion. This behavior can only be observed using domain sampling methods such as CGMD simulations, which are facilitated by the coarse-grained representation of the molecules.

In contrast, the trapping probability of HCB correlates more with the coordination numbers of HS beads, i.e., with the probability of HCB being found in the HS phase (Fig. 3), emphasizing the influence of the direct interaction of HCB with SOM on diffusion. The Pearson correlation coefficient between coordination numbers and HCB trapping probabilities is r = 0.97, indicating a strong positive relationship.

However, the factors that influence the diffusion of the pollutant molecules in SOM appear to be complex. For example, when comparing SOM III and IV, an interesting behavior emerges: SOM III induces a much higher trapping probability than SOM IV, although the diffusion of HCB is similar in both cases. This suggests that the low diffusion of HCB can be caused by different mechanisms, which would be either the direct trapping of the molecule (SOM III) or the slowdown of the molecules via the interactions with the HS molecules without trapping (SOM IV). The difference between HCB and water diffusion for the respective cases, which is higher for the SOM III case, should probably also be taken into account here. Furthermore, it is assumed that the high trapping probabilities (above 0.5) are connected to the strong interaction between HCB and the HS aromatic compounds. This is supported by a complementary analysis of ring–ring contacts, which showed that trapping events coincide with an increased frequency of π–π stacking configurations (see Section S8 in the SI).

It is important to note that the SOM models used in this work represent highly hydrated, diluted systems. As such, the results concerning trapping and diffusion should primarily be interpreted in the context of such hydration levels. Since water facilitates molecular mobility, lower-moisture systems would be expected to exhibit reduced diffusion rates and increased trapping probabilities, particularly for polar pollutants such as SAA. Therefore, the observed interactions and mobility patterns may differ under drier or more heterogeneous soil conditions.

Using the presented methodology, molecular-level insights into how pollutants interact with SOM functional groups and how SOM composition influences pollutant mobility can be revealed, which has further implications for remediation strategies. Identifying preferential binding sites of pollutants can guide the design of amendments that either immobilize contaminants or enhance their degradation.¹⁰⁶ Future studies should aim to systematically characterize these binding preferences and define the conditions under which different pollutants are stabilized or mobilized, thereby allowing predictions of how amendments may alter these processes.

3.6 Implications of the Martini representation

The coarse-grained simplification that underpins the Martini models accelerates simulations, but inevitably sacrifices chemical resolution and can lead to artifacts in the simulation of the molecules. Many of the shortcomings of the Martini 2 force field have been alleviated in Martini 3. For instance, in Martini 2, the models' interactions caused an overstabilization of biomolecular interactions for soluble¹⁰⁷ and membrane proteins,^108,109 and also for carbohydrates.^110–112 The consequence was a too stable self-association in solution and an exaggerated tendency of membrane proteins to aggregate. Martini 3 alleviated this “stickiness” through a global reparametrization of bead interactions,^69,112 yet overcompaction for intrinsically disordered or multidomain proteins was noted.¹¹³

A second artifact arose when very short bond distances or mixed bead sizes were employed: the resulting high density of interaction sites generated artificial desolvation barriers and overstabilized dimers.¹¹¹ Contrary to previous versions, Martini 3 introduces a “size–shaped” mapping together with additional new bead sizes, aimed to restore the correct molecular volume and to improve the oil–water partitioning of small molecules.^69,112 This type of mapping was applied to all molecules investigated in this study. However, although Martini 3 remedies many of the deficiencies identified in the earlier version, one should continue to critically examine whether the remaining approximations are acceptable for the physicochemical properties under investigation.

Coarse-graining eliminates certain configurational degrees of freedom and thus reduces the intrinsic molecular entropy. In the Martini framework, this entropy loss is empirically compensated by enhancing non-bonded interaction strengths so that resulting free energies match reference experimental or atomistic data. While this ensures accurate reproduction of free energy differences, it distorts the enthalpy–entropy balance, making coarse-grained interaction energies model-specific and not directly interpretable as physical enthalpies.^114–118 Therefore, the absolute values in Fig. 5 should be interpreted with caution and used in combination with non-energetic measures such as coordination numbers (Fig. S6 in the SI). A more rigorous treatment would involve calculating sorption free energies, similar to the approach of Petrov et al.,¹¹⁹ who investigated small organic compounds in LHA systems using atomistic simulations. This approach could also be applied to drier SOM systems, which exhibit slower solute diffusion as a result of reduced water activity.^120,121 While this was beyond the computational scope of the present study, such approaches could be explored in future work to quantify specific pollutant–SOM interactions in a statistically rigorous framework.

3.7 Comparison with literature

Experimental and computational studies have shown that HCB and SAA become trapped in SOM, including reports of irreversible sorption into SOM voids.^24,122,123 While the coarse-grained HS model used here cannot represent the full structural and chemical complexity of natural SOM, it reproduces key qualitative patterns of pollutant immobilization. This consistency provides a mesoscale confirmation of known behaviors and highlights the utility of coarse-grained approaches in bridging molecular insight with environmental timescales.

Ahmed et al.⁵¹ combined sorption experiments with all-atom MD simulations for simplified SOM models to show that HCB binds more strongly to SOM than SAA, attributed to HCB's hydrophobicity and aromaticity, which agrees with the calculated coordination numbers (Fig. 3) and interaction energies (Fig. 5). In a related density functional theory (DFT) study, Ahmed et al.⁵² demonstrated that water destabilizes SAA–SOM complexes and that SAA shows a strong affinity for ionic sites, consistent with the observations of the present study.

Furthermore, previous work indicates that SAA sorption onto SOM is strongly composition-dependent.^15,51 This behavior was not observed in the present study, probably because the high water content favored SAA partitioning into the aqueous phase. Future work on varying water content or pH with coarse-grained simulations may further clarify these differences. Regarding HCB, DFT calculations on HCB–SOM interactions revealed that it strongly binds to functional groups such as alkylated aromatics, phenol, and lignin.⁵⁰ In the present study, similar contributions from aromatic, alkane, and phenol groups were evident in the binding of HCB.

4 Conclusions

While pollutants continue to pose a persistent threat to terrestrial ecosystems, their molecular–scale interactions with SOM remain only partially understood, a gap this study helps to address. Using CGMD simulations, we investigated the mobility of SAA and HCB in SOM systems of varied composition.

Diffusion analysis revealed contrasting behavior: SAA followed closely the diffusion of water beads, whereas HCB displayed reduced diffusion except in systems with high carboxylate content. Based on the calculation of coordination numbers, the partitioning behavior of the pollutants indicated that SAA's distribution between water and HS phases was nearly composition-independent, while HCB accumulated in HS matrices containing low carboxyl and low carbonyl content, or high aromatic content.

Interaction-energy calculations revealed a preference of SAA to carboxylate groups and other hydrophilic functionalities. Aromatic moieties also contributed, displaying a site-specific interaction of SAA to HS. For HCB, interactions were dominated by contacts with aromatic moieties, alkane chains, as well as methoxy and acid-anhydride functional groups.

Analysis using a hidden Markov model resolved distinct trapping states for both pollutants, encompassing fixation within the SOM matrix and periods of free diffusion. The extent of SAA trapping was anticorrelated with water diffusion, underscoring the influence of solvent mobility on its behavior. In contrast, HCB trapping correlated with encounters with HS molecules, indicating direct interactions with the HS matrix as an immobilizing factor.

The coarse-grained approach employed here enables the investigation of dynamic processes over extended time and length scales, thereby complementing previous studies and integrating into the current research framework of pollutants in soil. The method offers broad potential for future applications, including studies of varied SOM compositions, the integration of other pre-parametrized pollutants from the Martini 3 framework, and the investigation of large-scale SOM structural organization.

Conflicts of interest

There are no conflicts to declare.

Data availability

All topology files required to reproduce the simulations from this study are available in Zenodo at DOI: https://doi.org/10.5281/zenodo.14749817

Data supporting this article have been included as part of the Supplementary Information (SI). See DOI: https://doi.org/10.1039/d5va00237k. The SI includes details on the pollutant parametrization and HS samples, additional results on coordination numbers and interaction energies, and analyses of π–π stacking configurations between pollutants and SOM.

References

1. R. R. Weil and N. C. Brady, The Nature and Properties of Soils, Pearson, Upper Saddle River, NJ, 15th edn, 2016.
2. B. D. Hudson, J. Soil Water Conserv., 1994, 49, 189–194.
3. W. Rawls, Y. A. Pachepsky, J. Ritchie, T. Sobecki and H. Bloodworth, Geoderma, 2003, 116, 61–76.
4. Food and Agriculture Organization of the United Nations, Plant Nutrition for Food Security, a Guide for Integrated Nutrient Management, Food & Agriculture Organization of the United Nations (FAO), Rome, Italy, 2006.
5. P. K. R. Nair, B. M. Kumar and V. D. Nair, in Soil Organic Matter (SOM) and Nutrient Cycling, Springer International Publishing, 2021, pp. 383–411.
6. B. Xu, G. Yang, A. Lehmann, S. Riedel and M. C. Rillig, Soil Ecol. Lett., 2022, 5, 108–117.
7. M. N. Hanif, N. Aijaz, K. Azam, M. Akhtar, W. A. Laftah, M. Babur, N. K. Abbood and I. B. Benitez, Int. J. Environ. Sci. Technol., 2024, 21, 10277–10318.
8. L. Ritter, K. R. Solomon, J. Forget, M. Stemeroff and C. O'Leary, Persistent Organic Pollutants: an Assessment Report on DDT, Aldrin, Dieldrin, Endrin, Chlordane, Heptachlor, Hexachlorobenzene, Mirex, Toxaphene, Polychlorinated Biphenyls, Dioxins and Furans, International Programme on Chemical Safety (IPCS), United Nations Environment Programme Technical Report PCS/95.39, 1995.
9. K. C. Jones and P. De Voogt, Environ. Pollut., 1999, 100, 209–221.
10. T. Luarte, V. A. Gómez-Aburto, I. Poblete-Castro, E. Castro-Nallar, N. Huneeus, M. Molina-Montenegro, C. Egas, G. Azcune, A. Pérez-Parada, R. Lohmann, P. Bohlin-Nizzetto, J. Dachs, S. Bengtson-Nash, G. Chiang, K. Pozo and C. J. Galbán-Malagón, Atmos. Chem. Phys., 2023, 23, 8103–8118.
11. J. L. Wilkinson, A. B. A. Boxall, D. W. Kolpin, K. M. Y. Leung, R. W. S. Lai, C. Galbán-Malagón, A. D. Adell, J. Mondon, M. Metian, R. A. Marchant, A. Bouzas-Monroy, A. Cuni-Sanchez, A. Coors, P. Carriquiriborde, M. Rojo, C. Gordon, M. Cara, M. Moermond, T. Luarte, V. Petrosyan, Y. Perikhanyan, C. S. Mahon, C. J. McGurk, T. Hofmann, T. Kormoker, V. Iniguez, J. Guzman-Otazo, J. L. Tavares, F. G. D. Figueiredo, M. T. P. Razzolini, V. Dougnon, G. Gbaguidi, O. Traoré, J. M. Blais, L. E. Kimpe, M. Wong, D. Wong, R. Ntchantcho, J. Pizarro, G.-G. Ying, C.-E. Chen, M. Páez, J. Martínez-Lara, J.-P. Otamonga, J. Poté, S. A. Ifo, P. Wilson, S. Echeverría-Sáenz, N. Udikovic-Kolic, M. Milakovic, D. Fatta-Kassinos, L. Ioannou-Ttofa, V. Belušová, J. Vymazal, M. Cárdenas-Bustamante, B. A. Kassa, J. Garric, A. Chaumot, P. Gibba, I. Kunchulia, S. Seidensticker, G. Lyberatos, H. P. Halldórsson, M. Melling, T. Shashidhar, M. Lamba, A. Nastiti, A. Supriatin, N. Pourang, A. Abedini, O. Abdullah, S. S. Gharbia, F. Pilla, B. Chefetz, T. Topaz, K. M. Yao, B. Aubakirova, R. Beisenova, L. Olaka, J. K. Mulu, P. Chatanga, V. Ntuli, N. T. Blama, S. Sherif, A. Z. Aris, L. J. Looi, M. Niang, S. T. Traore, R. Oldenkamp, O. Ogunbanwo, M. Ashfaq, M. Iqbal, Z. Abdeen, A. O'Dea, J. M. Morales-Saldaña, M. Custodio, H. de la Cruz, I. Navarrete, F. Carvalho, A. B. Gogra, B. M. Koroma, V. Cerkvenik-Flajs, M. Gombač, M. Thwala, K. Choi, H. Kang, J. L. C. Ladu, A. Rico, P. Amerasinghe, A. Sobek, G. Horlitz, A. K. Zenker, A. C. King, J.-J. Jiang, R. Kariuki, M. Tumbo, U. Tezel, T. T. Onay, J. B. Lejju, Y. Vystavna, Y. Vergeles, H. Heinzen, A. Pérez-Parada, D. B. Sims, M. Figy, D. Good and C. Teta, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2113947119.
12. B. Gworek, M. Kijeńska, J. Wrzosek and M. Graniewska, Water, Air, Soil Pollut., 2021, 232, 145.
13. W. Zheng and M. Guo, Curr. Pollut. Rep., 2021, 7, 510–523.
14. L. Wang, R. Govind and R. A. Dobbs, Environ. Sci. Technol., 1993, 27, 152–158.
15. S. Thiele-Bruhn, T. Seibicke, H. Schulten and P. Leinweber, J. Environ. Qual., 2004, 33, 1331–1342.
16. J. Gao and J. A. Pedersen, Environ. Sci. Technol., 2005, 39, 9509–9516.
17. A. A. MacKay and B. Canterbury, J. Environ. Qual., 2005, 34, 1964–1971.
18. Q. Wang, M. Guo and S. R. Yates, J. Agric. Food Chem., 2005, 54, 157–163.
19. T. L. ter Laak, W. A. Gebbink and J. Tolls, Environ. Toxicol. Chem., 2006, 25, 904–911.
20. S. T. Kurwadkar, C. D. Adams, M. T. Meyer and D. W. Kolpin, J. Agric. Food Chem., 2007, 55, 1370–1376.
21. W. Lertpaitoonpan, S. K. Ong and T. B. Moorman, Chemosphere, 2009, 76, 558–564.
22. R. A. Figueroa-Diva, D. Vasudevan and A. A. MacKay, Chemosphere, 2010, 79, 786–793.
23. A. Białk-Bielińska, J. Maszkowska, W. Mrozik, A. Bielawska, M. Kołodziejska, R. Palavinskas, P. Stepnowski and J. Kumirska, Chemosphere, 2012, 86, 1059–1065.
24. J. Schwarz, S. Thiele-Bruhn, K.-U. Eckhardt and H.-R. Schulten, ISRN Soil Sci., 2012, 2012, 1–10.
25. T. Müller, I. Rosendahl, A. Focks, J. Siemens, J. Klasmeier and M. Matthies, Environ. Pollut., 2013, 172, 180–185.
26. W. J. Farmer, M. S. Yang, J. Letey and W. F. Spencer, Soil Sci. Soc. Am. J., 1980, 44, 676–680.
27. G. Deane, Z. Chroneer and W. Lick, J. Environ. Eng., 1999, 125, 689–696.
28. S. N. Meijer, W. A. Ockenden, A. Sweetman, K. Breivik, J. O. Grimalt and K. C. Jones, Environ. Sci. Technol., 2003, 37, 667–672.
29. Y.-X. Chen, H.-L. Chen, Y.-T. Xu and M.-W. Shen, Environ. Int., 2004, 30, 31–37.
30. J. D. Kubicki, Environ. Sci. Technol., 2006, 40, 2298–2303.
31. D. Tunega, M. H. Gerzabek, G. Haberhauer, K. U. Totsche and H. Lischka, J. Colloid Interface Sci., 2009, 330, 244–249.
32. A. A. Ahmed, O. Kühn and P. Leinweber, Sci. Total Environ., 2012, 441, 151–158.
33. H. Zhang, Y. Wang, C. Sun, M. Yu, Y. Gao, T. Wang, J. Liu and G. Jiang, Environ. Sci. Technol., 2014, 48, 1525–1531.
34. L. Böhm, P. Grančič, E. Scholtzová, B. J. Heyde, R.-A. Düring, J. Siemens, M. H. Gerzabek and D. Tunega, Environ. Sci. Pollut. Res., 2022, 30, 36824–36837.
35. J. Zhang, Z. Zhou, L. Zeng, C. Wang, R. Han, X. Ren, W. Wang, M. Xiang, S. Chen and H. Li, Sci. Total Environ., 2024, 947, 174657.
36. S. Zeng, Z. Dai, B. Ma, R. A. Dahlgren and J. Xu, Earth Critical Zone, 2024, 1, 100002.
37. International Humic Substances Society, https://humic-substances.org.
38. J. Gerke, Agronomy, 2018, 8, 76.
39. M. Kleber and J. Lehmann, J. Environ. Qual., 2019, 48, 207–216.
40. M. H. Hayes and R. S. Swift, in Vindication of Humic Substances as a Key Component of Organic Matter in Soil and Water, Elsevier, 2020, pp. 1–37.
41. T. C. de Aguiar, D. F. de Oliveira Torchia, T. A. van Tol de Castro, O. C. H. Tavares, S. de Abreu Lopes, L. de Souza da Silva, R. N. Castro, R. L. L. Berbara, M. G. Pereira and A. C. García, Sci. Total Environ., 2022, 833, 155133.
42. D. C. Olk, P. R. Bloom, E. M. Perdue, D. M. McKnight, Y. Chen, A. Farenhorst, N. Senesi, Y. Chin, P. Schmitt-Kopplin, N. Hertkorn and M. Harir, J. Environ. Qual., 2019, 48, 217–232.
43. D. C. Olk, P. R. Bloom, M. De Nobili, Y. Chen, D. M. McKnight, M. J. M. Wells and J. Weber, J. Environ. Qual., 2019, 48, 1633–1643.
44. H. R. Schulten and M. Schnitzer, Naturwissenschaften, 1995, 82, 487–498.
45. H.-R. Schulten and M. Schnitzer, Soil Sci., 1997, 162, 115–130.
46. H.-R. Schulten and P. Leinweber, Biol. Fertil. Soils, 2000, 30, 399–432.
47. H.-R. Schulten, Soil Mineral-Organic Matter-Microorganism Interactions and Ecosystem Health, Dynamics, Mobility and Transformation of Pollutants and Nutrients, Elsevier, 2002, pp. 351–381.
48. A. J. A. Aquino, D. Tunega, G. Haberhauer, M. H. Gerzabek and H. Lischka, Eur. J. Soil Sci., 2007, 58, 889–899.
49. A. J. A. Aquino, D. Tunega, G. E. Schaumann, G. Haberhauer, M. H. Gerzabek and H. Lischka, J. Phys. Chem. C, 2009, 113, 16468–16475.
50. A. A. Ahmed, O. Kühn, S. G. Aziz, R. H. Hilal and P. Leinweber, Sci. Total Environ., 2014, 476–477, 98–106.
51. A. A. Ahmed, S. Thiele-Bruhn, S. G. Aziz, R. H. Hilal, S. A. Elroby, A. O. Al-Youbi, P. Leinweber and O. Kühn, Sci. Total Environ., 2015, 508, 276–287.
52. A. A. Ahmed, S. Thiele-Bruhn, P. Leinweber and O. Kühn, Sci. Total Environ., 2016, 559, 347–355.
53. P. Gros, A. Ahmed, O. Kühn and P. Leinweber, Sci. Total Environ., 2017, 586, 527–535.
54. H. Feng, H. Zhang, H. Cao, Y. Sun, A. Zhang and J. Fu, Environ. Sci. Technol., 2018, 52, 14228–14234.
55. Q. Xue, Z. Jiao, X. Liu, W. Pan, J. Fu and A. Zhang, Environ. Sci. Technol., 2023, 58, 1531–1540.
56. A. Sündermann, R. Solc, D. Tunega, G. Haberhauer, M. H. Gerzabek and C. Oostenbrink, J. Mol. Graphics Modell., 2015, 62, 253–261.
57. Y. Escalona, D. Petrov and C. Oostenbrink, J. Mol. Graphics Modell., 2021, 103, 107817.
58. M. H. Gerzabek, A. J. A. Aquino, Y. I. E. Balboa, E. Galicia-Andrés, P. Grančič, C. Oostenbrink, D. Petrov and D. Tunega, J. Plant Nutr. Soil Sci., 2022, 185, 44–59.
59. V. A. Harmandaris, N. P. Adhikari, N. F. A. van der Vegt and K. Kremer, Macromolecules, 2006, 39, 6708–6719.
60. W. G. Noid, J. Chem. Phys., 2013, 139, 090901.
61. S. Kmiecik, D. Gront, M. Kolinski, L. Wieteska, A. E. Dawid and A. Kolinski, Chem. Rev., 2016, 116, 7898–7936.
62. J. Jin, A. J. Pak, A. E. P. Durumeric, T. D. Loose and G. A. Voth, J. Chem. Theory Comput., 2022, 18, 5759–5791.
63. M. Orsi and J. W. Essex, PLoS One, 2011, 6, e28637.
64. L. Darré, M. R. Machado, A. F. Brandner, H. C. González, S. Ferreira and S. Pantano, J. Chem. Theory Comput., 2015, 11, 723–739.
65. M. R. Machado, E. E. Barrera, F. Klein, M. Sóñora, S. Silva and S. Pantano, J. Chem. Theory Comput., 2019, 15, 2719–2733.
66. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. de Vries, J. Phys. Chem. B, 2007, 111, 7812–7824.
67. L. F. Dettmann, O. Kühn and A. A. Ahmed, J. Chem. Theory Comput., 2024, 20, 5291–5305.
68. L. F. Dettmann, O. Kühn and A. A. Ahmed, J. Chem. Theory Comput., 2024, 20, 10684–10696.
69. P. C. T. Souza, R. Alessandri, J. Barnoud, S. Thallmair, I. Faustino, F. Grünewald, I. Patmanidis, H. Abdizadeh, B. M. H. Bruininks, T. A. Wassenaar, P. C. Kroon, J. Melcr, V. Nieto, V. Corradi, H. M. Khan, J. Domański, M. Javanainen, H. Martinez-Seara, N. Reuter, R. B. Best, I. Vattulainen, L. Monticelli, X. Periole, D. P. Tieleman, A. H. de Vries and S. J. Marrink, Nat. Methods, 2021, 18, 382–388.
70. R. Alessandri, J. Barnoud, A. S. Gertsen, I. Patmanidis, A. H. de Vries, P. C. T. Souza and S. J. Marrink, Adv. Theory Simul., 2021, 5, 2100391.
71. F. Grünewald, M. H. Punt, E. E. Jefferys, P. A. Vainikka, M. König, V. Virtanen, T. A. Meyer, W. Pezeshkian, A. J. Gormley, M. Karonen, M. S. P. Sansom, P. C. T. Souza and S. J. Marrink, J. Chem. Theory Comput., 2022, 18, 7555–7569.
72. C. Empereur-mot, K. B. Pedersen, R. Capelli, M. Crippa, C. Caruso, M. Perrone, P. C. T. Souza, S. J. Marrink and G. M. Pavan, J. Chem. Inf. Model., 2023, 63, 3827–3838.
73. L. Borges-Araújo, G. P. Pereira, M. Valério and P. C. Souza, Biochim. Biophys. Acta, Proteins Proteomics, 2024, 1872, 141014.
74. A. F. Brandner, I. P. S. Smith, S. J. Marrink, P. C. T. Souza and S. Khalid, J. Chem. Inf. Model., 2025, 65, 1537–1548.
75. F. Grünewald, R. Alessandri, P. C. Kroon, L. Monticelli, P. C. T. Souza and S. J. Marrink, Nat. Commun., 2022, 13, 68.
76. R. Shrestha, R. Alessandri, M. Vögele, C. Hilpert, P. Souza, S. J. Marrink and L. Monticelli, J. Chem. Theory Comput., 2025, 21, 9035–9053.
77. C. Empereur-Mot, L. Pesce, G. Doni, D. Bochicchio, R. Capelli, C. Perego and G. M. Pavan, ACS Omega, 2020, 5, 32823–32843.
78. C. H. Bennett, J. Comput. Phys., 1976, 22, 245–268.
79. Y. Escalona, D. Petrov, E. Galicia-Andrés and C. Oostenbrink, Agronomy, 2023, 13, 1044.
80. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
81. C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke and T. E. Oliphant, Nature, 2020, 585, 357–362.
82. J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95.
83. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.
84. R. Gowers, M. Linke, J. Barnoud, T. Reddy, M. Melo, S. Seyler, J. Domański, D. Dotson, S. Buchoux, I. Kenney and O. Beckstein, Proceedings of the Python in Science Conference, 2016.
85. F. Pérez and B. E. Granger, Comput. Sci. Eng., 2007, 9, 21–29.
86. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
87. R. Alessandri, S. Thallmair, C. G. Herrero, R. Mera-Adasme, S. J. Marrink and P. C. T. Souza, A Practical Guide to Recent Advances in Multiscale Modeling and Simulation of Biomolecules, AIP Publishing LLC Melville, New York, 2023, pp. 1–34.
88. C. Hansch, A. Leo and D. Hoekman, Exploring QSAR: Hydrophobic, Electronic, and Steric Constants, American Chemical Society, Washington, D.C., 1995, vol. 2, pp. 15–22.
89. G. Jung, M. Hanke and F. Schmid, J. Chem. Theory Comput., 2017, 13, 2481–2488.
90. V. Klippenstein, M. Tripathy, G. Jung, F. Schmid and N. F. A. van der Vegt, J. Phys. Chem. B, 2021, 125, 4931–4954.
91. T. Schilling, Phys. Rep., 2022, 972, 1–45.
92. M. S. Shell, J. Chem. Phys., 2008, 129, 144108.
93. M. S. Shell, J. Chem. Phys., 2012, 137, 084503.
94. Y. Rosenfeld, Phys. Rev. A, 1977, 15, 2545–2549.
95. J. C. Dyre, J. Chem. Phys., 2018, 149, 210901.
96. G. G. Rondina, M. C. Böhm and F. Müller-Plathe, J. Chem. Theory Comput., 2020, 16, 1431–1447.
97. J. Jin, K. S. Schweizer and G. A. Voth, J. Chem. Phys., 2023, 158, 034103.
98. P. K. Depa and J. K. Maranas, J. Chem. Phys., 2005, 123, 094901.
99. D. Fritz, K. Koschke, V. A. Harmandaris, N. F. A. van der Vegt and K. Kremer, Phys. Chem. Chem. Phys., 2011, 13, 10412.
100. M. K. Meinel and F. Müller-Plathe, J. Chem. Theory Comput., 2020, 16, 1411–1419.
101. M. K. Meinel and F. Müller-Plathe, J. Phys. Chem. B, 2022, 126, 3737–3747.
102. P. C. T. Souza, S. Thallmair, P. Conflitti, C. Ramírez-Palacios, R. Alessandri, S. Raniolo, V. Limongelli and S. J. Marrink, Nat. Commun., 2020, 11, 3714.
103. S. Redner, A Guide to First-Passage Processes, Cambridge University Press, 2001.
104. P. Danckwerts, Chem. Eng. Sci., 1953, 2, 1–13.
105. L. Rabiner, Proc. IEEE, 1989, 77, 257–286.
106. F. J. Larney and D. A. Angers, Can. J. Soil Sci., 2012, 92, 19–38.
107. A. C. Stark, C. T. Andrews and A. H. Elcock, J. Chem. Theory Comput., 2013, 9, 4176–4185.
108. M. Javanainen, H. Martinez-Seara and I. Vattulainen, PLoS One, 2017, 12, e0187936.
109. A. Majumder and J. E. Straub, J. Chem. Theory Comput., 2021, 17, 2513–2521.
110. P. S. Schmalhorst, F. Deluweit, R. Scherrers, C.-P. Heisenberg and M. Sikora, J. Chem. Theory Comput., 2017, 13, 5039–5053.
111. R. Alessandri, P. C. T. Souza, S. Thallmair, M. N. Melo, A. H. de Vries and S. J. Marrink, J. Chem. Theory Comput., 2019, 15, 5448–5460.
112. S. J. Marrink, L. Monticelli, M. N. Melo, R. Alessandri, D. P. Tieleman and P. C. T. Souza, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2022, 13, e1620.
113. F. E. Thomasen, T. Skaalum, A. Kumar, S. Srinivasan, S. Vanni and K. Lindorff-Larsen, Nat. Commun., 2024, 15, 6645.
114. N. J. H. Dunn, T. T. Foley and W. G. Noid, Acc. Chem. Res., 2016, 49, 2832–2840.
115. J. Jin, A. J. Pak and G. A. Voth, J. Phys. Chem. Lett., 2019, 10, 4549–4557.
116. K. M. Kidder, R. J. Szukalo and W. G. Noid, Eur. Phys. J. B, 2021, 94, 153.
117. Z. Jarin, J. Newhouse and G. A. Voth, J. Chem. Theory Comput., 2021, 17, 1170–1180.
118. T. D. Loose, P. G. Sahrmann, T. S. Qu and G. A. Voth, J. Chem. Theory Comput., 2024, 20, 9190–9208.
119. D. Petrov, D. Tunega, M. H. Gerzabek and C. Oostenbrink, Eur. J. Soil Sci., 2020, 71, 831–844.
120. V. P. Kalinitchenko, A. V. Swidsinski, A. P. Glinushkin, V. P. Meshalkin, S. V. Gudkov, T. M. Minkina, V. V. Chernenko, V. D. Rajput, S. S. Mandzhieva, S. N. Sushkova, A. A. Okolelova and A. A. Shestakova, Environ. Geochem. Health, 2023, 45, 8967–8987.
121. V. P. Kalinitchenko, A. V. Swidsinski, A. P. Glinushkin, M. R. Overcash, V. P. Meshalkin, T. M. Minkina, V. D. Rajput, V. V. Chernenko, S. S. Mandzhieva, S. N. Sushkova, A. E. Rykhlik, G. S. Larin, S. V. Gudkov, M. G. Baryshev, M. A. Sevostyanov and D. A. Makarenkov, Land Degrad. Dev., 2024, 35, 2367–2381.
122. A. Kan, W. Chen and M. Tomson, Environ. Pollut., 2000, 108, 81–89.
123. J.-j. Zhang, B. Wen, X.-q. Shan, S. Zhang and S. U. Khan, Environ. Pollut., 2007, 150, 234–242.

---