Toward Understanding Curcumin Delivery by Trimethyl Chitosan Polyethylene Glycol: Optimization of Polymer Concentration and Chain Length by Molecular Dynamics

The second main cause of death in the world and one of the major public health problems is cancer. Curcumin is a natural bioactive substance with good anti-cancerous effect. However, due to the low cellular uptake of curcumin anti-cancer drug, it is vital to exploit a noble formulation, which can contribute to a decrease in its hydrophobicity and enables the efficient therapeutic effect of curcumin. Biocompatibility and hydrophilicity of the polyethylene glycol cause it to be one of the most attractive drug carriers. Chitosan is also of great importance, considering its biocompatibility, and is used along with the drug-carrying polymers. In this study, for the first time, a combination of trimethyl chitosan and polyethylene glycol was employed to deliver curcumin. Herein, hydrophilicity, stability, and energy analysis of the systems have been investigated, from which it was found that the 60/40 is the optimum ratio concentration of chitosan to polyethylene glycol for Curcumin delivery. Another characteristic property of the hybrid drug delivery system was the PEG chain length, with its least magnitude being the optimal value. Results of the present molecular study give a practical insight into the curcumin drug delivery system and propose a novel hybrid carrier for efficient curcumin delivery, which can be further exploited to develop novel nanomedicine systems.


Introduction
Cancer occurs by a series of sequential mutations in genes, which alter the cell functions.
Furthermore, lack of tumor suppressor genes triggers the uncontrolled cell division. As presented in various researches, different molecular methods can be adopted to investigate the potency of gene expression and defective proteins which are also effective in determining the novel cancer biomarkers [1]. As reported, most of the tumor cells are more round than the normal cells, which can be attributed to the fact that they have a lower interaction energy with the extracellular matrix and the neighboring cells [2]. In order to improve the efficiency of the therapy and to minimize the severe side effects, delivering drug through a smart carrier has been extensively researched [3,4]. The great achievements in smart drug delivering systems could enhance the therapeutic index of anticancer drugs, through optimizing their pharmacokinetics and tissue distribution and to modify the targeted delivery to the action site [5]. The rational design and 3 optimization of drug delivery system cannot be fulfilled without no knowledge of the loading capacity and distribution of drug molecules [6].
Curcumin is a polyphenolic hydrophobic substance derived from the rhizome of Curcuma longa, which is typically provided in forms of tablets, capsules, creams, extracts, gels, nasal sprays, etc. by pharmaceutical companies [7]. The structure of curcumin has been broken down in a systematic style to identify the noticeable sites of structural modifications (Fig. 1). The sites can be broadly categorized as the aromatic group, the di-keto moiety, the active methylene site, and the carbon linker. Apart from these sites, monocarbonyl analogues have also been extensively studied to improve the anticancer activity [8]. Curcumin possesses various functions in pharmacology and biology such as anti-cancer, antioxidation, anti-inflammatory, anti-tumor and anti-HIV, antimicrobial, anti-phlogosis, and antilipidemic effects, which have enabled it to act against anorexia, diabetes, coughs, hepatic disorders, alzheimer disease and rheumatism. However, its hydrophobicity and instability in exposure to different media have limited its extensive applications [9]. To administer the therapeutic agents which have low oral bioavailability, developing an effective drug delivery 4 system is an effective strategy which can bring about magnificent consequences in medicine [10]. Kim et al. (2011) have examined the cellular uptake kinetics of the natural form of Curcumin and Curcumin− Polyethylene Glycol (PEG), where, the cytotoxic effect of proliferating the preadipocytes and antiadipogenic properties in distinguishing preadipocytes have also been studied. The Curcumin and Curcumin−PEG were found to be differently absorbed onto 3T3-L1 preadipocytes and adipocytes with a limited amount of Curcumin−PEG absorption in the cell.
The improved water solubility of Curcumin−PEG can be correlated with an increased cellular retention of Curcumin in 3T3-L1 cells, particularly in preadipocytes [11].
Chitosan (CS), a cationic polysaccharide with its many advantages including good biocompatibility, biodegradability, antibacterial activity, and nontoxicity has been extensively applied in drug delivery systems [12]. Due to the large quantities of primary amino groups in chitosan, its solubility in water can be controlled; this pH-responsive biopolymer can play interesting roles in the fields of controlled drug delivery [13]. Wan et al. (2012) have synthesized Curcumin/chitosan micro-particles containing curcumin using an ionic cross-linking method with a maximum concentration of 270.24 ng/ml. The encapsulation efficiency and the extent of drug loading by Cur/chitosan micro-particles were reported to be 85% and 33.5%, respectively [14]. Shelma et al. (2013) have synthesized amphiphilic lauroyl-sulfated chitosan which loaded curcumin (Cur/LSCS). They also showed that the serum concentration of curcumin by oral administration of curcumin containing Cur/LSCS was greater and persisted for a longer duration of time (7 days) than that observed through oral (6 hr) or intraperitoneal (24 hr) administration of free curcumin. Moreover, 5 pharmacokinetic studies showed that the bioavailability of curcumin from Cur/LSCS was 11.5times greater than that of the free curcumin [15].
In this work, for the first time, potential of targeted delivering of curcumin by smart hybrid carriers made of PEG and chitosan has been studied. To achieve the optimum and controlled release of curcumin in cancer cells, the PEG to chitosan ratio and PEG chain length should be optimized, which have not been considered and studied elsewhere. The present work was mainly focused on investigating the optimization conditions, stability, and hydrophilicity of the drug delivery system. The role of molecular dynamics simulation will be highlighted due to its prediction capabilities prior to experimental studies, which not only saves time, material, lab work, and costs, but also is promising to reduce the drug side effects on patients. Moreover, molecular dynamic (MD) simulation can be used to clarify the interactions between the drug and carrier at molecular level, which is difficult to be observed by experiments, and MD simulation can reveal the trajectory of the atoms at spatial and temporal scales, providing detailed information about the conformational changes. For this objective, the gyration radius, Rootmean-square deviation (RMSD) have been investigated, as obvious clues to predict the stability of the accumulated molecules and a guarantee for the efficiency of the drug delivery process. Furthermore, to indicate the increased hydrophilicity of the drug-carrier complex, hydrogen bonding has been investigated.

Molecular Dynamics
The chains were sequentially optimized using Avogadro and HyperChem software's UFF and OPLS-AA force field. In the next step, optimizations were carried out using Gaussian 09 software by considering the b3lyp function and 6-31 + G* basis set. To carry out the simulations 6 we used GROMACS 2019.5 (OPLS-AA force field) in the EM (10 kJ/mol/nm minimum force), NVT (500 ps), NPT (500 ps) and then the MD (100 ns) simulation in 2 fs time steps. All cases contained 6000 water molecules (SPC/E water model) (Fig. 2B) in boxes with 7×7×7 nm 3 .
The cutoff radius was adjusted to 1.4 nm for the van der Waals and Colomb interactions. We used Colomb energy algorithm and Particle Mesh Ewald (PME). The MD simulation step was carried out with isotropic Parrinello-Rahman algorithm at 1bar and with nose-hoover (velocityscaling algorithm in NVT and NPT) at 300K. The constraint algorithm was based on the Links algorithm that was only used for hydrogen bonds. The partial atomic charges of the structures are calculated using Gaussian software (pop=esp).
Newton's motion equations, for all system particles, were applied in Molecular dynamics that can provide a set of consecutive atomic positions to predict the upcoming moments on the basis of the present condition [16]. The molecular dynamic study was done in various stages. In the first stage, the initial configuration of the particles was obtained. In the second stage, calculation of neighbors lists was made within the force range of the targeted atom in the system. Finally, according to the applied force on each atom based on the configuration, initial conditions and acceleration of each particle, calculations were done through integral methods. Newton's motion equation were solved in short time intervals (1-10 fs) for an integrated calculation. In each step, calculation of the applied force on atoms with the present situation and velocities was done to reach the upcoming positions and velocities in the following step. A set of forces was considered as the atoms displace to new positions and the process will be repeated. In this respect, molecular dynamics simulation can be a beneficial path to describe the dynamic variables changes by passing time [17]. 7 The accuracy and validity of simulated results were determined by using Force Fields.

Force Fields
Interactive energies (potential energy) of inter-particle will be changed by changing their distances. Equation (1) presents the potential function; while the force function of each "i" atom in a N-atom system is achievable from equation (2), which is extracted from the potential function. These equations were concomitantly solved. In addition, the force can be allocated to time and the atomic position by using equation (3) [18, 19]. Simple potentials such as the hard sphere potential can be applied in the primary molecular simulation. In this model, it was assumed that the particles move in straight lines with a constant velocity. When the distance of spheres becomes equal to the sum of their radii, relatively elastic collisions will happen. Then, a new velocity will be calculated based on the principle of conservation of linear motion size. By using the hard sphere model, beneficial results will be achievable, although in atomic or molecular system simulations it is not ideal.
Based on the Van der Waals potential, as interatomic or intermolecular distances vary, their forces change. However, no force was considered among particles unless they collide with each other. Van der Waals potential was represented in Equation (4), where "σ" shows the potential good depth and "q" denotes the distance at which the potential becomes zero. Calculations can be used to specify these parameters by Fitting with laboratory data or exact quantum chemistry.
"r" illustrates the distance between two atoms, and their inter-atomic potential was shown by "V" [20]: Drug diffusion coefficient can be calculated using Equations (5) and (6). Mean-square displacement (MSD) was calculated to reveal the drug diffusion coefficient; Coordinates of atoms are also in view as "r" while "t" shows time. After MSD calculation, diffusion coefficient can be calculated for a three-dimensional system by using Einstein's relation (Equation (6)) [21].

System preparation
According to similar structures available in the OPLSaa force field, Charge and other relevant parameters of the nanostructured functional groups were defined. Calculation of non-bonding interactions (i.e., electrostatic and Van der Waals) was done through applying Lenard-Jones and Columbian potential models. The input structure was prepared by implementing the OPLSaa force field. All molecules were placed in the box and tip3p water model was utilized as a solvent to obtain molecules parameters (changed to script format). 50000 steps were employed to minimize the energy of all simulation systems. To omit Van der Waals interactions and create hydrogen bonds between water molecules and other species, the steepest descent method was applied. Nose-Hoover algorithm was used to increase the temperature of the system from 0 to 310 K in 100ps in a constant volume in the next stage. Moreover, temperature-coupling systems were considered as 0.5 ps. Then at a constant pressure in 200 ps, the system was balanced.
Parrinello-Rahman algorithm was employed to balance the system pressure. The temperature of 37ºC for 50 ns was used in the MD simulation. Cut-off distance was set to 1.2 and the electrostatic force calculation was conducted by employing particle mesh ewald (pme). The LINCS (linear constraint solver) algorithm was adopted for maintaining the bond lengths; while the bonds engaged in hydrogen atom was limited by applying SHAKE algorithm which will accelerate the calculations.

Concentration Optimization
The 3D images of polymers and drug molecules are displayed in Fig. 2. The green, red, and blue molecules exhibit chitosan, curcumin, and poly ethylene glycol, respectively. Due to the high hydrophobicity of curcumin, it will result in several problems which limit its extensive use. As an instance, curcumin molecules would be aggregated in the blood stream so that their transport through cellular membrane would not be possible. The more hydrophilicity of curcumin, the better distribution and solubility in water can be attained which are essential in developing an efficient smart nanomedicine.
Curcumin uptake in the digestive system is extremely low. In general, hydrophobicity of curcumin cause it to be less biocompatible. Therefore, a great deal of research has been attracted to improve the curcumin biocompatibility. The use of polymeric carriers is one of the best approaches to enhance the curcumin hydrophilicity. Among different types of polymers, the combination of chitosan polymer and polyethylene glycol can be the best multi-objective option for the absorption, transfer, and release of curcumin; because chitosan polymer is extremely biocompatible and of a natural origin. Furthermore, chitosan is sensitive to pH and can control the release of drug by varying the pH which can be exploited to develop smart drug carriers to 10 load drug in normal cells with neutral pH and release drug at the cancerous cells with an acidic pH. This is due to the fact that pH of cancer cells is quite different from the pH of Non-cancerous cell in the body [13].
Chitosan has hydrophilic functional groups, the use of chitosan can improve the hydrophilicity, biocompatibility and provides a controlled mechanism for curcumin release. Meanwhile, polyethylene glycol is also highly effective in carrying the drug and is one of the greatly hydrophilic polymers, which can solve the problems regarding the hydrophobicity of curcumin.
Given that, the simultaneous use of polyethylene glycol and chitosan can be totally effective in stimuli responsive and smart uptake, transfer, and release of curcumin. However, to optimize the curcumin carrier, the optimum concentrations ratio of chitosan to polyethylene glycol and PEG chain length must be determined.
To this end, the curcumin drug absorption was investigated at different concentration ratios of chitosan and polyethylene glycol. The three-dimensional Fig. 2 show the images before and after the simulation at each concentration. In the images prior to each simulation, it is apparent that the drug molecules and polymers were in a scattered form and in the simulation final images, the drug adsorption onto the polymer is obvious in Fig. 2.

Hydrophilicity
Hydrogen bonding between two atoms is defined as a receptor-acceptor pair with an angle between them being less than 30 degrees. The graph of the number of hydrogen bonds over time for all simulations is given in Fig. 3. It has been reported that the drug loading rate is greatly influenced by the hydrogen bond between the drug and the carrier [23,24]. Moreover, Hydrogen bonding is directly associated with the hydrophilicity. The higher the number of hydrogen bonds,

14
Gyration radius is a factor that enables us to study and analyzes the aggregation of molecules such as polymers, and resizing of biological macromolecules such as proteins over time. The initial radius of the aggregated molecules of chitosan, polyethylene glycol and curcumin and the final radius of the so-called molecules after simulation are presented in Table 1. Specifically, the smaller accumulation radius of the molecules leads to a better and stronger accumulation.
Furthermore, the better aggregation between chitosan and polyethylene glycol and drug molecules causes the complex to be more stable and be better for delivering curcumin.
According to Table 1, all systems except the system of 60/40 experienced an increase in the gyration radius.

Energy analysis
The binding energy for four complexes was determined by using the MM/PBSA method. The summery of van der waals, electrostatic, and total energies of 50/50, 60/40, 70/30, 80/20 ratios of chitosan-PEG content, over 100000ps is presented in Table 2. The negative electrostatic energy indicates a strong attraction between the drug and the carrier. It can be seen that for all samples, the electrostatic forces were the dominant interactions. The result is consistent with [26], whose energy profile revealed insights by describing the van der Waals (VdW) and electrostatic forces for the interactions of Curcumin and chitosan, with the electrostatic energy being much more prevalent than the VdW.
Moreover, the higher electrostatic energy at a neutral pH leads to a higher loading of the drug onto the nanocarrier surface at a given pH. It can be noted that for the drug delivery system Trend of the associated electrostatic energy for the chitosan to PEG ratio of 60/40 over 100000ps, is given in Fig. 6. Fig. 6. The electrostatic energy diagram of 60/40 chitosan to PEG ration over 100000ps. improved. However, a decrease in protein bioactivity with PEG conjugates was found. Their study reveals that there should be a tradeoff between two opposing impacts of PEGylation.

Optimization of PEG length chain
Actually, finding out the optimal PEG length can maximize the overall therapeutic efficacy of the drugs [28].
Therein, the effect of PEG chain at 5, 8, 10, 12 monomers was investigated in presence of 5 molecules of curcumin and 5 molecules of chitosan. The 3D images before and after the simulation at each concentration can be seen in Fig. 6. In the images prior to each simulation, it 20 is apparent that the drug molecules and polymers are scattered and in the simulation final images, the drug adsorption on the polymer is visible. 21 22 Fig. 7. The drug delivery systems before and after simulation at different PEG chain length Fig. 7 shows that the shortest PEG chain resulted in the highest amount of Curcumin (red molecules) loading onto the chitosan-PEG carrier, which is in complete agreement with [29], which have found out that the higher PEG chain length had smaller zeta potential value. The drug loading capacity will increase by augmenting the zeta potential [30].

Stability
The diagram of the gyration radius corresponding to the effect of PEG chain simulations is plotted in Fig. 8. The only system that experienced a reduction in gyration radius correspond to 5 mer PEG chain. This result indicates that the highest aggregation occurred in this polymer length which enabled these polymers to properly surround the drug and to form a more stable complex. To evaluate the variations in the residue positions, the RMSD analysis was performed. Fig. 9 shows the RMSD for four ratios of concentrations. At the beginning of simulation, instability of system was detected, then after 500 ps the slope approached zero and the system stability was kept constant, in other words, stable complexes were created. The magnified range that the slope was very near to zero (3×10 -6 ). This system reached the stable mode in a quick time, which can be due to the small size of the molecules and the fast and strong interactions of molecules with each other. In addition, the maximum and minimum RMSD values for each system are given in Fig. 9.
Comparing the fluctuations in all systems, it can be inferred that the sample with the fewest fluctuations was formed by the PEG length of 5 and 8 mer. The results of the sample with 5 mer is in a complete agreement with the results of gyration radius and H-bonding analysis, suggesting a more hydrophilic and stable nanosystem for curcumin drug delivery. However, for the system with 8 mer PEG chain, both the results of gyration radius and H-bonding were not promising. In another word, the system with 8 mer was both less hydrophilic and stable. Therefore PEGylation of curcumin is a convenient method to improve its solubility [32]. Comparison of the graphs at Fig. 10 reveal that for all the PEG chain lengths, the number of hydrogen bonds was favorable, and it contributed to the hydrophilicity of curcumin. Average number of hydrogen bonds for four systems was 0.713857, 0.394717107, 0.627781, and 0.330639113. It can be seen that the highest average number of hydrogen bonds was achieved for the system with 5 mer PEG. In another word, increasing the PEG chain length did not lead to a higher hydrophilicity of 27 drug delivery system. Therefore, the optimum length for PEG chain is 5mer. This finding confirms the gyration results. While increment of PEG chain length from to 5 to 10, leads to a less reduction in hydrophilicity of samples in comparison to increments from 5 to 8, or 5 to 12 mer.

Energy analysis
Since the hydrogen bond energy still originates from the dipole-dipole interaction of the donoracceptor groups, being added to the electrostatic potential, it can be inferred that the electrostatic interaction energy ascends as the number of hydrogen bonds increases [26]. Accordingly, for the chitosan/PEG carrier with 5-mer PEG, which possessed the highest average number of hydrogen bonds, the highest electrostatic interaction energy could be predicted. The binding energy for the four complexes was determined by using the MM/PBSA method. The summery of van der Waals, electrostatic, and total energy for 5 mer, 8 mer, 10 mer, and 12 mer of PEG in chitosan-PEG hybrid system, over 120000ps is presented in Table. 3. The negative electrostatic energy indicates a strong attraction between the drug and the carrier.
According to Table 3, a comparison of the energy data shows that the average energy for the 5mer polymer is more than other chain lengths of these polymers; this finding confirms the Hydrogen bonds analysis. The aggregation of 5-mer polymers has a stronger electrostatic attraction than those of the 8-mer and 12-mer. The van der Waals attraction is also higher in 5mer polymers. It is worth mentioning that the data must be normalized for comparison. In the normalized state, the energy superiority of the 5-mer polymers is much more obvious. As a result, it can be stated that the polymers aggregation in the 5-mer polymers is much more stable and has higher attraction energy than other states. Stronger attraction energy indicates a higher stability of polymer aggregation and the higher attraction energy, the more stable carrier for drug delivery.
The portion of electrostatic energy is higher than that of the VDW energy. The electrostatic attraction between curcumin and chitosan made this aggregation more stable. This is due to the negative charge of curcumin (caused by hydroxide functional groups); because these hydroxide groups interact with the chitosan amine groups and have electrostatic attraction. The amine groups have a positive charge, and the negative charge absorbs the hydroxide groups.
Electrostatic attraction between the amine and hydroxide groups is also present in the 8-mer and 12-mer polymers, but the important point is the spatial accessibility of these functional groups to each other. The 5-mer polymers have the highest spatial accessibility and the least steric hindrance for interaction between hydroxide groups and amine groups, so in this case the electrostatic attraction is maximum. 29 This same spatial accessibility cause the van der Waals energy to be even greater. Because arondro-aromatic interactions of van der Waals and other types of van der Waals interactions are greater in this case. The spatial accessibility of polymers and drugs to each other is greater in the 5-mer state and more concentrated and closer to each other. As a result, electrostatic and van der Waals interactions were greater in this case.

Validation test
To validate the research work done for the two polymer ( chitosan and PEG),we chose a relevant recent article and repeated the RMSD analysis for PEG/drug, which is shown in Fi.5.b in [33]. Mina Mahdavi et al. (2020), have repeated the simulation four times to show that there may be some variations in the result. It can be seen that our result is in good agreement with its counterpart in regard to the diagram trend and the initial and end points.
The final and average values of RMSD in the present work and [33] have been presented in Table.4. We used the plot digitizer to calculate the average RMSD values. It can be notified that the results are in excellent agreement.These measurements can provide proofs of this work and show consistency with previous works, suggesting that an equilibration of the drug adsorption on the nanocarriers.  Fig.11 The validation test of PEG/drug a) Simulation of the present work b) Simulation of [33]

Conclusion
Curcumin possesses great therapeutic properties, especially anticancer effects; however, it has a low aqueous solubility which is challengeous. Prior to efficiency assessment of the hybrid drug delivery system via clinical trials, MD simulations not only can provide valuable information at molecular level, but also reduce the risks of human trials. Curcumin drug absorption was investigated at different concentrations ratios of chitosan to polyethylene glycol, and the results showed that the 60/40 ratio is the optimum concentration ratio. The gyration radius analysis indicated that chitosan, polyethylene glycol and drug molecules at concentration ratio of 60/40 have the strongest aggregation, and that composition is the most stable. The RMSD analysis asserted that these small molecules had fast and strong interactions and more importantly, the fewest fluctuations were ascribed to the 60/40 composition, which validates the gyration radius results. The hydrogen bonds in the drug delivery systems exhibited the hydrophilicity. The energy analysis suggests that for the optimum composition of chitosan and PEG, the highest negative electrostatic energy was achieved, which confirms the highest adsorption for efficient drug delivery system. The second factor in the design of the hybrid drug delivery systems was 31 the effect of PEG chain length on the adsorption of the drug. The highest adsorption was observed for 5 mer chain due to spatial accessibility, and interestingly the only system that led to a reduction in gyration radius was 5 mer PEG chain, which highlights the optimum and stable drug delivery system. The RMSD results indicate that 5 and 8 mer had low fluctuations. The 5 mer sample is in complete agreement with the results of gyration radius and H-bonding, suggesting a more hydrophilic and stable hybrid for curcumin drug delivery. Actually, the highest average number of hydrogen bonds correspond to the system with 5 mer PEG. The energy results confirmed that aggregation in the 5-mer polymers was much more stable and exhibited higher attraction energy than the other states. This is due to the negative charge of hydroxide functional groups of curcumin interacting with the chitosan amine groups.