Molecular Dynamics Simulations of Plasma–Antifolate Drug Synergy in Cancer Therapy

2. Methods

2.1. Protein Modeling

The original hSLC19A1 conformation interacting with 5-MTHF and PMX was obtained from the RSCB Protein Data Bank (8GOE, 8GOF). hSLC19A1 consists of more than 400 residues and has a predominantly α-helical structure, with a distinct substrate-binding pocket in the centre spanning both ends of the phospholipid bilayer, as shown in Figure 1(c). In order to obtain the oxidized hSLC19A1 model, oxidative modification of specific residues was applied on VIENNA-PTM2.0[28,29,30]. It was reported that the sulphur-containing amino acids Met and Cys are highly reactive and susceptible to plasma oxidation, with Met being oxidized to (Met+O)H+ and Cys being oxidized to (Cys+3O)H+ in the presence of reactive particles, such as hydroxyl groups and hydrogen peroxide[31,32,33]. In order to investigate the specific oxidative modification sites of the plasma active particles, a solvent accessible surface area (SASA) analysis was carried out, residues with larger SASA have a higher probability of contact with the plasma and a higher probability of being oxidized. The results of the analyses are shown in Table 1, and it is clear that the SASA of Met119, 122 and Cys365, 396 is much lower than that of Met38, 130, 254 and Cys30, 33, suggesting that the former are buried inside the protein and have little access to the active particles. Therefore, modification of Cys30 and Cys33, Met38, 130, 254 residues should be considered as they have the highest probability of being oxidized.

2.2. Simulation Details

To investigate the effect of plasma on the transport of folate as well as pemetrexed by hSLC19A1 a total of four systems were created, two native systems as well as the oxidized proteins transporting 5-MTHF and pemetrexed, using 19A1 as an abbreviation for hSLC19A1, named 19A1-5MTHF, 19A1ox-5MTHF, 19A1- PMX, and 19A1ox-PMX, as shown in Table 2. The topologies of the ligand molecules were generated in the Automated Topology Builder[34,35] and introduced into the overall topology file. The protein was embedded in a phospholipid bilayer composed of 128 POPC molecules using the membed method, so that the protein was first scaled down to 0.1 times its original size in the xy-direction, and then gradually swelled to its normal size during the simulation in 1000 steps, so that the phospholipid was in close contact with the protein. A cubic box with dimensions of 7.15*7.15*8.95 nm³ was solvated by adding water molecules and NaCl, followed by energy-minimization using the conjugate gradient algorithm, positional constraints were then imposed on the hSLC19A1 backbone and ligand weight atoms, and constrained molecular dynamics simulations were carried out using a V-rescale thermal bath and a C-rescale pressure bath in the NPT ensemble, allowing the system to gradually reach equilibrium without strongly disturbing its original conformation. Subsequently, positional restraints were removed, the temperature was maintained at 298.15 K in order to facilitate comparison with the experimental binding data obtained afterwards, and the pressure was maintained at 1 atmosphere for a final equilibrium simulation of 100 ns. Electrostatic interactions were treated using the PME method. Periodic boundary conditions were applied in all three dimensions. The simulated time step is 2 fs and the LINCS algorithm is applied to maintain the rigidity of the molecular bonds. TIP3P water model was used to solvate all simulated systems. All simulations were performed using the GROMACS 2023.2 software package with the GROMOS54A7 force field containing non-standard residue descriptions provided on the PTM website, and a description of the phospholipid molecules was added to the force field. All molecular visualisations were done using VMD[36].

2.3. Free Energy Calculations

In order to calculate the effect of plasma oxidation on the binding energy between hSLC19A1 and the ligand and to analyse the contribution made by each residue to the binding energy, we also calculated the binding energy of all systems using the molecular mechanics Poisson-Boltzmann surface area(MMPBSA) method[37]. The MMPBSA method integrates three computational methods, the molecular mechanics force field (MM), the Poisson-Boltzmann equation, and the solvent-accessible surface area (SASA), is widely used for the calculation of interaction energies between biomolecules, which allows for a comparison of the relative energetic contribution of each residue to the total binding energy. Umbrella sampling was also conducted in order to show more precisely the effect of plasma oxidation on the binding of 5-MTHF and hSLC19A1. The protein and ligand were placed in a rectangular box elongated by 7.0 nm in the x-axis direction for pulling, and the constructed model was equilibrated in an NPT ensemble for 20 ns for Steered molecular dynamics simulations, where the ligand was slowly pulled along the x-axis at a rate of 0.01 nm/ps. Following the stretching simulations, a series of windows were extracted along the x-axis, followed by batch-independent simulations for each window to obtain a single free energy of dissociation, and finally the Weighted Histogram Analysis Method (WHAM) was applied to calculate the binding free energy of 5-MTHF/PMX dissociated from the hSLC19A1/hSLC19A1ox binding site.

3. Results and Discussion

3.1. Changes in Protein Structure

In this section, the structural changes of the protein after oxidation by plasma are analyzed. Firstly, the transmembrane channel radii of native and plasma-oxidized hSLC19A1 were calculated after equilibrium, as shown in Figure 2(a) and (b), to explain the effect of plasma oxidized modifications on the hSLC19A1 monomer structure and transmembrane transport function[38].

hSLC19A1 consists of 12 transmembrane domains, and each transmembrane domain is abbreviated as TM. In the green region of Figure 2(c), TM 1, TM 2, TM 7 and TM 8 are tightly coupled with each other, creating the narrowest constriction of the monomeric pore from the top, with an average radius of around 0.6 Å before oxidation, which is reduced to 0.5 Å in the 19A1ox-5MTHF system. The blue region is the main pocket for the binding of hSLC19A1 and ligand, with a total length of about 35 Å. At the end of the channel, TM 4, TM 5, as well as TM 10 and TM 11 dissociate to form the entrance to a large polar cavity within the transporter. The entire monomeric channel extends about 5.5 nm in length, longer than the thickness of the phospholipid bilayer (about 5 nm), and is therefore capable of transmembrane transport.

The red curve in Figure 2 (a) illustrates the channel radius of hSLC19A1_OX-5MTHF along the z-axis. Compared with the native system, the plasma-oxidized hSLC19A1 channel narrows significantly upon translocation of 5-MTHF, with the average channel radius narrowing from 2.384 Å in the native system to 2.173 Å. The most significant change in the radius of the substrate-binding pocket was observed between 3 nm and 6.2 nm along the z-axis, with the pocket length decreasing from 3.2 nm to approximately 2.8 nm. 5.9 nm from the z axis at the channel inlet radius decreased from 1.3 nm to 0.56 nm, suggesting that the outward conformation of hSLC19A1 tends to be closed, limiting substrate entry. Met130 is located at the bottom of the pocket where hSLC19A1 binds to the glutamate portion of 5-MTHF/PMX, and the distance between it and Arg373 on the channel’s opposite side, is able to reflect the radius of the protein’s channel at 5-6 nm. Figure 2(d) shows the variation of the distance between Met/Msx130 and Arg373 for the four systems; it can be seen that the distance between Msx130 and Arg373 in the 19A1ox-5MTHF system is significantly shortened, resulting in a reduction of the channel radius. The narrow channel would hinder folate uptake by cancer cells and may also adversely affect the binding affinity of the protein to folate. In contrast, there was no significant change in the radius of the pemetrexed substrate-binding pocket in Figure 2 (c), and the distance between Msx130 and Arg373 was not shortened in the 19A1ox-PMX system, suggesting that the oxidation of hSLC19A1 by CAP did not structurally affect pemetrexed transport while hindering folate uptake into the cancer cells, offering a mechanistic basis for plasma-conjugated antifolate therapies targeting cancer cells.

Secondary structure was analyzed to obtain the evolution in the secondary structure of each residue of hSLC19A1 and the proportion of the different secondary structures in it during the simulation of the four systems, as shown in Figure 3. A-Helix is the most common secondary structure and the main component of hSLC19A1. In 19A1-5MTHF, 19A1-PMX and 19A1ox-PMX, the content of A-Helix was stable at about 73%, indicating that the transport structure of the protein was not significantly affected. In contrast, the percentage of A-Helix in the 19A1ox-5MTHF system decreased by about 3% after 20 ns, suggesting that the oxidized hSLC19A1 may have a reduced ability to transport 5-MTHF. It was observed that residues near Cse30, Cse33, and Met38 in the 19A1_ox-5MTHF system tended to shift to T-turn, the appearance of which, as a short-range folding unit, may disrupt the continuity of the A-Helix, thereby affecting the aperture size of the transmembrane channel as well as the conformation of the substrate-binding pocket.

Futhermore, residues near Met130 shifted from A-Helix shifts to S bend and I 5-helix, suggesting a possible enhancement of flexibility in this region, which arises from plasma oxidation of the corresponding residues. Overall, the oxidative modification specificity of the plasma affects the helix stability of the 5-MTHF binding region.

Root mean square deviation (RMSD) is a widely used measure of protein structural deviation, and the system can be considered to have reached a steady state when the overall RMSD value stabilizes within a narrow range of values over time. Figure 4 (a) shows the variation of the protein backbone RMSD values of the four systems over the 100 ns simulation. The native system 19A1-5MTHF reached equilibrium quickly at the beginning of the simulation, indicating that the simulation parameters were reasonable, whereas after being oxidized by the plasma, the RMSD of 19A1ox-5MTHF was significantly higher than that of the native system, indicating that the protein structure was more unstable, which led to the weakening of the binding affinity to 5-MTHF. In contrast, the RMSD of 19A1_ox-PMX, although showing a rapid rise at 10 ns, stabilised to near 0.2 after 17 ns, and was lower than that of the native system at the later stage of the simulation. This suggests that the presence of PMX keeps the structure of the oxidized hSLC19A1 stable and is not overly disturbed by the oxidative modifications of the plasma.

The Cα root mean square fluctuation (RMSF) quantifies the atomic-level conformational flexibility of each residue in a molecular dynamics (MD) trajectory, and the part of the region with higher flexibility has more fluctuating protein motion and has larger RMSF values; while the region with lower flexibility is more stable. Figure 4 (b) shows the RMSF comparison between the oxidized system and the native system. From the figure, it can be seen that the RMSF of the protein in 19A1ox-5MTHF is obviously increased at residues Cys30 and 33. Cys30 and 33 are located on transmembrane domain 1a. Simulated snapshots in Figure 5 reveals that the transmembrane domain 1a of 19A1_ox is inwardly tightened, and the shift of the secondary structure of Cys33 from A-Helix to T-turn hinders the entry of folate into the binding pocket, which is in agreement with the results of the secondary structure analysis. In contrast, the change of protein flexibility in the 19A1_ox-PMX system mainly occurred near the Met130 residue, which is located on TM4, and it can be seen in Figure 5 that the increase of its flexibility led to the outward opening of the TM4 end, which was also the reason for the slight increase of the radius at the entrance of the channel, and the PMX-binding region remained structurally stable without observable conformational perturbations.

Overall, the oxidative modification of CAP narrows the radius of the channel of hSLC19A1 upon transporting 5-MTHF, disrupts its secondary structure and alters the conformational flexibility of critical residues, which impede 5-MTHF entry into cancer cells, thereby impairing one-carbon metabolism, and weakening their antioxidant defense systems. Meanwhile, pemetrexed (PMX) exhibits strong affinity binding to hSLC19A1, maintaining structural stability even under plasma-induced oxidative conditions, indicating that the oxidative modification of hSLC19A1 by plasma did not structurally affect PMX transport, and it could synergize with pemetrexed for antifolate therapy, which proposed a new strategy for cancer treatment.

3.2. Analysis of Changes in Protein-Ligand Interactions

In order to further investigate the effect of plasma oxidation of hSLC19A1 on its transport function, it is necessary to analyse the changes in hydrogen bonds and hydrophobic interactions and binding energies between 19A1 and 5-MTHF as well as PMX before and after plasma modification, for this purpose hydrogen bonds and hydrophobic interactions were analyzed. The binding free energy between 19A1 and the ligand was also calculated using the MMPBSA method as well as umbrella sampling, and molecular docking was also carried out to assess the binding affinity.

The phenomenon of hydrophobic groups aggregating close to each other to avoid water is known as the hydrophobic interaction and is the main driver of protein folding. In this paper, we used ligplot[39,40] to study the formation of hydrogen bonds and hydrophobic interactions between proteins and ligands after 100 ns simulation, as shown in Figure 6, and the probability and position of the appearance of each hydrogen bond within 100 ns were counted using VMD in Figure 7 and Table 3. The donor-acceptor distance was set to be 3.5 Å and the angular cut-off was 43.1°, corresponding to 30° for the GROMACS calculation of hydrogen bonding. As can be seen in Figure 6, in the native system, the pterin portion of 5-MTHF occupies the uppermost negatively charged pocket of the channel, forming hydrophobic contacts with Ile48 of TM 1, Leu72 of TM 2, and Tyr 126 of TM 4. The glutamate portion faces the inside of the cell and is close to the positively charged region of the protein. The γ-carboxylic acid group interacts with Arg133 of TM 4 and Arg373 of TM 10, the β-carboxylic acid group interacts with Tyr281/Tyr282 of TM 7 and Gln377 of TM 10, and the carbonyl group interacts with Arg133 of TM 4. The benzoyl group partially bridges the pterin- and glutamate-mediated interactions and forms hydrophobic interactions with Tyr126/Met130 of TM 4 and Val285/Tyr286 of TM 7. Upon oxidation, the hydrogen bonds formed by 5-MTHF with Tyr281 as well as the hydrophobic interactions were reduced, and Table 3 shows that the number of hydrogen bonds formed per frame in the simulations decreased from 0.9895 to 0.3628, and the number of hydrogen bonds formed with Arg133 decreased from 0.2724 to 0.0705, and the reduction of interactions with these two important residues is an important reason for the decrease in the binding of 19A1 to 5-MTHF. In the oxidized protein, the glutamate portion of 5-MTHF only formed hydrophobic interactions with two residues, Arg133 and Arg157, indicating that the binding ability of the glutamate portion of 5-MTHF to the protein was weakened.

The pyrimidine-pyrrole, benzoyl and glutamate portions of PMX form contacts with hSLC19A1 similar to those of 5-MTHF, highlighting a conserved substrate recognition mechanism for folate and antifolate substrates. However, in addition to these conserved functions, 5-MTHF and PMX still exhibit significant differences in their detailed conformations. Specifically, the planes of the pterin ring in 5-MTHF and the pyrimidine-pyrrole ring in PMX exhibit relative rotation and displacement, which in turn induces adjustments in the benzoyl and glutamate portions, leading to differences in hSLC19A1 binding to 5-MTHF and PMX before and after oxidation. After oxidation, the formation of hydrogen bonds between 19A1_OX and Tyr282, Tyr376, and Gln377 increased, with an average of one more hydrogen bond formed per frame, which made the binding of pemetrexed to hSLC19A1 somewhat more robust.

Molecular docking was performed based on the vina algorithm, as shown in Figure 8[41]. Vina score reflects the strength of affinity between the protein and the ligand, and at below -9 indicates that the two are able to bind stably. The vina scores of both native systems were below -9, indicating that both 5-MTHF and PMX are ligands with strong binding ability to hSLC19A1. PMX had a vina score as high as -9.5, which was stronger than 5-MTHF, and thus PMX was preferentially transported by hSLC19A1. Upon plasma oxidation, the predicted binding affinity of 19A1ox-5MTHF was reduced to -9.2, whereas the score of 19A1ox-PMX increased to -9.6, which is consistent with the results of the protein structural analyses and protein/ligand interaction analyses, indicating that the affinity of hSLC19A1 for PMX was significantly higher than that of 5-MTHF after plasma oxidation.

To further investigate the binding capacity between hSLC19A1 and 5-MTHF and pemetrexed under plasma action, we calculated the binding free energies broken down to individual residues using the molecular mechanics-Poisson Boltzmann surface area (MM-PBSA) method, which is widely used in the calculation of binding free energies. The overall protein-ligand binding free energy contributions to each branch are shown in Table 4. In the table, dG is the binding free energy, the MM energy is mainly the internal interaction energy, which contains two items, the Coulomb interaction (COU) and the van der Waals interaction (VDW), and the PBSA energy is used to describe the solvation free energy between the protein and the ligand, i.e., their contribution to the stability and solvation free energies in the solvent.

The binding energy of the native 19A1-5MTHF system was calculated by MMPBSA to be around -28 kcal/mol, which is ‘very strong’, and the binding energy of the native PMX system is slightly higher than that of the 5-MTHF system, which is consistent with the results of molecular docking and protein structure analysis. After plasma oxidation, the Coulomb interaction between 19A1_OX and 5-MTHF is slightly enhanced, and the van der Waals hydrophobic interaction does not change much. The repulsive force of the polar solvation energy was significantly increased, indicating that the disruption of the polar environment of the 5-MTHF binding pocket after plasma oxidation of 19A1 counteracted the enhancement of MM, and the final dG was significantly decreased to -16.866 KJ/mol, which is a 40 % reduction. In the PMX system, on the other hand, the repulsive force of the polar solvation energy was reduced, which was the main manifestation of the difference between the folate and pemetrexed transporter of hSLC19A1 after plasma oxidation, and the conformation of 5-MTHF and PMX showed relative rotation and shift, with the binding sites not being exactly the same as the root cause of this difference. In order to analyse the contribution of individual residues to the change in binding free energy in more depth, the residue energy contributions were plotted as shown in Figure 9, and the binding energies of several important residues are listed in Table 5.

In two 5MTHF systems, the binding free energy of Arg42 mutated from -12.945 kcal/mol to +2.088 kcal/mol, suggesting that plasma oxidation resulted in a shift from a strong binding contribution to a slight repulsion of this residue. This mainly stems from the significant enhancement of the polar solvation energy (PB) (+14.772 kJ/mol in 19A1-5MTHF and +37.547 kJ/mol in 19A1ox-5MTHF), which may be attributed to the stretching of the side chain of Arg42 after plasma oxidation, which partially exposes the originally buried charged group to the solvent, resulting in the need to overcome a stronger solvation The change in binding energy of Arg133 also followed the same trend as Arg42, and both residues were located near the negatively charged glutamate portion of 5-MTHF, suggesting that the positively charged pocket of the protein was disrupted in the oxidized system. Lys411 contributes as much as -5.756 kJ/mol to the native system, and is an important residue in hSLC19A1’s recognition of 5-MTHF, whose polar solvation energy also increased from 8.618 kJ/mol to 31.821 kJ/mol, shifting from strong binding to heavy repulsion.The binding energy of Tyr281 was reduced from -4.430 kJ/mol to -0.925 kJ/mol, in which the MM energy was reduced from -16.222 to -8.428 kJ/mol, and the coulombic interaction was significantly reduction and some reduction in van der Waals hydrophobic interactions, which is consistent with the analysis of the reduction in hydrogen bonding and hydrophobic interactions formed by Tyr281 as analyzed above. In contrast, in the two PMX systems, the binding energies of important binding residues such as Arg133, Tyr282, and Lys411 remain stable, although the mutation of Cys33 leads to an increase in repulsive force from 0.47 kJ/mol to 9.619 kJ/mol. The overall polar solvation energy of repulsion also increased, but the increase in MM energy resulted in no decrease in the overall binding energy. Differences in the intermolecular interaction energies of hSLC19A1 bound to 5-MTHF and PMX after oxidation by plasma, as well as contributions specific to each residue, were revealed by residue free energy decomposition.

Figure 10 quantifies the binding free energies of 5-MTHF as well as PMX to natural and oxidized hSLC19A1. It is clear from the figure that the binding free energies of 5-MTHF with natural and oxidized channels are 29.88 kcal/mol and 22.25 kcal/mol, respectively.This suggests that under plasma oxidation, the binding between 5-MTHF and hSLC19A1 is weakened, and the ability of the cancer cells to transport folate is reduced, thus hindering the one-carbon metabolic process of the cancer cells and destroying the antioxidant system of cancer cells. In contrast, the decrease in binding energy of PMX was extremely small, 34.94 kcal/mol before oxidizing the protein and 33.04 kcal/mol after oxidation. From the previous analysis, this is due to the rotational and fine-tuning relationship between the conformation of PMX and 5-MTHF. Since PMX would compete with folate for the binding site, the affinity of PMX after plasma oxidation of hSLC19A1 was much higher than that of 5-MTHF, and it would significantly reduce the folate uptake of the cancer cells when transporting PMX into the cancer cells, which indicated that the plasma could synergise with pemetrexed to work together to disrupt the antioxidant system of the cancer cells, and synergise with the plasma to further destroy the key intracellular components.

4. Conclusions

In this paper, the effect of Cold atmospheric plasma (CAP) oxidation on the uptake of 5-MTHF by cancer cells via hSLC19A1 was explored at the molecular level through classical molecular dynamics simulations, and the interference of cold atmospheric plasma on the transport of the antifolate chemotherapeutic drug pemetrexed while affecting the folate transport in cancer cells was considered. The results showed that plasma did not impede the uptake of the antifolate chemotherapeutic drug PMX by hSLC19A1 while disrupting its interaction with 5-MTHF. This opens new avenues to disrupt the antioxidant pathway in cancer cells.

We investigated both hSLC19A1 structural changes and changes in the interaction between hSLC191A1 and 5-MTHF/PMX. The results of the structural study showed that under plasma oxidation, the channel radius of hSLC19A1 upon transporting 5-MTHF narrowed at the bottom of the channel at 4.5 to 6 nm on the z axis, and the overall channel length was shortened; the stability of the secondary structure was reduced, and the secondary structures of Cys30,33 and Met130 shifted from A-Helix to structures such as T-turn and S bend, futhermore, increased flexibility led to inward tightening of the channel entrance; protein-ligand interaction analyses indicate that the hydrogen bonds formed between hSLC19A1 and 5-MTHF upon oxidation are reduced by an average of 0.71 per frame, followed by a significant reduction in the hydrophobic interaction formed by the glutamate portion of 5-MTHF; molecular docking showed that the predicted affinity of hSLC19A1 with 5-MTHF decreased after oxidative modification, and mmpbsa analysis showed that the repulsion of the free energy of non-polar solvation led to a 40 % decrease in the binding energy of hSLC19A1 with 5-MTHF; in contrast, plasma oxidation did not significantly affect the structure of hSLC19A1 when bound to PMX or the interaction between hSLC19A1 and PMX. This study demonstrates that plasma oxidation has no significant effect on the structure of hSLC19A1 when bounded with PMX and the interaction between hSLC19A1 and PMX. This study demonstrates that plasma can reduce the binding capacity of hSLC19A1 to 5-MTHF without disrupting the transport of pemetrexed. Disrupting the uptake of folate by cancer cells while combining with antifolate chemotherapeutic agents can help weaken the antioxidant defence system of cancer cells, increase the level of ROS in cancer cells, and promote the apoptosis of cancer cells.

Our simulations investigated the potential pathways by which plasma disrupts the antioxidant system of cancer cells and considered the possibility of synergistic treatment in combination with conventional antifolate chemotherapeutic agents, explored the intrinsic microscopic mechanisms, provided atomic-level insights into the mechanism of the integrated effect on cancer cells under plasma action, and provided theoretical references to the mechanism of optimal regulation of plasma-induced selective apoptosis of cancer cells.