Molecular Dynamics Simulation on the Diffusion of Flavor, O2 and H2O Molecules in LDPE Film

The diffusion of five flavor organic molecules, including D-limonene, myrcene, ethyl hexanoate, 2-nonanone, and linalool in low density polyethylene (LDPE) film were investigated by combined experimental and molecular dynamics (MD) simulation studies. The diffusion coefficients derived from the prediction model, experimental determination, and MD simulation were compared, and the related microscopic diffusion mechanism was investigated. The effects of the presence of the flavor organic molecules on the diffusion of O2 and H2O in polyethylene (PE) were also studied by MD simulation. Results show that: The diffusion of five flavor molecules in LDPE is well followed to Fick’s second law by the immersion experiment; MD simulation indicates the dual mode diffusion mechanism of the flavor molecules is in LDPE; the diffusion coefficients from MD simulation are close to those from the experimental determination, but are slightly larger than those values; the presence of the flavor organic molecules hinders the diffusion of O2 and H2O, which can be well explained from the fraction of free volume (FFV) and interaction energy calculation results derived from MD simulation.


Introduction
Polyethylene (PE) is a kind of thermoplastic, which has excellent chemical stability, heat sealing, toughness and high elasticity, as well as cold resistance and water vapor barrier [1,2]. Low density polyethylene (LDPE) is commonly used in food packaging because of its excellent thermal sealing and chemical stability [3]. Food flavor is an important quality characteristic of food [4]. It is reflected by volatile compounds in food, which determines the quality of food [5][6][7]. Food flavor molecules will be adsorbed by packaging materials after long-term storage, and the change of its quantity will lead to the change of food flavor, thus, affecting food sales and consumer complaints [8].
The adsorption of food flavor molecules in packaging materials mainly depends on the characteristics of flavor molecules, packaging materials, storage environment and other factors, which subsequently changes the flavor of food or the properties of polymer materials [9]. During the last decades, many researchers have studied the adsorption of flavor molecules into different packaging materials, such as PE [3,[10][11][12][13], polypropylene (PP) [14], polyethylene teraphthalate (PET) [5,15], polylactate (PLA) [13], etc. In packaging materials, the flavor adsorption amount of the rubber polymer (e.g., PE) is several thousand times higher than that of the glass polymer, such as PET [16]. The adsorbable flavor organic molecules of LDPE film have been extensively studied, including the alkenes, the esters, the aldehydes, the ketones, and the alcohols [3,[10][11][12][13]17,18]. The flavor of liquid food is easier to be lost than that of solid food [19]. The concentration range of flavor organic molecules in liquid food is ppm or less, which is within the threshold range of human taste [20].
The molecular dynamics simulation method is a good tool to illustrate the microscopic phenomenon of the molecule mass transfer in the packaging materials [21]. More detailed information and basic principles of the interaction between molecules and polymers can be obtained [22][23][24][25]. Börjesson et al. performed Monte Carlo (MC) and MD simulation to calculate solubility coefficients and diffusion coefficients of oxygen and water in PE to obtain a microscopic diffusion mechanism [26]. Yue et al. used MD method to study the diffusion of benzene, toluene and ethyl benzene in amorphous PE membrane, and found the simulation results had the same trend as the experimental results [27]. Wang et al. determined the diffusion coefficients of molecules with different molecular weights in amorphous PET by MD simulation [28]. Tao et al. used Grand canonical Monte Carlo (GCMC) and MD simulation to study the adsorption and diffusion of O 2 in PP and showed that the loading of O 2 increased and its diffusion coefficient in PP decreased with the increase in the polymerization degree of PP [29]. However, the adsorption of food flavor molecules in LDPE film was mainly studied by experiment [3][4][5]7], and there were few studies on the adsorption and diffusion of flavor substances in LDPE by MD simulation. It is also meaningful to analyze the diffusion behavior of flavor organic molecules in PE and their effects on the diffusion of O 2 and H 2 O in LDPE from the microscopic view.
In this work, the diffusion coefficients of food flavor organic molecules in LDPE film were studied by two-sided immersion experiment, Brandsch prediction model and MD simulation. The differences in diffusion coefficients obtained by the above three methods were compared. By calculating the interaction energies between the molecules and PE chain, the interaction energies between O 2 or H 2 O molecule and flavor molecules, and the FFV of PE cell, the related diffusion mechanism of diffusion molecules in LDPE film was further analyzed. Furthermore, the diffusion coefficients of O 2 and H 2 O in LDPE film before and after adsorption were calculated by the MD method, and the effects of flavor molecules on the diffusion of O 2 and H 2 O were investigated.

Film Adsorption Experiment
The contact process between liquid food and LDPE film was simulated by two-sided immersion method in the food simulation solution [10]. Five flavor reagents were mixed together into 10% (V/V) aqueous ethanol solution with a concentration of 500 mg/L for each flavor reagent to prepare a liquid simulator [32]. According to the contact ratio of actual food to material, a 5 cm × 3 cm film was placed into a 20 mL glass bottle with 20 mL of food simulator [33]. The bottle cap was filled with the polytetrafluoroethylene (PTFE) gasket. After the cap was sealed, the outer part was wound with a sealing film to ensure sealability. Storage temperature was set as 23 • C. After storage for a certain period of time, the film was taken out and repeatedly wiped with test paper until dry. The film was then cut into small pieces and placed in a 20 mL n-hexane sample bottle, and then sealed with the PTFE gasket and sealing film. Then the film was extracted in the ultrasonic water bath pot for 1 h. The ultrasonic wave frequency was 60 Hz, and the temperature was 23 • C. Sample bottles after extraction were stored at room temperature for 24 h. The amount of extraction sample was measured by Chromatography-Mass Spectrometer (GC-MS).

GC-MS Analysis Conditions
GCMS-QP2010 (Shimadzu, Kyoto, Japan) was used to measure the adsorption amount of flavor. The external standard method was used to analyze the measured flavor quantitatively. The GC-MS analysis conditions were as follows: The chromatographic column was DB-5 quartz capillary column (30 m × 0.25 mm × 0.25 µm), the detector was hydrogen flame ionization detector, and the injection volume was 1 µL. The inlet temperature was 250 • C, with a shunt ratio of 20:1, and the ion source temperature was 200 • C. The column temperature was programmed to an initial temperature of 40 • C, which was maintained for 3 min. Column was heated up to 150 • C at the rate of 4 • C/min, and then maintained at 150 • C for 1 min. Then it was heated up to 250 • C at a rate of 8 • C/min and maintained at 250 • C for 6 min.

Data Analysis
In order to better analyze the amount of flavor adsorbed by packaging film, the data measured by GC-MS need to be standardized. All the data were converted by Equation (1) to obtain the adsorption amount per unit mass of packaging film [3,10]: where M i,t is the adsorption amount of flavor i per unit mass of packaging film at time t, mg/g; C i is the concentration (mL/mL) of flavor i in the extract solution (n-hexane) detected by GC-MS; V E is the volume of the extract solution (n-hexane), mL; ρ i is the density (mg/mL) of flavor i and M P (g) is the quality of the adsorption packaging film.

Modeling
The molecular models of five flavor organic molecules and PE were constructed with the Materials Visualizer module of Materials Studio 8.0 (MS 8.0). Generally, setting a higher polymerization degree for PE can simulate a more realistic diffusion of molecules in PE [27]. However, considering the computational complexity, a polymerization degree of 160 was selected to build the PE chain (PE 160). After the competition of modeling, the above models were geometrically optimized by using Forcite module with the "Smart Minimizer" method. The amorphous cells containing the molecules and PE chains were built to simulate the composition of the sample prepared in this study (as listed in Table 1). In order to analyze the diffusion of O 2 and H 2 O in LDPE films, three O 2 or ten H 2 O were placed into the cell structure of PE 160. The simulation cells with the initial density of 0.923 g/cm 3 were constructed according to the actual density of LDPE film. Figure 1 shows an example of the simulation cells.

Simulations
The main procedures of simulation included geometry optimization, annealing and dynamic simulations, using the COMPASS II force field in all stages which is a newer version of the COMPASS force field and mainly used to simulate polymers, metals, metal oxides and small molecules [34]. The Coulomb force was evaluated by the Ewald summation method with an accuracy of 10 −3 kcal/mol, and the van der Waals force was computed by the atom-based summation method with a spline cutoff distance of 12.5 Å. The geometric optimization procedure was used to minimize the energy of the cell. The isothermal-isobaric ensemble (NPT) was selected to perform in the molecular dynamics simulations. Compared with other ensembles, NPT simulation is the most relaxed ensemble in the system [35]. The pressure of annealing simulation was set as 0.1 MPa controlled by Berendsen method. The initial temperature was 300 K, and the mid-cycle temperature was 600 K. Three annealing cycles were performed, and the total number of steps was 3.6 × 10 5 . After the annealing, the internal stress of cell decreased, and the unreasonable structure during cell construction was eliminated [36]. Then, the NPT MD simulation was performed for 1000 ps to reach balance and the other constant volume and temperature ensemble (NVT) MD simulation of 100 ns was performed to ensure the convergence of the MD simulation for diffusion analysis. It was found that the diffusion coefficients calculated by allowing only the diffusing species to move, but keeping fixed the rest of the system can be affected by significant errors [37,38]. Therefore, during all the simulations in this study, all the atoms, including the atoms of five flavor molecules, O2, H2O, and LDPE, were allowed to move (see Figure 2).

Simulations
The main procedures of simulation included geometry optimization, annealing and dynamic simulations, using the COMPASS II force field in all stages which is a newer version of the COMPASS force field and mainly used to simulate polymers, metals, metal oxides and small molecules [34]. The Coulomb force was evaluated by the Ewald summation method with an accuracy of 10 −3 kcal/mol, and the van der Waals force was computed by the atom-based summation method with a spline cutoff distance of 12.5 Å. The geometric optimization procedure was used to minimize the energy of the cell. The isothermal-isobaric ensemble (NPT) was selected to perform in the molecular dynamics simulations. Compared with other ensembles, NPT simulation is the most relaxed ensemble in the system [35]. The pressure of annealing simulation was set as 0.1 MPa controlled by Berendsen method. The initial temperature was 300 K, and the mid-cycle temperature was 600 K. Three annealing cycles were performed, and the total number of steps was 3.6 × 10 5 . After the annealing, the internal stress of cell decreased, and the unreasonable structure during cell construction was eliminated [36]. Then, the NPT MD simulation was performed for 1000 ps to reach balance and the other constant volume and temperature ensemble (NVT) MD simulation of 100 ns was performed to ensure the convergence of the MD simulation for diffusion analysis. It was found that the diffusion coefficients calculated by allowing only the diffusing species to move, but keeping fixed the rest of the system can be affected by significant errors [37,38]. Therefore, during all the simulations in this study, all the atoms, including the atoms of five flavor molecules, O 2 , H 2 O, and LDPE, were allowed to move (see Figure 2).

Interaction Energy between Diffusion Molecules and PE Chain
The interaction energy (Eint), indicating the intensity of interaction between the diffusion molecule and polymer, is derived according to the following equation:

Interaction Energy between Diffusion Molecules and PE Chain
The interaction energy (E int ), indicating the intensity of interaction between the diffusion molecule and polymer, is derived according to the following equation: where E total is the total energy of the system, E PE is the energy of PE chain, E diff,i is the energy of diffusion molecule i in the PE chain. A negative E int value is corresponding to stable interaction between the components [39,40]. More negative E int indicates a stronger interaction in the system. In this study, there was only physical bonding between diffusion molecules and polymer chains, and no chemical reaction and new chemical bonds involved. The non-bonding interactions between diffusion molecules and polymer chains, and the interatomic potential for the interactions among LDPE atoms and diffusion molecules were considered. The interaction energies of D-limonene, myrcene, ethyl hexanoate, 2-nonanone, linalool, O 2 and H 2 O with PE chain were calculated by MD simulation and were listed in Table 2. The stronger the interaction between diffusion molecules and polymer chains will result in the higher the energy barrier to be overcome and the smaller the corresponding diffusion coefficient. Therefore, if only take the effect of interaction energy on the diffusion into account, the diffusion coefficients of the diffusion molecules in PE should decrease in the following order: However, besides the interaction energy, molecular size and free volume of polymer cell are also important factors affecting the diffusion of molecules, which will be further discussed in a later section.

Calculation of D by Prediction Model
Diffusion coefficient is the most important analytical parameter in adsorption of flavor organic molecules by polymers. The classical prediction model of diffusion coefficient, Brandsch model (Equation (3)), was widely used to predict the D of molecules in polymers [41].
Brandsch model introduced two simplification factors, A P and τ, to describe the properties of polymer. Brandsch model was recognized as a general formula for characterizing the diffusion coefficient of molecule in plastic packaging film [42]: where D is the diffusion coefficient of the diffusion molecule in polymer (cm 2 /s), M r is the relative molecular weight of the diffusion molecule, A P is the characteristic parameter of polymer, and T is the absolute temperature (K). For different polymer matrix, the value of A P is different, which is mainly derived from a large number of experimental results. A P is usually a function of temperature, and it can be calculated from the below equation: where A P * is the adiabatic term of polymer, and τ is the deviation of diffusion activation energy from reference activation energy in polymer. For LDPE films, A P * equals 11.5, and τ equals 0 [43].

Experimental Determination of D
The adsorption of five flavor organic molecules in LDPE film at 23 • C was studied by film immersion experiment. As shown in Figure 3, for all five flavor organic molecules, adsorption equilibrium can be reached within 10 days. After reaching the adsorption equilibrium, the total amount of D-limonene adsorbed is the largest, i.e., 18.4212 mg/g, followed by myrcene, which is 7.5019 mg/g. For linalool, it has the lowest adsorption amount (3.2382 mg/g) among five flavor organic molecules. The adsorption equilibrium amounts of the organic molecules by LDPE films follow the order: D-limonene > myrcene > 2-nonanone > ethyl hexanoate > linalool. Although molecular weights of D-limonene and myrcene are the same, the adsorption amounts of the two molecules are quite different. fitted curve and the experimental scatter points is shown in Figure 3. Table 4 and Figure 3 show that the diffusion of five flavor organic molecules in LDPE well conforms to Fick's second law.   When the adsorption amount of the flavor organic molecule was known, the Fick's second law diffusion model can be used to calculate the diffusion coefficient [44,45]: where M i,t is the adsorption amount of flavor i adsorbed by films at t (s) time, mg/g; M i,e is the adsorption amount of flavor i at equilibrium time, mg/g; and L P is the thickness of films, cm. The experimental diffusion coefficient, D can be derived by fitting the Equation (5) with least square method. The diffusion coefficients of five flavor organic molecules in LDPE were fitted with the experimental data by using MATLAB Software and listed in Table 4. The comparison between the fitted curve and the experimental scatter points is shown in Figure 3. Table 4 and Figure 3 show that the diffusion of five flavor organic molecules in LDPE well conforms to Fick's second law.

Molecular Simulation Results
The diffusion coefficient in MD simulation was calculated by analyzing the mean square displacement (MSD) curve. Four methods can be used to analyze the MSD curve, i.e., Einstein method, Green-Kubo method, differential-zone variational method, and cluster analysis method [46].
Einstein method is a classical method, which correlates diffusion coefficient D with MSD of molecules [42]. The formula is as follows: where N is the number of diffusion molecule i in the system, t is the diffusion time, and r i (0) and r i (t) are the diffusion position vector of diffusion molecule i at the time of 0 and t, respectively. Figure 4 shows the MSD curves of five flavor organic molecules in PE derived from MD simulation. All MSD curves contain the linear part, unstable vibration part and the noisy part after a long simulation time. The non-linear and noisy parts usually encountered at the end of the curves may be attributed to the statistical error after longer simulation times. Reynier proposed that there were three diffusion modes of molecules in polymers: Crawling diffusion mode, jumping diffusion mode, and dual mode diffusion [47,48]. The flavor molecules (D-limonene, myrcene, ethyl hexanoate, 2-nonanone and linalool) contain both rigid and flexible groups; therefore, their diffusion may follow the dual mode diffusion (both the jumping and crawling and diffuse modes involved). Figure 4 shows that flavor organic molecules, because of their large size, cannot jump easily and diffuse in polymers like small gas molecules, and the linear part and unstable vibration part of their MSD curves show dual mode diffusion. In order to analyze molecular diffusion behavior more clearly, the MSD curves of five organic molecules are calculated by Log10 and fitted by segments (as shown in Figure 5). In Figure 5, m is the slope of the fitting straight line, which is anomalous (non-Einstein) behavior when m < 1, and Einstein behavior when m is close to 1 [49]. As can be seen from the figure, at the back end of the simulation time, five organic molecules transform from anomalous diffusion to Einstein diffusion. During this time, the diffusion coefficient can be calculated by Einstein Equation (6). The diffusion coefficients of five flavor molecules (D-limonene, myrcene, ethyl hexanoate, 2-nonanone and linalool) in PE are 9.31 × 10 −11 cm 2 /s, 7.90 × 10 −11 cm 2 /s, 6.90 × 10 −11 cm 2 /s, 1.56 × 10 −10 cm 2 /s, 1.72 × 10 −10 cm 2 /s, respectively. like small gas molecules, and the linear part and unstable vibration part of their MSD curves show dual mode diffusion. In order to analyze molecular diffusion behavior more clearly, the MSD curves of five organic molecules are calculated by Log10 and fitted by segments (as shown in Figure 5). In Figure 5, m is the slope of the fitting straight line, which is anomalous (non-Einstein) behavior when m < 1, and Einstein behavior when m is close to 1 [49]. As can be seen from the figure, at the back end of the simulation time, five organic molecules transform from anomalous diffusion to Einstein diffusion. During this time, the diffusion coefficient can be calculated by Einstein Equation (6). The diffusion coefficients of five flavor molecules (D-limonene, myrcene, ethyl hexanoate, 2-nonanone and linalool) in PE are 9.31 × 10 −11 cm 2 /s, 7.90 × 10 −11 cm 2 /s, 6.90 × 10 −11 cm 2 /s, 1.56 × 10 −10 cm 2 /s, 1.72 × 10 −10 cm 2 /s, respectively.  like small gas molecules, and the linear part and unstable vibration part of their MSD curves show dual mode diffusion. In order to analyze molecular diffusion behavior more clearly, the MSD curves of five organic molecules are calculated by Log10 and fitted by segments (as shown in Figure 5). In Figure 5, m is the slope of the fitting straight line, which is anomalous (non-Einstein) behavior when m < 1, and Einstein behavior when m is close to 1 [49]. As can be seen from the figure, at the back end of the simulation time, five organic molecules transform from anomalous diffusion to Einstein diffusion. During this time, the diffusion coefficient can be calculated by Einstein Equation (6). The diffusion coefficients of five flavor molecules (D-limonene, myrcene, ethyl hexanoate, 2-nonanone and linalool) in PE are 9.31 × 10 −11 cm 2 /s, 7.90 × 10 −11 cm 2 /s, 6.90 × 10 −11 cm 2 /s, 1.56 × 10 −10 cm 2 /s, 1.72 × 10 −10 cm 2 /s, respectively.

Comparison among Experiment, Prediction Model and MD Simulation Results
The differences among prediction model diffusion coefficient (Dpre), experimental diffusion coefficient (Dexp) and simulated diffusion coefficient (Dsim) were compared as listed in Table 5. The Dpre is higher than that obtained by experiments and MD simulation, which is similar to the results reported by Wang et al. [50] who investigated the diffusion of small molecules in PET. This is because Brandsch prediction model is based on the assumptions given under deteriorating conditions, and it only considers the properties of polymer substrates, the molecular weight of flavor molecules and the temperature, but does not take the interaction between molecules into account.
The diffusion coefficients obtained by MD simulation are close to those of experimental results, but are slightly larger than those values. This difference may be due to: (1) the setting of finite chain length for PE means that many relaxed free ends make the redistribution of free volume fraction easier and lead to the increase of diffusion coefficient. The simulation results reported by Wang et al. [50] showed a good agreement with that obtained from experiments, this may be due to the PET models built in their study had a polymerization degree of 100, which is close to that of actual PET. Typically, the actual LDPE has a polymerization degree of >1000. Building an LDPE model with such a high polymerization degree will lead to a very large model and an unbearable expense of calculation time in the subsequent MD simulation. (2) The amorphous PE cell was constructed for MD simulation, but the actual PE is an incomplete amorphous polymer which contains a certain crystalline phase. The crystallization zone will hinder the diffusion of flavor organic molecules. (3) MD simulation only considers the diffusion of flavor organic molecules in PE materials, and does not consider the distribution of flavor organic molecules from food simulation solution to PE material, which will lead to the deviation between MD simulation and experimental values. However, it can

Comparison among Experiment, Prediction Model and MD Simulation Results
The differences among prediction model diffusion coefficient (D pre ), experimental diffusion coefficient (D exp ) and simulated diffusion coefficient (D sim ) were compared as listed in Table 5. The D pre is higher than that obtained by experiments and MD simulation, which is similar to the results reported by Wang et al. [50] who investigated the diffusion of small molecules in PET. This is because Brandsch prediction model is based on the assumptions given under deteriorating conditions, and it only considers the properties of polymer substrates, the molecular weight of flavor molecules and the temperature, but does not take the interaction between molecules into account. The diffusion coefficients obtained by MD simulation are close to those of experimental results, but are slightly larger than those values. This difference may be due to: (1) the setting of finite chain length for PE means that many relaxed free ends make the redistribution of free volume fraction easier and lead to the increase of diffusion coefficient. The simulation results reported by Wang et al. [50] showed a good agreement with that obtained from experiments, this may be due to the PET models built in their study had a polymerization degree of 100, which is close to that of actual PET. Typically, the actual LDPE has a polymerization degree of >1000. Building an LDPE model with such a high polymerization degree will lead to a very large model and an unbearable expense of calculation time in the subsequent MD simulation. (2) The amorphous PE cell was constructed for MD simulation, but the actual PE is an incomplete amorphous polymer which contains a certain crystalline phase. The crystallization zone will hinder the diffusion of flavor organic molecules. (3) MD simulation only considers the diffusion of flavor organic molecules in PE materials, and does not consider the distribution of flavor organic molecules from food simulation solution to PE material, which will lead to the deviation between MD simulation and experimental values. However, it can be concluded that the simulation results are close to the experimental ones by taking the errors of these numbers (listed in Table 5) into account. In addition, the experimental process was more complex than the predictive model and MD simulation. Because the immersion of food simulation solution, i.e., ethanol and water, may be easier to first diffuse into the material than flavor molecules, thus, maybe hinder the diffusion of flavor molecules.

Fraction of Free Volume
The size and shape of a free volume play an important role in the diffusion behavior of diffusion molecules in polymers. Generally, the larger the free volume fraction of polymers results in the larger the diffusion coefficient. When calculating the free volume of the system, the hard sphere probe model in MS 8.0 was used to analyze the kinetic radius of the diffusion molecule. The free volume and occupied volume of the cell can be obtained when the molecular probe with a certain radius R P moves on the van der Waals surface. The ratio of free volume to simulated cell volume is defined as FFV [51]: where, V F is the free volume of the polymer, V O is the occupied volume of the polymer, and V S is the total volume of the polymer. The probe radii of 1.0 Å, 1.5 Å, 2.0 Å, 2.7 Å were selected to calculate the free volume of PE cells. Free volume morphologies with the different probe radius in PE cell are shown in Figure 6. Table 6 shows that with the increase of R P , the FFV of PE cell decreases. The free volume measured by 2.7 Å can meet the space requirement of jump diffusion for O 2 and H 2 O. It also means that for larger flavor organic molecules, there are fewer free holes to jump. The polymer chains can wriggle themselves, and when their wriggling creates enough jumping space, flavor organic molecules will jump from one hole to other holes.
where, VF is the free volume of the polymer, VO is the occupied volume of the polymer, and VS is the total volume of the polymer. The probe radii of 1.0 Å, 1.5 Å, 2.0 Å, 2.7 Å were selected to calculate the free volume of PE cells. Free volume morphologies with the different probe radius in PE cell are shown in Figure 6. The diffusion coefficients of O 2 and H 2 O before and after PE adsorption are listed in Table 7. Table 7 indicates that the diffusion of O 2 before and after PE adsorption is stronger than that of H 2 O. This is inconsistent with the results of interaction energy analysis (Section 4.1). Because the difference between the interaction energies of O 2 and H 2 O is relatively small, thus, effects of the interaction energies on the diffusion of both two molecules are not decisive. However, the size of O 2 molecule is smaller than that of H 2 O molecule; therefore, O 2 molecule is easier to diffuse in PE. The diffusion coefficients of O 2 and H 2 O in adsorbed PE are all smaller than those in pure PE. In order to better understand the influence of flavor organic molecules on the diffusion of O 2 and H 2 O in PE, the FFV of PE and the interaction energies of O 2 or H 2 O with flavor organic molecules were further analyzed.
The FFV of PE cells, before and after flavor adsorption, was calculated with a probe radius of 1.0 Å and listed in Table 8. The interaction energies of D-limonene, myrcene, ethyl hexanoate, 2-nonanone and linalool with O 2 or H 2 O were calculated by MD simulation and were listed in Table 9. Table 8 shows that the FFV of PE cell decreased with the loading of O 2 , while increased with the loading of H 2 O. This is because the alkane chain has strong hydrophobicity, and the free volume becomes larger after loading H 2 O. When flavor molecules were diffused in PE, the FFV of PE cell decreased, which led to the decrease in the diffusion coefficients of O 2 and H 2 O. long-term MD simulation, the diffusion of molecules in random motion is affected by the steric hindrance of PE chain, resulting in a decrease of diffusion coefficient. Compared with Figure 4, O2 or H2O compete with the flavor molecules for free holes in PE, which restricts the diffusion of O2, H2O and flavor molecules. After adsorption of flavor molecules, O2 and H2O diffuse faster than all flavors, and their MSD curves increase more rapidly. The reason is that O2 and H2O are smaller than flavor molecules and need less diffusion space. Moreover, the effect of flavor molecules on O2 is small; therefore, O2 diffusion is faster than H2O.  The diffusion coefficients of O2 and H2O before and after PE adsorption are listed in Table 7. Table 7 indicates that the diffusion of O2 before and after PE adsorption is stronger than that of H2O. This is inconsistent with the results of interaction energy analysis (Section 4.1). Because the difference between the interaction energies of O2 and H2O is relatively small, thus, effects of the interaction energies on the diffusion of both two molecules are not decisive. However, the size of O2 molecule is smaller than that of H2O molecule; therefore, O2 molecule is easier to diffuse in PE. The diffusion coefficients of O2 and H2O in adsorbed PE are all smaller than those in pure PE. In order to better understand the influence of flavor organic molecules on the diffusion of O2 and H2O in PE, the FFV of PE and the interaction energies of O2 or H2O with flavor organic molecules were further analyzed.
The FFV of PE cells, before and after flavor adsorption, was calculated with a probe radius of 1.0 Å and listed in Table 8. The interaction energies of D-limonene, myrcene, ethyl hexanoate, 2nonanone and linalool with O2 or H2O were calculated by MD simulation and were listed in Table 9.   In addition to the FFV, the interaction between O 2 or H 2 O molecule and flavor organic molecules also affects the diffusion coefficients of O 2 or H 2 O. Table 9 shows the interaction energies between O 2 or H 2 O molecule and five flavor organic molecules. All the interaction energies were negative, which indicates the presence of the flavor organic molecules will hinder the diffusion of O 2 and H 2 O molecules, and thus, lead to the decrease in the diffusion coefficients (Table 8). Besides, the interaction energies between H 2 O and five flavor organic molecules were stronger than those between O 2 and five flavor organic molecules, especially for ethyl hexanoate, 2-nonanone and linalool. The equilibrium adsorption configurations of H 2 O molecule with five flavor organic molecules were shown in Figure 9. The H atom in H 2 O molecule is positively charged and expected to absorb near to the negatively charged C atoms in D-limonene and myrcene (Figure 9a,b). O atom of higher electronegativity is contained in the ethyl hexanoate, 2-nonanone and linalool molecules, which can form the hydrogen bond with the H atom in H 2 O (Figure 9c-e). Therefore, the interaction between H 2 O and ethyl hexanoate, 2-nonanone and linalool molecules are stronger than those of D-limonene and myrcene. Therefore, in a word, the presence of the flavor organic molecules results in greater hindrance to the diffusion of H 2 O compared with the case of O 2 , which can further explain the greater decrease in the diffusion coefficient of H 2 O in the presence of flavor organic molecules (Table 7).  Table 8 shows that the FFV of PE cell decreased with the loading of O2, while increased with the loading of H2O. This is because the alkane chain has strong hydrophobicity, and the free volume becomes larger after loading H2O. When flavor molecules were diffused in PE, the FFV of PE cell decreased, which led to the decrease in the diffusion coefficients of O2 and H2O.
In addition to the FFV, the interaction between O2 or H2O molecule and flavor organic molecules also affects the diffusion coefficients of O2 or H2O. Table 9 shows the interaction energies between O2 or H2O molecule and five flavor organic molecules. All the interaction energies were negative, which indicates the presence of the flavor organic molecules will hinder the diffusion of O2 and H2O molecules, and thus, lead to the decrease in the diffusion coefficients (Table 8). Besides, the interaction energies between H2O and five flavor organic molecules were stronger than those between O2 and five flavor organic molecules, especially for ethyl hexanoate, 2-nonanone and linalool. The equilibrium adsorption configurations of H2O molecule with five flavor organic molecules were shown in Figure 9. The H atom in H2O molecule is positively charged and expected to absorb near to the negatively charged C atoms in D-limonene and myrcene (Figure 9a,b). O atom of higher electronegativity is contained in the ethyl hexanoate, 2-nonanone and linalool molecules, which can form the hydrogen bond with the H atom in H2O (Figure 9c−e). Therefore, the interaction between H2O and ethyl hexanoate, 2-nonanone and linalool molecules are stronger than those of D-limonene and myrcene. Therefore, in a word, the presence of the flavor organic molecules results in greater hindrance to the diffusion of H2O compared with the case of O2, which can further explain the greater decrease in the diffusion coefficient of H2O in the presence of flavor organic molecules (Table 7).

Conclusions
The diffusion of D-limonene, myrcene, ethyl hexanoate, 2-nonanone and linalool in LDPE were investigated by two-sided immersion experiment, Brandsch prediction model and MD simulation.

Conclusions
The diffusion of D-limonene, myrcene, ethyl hexanoate, 2-nonanone and linalool in LDPE were investigated by two-sided immersion experiment, Brandsch prediction model and MD simulation. The diffusion of O 2 and H 2 O in PE before and after the adsorption of the flavor substances were also studied. The following conclusions were drawn: (1) the immersion experiment showed the diffusion of five flavor molecules in LDPE well followed Fick's second law; (2) MD simulation indicated the dual mode diffusion mechanism of the flavor molecules in LDPE; (3) the diffusion coefficients from MD simulation are close to those from the experimental determination, but are slightly larger than those values; (4) the presence of the flavor organic molecules hindered the diffusion of O 2 and H 2 O, which can be well explained from FFV and interaction energy calculation results derived from MD simulation.