Effect of head group and lipid tail oxidation in the cell membrane revealed through integrated simulations and experiments

We report on multi-level atomistic simulations for the interaction of reactive oxygen species (ROS) with the head groups of the phospholipid bilayer, and the subsequent effect of head group and lipid tail oxidation on the structural and dynamic properties of the cell membrane. Our simulations are validated by experiments using a cold atmospheric plasma as external ROS source. We found that plasma treatment leads to a slight initial rise in membrane rigidity, followed by a strong and persistent increase in fluidity, indicating a drop in lipid order. The latter is also revealed by our simulations. This study is important for cancer treatment by therapies producing (extracellular) ROS, such as plasma treatment. These ROS will interact with the cell membrane, first oxidizing the head groups, followed by the lipid tails. A drop in lipid order might allow them to penetrate into the cell interior (e.g., through pores created due to oxidation of the lipid tails) and cause intracellular oxidative damage, eventually leading to cell death. This work in general elucidates the underlying mechanisms of ROS interaction with the cell membrane at the atomic level.

. Comparison between our calculated data and experimental or other computational data from literature, for pure DOPC and POPC bilayers.

DFT simulations
We also performed DFT calculations to obtain new parameters for the non-reactive force field (for the oxidized/modified bilayers), by carrying out a fitting of some standard potential energy functions to the DFT energy. For this purpose we used the ADF modeling suite package 17, 18 , employing the B3LYP functional 19 with the Slater-type basis set TZ2P 20 . The implemented new parameters to the force field are given in Table  S2. For the calculation of bond and angle (as well as improper dihedral) force constants, we restrained the bond lengths and angles at seven different values and then fitted a harmonic potential function to the energy profile. For the proper dihedral parameters, dihedral angles were restrained at 36 different values from 0 to 360°, and the corresponding dihedral function was fitted to the potential energy (see Table S2). Note that we performed the parametrization (or fitting) procedure only for those bonds or angles or dihedrals which are not included in the force field. Moreover, we used only those Lennard-Jones parameters which are present in the force field.   Fig. S1).

Bonds
Bond b0 (nm) kb (kJ.mol -1 .nm -4 ) functional type (see 21 ) C1-N4 (Fig. S1a) 0.148 5.5×10 6 2 C3-N4 (Fig. S1a) 0.148 5.5×10 6 2 C2-N4 (Fig. S1a) 0.127 2.03×10 7 2 Angles Angle θ0 (deg) kθ (kJ.mol -1 ) functional type (see 21 ) C1-N4-C2 (Fig. S1a) 122.2 1257.42 2 C3-N4-C2 (Fig. S1a) 122.2 1257.42 2 C1-N4-C3 (Fig. S1a) 115.6 1020.91 2 H1-O2-C3 (Fig. S1b) 104. It should be mentioned here that initially we performed the parametrization for both OX1 and OX2 oxidation products, since our reactive (DFTB) MD simulation results showed these two oxidation products to be most abundantly formed. However, the experimental results did not reveal OX1 as a prominent oxidation product, and therefore we did not use the parameters of this oxidation product in our further non-reactive MD simulations. Figure S2 illustrates the reaction energies of the two possible subsequent reactions of the head group oxidized with OX1, with a new OH radical (see also Fig. 2 in the main paper). A new OH radical can either abstract the H atom from the CH2 group (see initial structure in the red curve of Fig. S2), forming a stable double C=C bond and water molecule (see final structure in the red curve) or it can react with the C radical (see initial structure in the black curve), forming an alcohol group (see final structure in the black curve). As is clear from Fig. S2, the calculated reaction energies employing the above described DFT method revealed that the reaction in which the alcohol group is formed (see final structure in the black curve) is energetically slightly more favorable (cf. ΔE2 in Fig. S2). Moreover, both reactions do not have a barrier. Note that to calculate the reaction energies of these two systems we performed a constrained optimization at 25 different points. In other words, we restrained the distance between an O atom of the OH radical and an H atom of the CH2 residue connected to the phosphate group (see initial structure in red curve) and the distance between the O atom of the OH radical and the C radical (see initial configuration in black curve), at 25 different values and calculated the energy profiles.

Reactive (DFTB) MD results
As mentioned above, the reactions of the OH radicals with the head group of the PL were studied using reactive (DFTB) MD simulations. To get some limited statistics we performed 100 MD simulations. According to the simulation results, we observed six reaction mechanisms which are presented in Fig. S3. For clarity, we named these mechanisms as OX1, OX2, …, OX6. As is clear from Fig. S3, all of these reaction mechanisms are initiated by H-abstraction from (different parts of) the head groups, but few of them give rise to further bond breaking. Note that the OX6 reaction mechanism can take place in both lipid tails (beneath the ester groups, see Fig. S3) and the percentage given is the sum of the OX6 occurrences in both tails. A similar summation over the CH3 groups is also applied for OX1. We did not investigate the further reaction of a new OH radical with the modified PL (by OX1-X6) due to the computational cost of the DFTB method, except for OX1, where we performed further DFT calculations (see above) to find out the energetically most favorable reaction product. It turned out that the formation of an alcohol group was most favorable, therefore we assumed the same for the other modified PLs. This is indeed the most logical for these structures.
We studied the longer-term behavior of only OX2 (as well as ALD, see the main text) in our non-reactive MD simulations and we did not study the consequences of the other mechanisms. The reason is that the experimental results suggested only OX2 together with ALD to be the most prominent oxidation products.

Non-reactive MD results
The surface area per lipid and the bilayer thickness are plotted in Fig. S4 as a function of the concentration of the oxidized PLs, for the OX2 reaction mechanism on the POPC and DPPC PLBs. As is clear, similar trends were obtained for the POPC and DPPC bilayers as for the DOPC PLB (see section Nonreactive MD results in the main paper for a detailed explanation). The only difference is that the relative changes (increase and decrease) in the values are more pronounced in the case of DPPC (see Fig. S4c,d). This is probably due to the difference in the melting temperature of the different PLs. At the simulated temperature (i.e., 300 K), DPPC is in the gel phase (Tm ≈ 328 22 ), while POPC and DOPC (Tm ≈ 294 K and 274 K, respectively 22 ) are already in the liquid state (i.e., with more disordered lipids). Detachment of a lipid tail (see Fig. 5c in the main paper) has a stiffening effect on the PLB, which makes that the OX2 oxidized DPPC bilayers remain in the gel phase. Finally, it is clear from Fig. S4c,d that both the area per lipid and the bilayer thickness level off after 50% oxidation with OX2, which indicates that the lipids are already nearly perfectly structured in the DPPC system. Hence, both the area per lipid as well as the bilayer thickness cannot further decrease/increase.  Figure S5: Mirror plot of DOPC control (top) and plasma treated DOPC (bottom). Data generated from direct infusion, electrospray ionization / TripleTOF5600 high resolution mass spectrometry. Data normalized to base peak count (DOPC, ≈ 20,000/s). Only minor differences were detected between treated and untreated DOPC (see text). Peaks with higher m/z belong to alkali and earth alkali adducts of DOPC. In the low m/z range in source fragments and impurities/traces from buffers occur. Figure S6: Mirror plot of DMPC control (top) and plasma treated DMPC (bottom). Data generated from direct infusion, electrospray ionization / TripleTOF5600 high resolution mass spectrometer. Data normalized to base peak count (DMPC, ≈ 40,000/s). Only minor differences were detected between treated and untreated DMPC (see text). Peaks with higher m/z belong to alkali and earth alkali adducts of DMPC. In low m/z range, major peaks belong to doubly charged ions of DMPC + magnesium or calcium.