Possible Mechanism of Glucose Uptake Enhanced by Cold Atmospheric Plasma : Atomic Scale Simulations

Cold atmospheric plasma (CAP) has shown its potential in biomedical applications, such as wound healing, cancer treatment and bacterial disinfection. Recent experiments have provided evidence that CAP can also enhance the intracellular uptake of glucose molecules which is important in diabetes therapy. In this respect, it is essential to understand the underlying mechanisms of intracellular glucose uptake induced by CAP, which is still unclear. Hence, in this study we try to elucidate the possible mechanism of glucose uptake by cells by performing computer simulations. Specifically, we study the transport of glucose molecules through native and oxidized membranes. Our simulation results show that the free energy barrier for the permeation of glucose molecules across the membrane decreases upon increasing the degree of oxidized lipids in the membrane. This indicates that the glucose permeation rate into cells increases when the CAP oxidation level in the cell membrane is increased.


Introduction
Diabetes is a chronic disease related to an abnormal increase of the glucose level in the blood.Approximately 75% of glucose in the body is consumed by skeletal muscle cells that are stimulated by insulin [1,2].Failure of glucose uptake leads to the presence of a high level of sugar in the blood.This is due to one of two mechanisms or a combination of both, i.e., (a) disproportional production of insulin or (b) inadequate sensitivity of cells to insulin [3][4][5].Most of the currently available pharmaceuticals are inefficient in the treatment of diabetes, thus alternative insulin independent healing methods need to be found in order to increase glucose uptake by muscle cells.
Recently, Kumar et al. investigated the impact of cold atmospheric plasma (CAP) on the regulation of glucose homeostasis [6].The experimental results showed that the glucose uptake is significantly enhanced in skeletal muscle cells after plasma treatment, exhibiting the beneficial effects of CAP.Furthermore, higher levels of intracellular Ca 2+ and reactive oxygen species (ROS) in CAP treated cells were observed.An increase in intracellular ROS and Ca 2+ ions helps to increase the glucose uptake by skeletal muscle cells [7,8].Increases in the Ca 2+ ion concentration as well as the uptake of middle-sized, membrane-impermeable molecules after direct CAP exposure were also observed in [9].The authors attributed these effects to an increase in cell permeability caused by CAP-generated electric stimulation and the delivery of OH radicals into cells.Moreover, Vijayarangan et al. investigated CAP parameters and conditions for drug delivery across HeLa cells [10].They determined the efficient treatment time (i.e., low toxicity) and the range of frequencies with an optimal number of pulses as key parameters for the cell membrane permeability [10].In addition, Leduc et al. determined the maximum radius for macromolecules that are capable of ingressing into HeLa cells, applying a specifically designed CAP source [11].Hence, an increase in the cell membrane permeability might play an important role in the delivery of the abovementioned species into the cell.The underlying mechanisms, however, still remain unclear, and need more thorough investigation.Computer simulations might be a useful tool to gain insight into the atomic level processes.
In the context of plasma medicine, we have already performed several computational studies, by means of molecular dynamics (MD) simulations, using a phospholipid bilayer (PLB) as a model system for the plasma membrane.Specifically, we studied the effect of lipid oxidation on phosphatidylserine translocation across the plasma membrane, which plays a vital role in apoptosis signaling [12].Furthermore, we investigated the ROS oxidation of the head groups and lipid tails in the membrane [13], the permeation of ROS across oxidized and non-oxidized membranes, including the synergistic effect of plasma oxidation and electric field [14], and the hampering effect of cholesterol [15,16].In general, our investigations showed that oxidation of the membrane leads to an increase in its fluidity and permeability to ROS, thereby affecting the abovementioned processes.
In this study, we investigate glucose translocation across native and oxidized membranes in order to provide a possible explanation for the mechanism of glucose uptake observed in previous experiments.In particular, we perform MD simulations to calculate the free energy barriers for glucose transport through oxidized and non-oxidized membranes.Comparison of the latter shows the effect of CAP oxidation on the glucose transport across the membrane, which might explain the increased level of glucose uptake after the plasma treatment of cells, as observed experimentally.

Simulation Setup
We performed MD simulations in order to study the glucose translocation across both intact and oxidized PLBs.The PLB is considered to be a simple model system that represents the eukaryotic cell membrane, since it determines the thickness of the bilayer.A schematic representation of the intact PLB is given in Figure 1a.It consists of 128 phospholipids (PLs), covered by 6000 water molecules, organized in two lamellae (i.e., 64 PLs, with a corresponding water layer at the top, and 64 at the bottom).As the PL molecule, we used 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), as depicted in Figure 1b.Glucose molecules were randomly placed in the xy-plane of the upper side of the PLB, i.e., in the water phase, as well as in the lipid tail region, see Figure 1a (and see below for more details).
Simulations were carried out using the GROMACS package (version 5.1) [17], by applying the GROMOS (43A1-S3) force field [18].To study the effect of oxidized PLs, we used aldehyde oxidation products (see Figure 1b, DOPC-ALD), which are one of the major oxidation products [19].Note that CAP yields a cocktail of reactive species and thus, could possibly form a range of products upon oxidation of the PLB, but the formation of aldehyde groups (i.e., DOPC-ALD) was prominently observed in CAP-treated vesicles [13].These oxidation products were also used in our recent simulation studies on various properties of the PLB [13][14][15].The force field parameters for the aldehyde groups in the oxidized PLs were obtained from [20] and the parameters of glucose were based on [21].The Packmol package was employed to create initial configurations of the intact and oxidized PLB systems [22].Two aldehyde oxidation products were created from the non-oxidized (i.e., native) PLBs containing 128 PLs, by replacing 32 and 64 DOPC molecules with DOPC-ALD, corresponding to concentrations of 25% and 50%, respectively.The presence of 25% and 50% DOPC-ALD in the oxidized PLB does not necessarily correspond to experimental oxidation levels, but performing the calculations for a lower degree of oxidation would require excessive calculation times to give the same qualitative conclusions.Hence, this oxidation degree was high enough to observe the effect of oxidation within an acceptable calculation time, but low enough so that pore formation did not occur within the simulated time scale [15].Thus, in total, three model systems were studied in our simulations, i.e., native (0%) as well as aldehyde oxidized (25% and 50%) PLBs.For each system, we created four different structures (e.g., four native PLBs) extracted from the last 40 ns trajectory of the 200 ns equilibration run, with a time interval of 10 ns, in order to obtain an average free energy profile (FEP) of glucose transition across each system (see next section).All structures (in total 12, including the native PLB) are initially optimized using the steepest descent algorithm and then equilibrated for 200 ns by so-called NPT simulations (i.e., at a constant number of particles (N), pressure (P) and temperature (T)), at 310 K and 1 bar, employing the Nose-Hoover thermostat [23] with a coupling constant of 0.2 ps as well as the semi-isotropic Parrinello-Rahman barostat [24] with a compressibility and coupling constant of 4.5 × 10 −5 bar −1 and 0.1 ps, respectively.For the non-bonded interactions, a 1.1 nm cut-off was applied.The long range electrostatic interactions were described by the particle mesh Ewald (PME) method [25], using a 1.1 nm cut-off for the real-space interactions in combination with a 0.15 nm spaced grid for the reciprocal-space interactions.Subsequently, a series of umbrella sampling (US) simulations [26] were run for 20 ns, applying, again, the NPT ensemble (see next section), of which the last 10 ns was used for further analysis.In all simulations, a time step of 2 fs was used.Periodic boundary conditions were applied in all three directions.
Plasma 2018, 2, x FOR PEER REVIEW 3 of 7 For each system, we created four different structures (e.g., four native PLBs) extracted from the last 40 ns trajectory of the 200 ns equilibration run, with a time interval of 10 ns, in order to obtain an average free energy profile (FEP) of glucose transition across each system (see next section).All structures (in total 12, including the native PLB) are initially optimized using the steepest descent algorithm and then equilibrated for 200 ns by so-called NPT simulations (i.e., at a constant number of particles (N), pressure (P) and temperature (T)), at 310 K and 1 bar, employing the Nose-Hoover thermostat [23] with a coupling constant of 0.2 ps as well as the semi-isotropic Parrinello-Rahman barostat [24] with a compressibility and coupling constant of 4.5 × 10 −5 bar −1 and 0.1 ps, respectively.For the non-bonded interactions, a 1.1 nm cut-off was applied.The long range electrostatic interactions were described by the particle mesh Ewald (PME) method [25], using a 1.1 nm cut-off for the real-space interactions in combination with a 0.15 nm spaced grid for the reciprocal-space interactions.Subsequently, a series of umbrella sampling (US) simulations [26] were run for 20 ns, applying, again, the NPT ensemble (see next section), of which the last 10 ns was used for further analysis.In all simulations, a time step of 2 fs was used.Periodic boundary conditions were applied in all three directions.

Umbrella Sampling (US) Simulations
We performed US simulations in order to determine the free energy profiles (FEPs) of glucose translocation across the native, as well as the oxidized, PLBs.In total, we obtained 12 FEPs for all model systems (i.e., four for the native, four for the 25% oxidized, and four for the 50% oxidized PLBs, see previous section) and averaging was performed over four FEPs for each system.For the calculation of each energy profile, we extracted 32 to 36 windows (36 for the native, 34 for the 25% and 32 for the 50% oxidized PLB, due to a decrease of the bilayer thickness) along the z-axis, which were separated by 0.1 nm.These windows were obtained by pulling glucose molecules against the z-axis direction (see Figure 1a) for 500 ps, using a harmonic bias between these glucose molecules and the center of mass of the PLB, with a force constant of 2000 kJ•mol −1 •nm −2 and a pulling rate of 0.01 nm•ps −1 .Note that, in principle, slow pulling rates can be used in the pulling simulations.However, dragging the glucose from water phase into the center of the PLB requires a long computation time.Higher pulling rates can solve this issue, but they can lead to significant disturbances in the PLB.Thus, by using an appropriate (standard) pulling rate, we performed short

Umbrella Sampling (US) Simulations
We performed US simulations in order to determine the free energy profiles (FEPs) of glucose translocation across the native, as well as the oxidized, PLBs.In total, we obtained 12 FEPs for all model systems (i.e., four for the native, four for the 25% oxidized, and four for the 50% oxidized PLBs, see previous section) and averaging was performed over four FEPs for each system.For the calculation of each energy profile, we extracted 32 to 36 windows (36 for the native, 34 for the 25% and 32 for the 50% oxidized PLB, due to a decrease of the bilayer thickness) along the z-axis, which were separated by 0.1 nm.These windows were obtained by pulling glucose molecules against the z-axis direction (see Figure 1a) for 500 ps, using a harmonic bias between these glucose molecules and the center of mass of the PLB, with a force constant of 2000 kJ•mol −1 •nm −2 and a pulling rate of 0.01 nm•ps −1 .Note that, in principle, slow pulling rates can be used in the pulling simulations.However, dragging the glucose from water phase into the center of the PLB requires a long computation time.Higher pulling rates can solve this issue, but they can lead to significant disturbances in the PLB.Thus, by using an appropriate (standard) pulling rate, we performed short pulling simulations to save computation time with minimum perturbations on the PLB.One of the glucose molecules (i.e., the glucose in the upper water phase, see Figure 1a) was pulled until it reached the center of mass of the PLB, while the second one moved towards the lower water phase from the center of the bilayer.Thus, each US simulation involved two glucose molecules.In this way, we saved computational resources, thereby increasing the number of sampling points.Note that these two glucose molecules were separated from each other at least by 3 nm in the z direction; hence, there was no interaction between these two molecules.In principle, we could have used three glucose molecules in each US simulation.However, in order to cause minimal disturbance to the PLB system, we chose two glucose molecules instead of three.It should be mentioned that, in reality, adsorption or chemical reaction of glucose might take place in the PLB.However, these processes cannot be studied with conventional non-reactive MD simulations, due to the limitations in time and reactivity.Nevertheless, the US simulations can predict how often the glucose transport occurs across the PLB before and after oxidation, through calculation of the FEPs.
As mentioned above, we extracted 32 to 36 US windows from our 500 ps pulling simulation.Hence, 32 to 36 US simulations were performed to construct a single FEP.Each US simulation lasted 20 ns, and the last 10 ns were used for analysis, i.e., to collect the US histograms and calculate the FEPs.Note that the pulling simulation trajectory was used only to extract windows/frames for the further US simulation to obtain the FEPs.During the pulling simulations, the glucose dragged water molecules with it into the hydrophobic core of the bilayer, but these molecules escaped from the hydrophobic core within the initial 10 ns of US simulation.Thus, the last 10 ns of the US simulation was an adequate time for calculating the FEPs, as there were not any water defects or hydration layers in the hydrophobic part of the membrane.In each US simulation, the glucose molecules were able to freely travel in the xy-plane, but their movement in the z-direction was restricted by applying a harmonic bias with a spring constant of 2000 kJ The FEPs were constructed using a periodic version of the weighted histogram analysis method (WHAM) [27], as implemented in GROMACS.We analyzed our simulation systems to identify underestimated possible "hidden barriers" based on previous literature [28].We did not define any hidden barriers because in our US simulations, the data was sufficiently sampled.Indeed, in the four model systems used, the glucose was randomly positioned in the xy-plane to escape trapped metastable states and to allow estimation of the error bars that are associated with choosing the initial model systems.As noted above, the final energy profiles were obtained by averaging four independently-built FEPs for each system which differed from one another based on their starting structure to allow for some statistical variation.The uncertainties associated with the FEPs were obtained by calculating the standard deviations between these four FEPs for each system.Thus, in total, 410 US simulations were performed for the calculation of the FEPs.

Results and Discussion
The US MD simulations allowed us to elucidate the glucose translocation across the native and oxidized PLBs which gave insight into the possible mechanism of glucose ingress triggered by plasma oxidation of the cell membrane.Figure 2 illustrates the symmetric FEPs of glucose transfer across native (0%) as well as oxidized (25% and 50%) PLBs.
It is clear that, in all three cases, the ∆G started to drop when glucose entered the hydrophilic head group region from the water phase (see Figure 1b).This means that the head group of the PLB is energetically the most favourable region, showing the minimum energy for glucose transfer.Moreover, this energy minimum slightly shifted towards the centre of the bilayer in the case of oxidized PLBs.We know from our previous studies that when lipid tails are oxidized, this eventually leads to a drop of the bilayer thickness [13,14].In other words, the shortened tail and aldehyde products move towards the water phase, thereby increasing the area and fluidity of the PLB.Consequently, this results in a decrease in PLB thickness, e.g., the native DOPC bilayer has a thickness of around 3.89 nm, whereas, in the case of a 50% DOPC-ALD bilayer, the thickness decreases to 3.33 nm [14].Therefore, the free energy minima in Figure 2 are found at around |z| = 1.9 nm and |z| = 1.76 nm for native and 50% oxidized PLBs, respectively.
Continuing the motion of glucose towards the hydrophobic tail region leads to an increase of the free energy for translocation, representing the role of the membrane as a permeation barrier.In case of the native PLB, the free energy barrier for the permeation of glucose was 51 ± 3 kJ/mol, which is within the range of experimental results [29,30].On the other hand, for the 25% and 50% oxidized PLBs, this barrier for the translocation of glucose decreased to values of 37 ± 4 kJ/mol and 28 ± 4 kJ/mol, respectively (see Figure 2).Hence, the obtained results show that the free energy barrier for the transport of glucose molecules across the PLB decreases when the oxidation degree is increased.This, in turn, leads to an increase in the probability of glucose permeation to the cell interior.The obtained simulation results can be correlated with the experimental observations [6,9] as the plasma treatment of cells most probably rise to oxidation of the cell membrane, thereby increasing the glucose (or other middle-sized molecules, as well as Ca 2+ ) translocation rate.
Plasma 2018, 2, x FOR PEER REVIEW 5 of 7 Continuing the motion of glucose towards the hydrophobic tail region leads to an increase of the free energy for translocation, representing the role of the membrane as a permeation barrier.In case of the native PLB, the free energy barrier for the permeation of glucose was 51 ± 3 kJ/mol, which is within the range of experimental results [29,30].On the other hand, for the 25% and 50% oxidized PLBs, this barrier for the translocation of glucose decreased to values of 37 ± 4 kJ/mol and 28 ± 4 kJ/mol, respectively (see Figure 2).Hence, the obtained results show that the free energy barrier for the transport of glucose molecules across the PLB decreases when the oxidation degree is increased.This, in turn, leads to an increase in the probability of glucose permeation to the cell interior.The obtained simulation results can be correlated with the experimental observations [6,9] as the plasma treatment of cells most probably gives rise to oxidation of the cell membrane, thereby increasing the glucose (or other middle-sized molecules, as well as Ca 2+ ) translocation rate.Note that our simulations only provide one possible explanation for the increased level of glucose [6] in cells after plasma treatment, while other mechanisms might play a role as well [9][10][11].Therefore, further investigations should be performed to obtain a more complete picture of the plasma effect on the membrane permeability.This might be achieved, for instance, by silencing the glucose transporter proteins (or GLUTs) and measuring the concentration of intracellular glucose molecules after CAP treatment.The latter would clearly reveal the role of cell membrane permeability (induced by CAP oxidation of lipids) in transporting glucose molecules.Note that our simulations only provide one possible explanation for the increased level of glucose [6] in cells after plasma treatment, while other mechanisms might play a role as well [9][10][11].Therefore, further investigations should be performed to obtain a more complete picture of the plasma effect on the membrane permeability.This might be achieved, for instance, by silencing the glucose transporter proteins (or GLUTs) and measuring the concentration of intracellular glucose molecules after CAP treatment.The latter would clearly reveal the role of cell membrane permeability (induced by CAP oxidation of lipids) in transporting glucose molecules.

Figure 1 .
Figure 1.(a) Intact or native 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) phospholipid bilayer (PLB), together with two glucose molecules in the water and lipid tail regions.For the sake of clarity, the N and P atoms of DOPC are shown with larger beads.The bilayer center is indicated by the red dashed line; (b) Schematic illustration of native (DOPC) and oxidized (DOPC-ALD) PLs.The head group consists of choline, phosphate and glycerol, whereas the lipid tails are two fatty acid chains.

Figure 1 .
Figure 1.(a) Intact or native 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) phospholipid bilayer (PLB), together with two glucose molecules in the water and lipid tail regions.For the sake of clarity, the N and P atoms of DOPC are shown with larger beads.The bilayer center is indicated by the red dashed line; (b) Schematic illustration of native (DOPC) and oxidized (DOPC-ALD) PLs.The head group consists of choline, phosphate and glycerol, whereas the lipid tails are two fatty acid chains.

Figure 2 .
Figure 2. (a) Symmetric free energy profiles for the translocation of glucose across the native and oxidized PLBs.A PLB is schematically illustrated in the background, to indicate the position of the water layer, the head groups and the lipid tails.For clarity, the zoomed extrema of the profiles are shown in (b,c).Errors associated with the umbrella sampling (US) calculations are depicted in pale color.

Figure 2 .
Figure 2. (a) Symmetric free energy profiles for the translocation of glucose across the native and oxidized PLBs.A PLB is schematically illustrated in the background, to indicate the position of the water layer, the head groups and the lipid tails.For clarity, the zoomed extrema of the profiles are shown in (b,c).Errors associated with the umbrella sampling (US) calculations are depicted in pale color.