INTRODUCTION
Currently, cancer treatment has received great attention and considerable achievements have been obtained on the drug safety, drug delivery efficacy and survival rates so that some certain tumors have recently been cured [1]. Nevertheless, cancer has been remained the second cause of death worldwide and its global incidence continues to increase [2]. The World Health Organization (WHO) announced on 30th May 2017 that cancer resulted in 7.9 million deaths in 2007 which was about 13% of all deaths, and 8.2 million deaths worldwide in 2012. Also, it was expected that the annual number to be increased from 14.1 million in 2012 to 21.6 million in 2030 [3]. For this reason, numerous cancer treatment methods (such as chemotherapy) are carried out in order to decline the cancer related deaths [4].
The pharmacological activities of most drugs are decreased due to their rapid metabolisms, low bioavailability, poor bio-distributions in pathological parts or little intracellular penetrations which lead to body resistance against cancer treatments [5]. Encapsulation of drugs into nanocarriers can somewhat overcome these problems [6-10]. Nonetheless, among nanomedicines reached the market (particularly in oncology), few of them are in advanced phase III clinical trials [11]. One reason is uncontrolled drug release (called burst release) which is due to the rapid drug leakage from the carrier. Accordingly, various drug design methodologies have been established in order to improve existing conventional chemotherapy procedures by alleviating undesirable side effects that affect healthy tissues [12]. Such an important objective will be attained by targeted pharmaceutics delivery to tumor/cancer cells [13]. Sustained drug delivery attempts to prolong the inter-dose time for a drug in chronic diseases but preserves almost constant plasma medication concentration in the therapeutic range [14]. A sustained release formulation has a lot of potential clinical advantages like decreased side effects of therapeutics along with low toxic thresholds, lower cost and more patient compliance for frequent, hard or invasive administrations [15].
Polymers are commonly used as effective drug delivery systems because of their pharmacokinetic benefits [16-23]. Poly(ethylene glycol), PEG, is a valuable polymer which is widely utilized by the pharmaceutical companies for drug carrying purposes [24]. PEG is a hydrophilic polyether that is biocompatible, non-toxic, non-immunogenic, bio-inert, highly bioavailable and extremely soluble in water that has been approved by Food and Drug Administration (FDA) for usage in medicinal and all biomedical applications [25]. Polylactic acid (PLA) is a thermoplastic polymer that is very extensively utilized for the drug encapsulation in controlled delivery applications, fabrication of cardiovascular devices, biodegradable sutures and various implants used to fix fractures [26]. It is notable that advanced drug carriers that are based on biodegradable polyesters are mainly prepared from PLA and poly(lactic-co-glycolic acid) (PLGA) polymers that are biocompatible and have shown very low toxicity and immunogenicity [27,28].
Among chemotherapeutic anticancer drugs, nucleoside analogs are one of the most significant pharmaceutics [29]. Gemcitabine (2′,2′-difluoro-2′-deoxycytidine) (GEM) is applied as a first-line medication for metastatic and advanced pancreatic cancer [30]. GEM is a very hydrophilic FDA-approved anticancer drug for usage in a wide range of human solid tumors such as non-small cell lung, bladder, pancreatic, breast, colon, ovarian, hepatocellular and cervical cancers [31]. Its anticancer activity includes inhibiting the DNA replication which results in apoptosis of tumor cells [32]. Nevertheless, GEM has a high hydrophilicity and low molecular weight which is simply degraded by cytidine deaminase and leads to its poor cellular penetration [33]. Consequently, in order to enhance the GEM cellular uptake, numerous drug delivery systems have been designed and investigated [34].
Molecular dynamics (MD) simulations are frequently performed to get insight about physical properties of various systems [35-39]. The MD simulations can describe intermolecular interactions occurred between drug molecules and polymeric chains in polymeric nanocomposite drug carriers [40,41]. Nonetheless, computational simulations have not broadly been carried out on the transport of anticancer drugs by means of polymer nanocomposites as drug delivery platforms [42-44].
There are a lot of experimental works on drug delivery shuttles [45,46] but literature review indicates that the “PEG-PLA-GEM/CNT derivatives” systems designed in this work have not yet been investigated elsewhere. Similar papers report using PEG/single walled carbon nanotubes (SWCNT) [47] and PEG/functionalized multiwalled carbon nanotubes (MWCNT) for doxorubicin delivery [48], aqueous-core PEG-coated PLA nanocapsules for gemcitabine hydrochloride (GEM) delivery [49], PLA-mPEG copolymer/Fe3O4-SWCNT for docetaxel delivery [50] as well as PEGylated oxidized MWCNT modified with angiopep-2 for doxorubicin delivery [51]. Also, as carbon nanotubes and their derivatives are widely used as fillers in the polymeric nanocomposite drug delivery systems [52], it was decided to use these inorganic nanomaterials to evaluate their effects on the drug delivery efficacy of the designed systems.
Herein, MD simulations were accomplished at two temperatures (298.15 and 310.15 K) in both gas phase and water solvent to investigate the PEG-PLA polymeric nanocomposite systems filled with various nanotubes as inorganic fillers including carbon nanotube and its derivatives to be used for the delivery of the GEM anticancer drug molecules. The diffusion of GEM into the PEG-PLA and PEG-PLA-FA systems was examined using the MD simulations. Moreover, the interaction energies between diffused molecules and polymeric chains, FV, FFV, X-ray diffraction patterns, radius of gyration, solubility parameters and the diffusion coefficients were evaluated and compared with each other to find the most suitable drug carrier for the GEM anticancer drug. It was found that the PEG-PLA-GEM-FA-CN could be used at 310.15 K in water as the best system to deliver the GEM drug.

SIMULATION METHODOLOGY
The condensed-phase optimized molecular potentials (COMPASS) force field (FF) was employed in both of the minimization and dynamics steps. All of the atomic charges were determined by the force field assigned method. The COMPASS FF is based on the Polymer Consistent Force Field (PCFF) [53] which is the first FF that has been certified and parameterized through condensed-phase criteria, ab initio and empirical data existing for a number of molecules. In addition to the bonded and non-bonded interactions which are applied in other FFs, the cross coupling contacts are also considered by COMPASS; thus it can exactly estimate numerous properties for different materials predominantly for polymers [54].
All of particles were built by means of the molecule modelling tool in Materials Studio software [55]. The amorphous cells contained five PEG chains with twenty repeating units, one PLA chain with twenty repeating units, one filler nanoparticle (armchair-(6,6) CNT, CN, CP, SiN or SiP nanotube with 8.14 Ǻ diameter and 4.92 Ǻ length), three GEM drug molecules and one folic acid molecule in order to assess the impacts of these systems loaded with different types of fillers on the drug delivery efficacy. Also, 200 water (H2O) molecules were added into the amorphous cells to examine the effect of water on the drug delivery. The amorphous cells were run by the MD simulations at two temperatures (298.15 and 310.15 K) and 1 atm pressure.
The energy minimizations of all particles as well as amorphous cells were performed using the smart minimizer algorithm which combines conjugate gradient, steepest descent and Newton minimization algorithms. All minimizations were conducted in a cascading manner for 4,000,000 iterations with energy and force resolutions of 0.001 kcal/mol and 0.5 kcal/mol/Ǻ, respectively, until reaching the relaxation. Then, dynamics simulations were done to fully relax the simulation systems. The NVT simulation was initially performed for 30 ns (time step: 1 fs, total simulation time: 3×104 ps) and then, the NPT simulation was carried out for 30 ns (time step: 1 fs, total simulation time: 3×104 ps). In the course of the relaxation step, non-bond interactions (van der Waals and electrostatic potentials) were estimated by Ewald summation method. Afterward, the MD simulations with NVT ensemble were accomplished at the relaxed density for 5 ns (time step: 1 fs, total simulation time: 5×103 ps). The last 1 ns trajectories were used to extract data and analyze energetic, dynamical and structural features of cells. The velocity Verlet algorithm was applied to employ the classical motion equation with a 1 fs time step. The Berendsen barostat and thermostat procedures (with a 0.1 ps decay constant) were used in all simulations to preserve the temperature and pressure at the certain values during the simulations. All analyses were done using Forcite module analysis existing in the Materials Studio software. The cutoff for achieving RDFs was 20 Ǻ with 0.02 Ǻ interval.

RESULTS AND DISCUSSIONS
Equilibrium and relaxation of amorphous cells
In this study, the physicochemical and structural properties of PEG-PLA polymeric nanocomposites containing various nanotube fillers including carbon nanotube (CNT), carbon nitride (CN), carbon phosphide (CP), silicon nitride (SiN) and silicon phosphide (SiP) have been evaluated by the MD simulations at two temperatures (298.15 and 310.15 K) in both gas phase and water solvent. These systems have been designed as the drug delivery shuttles for the GEM anticancer drug. It should be noted that the temperature 310.15 K (37 °C) and H2O solution are used in order to investigate the drug delivery process in an environment similar to the human body. Moreover, folic acid (FA) has been added into the systems in order to examine its effect on their properties as well as drug delivery efficacies. A summary of the MD simulation data acquired at 1 atm and two temperatures (310.15 and 298.15 K) are presented in Tables 1 and S1, respectively, for the PEG-PLA-GEM and PEG-PLA-GEM-FA systems containing different nanotube fillers. Figs. 1 and S1 indicate the snapshots of the relaxed PEG-PLA-GEM and PEG-PLA-GEM-FA systems filled with diverse nanotubes at 310.15 and 298.15 K, respectively. All of the systems have reached relaxation at the two temperatures and this is evident from the negligible variations in their temperatures, root mean square displacement (RMSD) values as well as kinetic, potential and non-bond energies by the simulation time. It is found from Figs. 2 and S2 that the temperature and RMSD are almost constant for all of the PEG-PLA-GEM and PEG-PLA-GEM-FA drug delivery systems containing different nanotubes at both temperatures (298.15 and 310.15 K).
Figs. 3 and S3 exhibit the variations in the energies of the PEG-PLA-GEM and PEG-PLA-GEM-FA systems at 310.15 and 298.15 K, respectively. The negligible changes in the kinetic, potential and non-bond energies of the PEG-PLA-GEM and PEG-PLA-GEM-FA systems versus simulation time verify that all of the cells have been relaxed. Tables 2 and S2 give the average energies values for all of the systems examined at 310.15 and 298.15 K, respectively. It is found that at 298.15 K, addition of the CP nanotube into both kinds of PEG-PLA-GEM and PEG-PLA-GEM-FA cells affords negative non-bond energies of -15 and -110 kcal/mol for the PEG-PLA-GEM-CP and PEG-PLA-GEM-FA-CP cells, respectively. However, positive values are achieved for the potential and kinetic energies of all systems at 298.15 K.
At 310.15 K, all of the systems illustrate negative non-bond energies but positive kinetic energies. However, only positive potential energies are measured for the cell loaded by CNT and CN nanotubes while other systems have negative potential energies with the systems containing CP and SiP nanotubes display the most negative potential energies. It may be concluded that the cell showing the smallest energy can be the most appropriate system but in order to select the most desirable drug delivery system, other features of the polymeric nanocomposites must be evaluated and compared. Hence, the relationship between types of polymeric nanocomposite systems and their drug delivery capacity will be found in following sections in order to find the most suitable cell for application in the delivery of the GEM drug.
It is observed in Tables 1 and S1 that the cell size and therefore the cell volume are increased at both temperatures (298.15 and 310.15 K) by adding nanotubes into the PEG-PLA-GEM and PEG-PLA-GEM-FA cells. In addition, greater cell volumes are obtained for the cells containing CP, SiN and SiP fillers instead of the CNT which is related to the substitution of larger heteroatoms by some of the carbon atoms in the CNT structure. However, among the systems filled with nanotubes, the cell size and cell volume are the smallest for the cells composed of the CN nanotube. That is, such variations in cell dimensions can be attributed to smaller size of the CN but greater sizes of the CP, SiN and SiP fillers compared to that of the CNT. The volumes of the CNT, CN, CP, SiN and SiP particles are equal to 1086.44, 907.91, 1522.39, 1925.08 and 2165.32 Ǻ3, respectively, approving the CN and SiP occupy the smallest and greatest volumes, respectively. Fig. 4 demonstrates the structures and volumes (displayed by blue color) for the CNT, CN, CP, SiN and SiP species. It can be seen that the volumes of the nanotubes change in the order CN<CNT<CP<SiN<SiP.
Comparing the cell sizes and the cell volumes at the two temperatures of 298.15 and 310.15 K reveals that larger values are attained at higher temperature. This can be related to much greater kinetic energies of molecules at higher temperature leading to their greater movements, less interactions and larger intermolecular distances. Also, the presence of H2O molecules in the systems interacting with the functional groups of all species results in increasing the cell sizes and cell volumes. For example, the cell sizes and the cell volumes for the PEG-PLA-GEM system are 20.98 Å and 9234.57Å3 at 298.15 K but they are 25.19 Å and 15983.96 Å3 at 310.15 K. Similar results are obtained for all other cells.
The density values of the PEG-PLA-GEM and PEG-PLA-GEM-FA systems at 298.15 K are measured to be 1.39 and 1.36 g/cm3 but they are increased after addition of nanotubes into the cells. Moreover, the densities of the PEG-PLA-GEM nanocomposites are larger than those of their related PEG-PLA-GEM-FA cells which can be due to the greater cell volumes of the latter. Among PEG-PLA-GEM systems, the lowest and the highest densities have been achieved for the PEG-PLA-GEM and PEG-PLA-GEM-SiP cells, respectively. Similar results are observed for the PEG-PLA-GEM-FA nanocomposites. These findings can be ascribed to the variations in the total mass and volume of each cell so that the density indicates a reverse relationship with the volume.
A comparison of the density values for the PEG-PLA-GEM and PEG-PLA-GEM-FA systems containing different nanotubes at the two temperatures of 298.15 and 310.15 K confirms that all of the systems have smaller densities at higher temperature. This can be associated with the much larger cell volumes at higher temperature leading to their smaller densities because the density has a reverse relationship with the volume. For instance, the densities for the PEG-PLA-GEM system are 1.39 and 1.17 g/cm3 at 298.15 and 310.15 K, respectively. Similar results are obtained for all other cells.
The density values are in agreement with those reported in literature for the PEG and PLA, e.g., the densities attained for the PEG and PLA were 1.10 [56,57] and 1.27 g/cm3 [58], respectively. Consequently, the simulated densities measured for the PEG-PLA-GEM and PEG-PLA-GEM-FA nanocomposite systems are acceptable and this result can be considered as another sign of reaching cell equilibration.

Surface area and free volume
The surface areas of the PEG-PLA-GEM and PEG-PLA-GEM-FA nanocomposites are given in Tables 1 and S1 at 310.15 and 298.15 K, respectively, indicating the surface areas of these systems are increased by adding nanotubes. Also, the surface areas of the PEG-PLA-GEM-CN and PEG-PLA-GEM-FA-CN cells containing CN are the smallest among other nanocomposites. This result can be attributed to the smallest surface area of CN. The surface areas of the CNT, CN, CP, SiN and SiP species are 835.16, 717.37, 1129.42, 1202.16 and 1262.71 Ǻ2, respectively. Moreover, the surface areas of all systems are enhanced by raising the temperature from 298.15 to 310.15 K. This can be attributed to the increase in the intermolecular and interatomic distances at higher temperature due to higher kinetic energies of molecules and atoms as well as the presence of H2O molecules in the systems interacting with the functional groups of all species. Similar results have been reported in literature in which the surface area was increased upon raising the temperature [59-62].
The molecular diffusion in a system can be related to its free volume (FV) value which is obtained by sum of the static holes (produced by chains’ packing) and transient gaps (generated through thermally rearranged chains) in order to provide a low-resistance transport path for the diffusing molecules [43]. Thus, the FV can affect the molecular diffusion in polymeric nanocomposites so that an enhanced FV can facilitate the diffusion [44]. The FV values are achieved using the Connolly surface method in which a hard sphere with a radius of 0.1 nm measures the FVs.
Tables 1 and S1 present the FV values for the PEG-PLA-GEM and PEG-PLA-GEM-FA polymeric nanocomposites at 310.15 and 298.15 K, respectively. Besides, Figs. S4 and 5 illustrate the FVs in these cells by blue color at 298.15 and 310.15 K, respectively. It is found that the FV is increased by addition of nanotubes to the systems. For instance, at 298.15 K, the FVs of the PEG-PLA-GEM, PEG-PLA-GEM-CNT and PEG-PLA-GEM-CN are equal to 3511.85, 4459.98 and 5558.04 Å3, respectively. Moreover, the FV is enhanced through introducing the FA into the nanocomposites. For example, the FV of the PEG-PLA-GEM-FA at 298.15 K is 3905.96 Å3 which is greater than that of the PEG-PLA-GEM. Analogous results are obtained for other cells. This can be attributed to the enhanced intermolecular interactions occurred among the polymeric chains, drug molecules, nanotubes and the FA particles. In addition, the FV is the largest for the systems filled with the CN nanotube. For both of the PEG-PLA-GEM and PEG-PLA-GEM-FA cells composed of different nanotubes, the FV values change by the nanotube type in the order of CN>SiN>CP>SiP>CNT. This result confirms that the total number and strength of the intermolecular interactions in the cells strongly affect the FV values so that the CN can lead to the strongest intermolecular interactions.
The fractional free volume (FFV) is defined as the ratio of the FV to the total cell volume. It is observed that the FFV values vary similar to the FVs. For example, the FFV of PEG-PLA-GEM systems filled by various nanotubes change as CN>SiN>CP>SiP>CNT. Therefore, the highest FFV value of the PEG-PLA-GEM-FA nanocomposite comprised of the CN nanotube may display that the diffusion/transport of the GEM drug in this system happens the most simply and it may be predicted that this nanocomposite affords the greatest diffusion coefficient for the drug molecules. However, in order to determine the most appropriate drug delivery system, further properties of all cells (such as radius of gyration, X-ray diffraction pattern, radial distribution function, solubility parameter and mean square displacement) must be examined and compared with each other.
It is found from Tables 1 and S1 that the FV is enhanced but the FFV is decreased at higher temperature. This result can be attributed to the increase in the intermolecular distances by raising the temperature in water leading to the greater FV. On the other hand, as the total cell volume is also very much enhanced at higher temperature, the FFV is diminished and does not indicate an identical trend with the FV. For instance, the FV and FFV of the PEG-PLA-GEM are 3511.85 Å3 and 38.03% at 298.15 K but they are 4716.30 Å3 and 29.51 at 310.15 K, respectively. Similar results are obtained for all other cells. Moreover, addition of FA into the systems results in increasing both of the FV and FFV values. This is similar to the result observed at 298.15 K.

Radius of gyration (Rg)
The Rg evaluates the dimension of a polymeric chain that is calculated using the following formula where mi and ri stand for the mass and position vectors of atom i related to the molecular mass center [63]. The Rg2 designates the mean square of the distance between the polymeric chains and the mass center; thus a bigger Rg reveals a greater distance between the polymeric chain and its center of mass [59].

                                              
The gyration radii data measured at 310.15 and 298.15 K for all of the PEG-PLA-GEM and PEG-PLA-GEM-FA polymeric systems are demonstrated in Tables 1 and S1, respectively. Also, the Rg diagrams of the systems at 310.15 and 298.15 K are displayed in Figs. S5 and 6, respectively. It is obvious that, at the two temperatures, adding nanotubes into both kinds of the PEG-PLA-GEM and PEG-PLA-GEM-FA cells provides smaller Rg values. For instance, the Rg values measured at 298.15 K for the PEG-PLA-GEM cells containing CNT, CN, CP, SiN and SiP fillers are 3.13, 3.08, 2.84, 2.95, 2.89 and 3.01 Å, respectively. Furthermore, the Rg values in both PEG-PLA-GEM and PEG-PLA-GEM-FA systems filled by different nanotubes change in the order of CN<SiN<CP<SiPConsequently, it is established that the polymeric chains in the PEG-PLA-GEM-FA-CN system are in the closest vicinity due to the greatest intermolecular interactions with themselves and other species which have led to their utmost contacts. This result is also in agreement with the FV and FFV data so that the greatest values were measured for the PEG-PLA-GEM-FA-CN among all of the nanocomposite systems. Accordingly, it may be expected that the delivery of the GEM drug can be the most efficient using this system which has the highest FFV and thus provides the most accelerated drug diffusion. A comparison of these findings with the experimental data indicates that the Rg values are acceptable. For instance, the Mw dependent Rg values for different PLA chains were obtained about 1‒5 Å [64] and 5.0‒6.2 [65]. As well, the solution concentration dependent Rg values for the PEG chains were measured near 19.5‒21.0 Å [66].

The X-ray diffraction (XRD) patterns
The XRD intensity, I(Q), of a material is calculated using the following equation [39,40] where θ shows the scattering angle, λ is the X-ray wavelength and Q stands for the magnitude of the scattering angle, Q = 4πsinθ/λ . The indices j and k vary for all atoms existing in a structure/material. The fi/fj and rjk exhibit shape factor and mutual distance values, respectively.

The XRD patterns obtained at 310.15 and 298.15 K for all of the PEG-PLA-GEM and PEG-PLA-GEM-FA systems in the 2θ range of 0-120° are depicted in Figs. S6 and 7, respectively. The XRD pattern indicates the crystallinity of a material by comparing the intensities of peaks so that sharp peaks are observed for a crystalline compound whereas an amorphous structure illustrates broad signal(s). As the XRD patterns of all PEG-PLA-GEM and PEG-PLA-GEM-FA cells at both temperatures reveal one intense semi-broad peak at around 20° having some sharp shoulders plus two weak peaks at about 35 and 60°, they have semi-crystalline structures. These patterns approve the existence of PEG, PLA and nanotubes in all of the systems. Literature review reveals that the PEO powder displays a crystalline XRD pattern containing two very sharp peaks near 20 and 24° plus four weak peaks at about 26, 27, 35 and 40° [67]. The PLA powder demonstrates a broad peak with the highest intensity near 17° [17]. The CNT exhibits its characteristic XRD peaks at 28 and 45° [68,69].
The inter-chain distances calculated using the Bragg equation (λ=2dsinθ) at 310.15 and 298.15 K are given in Tables 1 and S1, respectively. It is observed that adding nanotubes into the cells increases the diffraction angle of the maximum peak in the XRD patterns; thus, the d-spacing values are decreased for the nanocomposite systems. For instance, the inter-chain distances measured at 298.15 K for the pure PEG-PLA-GEM cell and those filled by CNT, CN, SiN and SiP are CS-GP-CP are 0.48, 0.47, 0.43, 0.45 and 0.46 Å, respectively. Similar results are achieved for the PEG-PLA-GEM-FA systems indicating the d-spacing changes as CN<SiN<CP<SiP

Radial distribution function (RDF)
In order to get insight about the interaction distances within the PEG-PLA-GEM and PEG-PLA-GEM-FA polymeric nanocomposites, the RDF value which is also called pair correlation function, gA–B(r) is calculated. The RDF estimates the probability of distributing B atoms nearby A atoms. The gA–B(r) is obtained using the following equation, where nB shows the number of B atoms placed around A atoms in a spherical shell with Δr thickness, NB stands for the whole number of B atoms added to the amorphous cell and V exhibits the cell volume.

The RDF values reflect that weak, moderate or strong interactions are happened between the interacting species. Figs. S7-S17 present the RDF plots and Tables S3 and S4 afford the RDF data for all of the PEG-PLA-GEM and PEG-PLA-GEM-FA systems containing various nanotubes at 298.15 K. Figs. 8 and 9 illustrate the RDF plots and Tables 3 and 4 give the RDF data for the inter-atomic interactions occurred in the PEG-PLA-GEM and PEG-PLA-GEM-CNT systems in water at 310.15 K. Additionally, the RDF plots for other PEG-PLA-GEM and PEG-PLA-GEM-FA systems containing various nanotubes are illustrated in Figs. S18-S27. It should be noted that the H(NH2-GEM)-O(PEG) defines the distance between hydrogen atom of NH2 group on the GEM drug and the oxygen atom on the PEG chain. This definition is utilized for all other RDFs.
It is recognized that in the PEG-PLA-GEM system, the drug-PEG interaction is stronger than that of the drug-PLA because the H(NH2-GEM)-O(PEG) is smaller than the H(NH2-GEM)-O(PLA). This result is also observed in all of the PEG-PLA-GEM and PEG-PLA-GEM-FA cells containing nanotubes. Moreover, the interactions of GEM, PLA and PEG with the nanotubes vary in the order GEM-nanotube<PLA-nanotubePLA-nanotube>PEG-nanotube which is completely opposite to the result observed at 298.15 K. For example, the F(GEM)-C(CNT), O(PLA)-C(CNT) and O(PEG)-C(CNT) distances for the PEG-PLA-GEM-CNT cell at 310.15 K in water are 2.25, 2.73 and 2.99 Å, respectively. Moreover, some of the RDFs are very much decreased at higher temperature (310.15 K) and this may be related to the presence of H2O molecules in the systems interacting with the functional groups of all species.
Comparing the RDF values for all of the cells reveals that smaller values are achieved for the systems containing FA molecules. This confirms that addition of the FA results in stronger H-bond and electrostatic interactions among the functional groups of the PEG, PLA, GEM, FA and nanotubes. For example, in the PEG-PLA-GEM-CN and PEG-PLA-GEM-CN-FA systems at 298.15 K, the O(PEG)-N(CN), F(GEM)-N(CN) and O(PLA)-N(CN) distances are 2.91, 3.78, 3.45 and 2.87, 3.60, 3.70 Å, respectively, but they are 2.89, 2.15, 2.63 and 2.81, 2.17, 2.65 Å at 310.15 K, respectively. Moreover, among all PEG-PLA-GEM and PEG-PLA-GEM-FA systems composed of different nanotubes, the smallest RDFs are achieved when CN is added into the cells. Consequently, it may be anticipated that the PEG-PLA-GEM-CN-FA can the most effectively deliver the GEM molecules and it may be known as the most suitable drug carrier.

Solubility parameter
The solubility parameter (δ) numerically estimates the degree of interaction between materials, commonly for polymers, to indicate their solubility [70]. According to the Hildebrand-Scatchard solution theory, it can evaluate the permeation, swelling, compatibility and bulk/solution properties of polymers [66]. It is the square root of the cohesive energy density (ECED), see Eq. 1.

                                                                                     (1)

The ECED is the energy required to completely remove unit volume of molecules from their neighbors to infinite separation (an ideal gas); this equals the (ΔHvap–RT)/V (Eq. 2), where ΔHvap indicates the molar vaporization energy, V reveals the molar volume, R stands for the gas constant and T is the temperature.

                                                                                     (2)

The Hildebrand parameter is comprised of three Hansen parameters, i.e. dispersion, polar and hydrogen bond forces called δd, δp and δh solubility parameters, respectively (Eq. 3) [71].

                                                                                     (3)

The solubility parameters for the PEG-PLA-GEM-CN and PEG-PLA-GEM-CN-FA polymeric nanocomposites and their constituents at 310.15 and 298.15 K are afforded in Tables S5 and 5, respectively. It is distinguished that the solubility parameter of GEM drug is closer to those of the PEG, PLA and FA compared to the values measured for the nanotubes. This reveals that the GEM molecules are appropriately compatible and miscible with the PEG-PLA and PEG-PLA-FA blends. Also, among all nanotubes, the solubility parameter is the greatest for the CN. The δ values for the CNT, CN, CP, SiN and SiP at 298.15 K are 7.89, 10.70, 8.78, 9.26 and 8.03 (J.cm-3)0.5, respectively. The highest δ value for the CN can be related to the greatest electronegativity of N among C, N and P elements which can lead to the highest hydrogen bonding and polar (dipole) interactions. Similarly, the second greatest solubility parameter is obtained for the SiN. The smallest value for the SiP is due to its greatest volume which gives the lowest δ (see Eq. 2). Comparable results are achieved at 310.15 K for these species.
Comparing the solubility parameters measured at the two temperatures demonstrates that the solubility parameter of a single molecule is decreased or increased by rising the temperature from 298.15 to 310.15 K and this can be related to the increase/decrease of dispersion, polar or hydrogen bond forces, i.e. δd, δp and δh solubility parameters, see Eq. 3. Nevertheless, the solubility parameters of all PEG-PLA-GEM and PEG-PLA-GEM-FA systems are enhanced at 310.15 K in water that is mostly attributed to the presence of H2O molecules in the systems with very strong intermolecular polar and hydrogen bond interactions.
The solubility parameters for the PEG-PLA-GEM nanocomposites containing different nanotubes at both temperatures change in the order of CNT<SiP<CP<SiN

Mean square displacement (MSD) and diffusivity
The MSD plots for the PEG-PLA-GEM and PEG-PLA-GEM-FA cells at 310.15 and 298.15 K are drawn in Figs. S28 and 10, respectively, which estimate the diffusion of the GEM molecules inside the polymeric nanocomposites. The nearly linear lines in the MSD graphs illustrate constant drug diffusion in these systems during the MD simulations. As well, the drug diffusion coefficients have been achieved using the slopes of the MSD diagrams. Comparing the two kinds of the PEG-PLA-GEM and PEG-PLA-GEM-FA systems confirms that the MSD values are smaller for the latter which is more desirable in order to effectively deliver the drug in a sustained manner. Additionally, it is seen that the PEG-PLA-GEM-FA has the highest MSD for the GEM diffusion but among the PEG-PLA-GEM-FA nanocomposites containing different nanotubes, the MSD changes in the order of CN<SiN<CP<SiP<CNT; consequently, the PEG-PLA-GEM-FA-CN system can be employed to achieve the most controlled/slowest transport/diffusion of the GEM molecules. Comparing the MSD curves at the two temperatures of 298.15 and 310.15 K reveals that, as expected, the GEM diffusion has been enhanced at higher temperature. This can be related to the higher kinetic energies of the GEM molecules at higher temperature leading to greater diffusion coefficients.
The diffusion coefficients for the GEM molecules in both types of the PEG-PLA-GEM and PEG-PLA-GEM-FA polymeric systems at 310.15 and 298.15 K are gathered in Tables S6 and 6, respectively. These data are consistent with the MSD curves as lower diffusion coefficients are measured for the PEG-PLA-GEM-FA systems among those of the PEG-PLA-GEM cells. Furthermore, in both systems, the diffusion coefficients vary with the nanotube type in the order of CN<SiN<CP<SiPAs there are not any experimental data on the systems examined in this work, the drug diffusion coefficients experimentally measured for other similar drug delivery systems [72-75] are presented to get insight about the drug diffusion coefficient values in real systems, see Tables S7-S10. A comparison of the diffusion coefficients calculated for the GEM molecules inside the PEG-PLA-GEM and PEG-PLA-GEM-FA systems with the experimental values reported for other drugs in various systems in Tables S7-S10 indicates that the diffusion coefficients for the GEM drug change between 2.68×10-7 and 7.75×10-7 that are almost comparable to most of the data in Tables S7-S9 but larger than those of Table S10. This result confirms that the diffusion coefficients computed for the GEM drug inside the PEG-PLA-GEM and PEG-PLA-GEM-FA systems are acceptable and comparable to the experimental values measured in real systems. It is notable that the very small diffusion coefficients in Table S10 can be related to the strong interactions of progesterone as the diffusing agent with the PLGA particles.

CONCLUSIONS
The MD simulations were performed on the PEG-PLA-GEM and PEG-PLA-GEM-FA polymeric systems which were filled by CNT, CN, CP, SiN or SiP nanotubes in order to achieve the most suitable system for the delivery of anticancer drug GEM. Also, the effects of temperature and solvent were examined by running the simulations at two temperatures (298.15 and 310.15 K) in both gas phase and water solvent. The Rg values in both the PEG-PLA-GEM and the PEG-PLA-GEM-FA systems contained diverse nanotubes were changed in the order of CN<SiN<CP<SiPACKNOWLEDGEMENTS
The financial support of this work by the Research Office of Amirkabir University of Technology (Tehran Polytechnic), Tehran, Iran Author is gratefully appreciated. Also, author is grateful to High Performance Computing Cluster of Amirkabir University of Technology (Tehran Polytechnic), Tehran, Iran for affording the computational facilities (software and hardware) to accomplish the MD simulations.

CONFLICTS OF INTEREST
The authors declare that there are no conflicts of interest regarding the publication of this manuscript.