Unveiling the multiphase fate of 2,4-dinitrophenol on aerosols: interfacial hydration modulates competing oxidation pathways and predicted toxicity indicators

This study elucidates the mechanisms of selected transformation pathways of 2,4-dinitrophenol (2,4-DNP) using an integrated computational framework. In the gas phase, the calculated initial reactions with hydroxyl radicals (•OH) and ozone (O₃) are kinetically more favorable than the direct reaction with nitrogen dioxide radical (•NO₂), which is kinetically hindered. In the aqueous phase, solvation substantially modifies the calculated reaction kinetics, and the cumulative electron-withdrawing effect of the two nitro substituents lowers the O₃ reaction rate coefficient by approximately three orders of magnitude relative to mononitrocatechol. For the key radical intermediate generated by •OH addition, hydrogen atom abstraction (HAA) by O₂ is kinetically favored over the competing •NO₂ pathway based on the calculated barriers and an illustrative comparison of effective rates under representative atmospheric concentrations. At the gas-aqueous interface, molecular dynamics (MD) simulations show a non-monotonic dependence of 2,4-DNP adsorption on hydration in the idealized C60/water system, with adsorption strengthening at moderate hydration and weakening at high hydration. Finally, QSAR calculations predict that certain ozonolysis products exhibit higher predicted Ames mutagenicity and developmental toxicity indicators than the parent compound, indicating product-specific changes in predicted toxicity profiles during oxidation. These screening-level predictions require experimental validation. These findings provide molecular-level mechanistic insights into selected transformation pathways and interfacial interactions involving 2,4-DNP and motivate further experimental investigation of its transformation products and their toxicological properties.

1 Introduction

The complex physicochemical processes of atmospheric aerosols are central to understanding their roles in climate and environmental quality (Perumpully and Gautam, 2025; McNeill, 2017). Among them, secondary organic aerosols (SOAs) exert substantial influences on the global radiation balance and public health (Li et al., 2022a; Kroll and Seinfeld, 2008). Their formation originates from the oxidation, nitration, and subsequent reactions of volatile organic compounds (VOCs), which generate semi-volatile and low-volatility products (Fan et al., 2022; Li et al., 2021; Mahilang et al., 2021). Within the diverse range of precursors, nitrophenolic compounds (NPs) are important toxic organic pollutants, primarily emitted from incomplete combustion of fossil fuels and biomass, and also formed via photochemical oxidation of precursors such as polycyclic aromatic hydrocarbons (Li et al., 2022b, 2020; Salvador et al., 2021). Among them, 2,4-dinitrophenol (2,4-DNP) is listed as a priority pollutant by the U.S. Environmental Protection Agency due to its acute and chronic toxicity, potential genotoxicity, and bioaccumulation (Shukla et al., 2009; Yan et al., 2021). Although ambient concentrations of 2,4-DNP are typically in the ng m⁻³ range, its low volatility and favorable aqueous partitioning promote its presence in condensed phases such as aerosol water and cloud/fog droplets (Lüttke and Levsen, 1997). Consequently, it can be deposited onto terrestrial and aquatic ecosystems via dry and wet deposition, posing sustained threats to the environment and human health. A molecular-level understanding of its multiphase chemical transformation mechanisms is therefore essential for improving mechanistic descriptions of its environmental transformation and for informing future fate and risk assessments.

Over the past two decades, extensive research has addressed the atmospheric chemistry of nitrophenols, yet the progress has been uneven across different phases. In the gas phase, kinetic measurements have established rate constants for the •OH initiated oxidation of nitrophenols (Bejan et al., 2007), while photolysis and •NO₃ mediated pathways are recognized as important loss processes, particularly at night (Bejan et al., 2007; Ng et al., 2017). In the aqueous phase, reactions with •OH radicals have been systematically quantified; for 2,4-DNP, the second order rate constants are (1.76 ± 0.05) ×10⁹ M⁻¹ s⁻¹ for the neutral form and (2.33 ± 0.11) ×10⁹ M⁻¹ s⁻¹ for the phenolate (Albinet et al., 2010), with similar kinetic data available for nitrocatechol, nitroguaiacol, and other dinitrophenols (Hems and Abbatt, 2018). However, these studies have mainly focused on homogeneous reactions, leaving the gas-liquid and gas-solid interfaces – key regions where reactivity often differs from the bulk phases – far less understood. Recent work has shown that oxidation of phenolic compounds at interfaces yields products distinct from those formed in the bulk (Pillar-Little et al., 2014; Pillar-Little and Guzman, 2017; Rana and Guzman, 2020; Rana and Guzman, 2022b, a), underscoring the importance of elucidating the molecular level interactions that govern interfacial adsorption and reactivity. Nevertheless, the behavior of nitrophenols on carbonaceous aerosol surfaces remains largely unexplored.

This knowledge gap stems primarily from the extreme complexity of real carbonaceous particles, such as black carbon and soot. These particles, products of incomplete combustion, exhibit physicochemical properties that depend strongly on pyrolysis temperature, formation conditions, and subsequent atmospheric aging (Andreae and Gelencsér, 2006). Although their polyaromatic (graphite-like) surfaces can form π-π electron donor-acceptor complexes with nitroaromatic compounds, adsorption strength is significantly influenced by steric effects due to micropore filling and by fouling with natural organic matter that blocks pores and competes for sites (Cornelissen et al., 2005; Poster et al., 2007). Moreover, chemical heterogeneity (oxygen containing functional groups, defect sites, metal impurities) and morphological complexity (aggregates, pore structures) make it extremely difficult to decouple adsorption from partitioning contributions or to identify specific non-covalent interaction types at the molecular level. A well-defined model system is therefore necessary to isolate specific interactions from the confounding effects of surface heterogeneity. Fullerene (C60), with its well-defined, stable carbon cage structure and abundant π-electron cloud, provides a well-defined model system for studying non-covalent interactions at organic compound-particle interfaces (Athar and Patnaik, 2022; Liu and Lu, 2024). By using C60 as a surrogate for carbonaceous aerosol surfaces, one can effectively bypass the complexities of real particles and systematically examine the adsorption configurations, energetics, and their regulatory factors. While phenolic compounds may undergo further multiphase aging to form SOA and brown carbon (Pillar-Little and Guzman, 2018; Guzman et al., 2022; Al-Abadleh et al., 2022), the present study focuses on their initial molecular-scale oxidation and interfacial adsorption.

In the complex and dynamic real atmosphere, the transformation of nitrophenols is a dynamic competition among multiple reactive species. The key oxidants include •OH (dominant during daytime), O₃ (with its high background concentration), and •NO₂, which can accumulate at night in polluted regions and initiate non-photochemical nitration pathways (Delaria and Cohen, 2023; Pillar-Little et al., 2014; Yang et al., 2023b). Photolysis of 2,4-DNP has been shown to be a minor loss route compared to •OH oxidation (Albinet et al., 2010). •NO₃ radicals represent an additional important component of nighttime phenolic chemistry that is not explicitly included in the present computational framework (Rana and Guzman, 2022b; Ng et al., 2017). Recent experimental studies have shown that heterogeneous NO₃ oxidation of catechol can generate nitrocatechols, oligomeric products, and light-absorbing brown-carbon chromophores (Rana and Guzman, 2026b), while related studies of phenolic aldehydes further underscore the complexity of heterogeneous nighttime oxidation pathways (Rana and Guzman, 2026a). These studies underscore the importance of NO₃-driven chemistry in the nighttime transformation of phenolic compounds. The present study, however, is restricted to the reaction pathways initiated by O₃, •OH, and •NO₂; therefore, it does not provide a complete description of the nighttime oxidation chemistry of 2,4-DNP.

Such multi-pathway transformation processes generate a series of intermediates and final products with diverse chemical properties, which can fundamentally alter the ecological risk profile of nitroaromatic pollutants (Tiwari et al., 2019; Tam et al., 2021). Traditional experimental approaches often struggle to comprehensively identify all reactive intermediates and transformation products, and computational toxicology offers a high-throughput, low-cost alternative for assessing the toxicity of atmospheric transformation products (Tam et al., 2021; vom Eyser et al., 2013). For instance, for a set of 104 secondary organic aerosol compounds, reliable in silico predictions were obtained for 82 % of the species, whereas experimental data were available for only 13 % (Decesari et al., 2018). Importantly, atmospheric oxidation may produce compounds with toxicity exceeding that of their parents (vom Eyser et al., 2013), highlighting the need to integrate toxicity prediction into transformation pathway assessment. The application of quantitative structure activity relationship (QSAR) to nitrophenolic transformation products, however, remains limited, and the inherent uncertainties require careful interpretation.

Previous studies on phenolic oxidation have predominantly focused on systems with electron-donating substituents that enhance reactivity. Pillar-Little et al. (2015, 2014) and Pillar-Little and Guzman (2017) demonstrated that catechols and their electron-donating (–CH₃, –OCH₃) derivatives readily form carboxylic acids, quinones, and polyphenols at interfaces; Rana and Guzman (2020, 2022b, a) systematically investigated the interfacial oxidation of phenolic aldehydes; Witt et al. (2025) found that ozonolysis of non-nitrated phenolic aldehydes yields smaller, less toxic molecules such as oxalic and glyoxylic acids; and Rana et al. (2023) showed that catechol can be converted to 4-nitrocatechol (4NC) in aqueous microdroplets via multiphase •NO₂ $/$ O₃ chemistry. These studies have laid a foundation for understanding substituent effects, but 2,4-DNP – with two strongly electron-withdrawing nitro groups – may exhibit fundamentally different reactivity and mechanisms.

Despite these advances, a systematic, molecular-level understanding of the multiphase fate of 2,4-DNP remains lacking. Specifically, three core scientific questions need to be addressed: (i) how do the two nitro substituents quantitatively alter the kinetics and selectivity of reactions with O₃, •OH, and •NO₂ relative to mononitrophenols? (ii) how does interfacial hydration regulate the adsorption of 2,4-DNP on carbonaceous surfaces at the molecular level? (iii) does oxidative processing of 2,4-DNP detoxify the products, or could it (based on computational predictions) lead to products with higher toxicity?

To answer these questions, this study employs an integrated multiscale theoretical framework combining density functional theory (DFT), molecular dynamics (MD) simulations, and computational toxicology to investigate selected oxidation and nitration pathways of 2,4-DNP in gas and aqueous phases, its interfacial adsorption behavior on a carbonaceous aerosol surrogate under different hydration conditions, and predicted toxicity indicators of its transformation products. The results provide molecular-level insights into selected transformation pathways and interfacial interactions of 2,4-DNP. This work focuses exclusively on 2,4-DNP as a representative dinitrophenol, and all quantitative conclusions are strictly limited to this compound; however, the qualitative framework established here may serve as a reference for understanding related nitroaromatic systems.

2 Computational methods

2.1 Density functional theory

All density functional theory calculations in this study were performed using the Gaussian 16 software package (Frisch et al., 2016). Geometry optimization, frequency analysis, and single-point energy calculations were conducted for all structures along the reaction pathway, including reactants (Rs), transition states (TSs), intermediates (IMs), and products (Ps), employing the M06-2X functional (Zhao and Truhlar, 2008). To validate the reliability of the chosen basis set, the representative reaction 2,4-DNP + •OH → IM3 was selected, and reaction energies using multiple basis sets were systematically compared, as summarized in Table S1 of the Supplement. The reaction energies computed with different basis sets were in excellent agreement. Among the tested combinations, the M06-2X/6-311++G(3df,2p) level yields reaction energies closest to the high-accuracy CCSD(T) benchmark (Spiegel et al., 2024; Nagy, 2024). Considering both computational accuracy and efficiency, the composite method M06-2X/6-31+G(d,p)//6-311++G(3df,2p) was consistently applied in all subsequent calculations. Transition states were located using the TS optimization algorithm and further validated through vibrational frequency analysis and intrinsic reaction coordinate (IRC) calculations (Maeda et al., 2015). Frequency analysis confirmed that all stable structures exhibit no imaginary frequencies, while each transition state possesses exactly one imaginary frequency, consistent with the characteristics of a first-order saddle point, thereby validating the reaction pathway. The IRC trajectories for key transition states are shown in Fig. S1 in the Supplement. In addition, the solvation model based on density (SMD) (Kromann et al., 2018) was employed to analyze the electronic structures, reactive sites, and free energy changes of key molecules in solvated states.

To gain deeper insights into the ground-state electronic structure and the delocalization characteristics of potential reactive sites, orbital-weighted dual descriptor (OWDD) (Pino-Rios et al., 2019) analysis was conducted using the Multiwfn program (Lu and Chen, 2012). Furthermore, the independent gradient model (IGM) approach was employed to systematically examine both chemical bonding and noncovalent weak interactions within the system (Lefebvre et al., 2017; Lu and Chen, 2022).

2.2 Molecular dynamics simulations

All MD simulations were performed using the GROMACS 2018.8 software package (Van Der Spoel et al., 2005). The simulation system was constructed within a periodic boundary box measuring 30 Å ×30 Å ×60 Å. Each system contained 6 C60molecules and 15 2,4-DNP molecules, solvated with varying numbers of water molecules (100, 200, 300, and 400), using the explicit SPC/E water model. To construct an idealized gas-aqueous interfacial model, a 20 Å vacuum layer was introduced along the z-axis to create an explicit air-water interface.

The AMBER99SB-ILDN force field was employed for C60 and water molecules, whereas the generalized AMBER force field (GAFF) was used to parameterize 2,4-DNP (Shirts et al., 2003; Wang et al., 2004; Lindorff-Larsen et al., 2010). All topology files were generated with the interactive sobtop program (Lu, 2026). Prior to production runs, each system underwent energy minimization to eliminate unfavorable atomic overlaps. The systems were subsequently equilibrated for 500 ps in the NVT ensemble at 298.15 K with a 2 fs integration time step, followed by production MD simulations of 1 ns for data collection and analysis. All molecular visualizations and trajectory analyses were performed using the VMD software (Humphrey et al., 1996).

2.3 Kinetic calculations

The rate constants for the primary reactions were calculated using the thermodynamic formulation of the conventional transition state theory (TST) (Tuñón and Williams, 2019) with support from the KISthelp tool (Canneaux et al., 2014) and TST calculator (Lu, 2025), as represented by the following equation:

$$\begin{array}{}\text{(1)}& \mathrm{kTST}=\mathit{\sigma }\times k_{\mathrm{B}}T/h\times e-\mathrm{\Delta }G\ne /\left(k_{\mathrm{B}}T\right)\end{array}$$

where σ is the reaction path degeneracy, k_B is the Boltzmann constant, T is the temperature, and h is Planck's constant. ΔG^≠ represents the standard Gibbs free energy of activation for the considered reaction. For reaction path simplicity, the following equation is followed:

$$\begin{array}{}\text{(2)}& \mathit{\sigma }={\mathit{\sigma }}_{\mathrm{rot},\mathrm{R}}/{\mathit{\sigma }}_{\mathrm{rot},\mathrm{TS}}\end{array}$$

where σ_rot,R and σ_rot,TS are rotational symmetry numbers of reactants and transition states (for N molecule reactions, σ_rot,R is the product of the rotational symmetry numbers of all the reactants); σ is generally 1 for reactions of macromolecular systems.

2.4 Computational toxicology

The toxicological profiles of 2,4-DNP and its atmospheric transformation products were characterized through a hierarchical QSAR strategy (Young et al., 2008). Key toxicity endpoints, including Ames mutagenicity (Sushko et al., 2010) and developmental toxicity (Cassano et al., 2010), were predicted using the Toxicity Estimation Software Tool (T.E.S.T.) version 5.1.2 developed by the U.S. Environmental Protection Agency (EPA, 2020). All predictions were subjected to applicability domain verification based on descriptor leverage analysis, and only data for compounds falling within the chemical space defined by the training set were retained to ensure the statistical validity and reliability of subsequent comparative toxicological analyses (Benfenati et al., 2009).

3 Results and discussion

3.1 Reaction mechanisms of 2,4-DNP with O₃, •OH, and •NO₂ in the gas and aqueous phases

Establishing the stable molecular configuration and regional reactivity of 2,4-DNP is essential for understanding the initial reaction mechanisms with O₃ (an electrophile), •OH (a radical), and •NO₂ (which exhibits both electrophilic and radical characteristics). Such characterization provides critical theoretical insights into the distinct mechanistic origins of these reaction processes. To elucidate the fundamental reaction pathways and selectivity differences of 2,4-DNP with various reactants in multiphase environments, this study initiates a theoretical characterization of its ground state, including its electronic structure and reactive sites.

3.1.1 Characterization of the electronic structure of 2,4-DNP

Based on geometric optimization and electronic structure calculations, the stable conformation of 2,4-DNP and its reactive sites analysis are shown in Fig. 1. The global minimum energy conformation stabilized by intramolecular hydrogen bonding is presented in Fig. 1a, while the OWDD isosurface in Fig. 1b clearly delineates the nucleophilic and electrophilic reaction regions of the molecule, thus providing a robust theoretical foundation for elucidating the subsequent reaction mechanisms.

Figure 1 (a) The most stable configuration of 2,4-DNP, stabilized by intramolecular hydrogen bonding. (b) Isosurface of the OWDD for 2,4-DNP, where blue and red regions denote nucleophilic and electrophilic sites, respectively.

The global minimum energy conformation of 2,4-DNP is depicted in Fig. 1a. The exceptional stability of this structure arises primarily from the synergistic effects of intramolecular hydrogen bonding and π-electron delocalization. Structural analysis reveals a strong intramolecular hydrogen bond (O–H$\cdot \cdot \cdot$O₂N–) between the hydroxyl group and the ortho-positioned nitro group (position 2), with an H$\cdot \cdot \cdot$O distance of 1.75 Å, a value characteristic of strong intramolecular hydrogen bonding (1.6–1.9 Å) (Steiner, 2002). This interaction forms a stable, planar six-membered ring motif that benefits from both the hydrogen bond and the rigidity of the aromatic ring. Additionally, both the ortho- and para-nitro groups (at positions 2 and 4, respectively) adopt nearly coplanar orientations relative to the benzene ring, with dihedral angles of 178.30 and 179.67°, approaching the ideal 180° for maximum orbital overlap. This coplanarity enables effective π–π conjugation with the aromatic system and further enhances the overall molecular stability (Szatylowicz et al., 2017). Importantly, these two stabilizing factors are mutually reinforcing: the hydrogen bond locks the ortho-nitro group into the plane, facilitating conjugation, while the enhanced conjugation further polarizes the nitro group, strengthening the hydrogen bond.

As shown in Fig. 1b, the OWDD isosurface analysis (Δfω) reveals well-defined regional chemical features of the 2,4-DNP. Blue regions (Δfω<0) are localized around specific carbon atoms on the benzene ring and the hydroxyl oxygen atom, indicating electron-rich sites that are susceptible to electrophilic attack. These electron-rich sites arise from the electron-donating effect of the hydroxyl group, which concentrates electron density at the ortho and para positions relative to the nitro groups. In contrast, green regions (Δfω>0) are mainly located at the outer ends of the oxygen and nitrogen atoms in the nitro and hydroxyl groups, showing electron-poor sites that are susceptible to nucleophilic attack, which arise from the polarization of N–O bonds. This electronic structure feature provides a reliable basis for predicting its reactivity toward various atmospheric oxidants. As a strong electrophile, O₃ preferentially targets the electron-rich carbon atoms (blue regions), initiating electrophilic addition at the sites with the highest nucleophilic character. Meanwhile, •OH and •NO₂ can participate in both ring addition and hydrogen-abstraction reactions. However, the preferred HAA sites cannot be inferred solely from the ground-state Δfω distribution and are evaluated explicitly from the calculated reaction energetics below.

3.1.2 Primary reactions of 2,4-DNP with O₃, •OH, and •NO₂

Based on the analysis of nucleophilic and electrophilic sites in the electronic structure of 2,4-DNP, this section systematically elaborates the primary reaction pathways of 2,4-DNP with O₃, •OH, and •NO₂. The competition among different channels is elucidated through thermodynamic and kinetic calculations, thereby characterizing the selected transformation pathways examined in this study.

Figure 2 Proposed initial reaction mechanisms of 2,4-DNP with O₃, •OH, and •NO₂.

The coexistence of an electron-donating hydroxyl group and strong electron-withdrawing nitro groups in the 2,4-DNP induces a highly polarized and asymmetric distribution of π-electron density in the benzene ring (Li et al., 2023). This distinct electronic structure directly governs the pronounced regioselectivity of its reactions with O₃, •OH, and •NO₂.

As a typical 1,3-dipole, ozone can attack the electron-rich aromatic π-system via cycloaddition. Calculations indicate that among the six possible 1,3-dipolar cycloaddition pathways, the reaction preferentially occurs at the C5–C6 bond located between the hydroxyl and nitro groups, as indicated by the red region in Fig. 2. The reaction proceeds via transition state TS6, in which the terminal oxygen atoms of O₃ simultaneously approach the C5 and C6 atoms, accompanied by significant elongation of the original C5–C6 bond (optimized geometries are compiled in Table S2 in the Supplement), indicating irreversible disruption of aromatic stabilization. This step represents the most favorable reaction pathway. Subsequently, the resulting unstable primary ozonide (POZ) intermediate IM6 undergoes concerted dissociation via transition state TS6-17, leading to the formation of the vibrationally excited Criegee intermediate (CI) IM17, which is energetically the most favorable product.

The reaction of 2,4-DNP with •OH mainly proceeds through two competing mechanisms, namely ring addition and HAA, as shown in the blue region of Fig. 2. In the electrophilic addition pathway, •OH preferentially attacks the C3 site, which is located between the two nitro groups and is meta to the hydroxyl substituent. This process proceeds through transition state TS3', in which the •OH oxygen approaches C3 with significant shortening of the forming bond distance, accompanied by elongation of the adjacent C2–C3 and C3–C4 bonds and a geometry that facilitates effective overlap between the p-orbital of •OH and the π-system of the benzene ring, thereby ultimately yielding the hydroxy-cyclohexadienyl radical intermediate IM3'. In the HAA pathway, •OH primarily captures the H atom directly bonded to the C5 site. This reaction proceeds through transition state TS9', characterized by a nearly linear O–H–C5 arrangement, which is characteristic of a hydrogen abstraction transition state. This step produces the phenoxy radical intermediate IM9' and one molecule of H₂O. Notably, the OWDD analysis indicates that the carbon at C5 is an electron poor site, which is susceptible to nucleophilic attack. This appears contrary to the conventional expectation that an electrophilic radical such as •OH would preferentially abstract a hydrogen atom from an electron rich site. However, the HAA reaction of •OH is governed primarily by C–H bond dissociation energies (BDEs) rather than ground state electron density (Mayer, 2011). The preference for C5 abstraction is therefore more appropriately interpreted in terms of the calculated reaction energetics and transition-state structure rather than the ground-state electron-density distribution alone.

The reaction of •NO₂ with 2,4-DNP proceeds via two competing pathways, radical addition and HAA, as shown in the green region of Fig. 2. Although the reaction pathways of 2,4-DNP with •NO₂ and •OH both involve radical addition and hydrogen abstraction, significant differences are observed in their reactivity and pathway selectivity. In the radical addition pathway, •NO₂ preferentially attacks the C3 site via transition state TS3*, leading to the formation of the nitrocyclohexadienyl radical intermediate IM3*. This regioselectivity is similar to that calculated for •OH addition. However, unlike O₃ and •OH, whose reactivities sites show a clearer relationship with the Δfω distribution, •NO₂ has only weak electrophilic character, and its addition pathway cannot be interpreted solely in terms of electrophile–nucleophile interactions. Radical stabilization and transition-state effects may also contribute to the calculated regioselectivity. The HAA pathway exhibits different site selectivity. Whereas •OH preferentially abstracts the C5-bound hydrogen in the gas phase, the lowest-energy HAA channel identified for •NO₂ involves abstraction of the OH-bound hydrogen through transition state TS7*. Nevertheless, this preference should be distinguished from the overall kinetic feasibility of the pathway. The calculated gas-phase barrier for OH hydrogen abstraction by •NO₂ is 35.09 kcal mol⁻¹, substantially higher than the 18.48 kcal mol⁻¹ barrier for the competing C3-addition pathway. Thus, although OH of C1 is the preferred H-abstraction site among the investigated direct •NO₂ HAA channels, direct HAA by •NO₂ remains kinetically unfavorable for parent 2,4-DNP.

Figure 3 Thermodynamics, kinetics, and temperature dependence of the reactions of 2,4-DNP with O₃, •OH, and •NO₂. (a–c) Reaction energy barriers (ΔE^≠) and reaction energies (Δ_rE) for the key pathways (unit: kcal mol⁻¹). (d–f) Intrinsic temperature-dependent reaction rate constants for the major channels.

The combined thermodynamic and kinetic analysis reveals that the reaction of O₃ with 2,4-DNP is both thermodynamically favorable and kinetically accessible in the gas and aqueous phases. As shown in Fig. 3a, the most favorable gas-phase pathway, R6, has a low energy barrier of 10.12 kcal mol⁻¹ and is strongly exothermic (${\mathrm{\Delta }}_{\mathrm{r}}E=-\mathrm{35.56}$ kcal mol⁻¹). In the aqueous phase, although the barrier increases moderately to 11.62 kcal mol⁻¹, the reaction remains significantly exothermic. This indicates that solvation does not fundamentally alter the thermodynamic driving force, and the calculated barrier remains moderate within the present reaction framework. In contrast, the reaction of 2,4-DNP with •OH shows even more favorable calculated kinetics than the O₃ and direct •NO₂ reactions examined here. The addition pathway in the gas phase proceeds with an exceptionally low barrier of 2.84 kcal mol⁻¹, substantially lower than the O₃-initiated pathway, as presented in Fig. 3b. The HAA pathway also presents a low barrier of 4.70 kcal mol⁻¹. The similar barrier heights suggest that both pathways may be kinetically competitive at the molecular level. By comparison, the direct reaction of 2,4-DNP with •NO₂ is much less favorable. The overall process is endothermic and thermodynamically less favorable within the calculated energy framework, as illustrated in Fig. 3c. Among the two plausible direct reaction pathways, radical addition to the aromatic ring via the nitration pathway R3* has a calculated energy barrier of 18.48 kcal mol⁻¹, whereas HAA pathway has an even higher barrier of 35.09 mol⁻¹, making it far less accessible. Although the nitration channel is kinetically favored over HAA in the •NO₂ system, its energy barrier is still much higher than those of the corresponding •OH-initiated pathways and is also less thermodynamically favorable than both •OH- and O₃-mediated reactions.

Given these pronounced kinetic limitations, combined with the typically lower atmospheric concentration of •NO₂ relative to •OH and O₃, the calculations might suggest that •NO₂ reactivity is negligible (Shi et al., 2025; Chapleski et al., 2016). To distinguish the different roles of •NO₂ within the reaction framework examined here it is useful to distinguish three distinct levels of •NO₂ involvement in the atmospheric chemistry of 2,4-DNP, a distinction that has been emphasized in recent multiphase studies (Rana et al., 2023). Level (i) is the direct reaction of •NO₂ with the parent compound 2,4-DNP, which shows substantially lower calculated reactivity than the •OH pathway under the concentration scenarios considered here. Level (ii) is the reaction of •NO₂ with primary oxidation intermediates such as IM3' (Sect. 3.1.3), which exhibits a remarkably low barrier of 2.01 kcal mol⁻¹ for HAA, indicating that this indirect •NO₂ pathway is mechanistically feasible. However, its atmospheric contribution depends on competition with O₂ and other reactants, as discussed in Sect. 3.1.3. Level (iii) is the multiphase •NO₂ $/$ O₃ chemistry in aerosol water or microdroplets, where •NO₂ can participate through N₂O₄ uptake and hydrolysis, leading to nitration via NO${}_{\mathrm{2}}^{+}$ or NO⁺ species (Rana et al., 2023). This multi-level distinction clarifies that, although the direct reaction of •NO₂ with parent 2,4-DNP is kinetically unfavorable in the investigated scenarios, the indirect pathways involving •NO₂ and oxidized intermediates or multiphase processes represent additional mechanisms within the present reaction framework. These pathways should not be considered a complete description of nighttime reactive-nitrogen chemistry because direct NO₃ oxidation was not included in the present calculations.

This study further shows that the reaction behavior of 2,4-DNP in the aqueous phase differs markedly from that in gas phase and from that of simpler nitrophenol analogues. This difference is fundamentally attributed to the restructuring of molecular electronic distribution and reaction coordinates by the dynamic hydrogen-bonding network established between water and 2,4-DNP (Lyu et al., 2025). In contrast to the common trend in which solvation stabilizes polar transition states and thus lowers energy barriers in mononitrophenols, 2,4-DNP exhibits the distinctive features of increased aqueous-phase energy barriers and solvent-induced switching of regioselectivity.

The increase in the calculated aqueous-phase barriers can be attributed to changes in the balance between intramolecular and solute-solvent interactions upon solvation. The simultaneous presence of the hydroxyl group and nitro substituents produces pronounced electronic polarization, while the ortho-substitution pattern enables intramolecular hydrogen bonding (Albinet et al., 2010). These structural features provide significant ground-state stabilization in the gas phase. In aqueous solution, water molecules interact with the polar functional groups and compete with the intramolecular hydrogen-bonding network, thereby altering the relative stabilization of the reactant and transition-state structures (Yang et al., 2023a). Solvent reorganization may also contribute to the calculated free-energy barrier. Consequently, the overall barrier depends on the relative solvation of the reactant and transition state and can therefore increase for the pathways considered here. The change in site selectivity likewise arises from solvent-mediated electron redistribution and pathway switching. In the ring-addition pathway of •OH, the C3 position is more favorable in the gas phase because of the dominant conjugation effects. In the aqueous phase, however, water molecules enhance the local polarization near the hydroxyl-adjacent carbon C6, making another site more competitive for attack. For the •OH HAA pathways, the preferred abstraction site also changes upon solvation. In the gas phase, abstraction from C5 is favored among the investigated HAA channels, whereas in the aqueous phase the hydroxyl hydrogen becomes more favorable in the calculated pathway. This change suggests that solvation modifies the relative energetics of the competing H-abstraction sites, likely through differential stabilization of the reactant and transition-state structures. In contrast, the preferred HAA channel identified for •NO₂ involves abstraction from OH in both phases, although this pathway remains kinetically unfavorable because of its high calculated barriers.

Kinetic calculations further quantify the feasibility differences among various pathways under atmospheric conditions. To ensure physical significance and environmental interpretability, all rate constants were converted into the standard unit commonly used in atmospheric chemistry (cm³ molec.⁻¹ s⁻¹). Aqueous-phase data were converted from their original units of M⁻¹ s⁻¹ using the Avogadro constant and volume-unit conversion (1 M⁻¹ s${}^{-\mathrm{1}}=$ (1000/NA) cm³ molec.⁻¹ s${}^{-\mathrm{1}}\approx \mathrm{1.66}\times {\mathrm{10}}^{-\mathrm{21}}$ cm³ molec.⁻¹ s⁻¹).

Table 1 Calculated energy barriers (ΔE^≠, kcal mol⁻¹), reaction energies (Δ_rE, kcal mol⁻¹), and rate constants (k_tot, cm³ molec.⁻¹ s⁻¹) at 298.15 K for the representative reaction channels of 2,4-DNP with O₃, •OH, and •NO₂ in the gas and aqueous phases. “–” indicates that the rate constant was not calculated because the corresponding pathway is kinetically insignificant (Complete rate constant data for all channels are provided in Table S3). ktot at 298.15 K for each oxidant is the sum of all contributing channels for that oxidant.

As shown in Fig. 3d–f, the rate constants of all major channels increase with temperature over the range of 200–400 K, exhibiting Arrhenius-type temperature dependence. The total rate constants of the dominant pathways at 298.15 K are compiled in Table 1. Among the reactions examined, the •OH pathway shows the highest calculated gas-phase rate coefficient, with total gas-phase and aqueous-phase rate constants of $\mathrm{3.18}\times {\mathrm{10}}^{-\mathrm{15}}$ and $\mathrm{2.25}\times {\mathrm{10}}^{-\mathrm{17}}$ cm³ molec.⁻¹ s⁻¹, respectively. The O₃ pathway also shows a calculated aqueous-phase rate coefficient of $\mathrm{3.17}\times {\mathrm{10}}^{-\mathrm{17}}$ cm³ molec.⁻¹ s⁻¹, indicating appreciable intrinsic reactivity in the aqueous phase. In contrast, the direct •NO₂ pathway is kinetically unfavorable, with calculated gas-phase and aqueous-phase rate constants of $\mathrm{1.75}\times {\mathrm{10}}^{-\mathrm{28}}$ and $\mathrm{3.93}\times {\mathrm{10}}^{-\mathrm{28}}$ cm³ molec.⁻¹ s⁻¹, respectively. These bimolecular rate coefficients describe the intrinsic reactivity of the investigated pathways, but do not by themselves determine their atmospheric importance, which also depends on oxidant concentrations and phase partitioning. To examine whether the very low intrinsic reactivity of the direct •NO₂ pathway could be offset by nighttime •NO₂ accumulation, pseudo-first-order rate coefficients were estimated under illustrative scenarios based on literature-reported oxidant concentrations.

The literature review indicates that •OH concentrations typically reach $\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}\sim \mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}$ molec. cm⁻³, while •NO₂ is rapidly photolyzed during daytime. Under polluted nighttime conditions, however, •OH concentrations drop substantially. Geyer et al. (2003) reported a nocturnal •OH concentration of (1.85±0.82) ×10⁵ molec. cm⁻³ during the BERLIOZ campaign, while nighttime •OH concentrations of 10⁵∼10⁶ molec. cm⁻³ have been reported in UK urban surface environments (Khan et al., 2011). In contrast, •NO₂ can accumulate significantly at night, with nocturnal maxima reaching 8.43±0.4 ppb in Guangzhou (Qin et al., 2009). Using these literature-reported concentrations as illustrative nighttime scenarios, the pseudo-first-order rate constants for the •OH pathway is estimated to be approximately 10${}^{-\mathrm{10}}\sim {\mathrm{10}}^{-\mathrm{9}}$ s⁻¹, based on $k_{\mathrm{OH}}=\mathrm{3.18}\times {\mathrm{10}}^{-\mathrm{15}}$ cm³ molec.⁻¹ s⁻¹. In contrast, for the direct •NO₂ pathway, a bimolecular rate coefficient of $k_{{\mathrm{NO}}_{\mathrm{2}}}\approx {\mathrm{10}}^{-\mathrm{28}}$ cm³ molec.⁻¹ s⁻¹ combined with a nighttime •NO₂ concentration of several ppb, yields a pseudo-first-order rate coefficient on the order of 10⁻¹⁷ s⁻¹. Therefore, under the illustrative nighttime conditions considered here, the direct •NO₂ pathway remains more than seven orders of magnitude slower than the •OH pathway and is unlikely to compete effectively with •OH oxidation. Importantly, these estimates compare only the investigated bimolecular reaction channels and do not represent the overall atmospheric lifetime of 2,4-DNP, which would also depend on photolysis, deposition, phase partitioning, heterogeneous chemistry, and other loss processes. The aqueous-phase calculations further indicate appreciable intrinsic reactivity toward O₃. Accordingly, the relative contributions of the investigated pathways to the overall atmospheric loss of 2,4-DNP cannot be quantified without additional information on phase partitioning and more comprehensive atmospheric chemistry modeling. The possible involvement of •NO₂ in subsequent reactions with primary oxidation intermediates is discussed separately below.

The atmospheric contribution of aqueous-phase reactions additionally also depends on the fraction of 2,4-DNP present in the aqueous phase. The Henry's law constant for 2,4-DNP has been reported as $\mathrm{5.70}\times {\mathrm{10}}^{-\mathrm{8}}$ atm m³ mol⁻¹ at 5 °C (Lüttke and Levsen, 1997). Its acid dissociation constant (pKa ≈4.09 at 25 °C) further affects aqueous partitioning because the deprotonated form is more water-soluble than the neutral species (Schwarzenbach et al., 1988). Aerosol liquid water content also determines the aqueous volume available for multiphase reactions. Without quantitative information on these partitioning parameters under specific atmospheric conditions, the aqueous-phase rate constants should be interpreted as intrinsic reactivity indicators rather than direct measures of atmospheric importance.

Table 2 The experimental rate constant (cm³ molec.⁻¹ s⁻¹) values of various phenols with O₃ in the gas phase.

The calculated rate coefficients were further compared with available experimental and computational data for structurally related aromatic compounds, as summarized in Tables 2 and S4. Although direct experimental measurements for this specific reaction are currently unavailable, the theoretical rate constant for the mononitro analogue nitrocatechol ($\mathrm{1.48}\times {\mathrm{10}}^{-\mathrm{20}}$ cm³ molec.⁻¹ s⁻¹) (Wang et al., 2026) is approximately three orders of magnitude lower than that of catechol (Sun et al., 2022; Tomas et al., 2003), consistent with a strong deactivating effect of the nitro group, as an electron-withdrawing functional group on the ozone reactivity of the aromatic ring. Building on this trend, the calculated rate constant for the dinitro compound 2,4-DNP ($\mathrm{4.23}\times {\mathrm{10}}^{-\mathrm{23}}$ cm³ molec.⁻¹ s⁻¹) is reduced by an additional three orders of magnitude, quantitatively reflecting the cumulative electron-withdrawing effect of the two nitro substituents. This trend is consistent with the established understanding of substituent effects in phenolic oxidation: electron-donating groups such as –CH₃ and –OCH₃ enhance reactivity (Pillar-Little et al., 2014, 2015; Pillar-Little and Guzman, 2017; Rana and Guzman, 2020), while electron-withdrawing nitro groups deactivate the aromatic ring toward electrophilic attack.

Furthermore, the predominance of substituent electronic effects in governing reactivity is independently corroborated. As detailed in Table S4, a systematic comparison of the theoretical rate constants for the reaction of 2,4-DNP with •OH with data for structural analogues demonstrates that the strong electron-withdrawing character of the nitro groups likewise controls the kinetics of this pathway. This consistency suggests that the electronic influence of nitro substituent is general across different oxidants rather than being restricted to a single mechanism. Collectively, these comparisons support the qualitative role of electron-withdrawing nitro substituents in modulating the calculated reactivity of 2,4-DNP and related nitroaromatic compounds.

Figure 4 The subsequent reactivity mechanism of the key intermediate IM3' (unit: kcal mol⁻¹).

3.1.3 The subsequent reactivity mechanism of the key intermediate IM3'

Among the initiation pathways examined here, •OH addition generates the key radical intermediate IM3' via electrophilic addition. This highly reactive species can rapidly undergo secondary reactions with surrounding atmospheric constituents. The following calculations examine selected secondary reactions of IM3' with O₂ and •NO₂. These reactions were considered to evaluate the relative energetics of representative HAA and radical-addition pathways.

A systematic analysis was conducted for the reaction pathways of the intermediate IM3' with O₂ and •NO₂. For the reaction with O₂, as shown in Fig. 4a, IM3' mainly evolves through two competing pathways, namely HAA and radical addition. The HAA pathway contains four distinct subchannels, among which the RB (hydrogen abstraction at C3) channel is the most favorable both thermodynamically and kinetically. It is characterized by a remarkably low transition-state energy barrier of 2.86 kcal mol⁻¹ and a highly exothermic reaction energy of −65.36 kcal mol⁻¹. In this process, the O₂ molecule approaches the activated C3–H bond. The reaction is directed to the C3 position because this carbon is flanked by two strongly electron-withdrawing nitro groups at the adjacent C2 and C4 positions. This electronic effect substantially lowers the electron density and weakens the C3–H bond, making it the most favorable site for hydrogen abstraction. Key structural parameters in transition state TSB reveal a concerted mechanism: the distance between the attacking oxygen atom Oa of O₂ and the transferring hydrogen shortens significantly, while the breaking C3–H bond elongates accordingly, indicating synchronous bond breaking and bond forming. This synchronous bond-breaking and bond-forming process yields the phenoxy radical intermediate IMB and the hydroperoxyl radical (•HO₂). In contrast, the other HAA channels involving abstraction from the C5, C6, or hydroxyl positions have substantially higher calculated barriers and are therefore less favorable within the investigated reaction framework. In the competing radical addition pathway, the most favorable channel is RI, in which O₂ adds to the C4 position via the transition state TSI. This process has an energy barrier of 3.19 kcal mol⁻¹ and is exothermic by −37.52 kcal mol⁻¹. Because this barrier is close to that of the optimal HAA channel RB, the HAA and radical addition pathways for the reaction of IM3' with O₂ are expected to compete directly at the molecular level.

The reaction of IM3' with •NO₂ similarly follows a competitive mechanism involving HAA and radical addition (Fig. 4b). Within the HAA pathway, the RB* channel shows a clear advantage analogous to that observed in the O₂ system, with a low energy barrier of 2.01 kcal mol⁻¹ and a highly exothermic reaction energy of −59.89 kcal mol⁻¹. The mechanism of this process is highly consistent with that of the O₂-involved RB channel, in which the oxidant abstracts the hydrogen atom from the C3 position, yielding the identical phenoxy radical intermediate IMB, while producing HONO or •HO₂, respectively. In the radical addition pathway, the most feasible route is the RI* channel, with a barrier of 2.78 kcal mol⁻¹ and a reaction energy of −30.99 kcal mol⁻¹. Similar to the O₂ addition pathway, •NO₂ preferentially adds to the C4 position, generating the nitro-substituted cyclohexadienyl radical intermediate IMI*. Comparison of the lowest energy barriers for the two pathways reveals that, within the investigated IM3'+ •NO₂ reaction system, the HAA pathway has a slightly lower calculated barrier than the radical-addition pathway. Thus, HAA is the lower-barrier pathway among the •NO₂ reaction channels examined here, consistent with the corresponding barrier comparison for the O₂ reaction system. This corresponds to Level (ii) of •NO₂ involvement as distinguished in Sect. 3.1.2: the reaction of •NO₂ with a primary oxidation intermediate. The calculated low barrier for the •NO₂ + IM3' HAA pathway (2.01 kcal mol⁻¹) indicates that this indirect •NO₂ pathway is mechanistically feasible. However, its atmospheric contribution cannot be inferred from the energy barrier alone and depends on competition with O₂ and other reactants.

To illustrate this competition, effective pseudo-first-order rates (k × [oxidant]) were estimated using representative atmospheric concentrations and literature-based rate coefficients. The atmospheric concentration of O₂ was taken as 5.4×10¹⁸ molec. cm⁻³ (0.21 atm). For the reaction of IM3' with O₂ via HAA (calculated barrier =2.86 kcal mol⁻¹), rate constants for analogous radical + O₂ reactions are typically in the range of 10${}^{-\mathrm{12}}\sim {\mathrm{10}}^{-\mathrm{11}}$ cm³ molec.⁻¹ s⁻¹ (Miyoshi et al., 1990; Morgan et al., 1982). Using a conservative estimate of $k_{{\mathrm{O}}_{\mathrm{2}}}\approx {\mathrm{10}}^{-\mathrm{12}}$ cm³ molec.⁻¹ s⁻¹, gives an effective pseudo-first-order rate coefficient of $k_{{\mathrm{O}}_{\mathrm{2}}}^{\prime }\approx {\mathrm{10}}^{-\mathrm{12}}\times \mathrm{5.4}\times {\mathrm{10}}^{\mathrm{18}}\approx \mathrm{5.4}\times {\mathrm{10}}^{\mathrm{6}}$ s⁻¹, corresponding to an illustrative reaction timescale of approximately 0.2 µs. This estimate is consistent with the generally rapid reactions of organic radicals with O₂ (Vereecken and Francisco, 2012; Calvert et al., 2015).

For the reaction of IM3' with •NO₂ via HAA (barrier =2.01 kcal mol⁻¹), assuming a representative nighttime •NO₂ concentration of 10¹⁰∼10¹¹ molec. cm⁻³ and a literature-based rate coefficient of $k_{{\mathrm{NO}}_{\mathrm{2}}}\approx {\mathrm{10}}^{-\mathrm{11}}$ cm³ molec.⁻¹ s⁻¹ gives $k_{{\mathrm{NO}}_{\mathrm{2}}}^{\prime }\approx {\mathrm{10}}^{-\mathrm{1}}\sim \mathrm{1}$ s⁻¹, corresponding to an illustrative reaction timescale of approximately 1–10 s. Under these assumptions, the effective O₂ reaction rate exceeds that of the investigated •NO₂ pathway by more than six orders of magnitude. Thus, although the low calculated barrier demonstrates the mechanistic feasibility of the indirect •NO₂ pathway, this illustrative effective-rate comparison suggests that reaction with O₂ would strongly outcompete the investigated •NO₂ pathway under the conditions considered here. Quantifying the atmospheric branching of IM3' would require reaction-specific rate coefficients and a more complete treatment of competing loss processes.

3.2 Molecular dynamics

Although the homogeneous oxidation mechanisms of 2,4-DNP in the gas and aqueous phases have been clarified, its heterogeneous interfacial behavior on aerosol particle surfaces remains to be elucidated. C60 provides a well-defined model system for probing the non-covalent interactions between pollutants and particulate matter, owing to its well-defined carbon framework and surface curvature (Athar and Patnaik, 2022). The results obtained with this surrogate provide molecular-level insights into specific interaction mechanisms, but should not be directly extrapolated to real atmospheric aerosols without considering the additional complexity of real particle surfaces. In this study, the IGM was employed to qualitatively and quantitatively analyze the weak interactions within the adsorption complexes (Lu and Chen, 2022), thereby revealing the interactions that govern the interfacial behavior at the molecular level.

The intermolecular interactions within the 2,4-DNP-C60 aggregate were characterized utilizing the IGM by defining each monomer as an individual fragment. The intermolecular interaction pattern at the 2,4-DNP-C60 interface is elucidated jointly by the three-dimensional IGM isosurface map (Fig. 5a) and the two-dimensional density map (Fig. 5b), which exhibit excellent consistency. According to the principles of IGM analysis, the green and light-green isosurface regions, where sign(λ2) ≈0, correspond to van der Waals (vdW) interactions, indicative of diffuse and non-directional weak attractions. In contrast, the blue regions, characterized by sign(λ2) <0, indicate directional interactions such as hydrogen bonds or strong electrostatic attractions, which are associated with more localized electron-density distributions.

Figure 5 IGM analysis of intermolecular interactions in the 2,4-DNP-C60 aggregate. (a) Three-dimensional (3D) IGM isosurface map. (b) Two-dimensional (2D) IGM density map.

Specifically, the extensive green/light-green isosurface between the C60 cage and the benzene ring of 2,4-DNP corresponds closely to the continuous blue-to-light-blue signal band in the 2D density map (δg_inter≈0.10–0.30), confirming that vdW interactions, mainly dispersion and induction forces arising from the π-electron clouds of the benzene ring and fullerene surface, dominate the interfacial binding between 2,4-DNP and C60. Meanwhile, distinct blue regions are observed near the polar functional groups of 2,4-DNP, namely the nitro and hydroxyl moieties. These features manifest as localized yellow-green hotspots in the 2D density map (δg_inter≈0.60–0.80), indicating stronger polar interactions. These interactions can be attributed to specific directional forces, including a O–H$\cdot \cdot \cdot \mathit{\pi }$ hydrogen bond involving the hydrogen atom of the –OH group and the delocalized π-electron cloud of C60, as well as a dipole-π interaction between the oxygen atoms of the –NO₂ groups and the fullerene surface. These synergistic polar interactions enhance both the binding selectivity and the overall stability of the 2,4-DNP-C60 interface.

Further analysis indicates that the 2,4-DNP molecule preferentially adsorbs with its benzene ring approximately parallel to the C60 surface. In this configuration, the nonpolar benzene ring maximizes the dispersion interaction area by closely approaching the carbon cage via vdW forces, while the polar –OH and –NO₂ groups remain oriented outward. This arrangement not only avoids steric repulsion with the C60 surface corresponding to strong steric hindrance, but also facilitates hydrogen-bond-network formation between the polar groups of neighboring 2,4-DNP molecules. In the 2D IGM cross-sectional map, this orientation is reflected by a clear spatial separation between the polar hotspots of the 2,4-DNP and the van der Waals signal band associated with C60, with the hotspot positions precisely corresponding to the atomic coordinates of the polar functional groups.

Figure 6 (a) Molecular dynamics simulation framework construction (MOL: 2,4-DNP; COL: C60; SOL: water). (b) z-axis number density profile of molecules (units: nm⁻³). (c) Hydration-dependent adsorption energies of 2,4-DNP on C60 (units: kJ mol⁻¹).

Building upon the fundamental interaction mechanisms revealed by IGM in the anhydrous system, molecular dynamics simulations were further carried out to investigate how hydration, as a key environmental factor, dynamically regulates the adsorption behavior of 2,4-DNP on C60 under atmospherically relevant conditions. Figure 6a illustrates the initial simulation setup, which consists of 15 molecules of 2,4-DNP and 6 molecules of C60 placed in a three-layer simulation box measuring $\mathrm{30}\times \mathrm{30}\times \mathrm{60}$ ų, with a gas-aqueous-gas sandwich structure. Different hydration levels were modeled by varying the number of water molecules from SOL100 to SOL400. The system then underwent energy minimization and molecular dynamics relaxation to obtain an equilibrated structure. In the aqueous phase, the distribution of 2,4-DNP and C60 gradually shifted toward the center of the water layer as the number of water molecules increased, indicating that the introduction of water significantly altered the spatial distribution pattern of the solutes.

The number-density distribution profiles along the z-axis in Fig. 6b further reveal the influence of hydration on the axial distribution and aggregation behavior of the solute molecules. Combined with the strong intermolecular interactions revealed by IGM in the anhydrous system, it can be inferred that under low hydration conditions, 2,4-DNP and C60 tend to form compact aggregates because of their relatively strong mutual attraction. As the number of water molecules increases from SOL200 to SOL400, the density peak shifts toward larger radial distances (4.0–4.5 nm) and its height gradually rises. This indicates that water molecules disrupt the close aggregation tendency between 2,4-DNP and C60 through competitive solvation, promoting the outward diffusion of solute molecules along the z-axis and thereby significantly enhancing the overall dispersity of the system.

Based on the temporal evolution of the interaction energies shown in Fig. 6c, the dynamic regulatory role of hydration on intermolecular adsorption behavior can be further interpreted from an energetic perspective. In the low-hydration system SOL100, the interaction energy remains stable at approximately −850 kJ mol⁻¹. Its relatively large negative value is consistent with the compact aggregated state shown in Fig. 6b, indicating strong direct adsorption between 2,4-DNP and C60 under these conditions. As the number of water molecules increases to SOL200, the absolute value of the interaction energy increases further, suggesting enhanced intermolecular attraction. This change may originate from moderate hydration, where water molecules may promote specific orientations and stabilize solute interactions through effects such as dielectric screening or indirect bridging (Lyu et al., 2025). However, when the hydration level further increases to SOL400, the absolute value of the interaction energy decreases, implying a significant weakening of intermolecular adsorption. This indicates that under high hydration conditions, water molecules suppress the direct interaction between 2,4-DNP and C60, thereby driving the system toward a more dispersed and stable solvated state. This non-monotonic energy trend, which first strengthens and then weakens, is consistent with the morphological evolution of solute distribution shown in Fig. 6b, namely a transition from an initially aggregated state to a more dispersed one. Together, these results reveal the complex regulatory role of hydration in interfacial adsorption behavior.

To assess statistical uncertainty arising from initial conditions, four independent replica simulations were performed for each hydration level by assigning distinct initial atomic velocities, thereby generating statistically independent trajectories in phase space. Error analysis demonstrates that the main conclusions regarding the adsorption configuration and distribution of 2,4-DNP at the gas-aqueous interface remain consistent across all replicas, underscoring the robustness and reproducibility of the results.

The MD simulations presented here are restricted to molecular-scale adsorption behavior in the idealized C60/water/2,4-DNP system. The results provide insights into how hydration modulates specific non-covalent interactions, but they do not account for key aerosol properties that govern humidity-dependent behavior in real atmospheric particles, including liquid water content, ionic strength, phase behavior, viscosity, and particle morphology. These properties vary with relative humidity, particle composition, and atmospheric conditions, thereby influencing particle physical state and water uptake of aerosol particles. The simulations therefore represent a simplified model system. The MD results should be interpreted as molecular-scale insights into adsorption behavior rather than as direct predictions of atmospheric persistence or humidity-dependent partitioning. The present simulations do not quantitatively constrain atmospheric residence time, gas-particle partitioning, or long-range transport. These quantities require additional partitioning calculations and atmospheric chemistry and transport modeling.

3.3 Predicted toxicity indicators

The reaction of ozone with aromatic rings follows a relatively well-defined mechanism and produces a comparatively clear product spectrum, making it a comparatively clear model for evaluating the directional impact of atmospheric oxidation on the toxicity of 2,4-DNP. By focusing on the O₃ pathway, this study provides a defined set of reactions for examining predicted changes in toxicity indicators along this specific transformation products. The conclusions drawn here may provide a useful reference for future investigations of the more complex atmospheric oxidation networks involving radicals such as •OH and •NO₂. Accordingly, the key intermediates and final products (IM17, IM25, IM44, IM54, P1–P4) identified in Fig. S2, were predicted using the T.E.S.T. tool (EPA, 2020) to predict Ames mutagenicity and developmental toxicity indicators and to examine changes in predicted toxicity indicators along this specific oxidation process. A conceptual diagram integrating the investigated DFT reaction pathways, MD-derived interfacial adsorption behavior, and QSAR toxicity predictions is provided in Fig. S3 to summarize the scope of the present computational framework. The framework considers gas-phase oxidation, bulk aqueous-phase oxidation, and interfacial adsorption behavior, with gas-aqueous partitioning influenced by Henry's law behavior and acid-base speciation. Oxidation pathways calculated by DFT generate key intermediates IM3' and IM17. These intermediates further transform into a series of products, and the Ames mutagenicity and developmental toxicity indicators of selected transformation products were evaluated using QSAR models.

Figure 7 Toxicity prediction of 2,4-DNP and its transformation products. (a) Ames mutagenicity. (b) Developmental toxicity.

Differences in the predicted Ames mutagenicity and developmental toxicity indicators between 2,4-DNP and the intermediates (IMs series) and final products (Ps series) generated from its reaction with ozone are intuitively illustrated in Fig. 7. Regarding Ames mutagenicity in Fig. 7a, the parent compound 2,4-DNP exhibits a predicted value of 0.13 that falls below the model classification threshold. The specific data are presented in Table S5. In contrast, the final products show higher predicted Ames mutagenicity indicators in the QSAR model. Notably, the predicted value for P1 (0.53) exceeds the model classification threshold. Within the QSAR framework, these results indicate that ozonolysis does not uniformly decrease the predicted Ames mutagenicity indicators of the modeled products.

For developmental toxicity (Fig. 7b), the parent compound has a predicted value of 0.55, whereas the modeled ozonolysis products generally show higher predicted values, with IM54 reaching 0.94. P2 and P4 have predicted values of 0.56, slightly higher than that of the parent compound. The differences between the Ames mutagenicity and developmental toxicity predictions further illustrate that the two QSAR endpoints do not necessarily vary in parallel across the modeled transformation products. Overall, the QSAR analysis predicts that some ozonolysis products exhibit higher predicted indicators of Ames mutagenicity and/or developmental toxicity than the parent 2,4-DNP. These compound-level predictions describe changes in modeled toxicity indicators and should not be interpreted as demonstrating increased environmental or health risk, which would additionally require information on product yields, environmental concentrations, exposure, and experimentally validated toxicity.

To further assess the reliability of the QSAR predictions, an independent cross-validation was performed using Toxtree (version 3.1.0) (Patlewicz et al., 2008). For Ames mutagenicity, Toxtree's Benigni/Bossa rule-base was employed, which uses rule-based structural alerts that are algorithmically distinct from the statistical consensus method used in T.E.S.T (Benigni et al., 2007). The cross-validation results (Table S5) show consistent predictions for the majority of compounds, with deviations observed only for compounds with probability scores close to the classification threshold. For developmental toxicity, Toxtree does not provide a dedicated prediction model. Instead, we employed its Cramer decision tree to perform a toxicological concern threshold (TTC)-based classification. The Cramer decision tree classifies compounds into three classes based on their chemical structure: Class I (low toxicity concern), Class II (moderate toxicity concern), and Class III (high toxicity concern) (Cramer et al., 1976). Class III compounds are assigned the lowest TTC value (90 µg per person per day) (Munro et al., 1996), indicating the highest potential for adverse health effects, including developmental toxicity. As shown in Table S5, all 2,4-DNP and its ozonolysis products were classified as Cramer Class III, consistent with the T.E.S.T. predictions that all compounds exhibit developmental toxicity potential. This approach follows established practices in computational toxicology for the developmental toxicity screening of nitro-organic compounds (Yong-Di et al., 2022).

The toxicity predictions presented here are screening level assessments based on QSAR models. They are intended for hazard prioritization and hypothesis generation, not as definitive toxicological conclusions. Several factors limit their predictive reliability. First, QSAR model performance is constrained by the chemical space of training sets and the defined applicability domains. For structurally distinct species such as 2,4-DNP and its various derivatives, the predicted confidence levels are lower because these models cannot fully elucidate specific toxicological modes of action. As highlighted by the Ames/QSAR International Challenge Projects, even the best performing models achieve a balanced accuracy of approximately 76.9 % at 100 % coverage (Uesawa, 2024). Second, categorical predictions with probability scores close to the classification threshold are inherently unstable, suggesting that the “not detected” Ames mutagenicity signals for intermediates such as IM17, IM25, IM44, and IM54 may reflect model limitations rather than a definitive absence of toxicity. Furthermore, these computational frameworks necessarily simplify complex biological endpoints such as developmental toxicity. While a quantitative weight of evidence (qWoE) approach has been shown to achieve accuracy comparable to experimental Ames tests (Tintó-Moliner and Martin, 2020), such approaches are beyond the scope of this screening level study. Overall, the findings should be interpreted as hypothesis-generating (vom Eyser et al., 2013). Experimental validation through in vitro Ames tests and zebrafish embryo development toxicity assays would be valuable for evaluating the QSAR-predicted toxicity indicators of selected transformation products.

4 Conclusions

This study systematically investigated selected chemical transformation pathways, molecular-scale interfacial behavior, and predicted toxicity changes of 2,4-DNP using an integrated computational framework combining DFT, MD simulations, and QSAR models. Among the direct gas-phase reactions examined here, the calculated rate coefficient for the •OH pathway is higher than those for O₃ and •NO₂, while the aqueous-phase calculations indicate appreciable intrinsic reactivity toward O₃; the direct reaction of •NO₂ with the parent compound is kinetically unfavorable (gas-phase rate constant on the order of 10⁻²⁸ cm³ molec.⁻¹ s⁻¹), although reactions of •NO₂ with oxidation intermediates (barrier as low as 2.01 kcal mol⁻¹) and multiphase •NO₂ $/$ O₃ chemistry represent additional mechanistic pathways that warrant further investigation; solvation effects and the electron-withdrawing nitro groups substantially affect the calculated reaction kinetics; within the investigated O₂ and •NO₂ reaction systems, HAA exhibits the lowest calculated barrier among the examined secondary pathways; MD simulations show a non-monotonic hydration dependence of 2,4-DNP adsorption in the idealized C60/water system, with changes in hydration altering the molecular adsorption configurations and interaction strengths. These MD results represent molecular-scale behavior within the idealized model system and do not provide quantitative constraints on atmospheric residence time, gas-particle partitioning, or long-range transport. QSAR predictions indicate that certain ozonolysis products exhibit higher predicted Ames mutagenicity and developmental toxicity than the parent compound. These results represent compound-level, screening-level predictions and do not establish increased environmental risk, which would require information on product yields, exposure, and experimentally validated toxicity. While the quantitative conclusions apply specifically to 2,4-DNP, the mechanistic relationships examined here among substituent electronic effects, reaction kinetics, interfacial behavior, and predicted toxicity indicators may provide a basis for future investigations of related nitroaromatic compounds.

Several limitations should be noted. The aerosol surrogate (C60), while useful for isolating specific non-covalent interactions, does not capture the full chemical complexity of real atmospheric particles, including oxygen-containing functional groups, surface defects, heteroatoms, and transition metals that may influence reactivity and product formation. Consequently, the present model may emphasize π-π interactions while underrepresenting polar and chemically heterogeneous surface interactions. The present 1 ns MD simulations are limited in their ability to capture slower dynamical processes, particularly those occurring on more complex aerosol surfaces; longer simulations and more advanced force-fields treatments would therefore be valuable in future studies. The QSAR toxicity predictions are screening-level estimates requiring experimental validation. Moreover, quantitative determination of the relative contributions of the investigated pathways to the atmospheric loss of 2,4-DNP would require atmospheric chemistry modeling that incorporates temporally resolved oxidant concentrations, gas-particle and gas-aqueous partitioning, aerosol liquid water content, and additional loss processes. Additionally, downstream multiphase aging processes such as oligomerization or metal-catalyzed reactions, which may influence the fate of phenolic compounds in real atmospheric particles, were not addressed in this study and would require complementary experimental and modeling approaches.

Future work should therefore focus on experimental validation of the calculated reaction kinetics and QSAR predictions, more realistic representations of carbonaceous aerosol surfaces, and atmospheric chemistry modeling that explicitly incorporates phase partitioning and relevant loss processes. Within these limitations, the present results provide molecular-level mechanistic insight that can guide future experimental and modeling studies of nitrophenol multiphase chemistry.