Molecular Dynamics Simulations on Polymeric Nanocomposite Membranes Designed to Deliver Pipobromane Anticancer Drug

INTRODUCTION
Currently, chemotherapy has become a common cancer treatment by means of cytotoxic drugs [1–4]. However, traditional anticancer molecules may lead to severe systemic toxicity by distributing in the whole body through the blood circulation system [5–7]. This causes several adverse side effects including despairing the immune system, muscle pain and organ damage. Thus, chemotherapy treatment can be painful which may result in its non-popularity [8]. Furthermore, most of anticancer drugs are hydrophobic with poor water solubility and this can hinder their widespread usage [9]. Consequently, it is crucial to develop targeted and controlled drug delivery systems (DDSs) for various drugs [10-12] and especially for the hydrophobic therapeutic anticancer drugs [13,14]. Nanocomposite DDSs materials composed of polymeric chains, inorganic nanoparticles and drug biomolecules have well been established as drug delivery vehicles for the anticancer drugs [15–17]. They have exhibited promising physicochemical properties which could be served as favorable candidates for novel DDSs.
Natural cationic polymers such as chitosan are interesting compounds for therapeutic applications since they are commonly non-toxic, biocompatible, derived from renewable resources and biodegradable with low immunogenicity [18]. The polymeric nanocarriers should be soluble at physiological conditions in order to be successfully employed in drug delivery purposes. Chitosan (CS) is of a great attention in producing novel DDSs [19] due to it can form complexes with anticancer drugs antibiotics and therapeutic proteins which is a desired property for its benign usage as a pharmaceutical drug delivery tool. Thus, chitosan micro and nanoparticles are worldwide investigated as efficient drug delivery vehicles [20-22]. It is noteworthy that these systems are applied in controlled and targeted release of practically all families of bioactive molecules. Polyethylene oxide (PEO) is a biocompatible polymer, which is a non-ionic, water soluble polymer, with good binding and film forming properties. It retards the release rate of drugs and hence it is commonly used in pharmaceutical formulations such as controlled release dosage forms, hot-melt technology and mucoadhesive dosage forms. Thus, it plays a central role in designing novel drug delivery systems for both highly and poorly soluble drugs [23]. Polylactic acid or polylactide (PLA) is a biodegradable, biocompatible and thermoplastic aliphatic polyester derived from renewable resources, such as corn starch, tapioca roots, chips, starch, or sugarcane [24-26]. Being able to degrade into innocuous lactic acid, PLA is used as medical implants in the form of anchors, screws, plates, pins, rods, and as a mesh [27]. Depending on the exact type used, it breaks down inside the body within six months to 2 years. This gradual degradation is desirable for a support structure, since it slowly transfers the load to the body (e.g. the bone) as that area heals. The strength characteristics of PLA implants have been well examined [28].
Pipobromane (PIP) is a cytoreductive agent predominantly used for the treatment of polycythemia vera (PV) and essential thrombocythemia (ET) [29]. Its ability to inhibit myeloproliferation was documented in the 1960s [30] which made this drug of specific attention in treatment of disorders such as PV and ET. Indeed, ET is a myeloproliferative disorder which is characterized by a high and persistent thrombocytosis that normally exceeds 600×109/L. It is recommended to use platelet-lowering agents to patients having a risk factor for thrombosis or haemorrhage [31]. When these patients remain untreated, they may be encountered with about 45% thromboembolic complication [32]. Therapeutic cytoreductive agents utilized in ET are pipobroman [33], hydroxyurea [34], busulphan [35] and interferon alpha [36]. Pipobroman produced by Abbott laboratories and sold in Europe with the commercial name Vercyte (25 mg tablets), is a bromide derivative of piperazine, 1,4-bis(3-bromopropionyl) piperazine. Its chemical structure is similar to those of the alkylating agents however it seems also to act as a metabolic rival of pyrimidine bases. PIP is the most inexpensive platelet lowering compound accessible for the ET treatment. PIP was attracted by its effectiveness in dropping the platelet values and in preserving thrombocytosis in rather low doses. It was validated that PIP was effective and well tolerated with a little risk of thrombosis, myelofibrosis and leukaemia in a great number of patients with polycythaemia vera [37]. PIP showed high efficacy in controlling the thrombocytosis. This result was similar to the information achieved in PV with hydroxyurea [38].
Nowadays, molecular dynamics simulations (MD) have found enormous attraction in pharmaceutical investigations due to the MD studies offer insights into the intermolecular interactions occurring in drug-loaded polymeric matrixes [39]. Even though the polymeric nanocomposite DDSs have experimentally been examined, limited studies have been accomplished on drug encapsulation at the molecular or mesoscopic levels and little information has been acquired in this area. Computer simulations afford an alternative approach to model the diffusion behaviors of drug molecules within the polymeric DDSs and to gain insights into the drug delivery in an aqueous media [40]. Indeed, MD simulations are widely accomplished to understand physical phenomena occurring in complex membranes at a molecular level [41-43]. The MD simulations can investigate the intermolecular interactions of polymeric carriers and drug molecules in polymeric nanocomposite drug delivery vehicles [44]. It is noteworthy that few theoretical simulations have so far been accomplished on the delivery of anticancer drugs by means of polymeric nanocomposites as DDSs [45,46].
Concerning the above-mentioned topics, it was found intriguing for us to fulfill the MD simulations to deliver the anticancer drug PIP molecules by means of three silica filled polymeric composites based on CS, PLA and PEG polymeric matrixes. The novelty of this work is using MD simulations to envisage the drug delivery efficacy of polymeric composites as most of the researches accomplished till now are experimental ones and there are a few theoretical studies on the polymeric composites as DDSs. Herein, the diffusion performance of PIP drug to the CS, PLA and PEG matrixes including SiO2 nanoparticles (NPs) was investigated using the MD simulations. Notably, the silica NPs was loaded into these systems due to they not only could enhance the drug delivery effectiveness of the polymer composite systems but also they have esteemed characteristics including non-toxicity and biocompatibility. Furthermore, the interaction energies between diffusing drug molecules and polymeric chains, free volume (FV), FFV and the diffusion coefficients were assessed and compared with each other with the purpose of determining which nanocomposite DDS would be the most apposite one to deliver PIP molecules.

SIMULATION METHODOLOGY
The COMPASS force field (FF) which is condensed-phase optimized molecular potentials based on the Polymer Consistent Force Field (PCFF) [47] which is the first FF validated and parameterized with condensed-phase criteria, ab initio and empirical data presented for various molecules. The COMPASS takes into account cross coupling contacts in addition to the bonded and non-bonded interactions applied in other FFs [48]. Charge equilibration approach was selected to calculate partial atomic charges existing on molecules [49]. Also, summation of valence bond, non-bond and cross-term interactions provided the total energy for a system. The charge-equilibration method is used to achieve partial charges existing on atoms. The Columbic long-range contacts are acquired by the Ewald method [50].
In a system, summation of valence (bond), non-bond and cross-term interactions gives total energy for a system which is written as:

Etotal = Evalence + Enon-bond + Ecross-term

The Evalence represents the potential energy of the system that includes bond stretching (Ebond), bending of valence angle (Eangle), torsion of dihedral angle (Etorsion) plus inversion (or out of plane interaction) represents as Einversion or Eoop. Modern FFs such as COMPASS comprise Urey-Bradley term (EUB) which denotes interactions between atomic pairs within 1–3 configurations (for atoms binding to an atom):

Evalence = Etorsion + Eangle + Ebond + EUB + Eoop
Cross terms occurring in Ecrossterm enhance the FF accuracy employing the correction factors onto valence energy in order to show the interdependence between diverse valence terms. As an example, the term Ebond-bond considers the stretching-stretching interactions between two adjacent bonds. Likewise, the COMPASS FF is comprised of bending-bending, stretching-bending, bending-torsion stretching-torsion and bending-bending-torsion terms. The term Enon-bond specifies interactions between non-bonded nuclei (or secondary interactions) and it is composed of three Coulomb electrostatic (ECoulomb), van der Waals (EvdW) as well as hydrogen bond (EH-bond) energies:

Enon-bond = ECoulomb + EvdW + EH-bond

All of species were constructed by means of the molecule modelling tool. All of amorphous cells were incorporated with 500 H2O molecules, 10 pipobromane drug molecules, 4 silica NPs (cristobalite-high) and 20 polymeric (chitosan, polyethylene glycol or polylactic acid) chains having twenty repeating units in order to assess the influence of on these polymers on the drug delivery efficacy. The three-dimensional cells were optimized by running the MD simulations at standard temperature (298.15 K) and pressure (1 atm). Also, the original amorphous cells were produced at a very low density of 0.25 g cm–3 to create the cells having periodic boundary conditions using the amorphous cell module which exists in Materials Studio software (version 4.3) [51]. This initial density was selected to highly equilibrate the prepared cells. The chemical structures of CS, PLA, PEG and PIP compounds are given in Fig. 1. Furthermore, complete compositions of the cells along with the simulation results are provided in Table 1.
The energy minimizations of all configurations in amorphous cells were primarily accomplished by the smart minimizer algorithm which is a combination of conjugate gradient, steepest descent and Newton minimization algorithms. The minimizations were done in a cascading manner for 2,000,000 iterations to reach relaxation and the convergence criteria were chosen as 1000 and 0.1 kcal/(mol.Å), respectively. Subsequently, dynamics simulations were carried out to entirely relax the simulation systems. Therefore, the NVT simulation was initially performed at 298.15 K for 2 ns and after that, the NPT simulation was done at 298.15 K and 1 atm for 2 ns. This methodology intended to completely relax the cells. Throughout the relaxation steps, the non-bond interactions (i.e., van der Waals and electrostatic potentials) were calculated by Ewald summation method. After that, MD simulation was accomplished by the NVT ensemble at 298.15 K at the relaxed density for 3 ns by applying cell multipole method (CMM). At last, the final 1 ns trajectories were utilized to analyze the structural, dynamical and energetic characteristics [52]. The velocity Verlet algorithm was chosen to unravel the classical motion equation with a 1 fs time step [53]. The Berendsen barostat and thermostat procedures with a 0.1 ps decay constant were applied for all of the NPT and NVT simulations to maintain the selected pressure and temperature during the MD runs [54].

RESULTS AND DISCUSSION
Equilibrium and relaxation of cells
In order to study and compare various features of CS, PLA and PEG polymeric systems composed of silica NPs, water and PIP drug molecules, their relaxations and equilibrations were performed. Hence, the relaxation properties such as temperature along with potential and non-bond energies were projected during the MD simulations. Fig. 2 exhibits the relaxed configurations of the CS, PLA and PEG polymeric composites achieved after 1 ns production run. The presence of the SiO2 NPs within the polymer composites loaded by the PIP drug as well as H2O molecules is obvious in Fig. 2. The water molecules have been added to the polymeric matrixes to simulate the drug delivery in the aqueous medium which is similar to the human body. The computed density and dimension (cell size) values for the amorphous simulation cells including polymer composites are provided in Table 1. Apparently, the greatest cell size and consequently the greatest cell volume is obtained for the CS, then PLA and the smallest value is acquired for the PEG system. The changes in cell dimension could be accredited to the greatest volume of CS polymeric chains relative to those of the PLA and PEG chains which may enhance the inter-chain spaces. The volumes measured for the CS, PLA and PEG chains are equal to 2905.90, 1329.80 and 905.18 Ǻ3, respectively, which approves the CS polymer occupies the greatest volume among the three nanocomposite systems. The density of the CS matrix is 1.25 g/cm3 whereas it is dropped to 1.16 g/cm3 in PEG and 1.03 g/cm3 in PLA. This is in agreement with the theoretically and experimentally values stated for the CS, PLA and PEG composites. For instance, the densities for the CS [55], PLA [55] and PEG [56,57] were measured near 1.42, 1.27 and 1.10 g/cm3, respectively. Therefore, the simulated densities for these polymeric nanocomposites appear to be reasonable. This can also be accounted for another sign of attaining cell equilibrium.
Fig. 3 presents deviations in the non-bond and potential energies as well as temperature versus simulation time for all of the completely optimized CS, PLA and PEG polymeric systems. It is observed that the changes are very much little and the systems have reached the relaxation states. Additionally, the non-bond energies of CS, PLA and PEG DDSs are smaller than their related potential energies by ~5500, 3850 and 1550 kcal/mol, respectively. Comparing the potential and non-bond energies of CS, PLA and PEG composites discloses that PEG system provides the most negative energies, then PLA and at last the most positive values are measured for the CS system. The potential energies of the three CS, PLA and PEG nanocomposites are about 19500, 3700 and –4600 kcal/mol while their corresponding non-bond energies are approximately 14000, –150 and –6150 kcal/mol, respectively. Thus, it may be established that the PEG DDS with the most negative energies is the most apposite system. Nonetheless, in order to choose the most appropriate DDS with the most effectual drug delivery, other characteristics of the CS, PLA and PEG polymeric composites must be estimated and compared with each other. As a result, the relationships between polymeric nanocomposite type and PIP delivery efficacy will be analyzed in following sections to be able to state which system will be the most suitable candidate to be utilized for the delivery of PIP anticancer drug.

Surface area and free volume
Two different phases exist in a polymeric matrix that are solid phase polymeric chains plus other species as well as an empty FV space known [58]. The FV is obtained by the summation of the static holes produced by chains packing or transient gaps made by thermally rearrangement of the chains; thus a low-resistance transferring route is formed for the diffusing molecules [59]. It is recognized that diffusion inside a polymeric matrix depends on its FV and morphology. The FV has a direct effect on polymer diffusivity and enhancement in the FV can expedite the molecular diffusion [60]. In order to obtain the FV values for the CS, PLA and PEG polymeric nanocomposites, the Connolly surface method was applied. In the Connolly surface method, the probe fragment is modeled by a hard sphere having a radius equal to 0.1 nm. The FV is acquired by means of the hard sphere probe which is noticeably depends on the probe size.
The surface areas of the CS, PLA and PEG polymeric systems are given in Table 1. It is observed that the highest surface area belongs to the CS system (54524.24 Å2), then PLA (35041.47 Å2) and the least surface area is measured for the PEG matrix (26816.46 Å2). This may be correlated to the maximum volume as well as surface area of the CS chains among those of the PLA and PEG polymers. Fig. 4 exhibits the chemical structures and surface areas (indicated by blue color) for the CS, PLA and PEG chains. The surface areas measured for the CS, PLA and PEG chains are equal to 2266.26, 1031.68 and 727.62 Ǻ2, respectively, which supports the CS polymer affords the utmost surface area among the three nanocomposite systems.
Table 1 presents the FV values for the simulation cells composed of CS, PEG and PLA DDSs. Also, Fig. 5 demonstrates the FVs by blue color for these polymeric composites. It is found that the FV is the greatest for the CS nanocomposite (260786.68 Å3) and it is decreased to 92242.32 Å3 in PLA system and the lowest value is achieved for the PEG (91440.05 Å3). The fractional free volume (FFV) is equal to the ratio of FV to the total cell volume. It is found that CS indicates the maximum FFV (77.232%) whereas PLA nanocomposite has the smallest FFV (65.804%). High FFV values of all DDSs reveal that the diffusion and transport of the PIP drug in these cells would occur very simply and it is estimated that high diffusion coefficients will be achieved for these systems. The reason for the decreased FV and FFV values can be associated to the restrained segmental movement of polymeric chains near the silica–polymer interface that makes a number of nanosized holes as suitable spaces to diffuse/transport the PIP drug molecules. As the maximum FV and the FFV values belong to the CS among all of the composites and the smallest values are obtained for the PEG system, it may be expected that the most effective PIP delivery will be done using the CS system. Nonetheless, in order to determine the most appropriate DDS to be employed for the PIP drug delivery, additional parameters should be assessed and compared with each other.
An upsurge in the FV and FFV values typically give rise to more efficient penetrations inside the polymeric matrixes. As an example, the FV and alcohol transport features in the PDMS membranes were evaluated and it was shown that raising the temperature enhanced polymeric chain motion and an enlarged the FV which helped the feed passage in the membrane [59]. Another group studied three polysulfone membranes with various FVs using different solvents and the FV influence on the CO2 sorption/permeation was assessed and it was found that permeability and solubility coefficient were both increased by the FV enhancement [61]. The sorption, swelling and FV of polybenzimidazole (PBI) and PBI/zeolitic imidazolate framework (ZIF-8) nanocomposite were tested for pervaporation application and the positron annihilation lifetime spectroscopy established that great pervaporation permeability of PBI/ZIF-8 nanocomposite was related to the high FFV produced by big cavities of the ZIF-8 species [62].

Radius of gyration
The radius of gyration (Rg) defines the dimension of a polymeric chain [63]. The Rg2 is the mean square of the distance between the beads and the mass center. This means that a larger Rg indicates a greater distance between the polymeric chain and its mass center [63]. The Rg values and graphs for all of the CS, PLA and PEG polymeric systems are in Table 1 and Fig. 6a, respectively. It is realized that the CS provides the smallest Rg (1.85 Å) while the PEG has the greatest Rg value (2.96 Å). The Rg for the PLA is equal to 2.28 Å. These data specify that the CS chains have the closest contacts and intermolecular interactions with each other that cause their utmost proximity. This can bring about the highest FV and FFV values for the CS nanocomposite DDS among the three systems. As a result, it may be decided that the CS matrix with the highest FV will the most efficiently deliver the PIP drug and leads to the most facilitated carriage of the PIP drug molecules inside the CS drug delivery system. As a comparison with the experimental data, it is found that the Rg of CS chains in the unperturbed state was equal to 5.25 Å [64]. The Rg for the CS, in another work, was reported equal to 2.7 Å [65]. It is notable that the Rg depends on the degree of polymerization (DP). It was indicated that the Rg was changed in the range of 3.44 to 8.25 Å for the CS polymer with DP changing from 646 to 2100 [66]. Similarly, the Mw dependent Rg values for diverse PLA chains were measured around 5.0‒6.2 [67] and 1‒5 Å [68]. The solution concentration dependent Rg values for the PEG were attained close to 19.5‒21.0 Å [69].
The X-ray diffraction (XRD) pattern
The scattering analysis from the Forcite module of the simulation software was employed to produce the XRD patterns of all CS, PEG and PLA systems. The diffraction angle was altered from 0° to 120° with a 0.05° step size. The peak with maximum intensity in the XRD pattern was used to calculate the d-spacing values that are illustrative of inter-segmental distances between polymeric chains.
The XRD patterns of the CS, PEG and PLA polymeric nanocomposite systems are displayed in Fig. 6b. The crystallinity of a compound can be estimated using these diagrams by comparing the background intensity to the intensities of sharp signals. Commonly, a crystalline material reveals some sharp peaks while an amorphous compound (like glasses and liquids) illustrates a broad signal. It is obvious from the expanded XRD patterns in Fig. 6c that there is one peak near 20° which is indicative of the CS, PLA, PEG, silica NPs and H2O molecules presence in all of the nanocomposite systems. Fig. 6d exhibits the XRD patterns obtained for the isolated CS, PLA, PEG, silica NPs and H2O molecules which reveal a broad peak for each species at about 20, 17, 19, 19.5 and 23°, respectively. The XRD pattern of the PIP drug in Fig. 6d presents some broad peaks located at around 18, 21, 24 and 26°. As all of the cell components for the CS, PEG and PLA systems have broad peaks at about 20°, it can be decided that the appearance of a broad peak in the XRD pattern of each CS, PEG and PLA drug delivery system is acceptable.
Similar results were described in literature in which the amorphous SiO2 showed a characteristic XRD peak at 2θ≈23° [70]. The CS [71] and PEG [72,73] exhibited their corresponding XRD peaks at 2θ≈20°. The PLA indicated a broad peak at about 17° [74]. The reason for the appearance of the XRD peak at around 20° in our PLA system can be attributed to the occurrence of SiO2 NPs, PIP drug as well as water molecules within the polymeric nanocomposite matrix.
It is observed that the intensity of the peak at 2θ≈20° is the highest for the CS but it decreases in case of PEG and the lowest intensity is measured for the PLA system. This may be related to the highest crystallinity of CS while decreased crystallinity of the PEG and PLA systems. Further, the position of the XRD peak at 2θ≈20° moves to smaller 2θ value from CS nanocomposite to the PEG and PLA which can be described by the smallest inter-chain distance in CS but the greatest one in PLA.
The Bragg equation, (λ=2dsinθ), was used to calculate the inter-chain distances. Table 1 affords the 2θ values of the maximum XRD peaks and their corresponding d-spacing values for the nanocomposites. It is found that the inter-chain distances for the CS, PEG and PLA systems are obtained equal to 1.025, 1.606 and 2.605 Å, respectively. The decreased d-spacing may be accounted for greater intermolecular interactions which are due to the formation of the strong as well as weak hydrogen bonds among the polymeric chains, SiO2 NPs, PIP drug as well as H2O molecules. These data are in well agreement with the FV and FFV values so that the CS illuminates the largest FV and FFV among all of the DDSs, which is related to the greatest inter-particle interactions.

Radial distribution function (RDF)
The interaction distances of the CS, PEG and PLA polymeric nanocomposites were measured by the RDF values. RDF is also famous as pair correlation function, gA–B(r), indicating the probability of distributing B atoms near A reference atoms. Different RDF diagrams demonstrating the inter-atomic interactions happened in the CS, PLA and PEG drug delivery systems are provided in Figs. 7-9. Various RDF values reflect weak, moderate and strong interactions occurred between the species. It is noteworthy that the H(CS)–C(C=O,PIP) exhibits the distance between H atom of CS and C atom of the carbonyl group of the PIP drug molecule. Similarly, this type of definition is applied for other kinds of RDF distances offered in Figs. 7-9. Additionally, the RDF values for the intermolecular interactions took place between polymeric chains, silica NPs, PIP drug and H2O molecules in the polymeric nanocomposite systems are given in Table 2.
In the CS DDS, the silica-polymer interactions are very stronger than those of the silica-drug interactions due to the silica-polymer intermolecular distances are smaller. For example, the H(CS)-O(SiO2) and N(CS)-Si(SiO2) distances are 2.53 and 3.21 Å that can be considered as the silica-polymer interactions but the Si(SiO2)-O(C=O,PIP), Si(SiO2)-N(PIP), Si(SiO2)-Br(PIP) and O(SiO2)-C(C=O,PIP) distances are measured equal to 5.37, 3.01, 5.37 and 2.67 Å which are known as the silica-drug interactions. Comparable results are attained in case of the PEG system while a reverse trend is achieved for the PLA in which the silica-polymer interactions are weaker than their corresponding silica-drug interactions (see Table 2). Besides, in CS and PEG systems, the drug-polymer interactions are stronger (smaller RDFs) than those of their related drug-silica interactions but an opposite observation is achieved for the PLA for which the drug-polymer interactions are weaker than the drug-silica interactions. The H2O molecules can create hydrogen bonds with all of the species existing inside the polymeric DDSs. In the CS system, the water molecules have a stronger interaction with the CS chains compared with silica and PIP particles. This is understood from smaller O(W)-H(CS) RDF (2.07 Å) compared to those of the O(W)-C(PIP) and O(W)-Si(SiO2) RDFs (2.69 and 2.71 Å). The water-silica, water-polymer and water-drug RDFs are equal in both of the PEG and PLA systems that are 2.65 and 2.73 Å, respectively. These RDFs confirm that the H2O molecules form strong hydrogen bonding interactions with the CS polymeric chains, weaker electrostatic interactions with the PEG and the weakest interactions with the PLA.
The drug-polymer interactions in the CS DDS display that the oxygen atoms of carbonyl groups on the PIP molecules form very strong C=O…H-N hydrogen bonding interactions with the amino hydrogen atoms of CS chains, H(CS)-O(C=O,PIP)=2.09 Å; this could propose that the PIP molecules would diffuse the most slowly within the CS system. This finding may validate that the PIP drug molecules have the strongest interactions with the CS chains among other polymers. As a result, it can be concluded from the RDF data that the PIP drug molecules will show the most controlled diffusion/transport in the CS nanocomposite DDS.

Mean square displacement (MSD) and diffusivity
The PIP molecular diffusion into the CS, PLA and PEG polymeric nanocomposite systems was estimated from their MSD diagrams that are presented in Fig. 10. Also, the diffusion coefficient (D) for the PIP molecules was obtained using the slopes of these curves. The nearly linear lines in the MSD plots deduce constant diffusion of PIP drug molecules in the CS, PLA and PEG DDSs during the MD simulations. As well, it is noticeable that the CS provides the smallest MSD for the PIP diffusion while PLA shows the highest MSD. Therefore, in order to have a proficient PIP delivery, a system with the smallest MSD must be nominated. This leads in the most controlled and slowest passage/diffusion of PIP molecules within the CS DDS but the fastest transfer in the PLA nanocomposite. It is found from Fig. 10 that the CS is the most effective and suitable nanocomposite system for the PIP drug delivery. Table 3 affords the diffusion coefficients for the PIP molecules in the three CS, PLA and PEG polymeric nanocomposite systems. These data are well consistent with the MSD curves so that the highest diffusion coefficient (0.0183×10–4 cm2/s) is obtained for the PLA whereas the lowest diffusion coefficient is attained for the CS (0.0154×10–4 cm2/s). The PEG system illustrates a diffusion coefficient equal to 0.0163×10–4 cm2/s which is between the values measured for the CS and PLA systems. Thus, the diffusion of PIP in the CS matrix is the slowest which permits the most efficient drug delivery. This result supports that the delivery rate of the PIP drug using the three examined systems varies in the order of CS<PEG<PLA.

CONCLUSIONS
The MD simulations were performed on the CS, PLA and PEG polymeric nanocomposite systems containing silica NPs and water molecules in order to attain the best DDS with the most proficient PIP drug delivery in the aqueous medium which is similar to the human body. It was found that CS had the greatest FFV but PLA exhibited the least FFV which reflected the diffusion and transport of the PIP drug in the CS would happen the most simply. The CS system showed the smallest Rg value while the PLA DDS had the greatest Rg confirming the CS chains had the closest contacts and intermolecular interactions with each other which led to their utmost vicinity that could cause the greatest FV and FFV values for the CS DDS among the three nanocomposites. The XRD patterns of all DDSs displayed one peak at about 20° approving the presence of CS, PLA, PEG and silica NPs in all of the nanocomposites. It was deduced from the RDF values that the oxygen atoms of carbonyl groups existing on the PIP molecules created very strong C=O…H-N hydrogen bonds with the amino hydrogen atoms of CS chains which could suggest the PIP molecules would diffuse the slowest and in a controlled manner inside the CS system. The highest diffusion coefficient was obtained for the PLA but the lowest D was achieved for the CS. The PEG DDS illustrated a D value of 0.0163×10–4 cm2/s which was between the values obtained for the CS and PLA justifying the PIP diffusion in the CS was the slowest which certified the most controlled and effectual drug delivery.

ACKNOWLEDGEMENTS
Authors would like to appreciatively express their thanks to the 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.

CONFLICT OF INTEREST
The authors declare that there is no conflict of interests regarding the publication of this manuscript.