Unraveling the Transport Properties of RONS across Nitro-Oxidized Membranes

The potential of cold atmospheric plasma (CAP) in biomedical applications has received significant interest, due to its ability to generate reactive oxygen and nitrogen species (RONS). Upon exposure to living cells, CAP triggers alterations in various cellular components, such as the cell membrane. However, the permeation of RONS across nitrated and oxidized membranes remains understudied. To address this gap, we conducted molecular dynamics simulations, to investigate the permeation capabilities of RONS across modified cell membranes. This computational study investigated the translocation processes of less hydrophilic and hydrophilic RONS across the phospholipid bilayer (PLB), with various degrees of oxidation and nitration, and elucidated the impact of RONS on PLB permeability. The simulation results showed that less hydrophilic species, i.e., NO, NO2, N2O4, and O3, have a higher penetration ability through nitro-oxidized PLB compared to hydrophilic RONS, i.e., HNO3, s-cis-HONO, s-trans-HONO, H2O2, HO2, and OH. In particular, nitro-oxidation of PLB, induced by, e.g., cold atmospheric plasma, has minimal impact on the penetration of free energy barriers of less hydrophilic species, while it lowers these barriers for hydrophilic RONS, thereby enhancing their translocation across nitro-oxidized PLB. This research contributes to a better understanding of the translocation abilities of RONS in the field of plasma biomedical applications and highlights the need for further analysis of their role in intracellular signaling pathways.


Introduction
A cold atmospheric plasma (CAP) is a non-equilibrium plasma operating at atmospheric pressure and room temperature [1][2][3]. It is generated by applying a high-voltage a better understanding of the effect of CAP in plasma-cancer treatment and in plasmaagriculture applications.
Here, we investigated the impact of a combination of lipid nitration and oxidation (nitro-oxidation) on the phospholipid membrane and its permeability by translocating RONS across the membrane. For this purpose, we applied molecular dynamics (MD) simulations, together with umbrella sampling (US) to reveal the permeation of RONS through various nitro-oxidized damaged membranes.

Simulation Setup
We conducted MD simulations to investigate the translocation of RONS, namely hydroxyl radical (OH), hydroperoxyl radical (HO 2 ), hydrogen peroxide (H 2 O 2 ), ozone (O 3 ), nitric oxide (NO), nitrogen dioxide (NO 2 ), dinitrogen tetroxide (N 2 O 4 ), nitric acid (HNO 3 ), s-cis-nitrous acid (s-cis-HONO), and s-trans-nitrous acid (s-trans-HONO), across both native and nitro-oxidized membranes. As a model system, we used the phospholipid bilayer (PLB), which served as a representation of the cell membrane. The area per lipid in simulations can vary depending on the conditions, and it has been reported to range from 65 to 70 Å 2 in various studies [40,41]. These model systems typically include 40 to 65 water molecules per lipid, to ensure proper hydration. The experimentally measured area per lipid for DOPC was 67.4 ± 1.0 Å 2 [42], which falls within the range of theoretical calculations. In our model system, we included 62.5 water molecules per lipid, and the calculated area per lipid was found to be 65 Å 2 (cf. below Table 1 and Section 3), which is consistent with the experimental results. Our PLB consisted of 200 phospholipids (PLs) that were equally distributed in both layers of the membrane, in addition to 12,500 water molecules surrounding them from above and below (see Figure 1a). It is important to note that we used a higher number of water molecules per lipid in our model system, to facilitate umbrella sampling simulations. Moreover, inserting six of one of the RONS (e.g., 6 OH radicals) along the z-axis in our system required a larger simulation box. This was necessary to prevent interactions with periodic images. As a result, we included more water molecules, to ensure an appropriate distance between the RONS and the periodic boundary conditions. We considered native PLB consisting of 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) lipid molecules, as well as treated PLB consisting of both oxidized and nitrated forms of DOPC (see Figure 1b). As is clear from figure, the oxidation of DOPC led to the truncation of its fatty acid chain and the formation of aldehyde groups, while the nitration of DOPC led to the formation of a nitro group (-NO 2 ) on one of its fatty acid tails [43]. The modifications possess three variations of the treated PLB, each of which differed in its composition (see Table 1). We selected the above (native and modified) PLs based on the following reasons: (i) DOPC lipid is a fundamental constituent of the plasma membrane in both the extracellular and cytoplasmic leaflets [44], (ii) aldehyde-oxidized lipid (DOPC-ALD) is one of the frequently observed oxidation products [45,46] and the main one responsible for increasing membrane permeability [47][48][49], and (iii) nitro PL (DOPC-NIT) is one of the main nitrated lipids observed during PL nitration [34,35]. It should be noted that the cellular membrane is composed of a variety of membrane constituents, such as PLs, proteins, sterols, and so on. In practice, the computational resources required for simulating the realistic membrane composition are quite expensive, since even the most elementary plasma membrane of a red blood cell is composed of more than 150 lipid species [50]. In this respect, we considered the principal lipid constituents in our simulations. Hence, in our model system, we chose DOPC as the primary component of the PLB. This particular lipid molecule is ubiquitous, being present in the endoplasmic reticulum, mitochondrion membrane, and liver cell plasma membrane at levels of approximately 40%, 44%, and 24%, respectively [50]. The initial configurations of the native and modified PLBs were generated using the Packmol software package [51]. All simulations were carried out using the GPU accelerated version of GROMACS-2021.4 software package [52,53]. We employed the GROMOS (53A6) united atom force field to simulate the behavior of native PLs [54]. To account for the presence of the aldehyde and nitro functional groups, we used the force field parameters developed for oxi-and nitro-lipids [35,55]. Additionally, we obtained the GROMOS-type force field parameters for the RONS investigated in this study from [56][57][58]. The choice of these parameters was based on their demonstrated ability to produce accurate water-to-alkane partition coefficients for ROS and RNS, as these coefficients are crucial for accurate modeling of transmembrane permeation. The force field parameters we utilized were further complemented by the parameters for NO particles that we had previously developed using the same methodologies as described in [56,57].  [54]. To account for the presence of the aldehyde and nitro functional groups, we used the force field parameters developed for oxi-and nitro-lipids [35,55]. Additionally, we obtained the GROMOS-type force field parameters for the RONS investigated in this study from [56][57][58]. The choice of these parameters was based on their demonstrated ability to produce accurate water-to-alkane partition coefficients for ROS and RNS, as these coefficients are crucial for accurate modeling of transmembrane permeation. The force field parameters we utilized were further complemented by the parameters for NO particles that we had previously developed using the same methodologies as described in [56,57].  The model membranes were initially minimized employing the steepest descent algorithm prior to a 100 ps equilibration in the isothermal-isobaric (NPT) ensemble at 1 bar and 310 K. These equilibration runs utilized a semi-isotropic Parrinello-Rahman barostat with a compressibility and coupling constant of 4.5 × 10 −5 bar −1 and 5.0 ps, respectively, along with a Nose-Hoover thermostat with a coupling constant of 0.5 ps. A cut-off radius of 1.0 nm was used for non-bonded (Lennard-Jones) and electrostatic (Coulomb) interactions. The particle mesh Ewald method was employed to handle electrostatic interactions, with a real space cut-off of 1.0 nm, coupled with a 0.16 nm spaced-grid for the reciprocalspace interactions. All simulations were accomplished using a time step of 2 fs, and periodic boundary conditions were implemented in all three dimensions. Furthermore, the The model membranes were initially minimized employing the steepest descent algorithm prior to a 100 ps equilibration in the isothermal-isobaric (NPT) ensemble at 1 bar and 310 K. These equilibration runs utilized a semi-isotropic Parrinello-Rahman barostat with a compressibility and coupling constant of 4.5 × 10 −5 bar −1 and 5.0 ps, respectively, along with a Nose-Hoover thermostat with a coupling constant of 0.5 ps. A cut-off radius of 1.0 nm was used for non-bonded (Lennard-Jones) and electrostatic (Coulomb) interactions. The particle mesh Ewald method was employed to handle electrostatic interactions, with a real space cut-off of 1.0 nm, coupled with a 0.16 nm spaced-grid for the reciprocalspace interactions. All simulations were accomplished using a time step of 2 fs, and periodic boundary conditions were implemented in all three dimensions. Furthermore, the production run was performed for 500 ns which was a sufficiently long to achieve an equilibrated state for our four model membranes.
To compute the free energy profiles (FEPs) associated with the translocation of ROS and RNS across the membranes, we utilized the US method [59,60]. The initial membrane structures for US simulations were chosen from the final 30 ns of the equilibration trajectory (of the production run) by extracting frames with 10 ns intervals (i.e., 470, 480, 490, and 500 ns). Thus, four replicas were used in our US simulations. In order to efficiently use the computational resources, the US simulations were carried out as follows: as mentioned above, four membrane structures were selected from the last 30 ns and, in each of the US simulations, six permeants were positioned along the membrane normal, while maintaining a distance of 1.1 nm between them (see Figure 2a). Moreover, to gain more statistics, we inserted 18 more permeants of the above species, varying their position in the xy-plane. These species were located sufficiently far from each other, i.e., beyond the cut-off radius, and they were allowed to move in the xy-plane, within the range of the applied flatbottomed position restraint with a force constants of 500 kJ·mol −1 nm −2 and a radius of 0.5 nm (see Figure 2b). To clarify, the inserted species in our simulations did not interact with each other, as they were located beyond the specified cut-off radii. Despite traveling along the edge of the flat bottom restrained potential, the minimum distance between them in the xy-plane was approximately 3 nm (see Figure 2b). This ensured that the species remained sufficiently separated and prevented any unintended interactions during the simulation. To maintain the positions of permeants relative to the membrane center, their center-of-mass motion was also restrained along the z-axis using the harmonic bias with a force constant of 2000 kJ·mol −1 nm −2 . The permeants were allowed to move freely in the xy-plane. Each system underwent a US simulation at NPT, consisting of 2 ns of equilibration followed by 4 ns of sampling (i.e., a total of 6 ns US simulation). The histograms obtained from the US simulations were collected and analyzed using the weighted histogram analysis method to calculate the associated FEP for each system [61]. We conducted eleven individual US simulations, in order to construct a single FEP. As a result, (6 × 11) × 4 = 264 US data points were collected along the membrane normal, each separated by 1.1 Å, for four FEP (i.e., 66 for each FEP). Thus, we obtained 16 individual FEPs (i.e., 4 permeants in each US run × 4 replicas) for each RONS (i.e., OH, HO 2 , H 2 O 2 , O 3 , NO, NO 2 , N 2 O 4 HNO 3 , s-cis-HONO, and s-trans-HONO). In this way, the final FEP of each particle was obtained by averaging 16 individual FEPs. Since we used ten RONS and four PLB systems (see Table 1), a total of 11 US runs × 10 RONS × 4 PLB structures × 4 replicas = 1760 US simulations were performed, each lasting for 6 ns. Hence, a total of 10.56 µs of US simulations were carried out in this study. It is important to note that the data obtained for the FEP calculation heavily relied on the lipid composition surrounding the randomly placed RONS in the PLB. In order to obtain accurate free energy profiles (and corresponding barriers), a good statistical analysis is required. Additionally, to evaluate the convergence of the free energy results, it was necessary to examine the symmetry of the FEPs of each RONS to the center of the PLB. Our simulation results revealed that the full FEPs obtained in this study were very similar to the symmetrical ones, indicating the convergence of the simulations (see Appendix A, Figure A1 and cf the FEPs of HO 2 (hydrophilic) and NO (less hydrophilic) across the nit25-ox25 system). It is worth noting that in reality, chemical reactions involving RONS may occur in the PLB. However, conventional non-reactive MD simulations cannot fully capture these processes, due to the limitations of the potential used. The electronic degrees of freedom required to describe chemical reactions are not explicitly taken into account in classical MD simulations. Owing to these limitations, we simulated membranes after oxidation/nitration. Nevertheless, US simulations can provide valuable information on the RONS permeation rate across the PLB before and after modification, as well as the most probable accumulation regions for RONS. These results allow for the study of position-dependent specific interactions in the membrane, which may have implications for biomedical applications.

Results and Discussion
In this study, our goal was to investigate the mechanisms of the permeation of RONS, i.e., OH, HO2, H2O2, O3, NO, NO2, N2O4, HNO3, s-cis-HONO, and s-trans-HONO through native and modified PLBs. Figure 3 shows the FEPs of RONS across both native and modified PLBs.

Results and Discussion
In this study, our goal was to investigate the mechanisms of the permeation of RONS, i.e., OH, HO 2 , H 2 O 2 , O 3 , NO, NO 2 , N 2 O 4 , HNO 3 , s-cis-HONO, and s-trans-HONO through native and modified PLBs. Figure 3 shows the FEPs of RONS across both native and modified PLBs.    The associated free energy barriers (∆G) are given in Table 2. In reference to hydrophilic species, each energy barrier was quantified using the difference of the minimum and maximum of the FEPs. Conversely, for the case of less hydrophilic species [62][63][64], the free energy barrier was observed within the headgroup region of PLB. When the hydrophilic RONS (i.e., HNO 3 , s-cis-HONO, s-trans-HONO, H 2 O 2 , HO 2 , and OH) are transferred from the aqueous phase to the interior of the membrane, the FEP initially decreases, attaining a minimum at the headgroup region, and subsequently rises as the RONS are become positioned deeper into the native PLB (see Figure 3a). As is clear, the free energy barrier was obtained at the core of the PLB in all cases. Adsorption at the headgroup region was considerably more robust for H 2 O 2 and HO 2 , with H 2 O 2 displaying the highest permeation barrier among all hydrophilic RONS. A previous simulation study found specific interactions between ROS and the membrane, which could be responsible for these trends [56]. It has been postulated that the adsorption of RONS onto the surface of PLBs is principally driven by two factors: (i) hydrogen bonding interactions between RONS and the carbonyl ester groups of the PL, and (ii) dispersion interactions with the headgroup region. Although all the hydrophilic RONS investigated in this study are capable of acting as hydrogen bond donors, H 2 O 2 and HO 2 possess an additional oxygen atom, which results in stronger dispersion interactions in comparison to the smaller OH radical. This may explain the differences in the adsorption tendencies of the hydrophilic RONS. Additionally, it is known that H 2 O 2 has the ability to form twice as many hydrogen bonds in water as compared to OH and HO 2 [56]. This may elucidate the reason behind the substantial permeation barrier observed for H 2 O 2 . It is noteworthy how fundamental principles of chemical bonding can be employed to elucidate distinctions in the comportment of s-ciss-trans isomeric forms of nitrous acid. In the case of s-trans-HNO 2 , the alignment of discrete bond dipoles occurs in a harmonious manner, fostering a propitious interplay among them. Such felicitously aligned bond dipoles leads to the formation of a discernible molecular dipole, ultimately resulting in favorable hydration. On the contrary, the s-cis conformation gives rise to an adverse alignment of bond dipoles, rendering it less polar when juxtaposed with its s-trans counterpart. Our straightforward qualitative explanation coheres with experimentally determined dipole moments: 1.85 D [65] and 1.42 D [66] for s-trans and s-cis isomers, respectively. Furthermore, the dissimilarities in the chemical structures of these isomers also contributed to the contrasting topological polar surfaces exhibited by them, with the s-trans form boasting a larger polar surface in comparison to the s-cis form. These electronic and structural disparities were evinced in the ∆G of these isomers: the s-trans isomer underwent more favorable hydration and exhibited a larger ∆G in the headgroup region, while the s-cis isomer experienced less favorable hydration and confronted a diminished ∆G barrier in the headgroup region (see Table 2 and Figure 3a-d). In the event that the PLB underwent nitro-oxidation, the permeation pattern exhibited similar qualitative tendencies, albeit with comparatively less free energy barriers (see Figure 3b-d). The oxidation of PL produced functional groups and fragments that enhanced the hydrophilicity of the membrane core, thereby increasing the permeability of the PLB to RONS [67]. Whilst the permeation barrier remained comparatively high for H 2 O 2 , it became significantly reduced for OH in all cases of nitro-oxidized PLB. Nevertheless, the nitro-oxidation of the PLB reduced the translocation free energy barriers of all hydrophilic species compared to native PLB. Notably, with regards to HO 2 , nitro-oxidation of the PLB apparently induced a modification in its partition behavior. The FEPs showed that HO 2 may manifest a preference for the PLB core, as nitro-oxidation occurs. This is a particularly important effect given that HO 2 is the protonated form of the biologically significant superoxide radical (O − 2 ). These radicals play a pivotal role in biological systems, acting as essential mediators in redox signaling pathways. Their controlled production and regulation are crucial for maintaining cellular homeostasis, as they participate in various physiological processes, such as cell signaling, immune response modulation, and regulation of vascular tone. Nevertheless, we would like to emphasize that this conclusion is still hypothetical. The HO 2 model was initially parametrized to replicate its experimentally determined hydration free energy [56]. The consideration of less hydrophilic solvation was precluded, due to a lack of available experimental data. To draw a more definitive conclusion, it would be necessary to assess the HO 2 model's ability to describe its solvation free energy in non-polar solvents.
Overall, the less hydrophilic RONS exhibit substantially lower permeation barriers than their hydrophilic counterparts (see Figure 3e-h). They do not participate in hydrogen bonding interactions in the aqueous layer, which facilitates their diffusion from the aqueous phase to the core of the PLB. The maximum energy of all less hydrophilic RONS is located near the water-lipid interface, specifically in the headgroup region, which is mainly attributed to the polarity of the headgroup region. The values of maximum energy for all less hydrophilic RONS in all nitro-oxidized PLB systems are similar, having a slightly higher barrier in the case of the native PLB (Figure 3e-h). Additionally, the FEPs of the less hydrophilic RONS exhibit their minima at the center of the PLB, indicating the preference of these species (especially, NO, NO 2 , and O 3 ) to accumulate in the PLB center. The FEPs imply that NO, being virtually nonpolar, has a high tendency to accumulate inside the PLB, and triatomic species such as NO 2 and O 3 also represents similar FEPs, which conserve a residual dipole moment [57]. N 2 O 4 exhibits a behavior distinctive from other hydrophilic species, as it does not encounter any energy barrier in the headgroup region; rather, a minimum free energy is observed in close proximity to this area. When it penetrates deeper into the PLB core, the free energy gradually rises, but it is still lower than that in bulk water. The barrier near the center of the native PLB (i.e., at a distance of about 0.7 nm) is caused by a double bond present in the oleoyl tails, which makes the lipid chains bent. This barrier also represents similar changes for the other less hydrophilic species (see NO 2 and O 3 in Figure 3e). These findings imply that N 2 O 4 has strong interaction with the membrane surface but partitions almost equally between the aqueous phase and the PLB core. The FEP of N 2 O 4 exhibits characteristics of both hydrophilic and less hydrophilic species. The geometry of N 2 O 4 is such that it has no dipole moment, but the partial charges of N and O atoms are relatively significant (+0.584e and −0.292e, respectively) [57]. It is possible for N 2 O 4 to exhibit a strong quadrupole moment, which may explain its tendency to accumulate around the polar headgroup region, despite being less hydrated in this area compared to the aqueous phase. Despite undergoing nitro-oxidation, the permeation FEPs of less hydrophilic RONS remain largely unaffected, except for their smoothing and the near disappearance of the barrier near the headgroup (see Figure 3e-h). The local free energy barriers present at approximately |0.7| nm also disappears (see Figure 3f-h). This phenomenon can be attributed to the cleavage of the lipid tails (in the case of oxidation) and the consequent formation of aldehyde groups (see Figure 1b). The elimination of the localized free energy barriers can be attributed to the fluidity of the membrane. This fluidity is caused by an increase in the surface area of the PLB (area per lipid) and a decrease in the membrane thickness due to disordered lipid tails (see Table 3) [67]. Nevertheless, less hydrophilic species, except for N 2 O 4 , still tend to accumulate at the center of the PLB. The findings from our simulations are in accordance with previous experimental studies on PLB permeability. Specifically, our simulation results demonstrate that less hydrophilic species such as NO exhibit a permeability that is 3-6 orders of magnitude higher than that of hydrophilic RONS such as H 2 O 2 [68,69]. Thus, the incorporation of membrane-embedded aquaporin channels or pores induced by powerful electric fields [67] becomes crucial, to enhance the transport of hydrophilic RONS into and out of the cell. In fact, our recent findings indicate that the permeability of H 2 O 2 through AQP is two orders of magnitude greater than through the PLB [32]. Indeed, for less hydrophilic RONS, transmembrane permeation can occur, even without the presence of channels and pores. The literature has also demonstrated this, where a higher permeability coefficient was reported [70], ranging from 18 to 93 cm s −1 for NO [69,71,72], 5 cm s −1 for NO 2 [72], and 12 to 157 cm s −1 for O 2 [69,71,[73][74][75], whereas for H 2 O 2 , the permeability coefficient varies between 4 × 10 −4 and 1.2 × 10 −2 cm s −1 [70].
As previously discussed, the cell membrane is inherently intricate and composed of various membrane constituents. As a result, these constituents may also impact the penetration of RONS into the cytoplasm. It is evident that examining the impact of all these constituents is unfeasible. In a previous investigation, the impact of cholesterol on the FEPs of ROS (i.e., OH, HO, H 2 O 2 and O 2 ) was explored [76]. According to this study, the presence of cholesterol raised the lipid order, leading to an increase in the free energy barrier for ROS penetration. Thus, the membrane system which included cholesterol may influence the FEPs of the studied RONS, by raising their free energy barriers.

Conclusions
This research is of significant importance for the field of plasma medicine and marks a significant stride towards comprehending the diverse translocation abilities of various RONS (produced by CAP) across the native and nitro-oxidized cell membrane.
Overall, we may infer that less hydrophilic species, namely NO, NO 2 , N 2 O 4 , and O 3 , exhibit a higher penetration ability across both native and nitro-oxidized PLBs, in comparison to their hydrophilic counterparts, such as HNO 3 , s-cis-HONO, s-trans-HONO, H 2 O 2 , HO 2 , and OH. The nitro-oxidation of the PLB induced by processes such as CAP does not have a significant impact on the free energy barriers of less hydrophilic species. As is clear, the free energy barrier is slightly higher when there is 25% nitrated lipids + 25% oxidized lipids (especially for HNO 3 , s-trans-HONO, and HO 2 ). This result is in agreement with the previous results of the study, where a mixture between 50% oxidized lipids and 50% nitrated lipids presented a higher free energy barrier for water permeation than 100% nitrated membranes [38]. However, this reduces the barriers of hydrophilic RONS, thereby increasing their probability of translocation across the nitro-oxidized PLB.
These computational studies shed light on the mechanisms underlying the action of RONS on native and modified membranes and provide valuable information about processes at the molecular level. Further detailed analyses are necessary, to unravel the synergistic function of these species in intracellular signaling pathways. This, in turn, will help in better understanding the behavior of biological systems, as well as provide a basis for future research and the development of new treatment options and targeted technologies.  100% nitrated membranes [38]. However, this reduces the barriers of hydrophilic RONS, thereby increasing their probability of translocation across the nitro-oxidized PLB. These computational studies shed light on the mechanisms underlying the action of RONS on native and modified membranes and provide valuable information about processes at the molecular level. Further detailed analyses are necessary, to unravel the synergistic function of these species in intracellular signaling pathways. This, in turn, will help in better understanding the behavior of biological systems, as well as provide a basis for future research and the development of new treatment options and targeted technologies.