Molecular dynamics simulation of salty ions diffusion in phospholipid bilayer

Abstract Six molecular simulation systems were established to investigate the diffusion property of salty ions in phospholipid bilayers at different concentrations and intensities of electric field. The results show that the electric field intensity is inversely proportional to the value of mean square displacement (MSD) due to the strong disturbance of charged ions by molecules in the phospholipid bilayers. The strength of the electric field increases the free volume fraction (FFV) of ions. However, the existence of electric field makes the surface charges of phospholipid bilayer redistribute and accumulate, and the increase of electrostatic potential intensifies the disturbance of carbon chain to salty ions. This study reveals the diffusion mechanism of ions on the surface of muscle and skin at the microscale, which can provide guidance for muscle energy storage.


Introduction
The migration of salty ions in the human body is of great significance to the lasting movement of muscles and the good operation of various organs in the human body (Rutkove, Hyeuknam, Maria et al., 2018a;Suleman & Riaz, 2020). In the case of strenuous exercise, ions in the human body will penetrate the skin surface through sweat, resulting in a large loss of salty ions (Bazant et al., 2009;Bhuiyan & Outhwaite, 2009;Pirtini & Herman, 2010). Under the condition of high-intensity exercise, the body surface temperature of athletes increases with the continuous rise of body temperature, which stimulates the relaxation of body surface pores and accelerates the loss of electrolyte in the body (Aaaf et al., 2020;Halter et al., 2015;Yan & Ch, 2018;Zhou & Herman, 2018). This will directly cause the functional failure of human organs in the case of extreme water shortage and electrolyte, which will affect the athletes' performance in the competition (Andreassen et al., 2009;Ekstrand, Healy, Waldén et al., 2012a;Kwon et al., 2017).
Based on this, the diffusion of electrolyte needs a lot of research under micro-conditions to reveal the mechanism of salty ion loss. Many researchers have done a lot of research on muscle injury and muscle imaging through various methods (Ekstrand, Healy, Waldén et al., 2012b;Rutkove, Hyeuknam, Maria et al., 2018b;Slavotinek, 2010). By studying sciatic nerve regeneration and functional recovery, Fu et al. (Tf et al., 2020) found that after nerve injury and immediate repair (IR), electrical stimulation of muscle accelerates axon regeneration and functional recovery by promoting autophagy flow of distal nerve segments, but its mechanism is unclear. Cai et al. (Cai et al., 2020) found that fibroblast growth factor (FGF6) can reduce skeletal muscle atrophy and promote the transformation from slow muscle to fast muscle fiber, so as to promote the functional recovery of regenerated skeletal muscle after innervation.
As a powerful tool to study and observe the microscopic mechanism among atoms and molecules, molecular dynamics (MD) simulation is widely used in the fields of aerospace, advanced materials and biomedicine Hu et al., 2018;Li et al., 2021;Zhang et al., 2019). So far, researchers have clearly proved the dynamics of multistep nucleation process through MD simulation, including the characteristic metastable clusters before supercritical smectic nucleation, which cannot be explained by classical nucleation theory (Li et al., 2017;Stock et al., 2017;Takahashi et al., 2021). Zachary (Rollins et al., 2020) employed directed MD simulations, the results show that mutations in the peptide of the pMHC binding groove, which is known to discretely change the response of T cells to antigens, changed the MHC conformation at equilibrium. Nafiseh (Farhadian et al., 2022) MD simulation is applied to investigate drug transport in both pure state and conjugated with neutral gold nanoparticle (AuNP) as a drug carrier inside dipalmitoylphosphatidylcholine (DPPC) membrane, results demonstrate that the interaction between drug molecules and DPPC changes after drug conjugating with AuNP, and solving the problems in nanodrug translation. Terence et al.  investigated the effects of water-based substitutional defects in zeolitic imidazolate frameworks (ZIF)-8 membranes on their reverse osmosis (RO) desalination performance with MD simulation, the results show that ion adsorption on the membranes occurs at either the nitrogen atoms or the defect sites and the presence of linker defects increases membrane hydrophilicity.
Although researchers have made great progress in the fields of muscle injury rehabilitation and tumor treatment, few people hope to prevent muscle damage by observing ion diffusion at the micro-scale, so it is difficult to prevent muscle damage from the perspective of body fluid electrolyte diffusion mechanism. However, the migration characteristics of microions have important research value for energy storage and prevention of electrolyte loss. In this study, the ion diffusion ability under three different electric field intensities of 10 V/m, 15 V/m and 20 V/m is discussed through the mixing model of component phospholipid bilayer and electrolyte charged ions. In addition, three simulation systems with different ion concentrations were established to study the effect of electric field on ion diffusion in body fluid.

Model of the computing system
The phospholipid bilayer structure and salty ion components are established, as shown in Figure 1. All simulation systems are performed by Materials studio (MS; Liu & Song, 2021;song et al., 2020).

Figure 1. Phospholipid bilayer and body fluid composition.
The size of the simulation system is 58 Å × 58 Å × 58 Å, the periodic boundary conditions are applied to the X, Y and Z directions, respectively. Furthermore, the temperature was controlled by deploying Nosé-Hoover thermostat method. The interaction is described by the force field type of COMPASS (Sun, 1998). In order to explore the effect of concentration on ion diffusion ability, three systems with different ion concentrations (P-1, P-2 and P-3) are constructed, as depicted in Table 1. Besides, the runtime of each simulation system is 100 ps and three times repeated to increase the reliability of calculation.

Mean square displacement and radial distribution function
In order to measure the diffusion property of ions in electrolyte, mean square displacement (MSD) is always utilized by researchers to investigate the diffusion performance of different ions in microscale. In this section, all the detailed information of MSD under multiple electric field and different concentrations is depicted in Figures. 2 and 3. The temperature of each simulation system is set in 300 K in order to eliminate the effect of temperature. MSD value of K + ions is the highest, while that of the SO 4 2ions is the smallest under different electric field, as shown in Figure 2. Besides, Figure 2 (a-e) indicates that K + , Mg 2+ , SO 4 2ions show the strongest diffusion performance under the electric field of 10 V/m, which is an interesting phenomenon. Additionally, almost all the MSD values of ions are the lowest under 20 V/m in the simulation systems, this indicates the increase of intensity of electric field that induces the reduce of MSD value of each ion in the simulation systems. However, different from the other ions in the simulation systems, the change of electric field has little impact on the diffusion of HCO 3 − , which indicates that the electric field has the greatest impact on the diffusion of K + and Mg 2+ ions. Figure 3 shows the MSD curve of salty ions at different concentrations; it can be seen that the MSD value decreases with the rise of concentration, and HCO 3 − has the strongest diffusion performance at the concentration of P-1, while the diffusion capacity of Cl − ion at this concentration is the smallest. Among them, the MSD value of HCO 3 − ion at P-1 concentration is the highest, which means that it has the strongest diffusion capacity, as shown in Figure 3(f). Followed by K + and Na + ions, their MSD curves at different concentrations show consistency, which has a great correlation with their atom size and the number of charges, as shown in Figure 3(a-b). It is worth noting that the MSD curve trend of Na + ions in P-1 and P-2 systems is consistent, which indicates that the diffusion performance of Na + ions is not affected under the condition of high concentration. In addition, the electric field has the least effect on the diffusion performance of Cl − and SO 4 2ions, and the MSD of both at each concentration does not exceed 100, as shown in Figure 3(d-e). Different from HCO 3 − , the MSD values of SO 4 2at the concentration of P-1 and P-2 are higher than that at the concentration of P-3.
Radial distribution function (RDF) is a key parameter to reveal the interatomic and intraatomic interaction. The RDF of each ion and phospholipid molecule is depicted in Figure 4 and Figure 5. As we can see, from Figure 4, the RDF value of Cl 1and phospholipid molecule is the highest, followed by the value of HCO 3 − and phospholipid molecule, and this indicates a strong interaction energy exists between Cl 1and phospholipid molecule under three different concentrations, which can explain the reason for the low level of Cl − ion's MSD throughout the whole process of simulation. Figure 5 also shows the RDF of each ion and phospholipid molecule under three different intensities of electric field, the RDF value of Cl 1and phospholipid molecule still remain the highest, and this proves the strong interaction between Cl 1and phospholipid molecule has no relationship with the intensity of electric field.

Free volume
Free volume fraction (FFV) is used to characterize the free diffusion space of solvent particles in liquid, which is an important method to explore the diffusion mechanism of particles in liquid. According to the free volume theory of Fox and Flory, the calculation formula of FFV is given in Eq. (1): where V 0 is the volume occupied by molecules, and V F is the free volume not occupied by molecules. Figure 6 shows the FFV of the model when the temperature reaches 300 K, in which the blue part represents the free volume and the gray area represents the occupied volume. The comparison shows that the size and position of the free volume of the hybrid model have changed significantly before and after the application of electric field. Table 2 shows the FFV of the model under different electric field intensity. FFV of the system increases with the increase of electric field intensity, which indicates that the increase of electric field makes the diffusion space of charged ions in phospholipid bilayer bigger.    The cross-section of free volume for the simulation systems under multiple intensity of electric field is shown in Figure 7 to reveal the free volume distribution inside the phospholipid bilayer, and area covered by red color represents the phospholipid molecules. Two sides of the phospholipid layers covered by the color of blue, indicating that the position around phospholipid molecules is not available for the diffusion of ions. Additionally, the FFV is inversely proportional to the electric field strength, which is consistent with the data in Table 2.

Volume charge distribution
The volume charge distribution has an important influence on the ion migration in liquid electrolyte. We extracted the charge molecules on the particle surface in the X-axis direction to explore the effect of phospholipid bilayer on the ion displacement under electric field. It can be seen from Figure 8 that compared with the volume charge distributions of the systems under the electric field of 10 V/m, 15 V/m and 20 V/m, the volume charge distributions under 0 V/m are much more uniform. Besides, the volume charge distributions under 10 V/m and 15 V/m electric fields are relatively similar. After experiencing the first peak within 0-1 nm, the second peak under 10 V/m electric field has a large jump. On the contrary, the second peak at 15 V/m is much smaller. The amount of charge under 20 V/m electric field further decreases and approaches 0, which means that the aggregation of negatively charged ions under high voltage intensifies in the range of 0-2 nm. Interestingly, the amount of charge in the range of 2.5-3.5 nm increases sharply, which indicates that a large amount of positive charge in the range of 2.0-2.5 nm migrates to the range of 2.5-3.5 nm. This is because the position of 3 nm is in the channel between the phospholipid bilayers, and ions move rapidly in this channel under the influence of electric field.
Electrostatic potential is an important method to explore the charge distribution on the surface of molecules and the resulting potential. For phospholipid molecules, the electrostatic potential distribution and finite volume of its head group will affect the internal ion diffusion (Lebar et al., 2016). The electrostatic potential (ESP) of SO 4 2-, HCO 3 − and phospholipid molecule is shown in Figure 9. The ESP distribution of SO 4 2is highly symmetrical, and the negatively charged part occupies the most surface area of molecules, which means that SO 4 2is more likely to be attracted by the positively charged on the left half part of phospholipid molecules. However, although most of the molecular structures of HCO 3 − ions show weak electronegativity (0.29 kcal/mol), the surface of its H atom shows very strong electronegativity (29.8 kcal/mol), which means that HCO 3 − ions are more easily limited by the two carbon chains of phospholipid molecules (−0.3 kcal/mol), as shown in Figure 9(c). Additionally, the head group region shows great positive polarity (29.79 kcal/mol), which makes it easier to attract negatively charged particles.
The whole diffusion and interacting process of salty ions in the phospholipid bilayers is depicted in Figure 10. The inhomogeneous distribution of the surface electrostatic potential makes the charge on the surface of molecules in the phospholipid bilayer redistribute under the electric field, and the carbon chains obtain additional power due to the charge transfer, as shown in Figure 10 (a). Meanwhile, the interacting force between carbon chains and salty ions enhanced under the high electric field intensity, resulting in greater displacement of salty ions, as shown in Figure 10(c). The channel between two layers widened under the high intensity of electric field, leading to the increase of free volume in the simulation systems, as shown in Figure 10

Conclusion
MD simulation is utilized to investigate the diffusion process of salty ions at multiple concentration in phospholipid bilayer, and the impact of different electric field intensity on the ion diffusion is also studied. The results show that the diffusion capacity of ions decreases with the increase of the electric field intensity, and the main reason is that the charged ions have a strong combination with the molecules in the phospholipid bilayer, and the higher electric field intensity makes the charge inside the system accumulate on the surface of molecules in phospholipid bilayer, which leads to the enhanced connection between ion and phospholipid molecules. Besides, K + ion has the strongest diffusion capacity due to the weak interaction between K + ions and phospholipid molecules, while the displacement performance of Cl 1is the worst because of the strong interaction between phospholipid molecules. Meanwhile, the increase of ion concentration leads to the decrease of ion diffusion ability. Additionally, the free volume is inversely proportional to the electric field strength, which is also an important cause for the limit of ion's diffusion in phospholipid bilayer. Furthermore, through the ESP calculation, the carbon chain of phospholipid molecules is much easier to connect negatively charged ion. This paper will provide a reference for the fluidity of electrolytes in human body.