J Mol Model (2015) 21: 165 DOI 10.1007/s00894-015-2705-2

ORIGINAL PAPER

Molecular dynamics simulation study of chitosan and gemcitabine as a drug delivery system Fariba Razmimanesh 1 & Sepideh Amjad-Iranagh 1 & Hamid Modarress 1

Received: 14 March 2015 / Accepted: 17 May 2015 / Published online: 6 June 2015 # Springer-Verlag Berlin Heidelberg 2015

Abstract By using molecular dynamics (MD) simulation, biodegradable biopolymer chitosan as a carrier for the drug gemcitabine was investigated and the effect of three initial drug concentrations (10, 40, and 80 %) on its loading efficiency was studied. Then water was added to the systems of drug and biopolymer and the effects of water on the interactions of drug and chitosan and on the drug loading efficiency were examined. From the results it was found that the maximum loading of the drug occurred at 40 % of the drug concentration. The radial distribution function calculations indicated that in the absence of water molecules, the drug molecules were located at shorter distance from chitosan and the loading efficiency of the drug in these systems was higher. Keywords Biopolymer . Chitosan . Gemcitabine . Drug delivery . Molecular dynamics simulation

Introduction Polymers are used widely in drug delivery systems, as they can encapsulate the drug and deliver them to the target cells and then release them gradually over a long period of time, whereby the repetitive dosing of the drug would be avoided. Polymeric drug delivery systems are able to improve therapeutic effects on the body by increasing efficiency of the drugs or reducing their side effects. Several polymeric drug delivery

* Hamid Modarress [email protected] 1

Department of Chemical Engineering, Amirkabir University of Technology (Tehran Polytechnic), Hafez Avenue, Tehran, Iran

systems such as nanoparticles [1, 2], micelles [3–6], hydrogels [7–12], and dendrimers [13] have been studied extensively. Usually these systems, include a combination of biocompatible and biodegradable polymer and an active pharmaceutical agent that is dispersed or covalently bound to the polymer. Usually penetration of drug through biopolymer is the cause of drug release and due to biodegradability of the biopolymers, they are highly considered as favorable carriers because they will be eliminated after the drug release [14]. Natural polymers such as polysaccharides [15], polypeptides, and phospholipids are commonly used as structural blocks of carriers [14]. Among various polymers, the biopolymer chitosan has excellent properties such as, film forming ability, mechanical tension, high permeability to water, high biocompatibility, non-toxic, chemical modification susceptibility, biodegradability [14–19], hydrophilicity [20, 21], antibacterial, haemostatic, fungistatic, and antitumoral [22, 23]. The hydrophilic nanoparticles can remain for long period of time in the circulating blood and can deliver the drug more effectively to the patient tissue [20]. Chitosan has the structure of α (1→4)-linked 2-amino-2deoxy-ß-D-glucopyranose and as a linear polymer is the main derivative of chitin [16, 17] (Fig. 1a). Gemcitabine is an anticancer drug which is effective in the treatment of a wide variety of tumors, where, the cytotoxic effects of gemcitabine are utilized by incorporation into DNA and in this case, it is important to obtain the best efficiency by using low doses of gemcitabine loaded on the polymeric nanoparticles as a delivery system to be targeted to the tumor sites [24–26]. Major problems of drug delivery systems are drug loss and unentrapped drug which require application of special techniques to recover the wasted drug. In some cases application of these techniques become very expensive due to high-cost of drugs [27].

165 Page 2 of 14

The stability of drug in human body is important. Stability is closely related to the physicochemical properties of the microenvironment in the delivery system, and to the interactions between the drug and the polymer [27]. Wei et al. [28] used salicylic acid-grafted chitosan oligosaccharide nanoparticles for delivery of the drug paclitaxel. They found that the drug concentration and the molecular weight of the polymer affected the drug loading efficiency. Their results indicated that increasing concentration of paclitaxel at constant polymer concentration, reduced the drug loading efficiency. Couvreur et al. [26] studied the effect of an anticancer drug, gemcitabine, concentration on the drug loading in chitosan combined with magnetite, and reported that the drug loading increased with increasing the feed ratios of the drug. High-priced drugs are accessible only in small quantities, therefore a serious limitation should be imposed on using drugs in standard experimental methods to optimize the drug delivery system. Molecular dynamics (MD) simulation can be used to study the structure and key properties of materials such as their stability, diffusion, and interactions between their molecules. MD simulation, without paying a high price for experimental methods, can be used to optimize drug delivery systems and subsequently to study the properties of the drug and its interaction with other molecules [27]. In addition, MD simulations can be applied to investigate the behavior of complex materials in molecular scale in various branches of sciences especially nanotechnology and biology [24]. Attempts have been made to study drug delivery systems by using MD simulations: Subashini et al. [27] investigated the interaction between drugs (doxorubicin, silymarin, and gliclazide) and six polymers (chitosan, alginic acid, sodium alginate, Gantrez AN119, Eudragit L100, and Eudragit RSPO) and they determined the best polymer for each drug by comparing the results obtained for the interaction energy and the percentage encapsulation efficiency. Macháčková et al. [29] studied six biodegradable polymers (chitosan, L-polylactide, D-polylactide, polyglycolic acid, polyethylene glycol, and cellulose) as carriers for cyclosporine A and they found a relationship between mixing energy, chain length, and chain flexibility of the studied polymer/drug delivery systems. Paclitaxel encapsulating within the salicylic acid-grafted chitosan oligosaccaride was studied by Wang et al. [30]. They examined the extent of interactions between paclitaxel and salicylic acid as a hydrophobic core and chitosan oligosaccaride as a hydrophilic shell in the drug delivery system. In this study, MD simulation is employed to investigate gemcitabine-loaded into chitosan polymer as a drug delivery system. In most simulation methods, the interaction between constituent molecules of the system have been studied in absence of water molecules but in real systems water molecules form the medium of drug delivery systems and due to their

J Mol Model (2015) 21: 165

polarity, they have an important effect on functioning of the systems. Therefore, the effect of feed drug concentration on loading efficiency in two cases; in absence of water molecules and in presence of water molecules were examined. The effect of water on interaction of the drug with the polymer and on loading amount of the drug were evaluated.

Simulation method MD simulation provides a clear understanding of molecular interactions. To investigate the effect of initial drug concentration on the loading efficiency of the drug, understanding the interactions between the components of the studied system is essential. Previous simulation studies indicated that condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS) force field effectively and accurately represents the molecular interactions [31–33]. Therefore, in this work the COMPASS force field was used to represent all the interactions between the polymer, as the carrier, and the drug in the studied drug delivery system. All simulations were performed by using Materials Studio 4.3 (Accelrys) at ambient temperature (298 K). The COMPASS force field is based on the polymer consistent force field (PCFF) model and has been widely used to optimize and predict the structural, conformational, and thermophysical properties of molecules in isolated or condensed phases including polymers [34, 35]. Ignoring the electronic motions, the total potential energy in the COMPASS force field consists of stretching of bond, bending angle, torsion, and nonbonded interaction terms. The functional form of COMPASS force field terms are given in reference [36]. Configuration of the polymer, chitosan Number of repeating units in a polymers for MD simulation can be obtained by considering the stability of solubility parameter (δ). δ is calculated according to the following equation [21, 33, 36]. δ¼

pffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi CED ¼ Ecoh =V

ð1Þ

where the cohesive energy density (CED) is evaluated as the molar cohesive energy (Ecoh) devided by molar volume (V) of the polymer. Polymers were first subjected to energy minimization and then the MD simulation was performed to find the solubility parameter. The evaluated solubility parameter were plotted versus number of repeating units of the polymer where a stable value of solubility parameter

J Mol Model (2015) 21: 165

Page 3 of 14 165

Fig. 1 a chemical structure of chitosan ((C6 H13 N O5)n) [18] b spatial structure of chitosan

(δ) indicated that the number of repeating units were proper for the simulated polymer [21, 33, 36]. For example, for the number of repeating units of chitosan equal to 20, a constant value of δ close to its experimental value, 10.6±0.8 (cal/cm3)0.5 was obtained [36]. In this study to confirm the accuracy of the simulation procedure, the simulations were run in a simulation box containing four chitosan chains where each chain had 20

Fig. 2 a chemical structure of gemcitabine (C9 H11 F2 N3 O4) (http://www.rcsb.org/pdb/ligand/ gemcitabine) b spatial structure of gemcitabine

monomers in its isotactic stereochemical structure with molecular weight of 3241.15 g mol-1. The density of this polymer determined experimentally is 0.67 g/cm3 [33, 36] at 298 K. For having energy minimization the polymer chain represented in Fig. 1b was subjected to the geometry optimization by the conjugate gradient [37] method. Then, the simulation box was minimized for energy by using subsequently the steepest descent [38] and conjugate gradient methods with the

165 Page 4 of 14

Fig. 3 Configuration of simulation boxes with 40 % of the drug concentration. a in absence of water. b in presence of water. (drug molecules are in red color)

convergence criterion of 1×10−5 kcal mol−1 Å−1. After this minimization step, the equilibration step was run for 200 ps by using NPT ensemble and the density and solubility parameter of polymer were obtained as 0.74 g/cm3 and 10.2 (cal/cm3)0.5 which were respectively close to the experimental values [33, 36] 0.67 g/cm3 and 10.6±0.8 (cal/cm3)0.5 . In the next step the drug, gemcitabine, and water molecule, were optimized by the same procedure as used for the polymer chain. The molecular structure of gemcitabine is shown in Fig. 2.

J Mol Model (2015) 21: 165

molecules. The drug concentrations on weight percent were respectively; 10, 40, and 80 % with the initial density of 0.3 g/cm3 (Fig. 3a). These simulation boxes are represented as B1, B2, and B3, respectively, for the drug concentrations of 10, 40, and 80 %. Also, three simulation boxes were constructed containing the same number of drug molecules and polymer chains as previous boxes but in presence of water molecules. The simulation boxes containing water molecules are represented as BW1, BW2, and BW3, which are respectively, for the drug concentrations of 10, 40, and 80 %. These three concentrations cover a wide range, for the amount of the drug in the simulation boxes and it is expected that the obtained results will show the effect of drug concentration on the desired behavior of chitosan polymer chain in the drug loading. We believe that these three concentrations sufficiently represent the simulation systems, as the final results would confirm this opinion. Then proceeding with the simulation process in this manner, would reduce the computational time and keeps the simulation procedure less expensive. The appropriate number of water molecules in these simulation boxes were evaluated by the Material studio module for maximum drug concentration (80 %) at 298 K and at constant pressure of 101.325 kPa. The results of calculation showed that the required number of water molecules to achieve a stable homogeneous equilibrium configuration for the simulation box, was 500. This number of water molecules was used throughout all simulations. The initial density of the simulation boxes were fixed at 0.5 g/ cm3 (Fig. 3b). The characteristics of the simulation boxes are summarized in Table 1. After constructing each simulation box, the periodic boundary condition was used and then to obtain proper blending of the components, each simulation box was examined for space filling regularly. If the three components were not well mixed in the initial configuration, the simulation box construction was rejected and a new construction was examined. This procedure was repeated for ten independently generated structures and then a proper structure was selected. The initially constructed atomistic configuration was subsequently optimized and then minimized to eliminate the undesirable contacts, including overlapping and close contacts, and to achieve this, a 10,000step energy minimization was adopted with the steepest descent and conjugate gradient methods. Molecular dynamics simulation

Simulation box In the present work, three chains of chitosan (20 monomer units in each chain) were put in each simulation box containing 4, 15, and 30 gemcitabine drug

MD simulations under constant temperature and pressure (NPT ensemble) have been performed for each configuration. Andersen thermostat [39] and Berendsen barostat [40] were respectively used for controlling the

J Mol Model (2015) 21: 165

Page 5 of 14 165

Table 1 Simulation boxes characteristics; drug (gemcitabine) concentration, number of drug molecules, number of water molecules, number of chains of polymer (chitosan), repeating units of polymer chain, number of atoms and box sizes (X, Y, and Z (Å)) Simulation box

Drug conc

Drug molecules

Water molecules

Chains

Rep. units

atoms

X

Y

Z

B1 B2 B3 BW1

10 40 80 10

% % % %

4 15 30 4

0 0 0 500

3 3 3 3

20 20 20 20

1445 1764 2199 2945

39.07 42.30 46.03 40.35

39.07 42.30 46.03 40.35

39.07 42.30 46.03 40.35

BW2 BW3

40 % 80 %

15 30

500 500

3 3

20 20

3264 3699

42.23 44.55

42.23 44.55

42.23 44.55

mass of drug Drug concentration ¼ mass of polymer  100

temperature (298 K) and pressure (101.325 kPa). A group-based cutoff of 15.5 Å with a splinter width of 5 Å was applied to calculate the non-bonded

interactions. The simulations were continued until the total energy and density of the system was stabilized at a constant value. The dynamic trajectory length of 800 ps was found to be sufficient for the system equilibration. The results were stored every 1 ps. Two snapshots of the simulation boxes B2 and BW2 are shown in Fig. 4.

Results and discussion To evaluate the effect of water on the gemcitabine loading and on the interaction strength between gemcitabine and chitosan, two simulation systems; one in absence of water molecules and one in presence of water molecules were designed. The dimensions of simulation boxes after performing the molecular dynamics simulation are presented in Table 2. As shown in Figs. 5 and 6, the energy and density have reached the stable state in the given simulation time (800 ps) and the equilibrium have been achieved in the course of simulation. Drug loading efficiency The loading efficiency is an important parameter which represents encapsulation capacity of the polymer as a drug carrier. The higher value of this parameter is an indication of the drug delivery system efficiency. The experimental loading efficiency (ELE) of drug is defined as [28]:

ELE ¼ ðdrug loaded on the carrierÞ=ðdrug initially fed Þ ð2Þ

Fig. 4 Simulation boxes with 40 % of the drug concentration after performing molecular dynamics simulation. a in absence of water. b in presence of water (drug molecules are in red color)

However, in molecular simulations, the simulated loading efficiency (SLE) of drug is conveniently defined as:

165 Page 6 of 14

J Mol Model (2015) 21: 165

Fig. 5 Energy profile of simulation box with 80 % of the drug concentration in presence of water

SLE ¼ ðnumber of drug molecules placed at predefined certain distances from polymer chainsÞ=ðtotal number of drug moleculesÞ

ð3Þ

To estimate the SLE of the studied simulation boxes, the number of drug molecules located at certain distances from chitosan chains were determined. As seen from Tables 2 and 3, the simulation boxes B2 and BW2 with 40 % of the drug concentration, have higher SLE, than other simulation boxes at the two predefined radii (approximately 3 and 4 Å).

By comparing the values of SLE in the Tables 3 and 4, it is found that in the simulation boxes B2 and BW2, the drug molecules are located at a closer distance from chitosan chains compared with the simulation box B1 and BW1. It seems that this behavior is controlled by the interaction of the drug molecules with the active sites on the chitosan chain, where the molecules of the drug are attracted toward the chitosan chains to occupy these sites. As seen from the results in Tables 3 and 4, at 80 % drug concentration, the SLE has decreased. This can be explored further by considering Eq. 3, where the number of drug molecules in the denominator has increased significantly, whereas, the numerator had a slight increase due to occupation of active sites in the chitosan chain. This point becomes evident by considering Fig. 7 which indicates that at 40 % of the drug concentration, the drug molecules are uniformly distributed closely along the chitosan chain, whereas, in 80 % of drug concentration, the drug molecules are mostly accumulated at farther distance from the chains. Therefore, it can be suggested that 40 % of the drug concentration is more beneficial to be used for efficient drug delivery purposes.

Table 2 The dimensions of simulation boxes after performing the molecular dynamics simulation

Fig. 6 Density profile of simulation boxes with 10, 40, and 80 % of the drug concentration. a in absence of water. b in presence of water

Box sizes (Å)

B1

B2

B3

BW1

BW2

BW3

X Y Z

24.52 24.52 24.52

27.48 27.48 27.48

27.70 27.70 27.70

30.05 30.05 30.05

31.10 31.10 31.10

32.60 32.60 32.60

J Mol Model (2015) 21: 165

Page 7 of 14 165

Table 3 Simulated loading efficiency (the number of drug molecules located at certain distances from chitosan chains Simulation box r (Å)

B1

B2

B3

r~3 (Å) r~4 (Å)

0.25 (1) 0.25 (1)

0.27 (4) 0.33 (5)

0.1 (3) 0.17 (5)

studied simulation boxes, where the subscripts i and j represent pair atoms in two molecules. The gi−j(r) are represented as follows: g Oc−Og ðrÞ : g N c−Og ðrÞ : g Og−Ow ðrÞ :

Radial distribution function (RDF), coordination number Radial distribution function (RDF) and coordination number (n) are two factors that can be used to study the position and distribution of specified atoms and their interaction with other molecules in the studied systems. The RDF represents the probability of finding an atom at a spherical shell of a certain thickness at a distance (r) from the reference atom which would generally be located at the origin. The RDF, g(r), can be expressed by the following equation [33]: gAB ð⋅ r⋅ Þ¼ < nAB ð⋅ r⋅ Þ > =4πr2 Δr ρAB

ð4Þ

whererepresents the average number of atom pairs between r and r+Δr and ρ is the number density of atom pairs A and B. Also, g(r) gives a measure of spatial conformation of atoms about the focal atom and demonstrates long-range order of the structure. Therefore, g(r) can be used to distinguish between amorphous and crystalline structures of various materials such as polymers. RDF can provide an insight as to how the atoms are located in an amorphous structure [33]. The position of the first peak obtained in the plot of g(r) versus r, represents the closest distance between the two atoms, and the probability that an atom appears at this distance would be demonstrated by the first peak height in this plot [36]. As shown in Fig. 1a, chitosan has one amine group (-NH2) and two hydroxyl groups (-OH) in each of its repeating units. Gemcitabine has one amine group (-NH2) and two hydroxyl groups (-OH) and one oxygen with double bond (-C=O) and fluorine atoms (F) that can participate in formation of hydrogen bond with hydrogen atoms. According to the following representations, the gi−j(r), radial distribution function (RDF), were calculated for the

Table 4 Simulated loading efficiency (the number of drug molecules located at certain distances from chitosan chains) Simulation box r (Å)

BW1

BW2

BW3

r~3 (Å) r~4 (Å)

0 (0) 0.25 (1)

0.13 (2) 0.27 (4)

0.06 (2) 0.1 (3)

g N g−Ow ðrÞ : g Oc−Ow ðrÞ : g N c−Ow ðrÞ :

ðOc Þ; oxygen atoms of hydroxyl groups of chitosan Og ; oxygen atoms of hydroxyl groups of gemcitabine ðNc Þ ; nitrogen atoms of amine groups of chitosan Og ; oxygen atoms of hydroxyl groups of gemcitabine Og ; oxygen atoms of hydroxyl groups of gemcitabine oxygen atoms of water ðOwÞ ; Ng ; nitrogen atoms of amine groups of gemcitabine ðOw Þ ; oxygen atoms of water ðOc Þ; oxygen atoms of hydroxyl groups of chitosan ðOw Þ ; oxygen atoms of water ðNc Þ; nitrogen atoms of amine groups of chitosan ðOw Þ ; oxygen atoms of water

The RDF analyses were performed at the simulation intervals where the simulation results showed a stable behavior and then, the last hundred ps of the simulation time was used for calculating RDF.

Simulation systems in absence of water molecules It is well known that when two hydroxyl groups in two molecules form hydrogen bond, the distance between the oxygen atoms of the two hydroxyl groups will be about 2.8 Å [36]. This distance is the same for hydrogen bond formation between atom pairs, nitrogen of an amine group in a molecule and oxygen of a hydroxyl group in another molecule [36]. Figure 8a shows the gOc−Og(r) for atom pairs Oc and Og and indicates that the first peaks are located at 2.6 Å. Therefore, hydroxyl groups of chitosan (OH)c interact strongly with hydroxyl groups of gemcitabine (OH)g by hydrogen bond formation. The RDF plot in Fig. 9a indicates that the height of the peak at 40 % drug concentration is higher than that of 10 % concentration. This, as is quite obvious, means that at higher concentrations more drug molecules are available to form hydrogen bond with chitosan. On the other hand, it is seen that at 80 % drug concentration, the height of the RDF peak has decreased. The interpretation for this behavior can be presented by considering the definition of RDF,(g(r)), as expressed by the following equation [41]: gðrÞ∼expðuAB =kT 0 Þ

ð5Þ

The above equation shows, that g(r) is expressed as (NAB/N) (the ratio of the number of molecular pairs A and B involved in the interactions to form hydrogen bond, over N the total number of molecules present in the system. Also g(r) is related to uAB, the interaction potential between molecules A and B. The results obtained for g(r), by considering Eq. 5, can be explained as follows; at about 40 % of the drug

165 Page 8 of 14

J Mol Model (2015) 21: 165

Systems containing water molecules Figure 9a represents the calculated values of gOc−Og(r) versus (r), and the first peak of RDF for 10 % of the drug concentration appears at a distance of 3.15 Å which is higher than the expected distance of 2.8 Å for hydrogen bond formation. The reason for this difference can be explained by considering competition of water and drug molecules for hydrogen bond formation with chitosan where, in this case, smaller water molecules are more available and are more susceptible, than drug molecules, to form hydrogen bond with chitosan. Therefore, water molecules get closer to the chitosan chain and push the drug molecules further apart. The first sharp peaks which belong to 40 and 80 % of the drug concentrations (in the boxes BW2 and BW3) are located at about 2.8 Å and the peak for the 40 % of the drug concentration has higher intensity. Again Fig. 9b (gNc−Og(r)) shows that; the first peak for 10 % of the drug in the simulation box BW1 appears at 4.95 Å, but for 40 and 80 % of the drug concentrations, in the simulation boxes BW2 and BW3, the first peaks appear at 2.8 Å, and the first

Fig. 7 a a cell with 40 % of the drug concentration showing distribution (b) a cell with 80 % of the drug concentration showing accumulation

concentration all the active sites for hydrogen bond formation have been consumed and no more active sites are available to form hydrogen bond. Therefore, while the total number of molecules (N) has increased significantly, due to higher concentration of the drug, NAB has slightly increased or remained the same as that of 40 % drug concentration. However, as shown in Fig. 8b, gNc−Og(r) for atom pairs Og and Nc, the first peaks which corresponds to 40 and 80 % of the drug concentrations are located at approximately 2.8 Å and this confirms conformation of hydrogen bonds between (-OH) group of the drug and (-NH2) group of chitosan. Also in Fig. 8b, the first peak which corresponds to 10 % of the drug concentration is located at 3.15 Å and, as seen in this figure, the height of the peak for 40 % of the drug concentration is higher than those of 10 and 80 %. These results, again substantiate our results obtained for the drug loading efficiency; that is, the 40 % of the drug (gemcitabine) concentration, is the optimum concentration to have the highest performance for chitosan, as a drug carrier, for this drug delivery system.

Fig. 8 a the RDF, gOc − Og(r), for oxygen atom of hydroxyl group of chitosan and oxygen atom of hydroxyl group of gemcitabine. b the RDF, gNc − Og(r), for nitrogen atom of amine group of chitosan and oxygen atom of hydroxyl group of gemcitabine (in absence of water)

J Mol Model (2015) 21: 165

Page 9 of 14 165

drug molecules can be loaded on the chitosan and then the unloaded drug molecules will be surrounded by water molecules.

Chitosan and water In Figs. 11a and b, gOc−Ow(r) and gNc−Ow(r) are presented and they demonstrate that the increase in the drug concentration has no effect on the interaction of chitosan and water molecules. That is, the water molecules move toward chitosan and take position at close distance from the functional groups of chitosan to form hydrogen bond. An average coordination number representing the number of B atoms surrounding an A atom can be obtained by evaluating the integral in the following equation [42]: NB hnA−B i ¼ V

ZRmin 4πr2 g ðrÞdr

ð6Þ

0

where NB is the total number of B atoms and V is the volume of the simulation box. Rmin is the distance of the first minimum

Fig. 9 a the RDF, gOc − Og(r), for oxygen atom of hydroxyl group of chitosan and oxygen atom of hydroxyl group of gemcitabine. b the RDF, gNc − Og(r), for nitrogen atom of amine group of chitosan and oxygen atom of hydroxyl group of gemcitabine (in presence of water)

peak for BW2 has higher intensity which is consistent with our previous results and the above mentioned interpretations. Gemcitabine and water As shown in Fig. 10a, the calculated RDF for oxygen atom pairs in gemcitabine and water, gOg−Ow(r), have an intense peak at 2.8 Å which corresponds to the hydrogen bond formation between (-OH)g and (-OH)w groups. The sharp peak which belongs to 80 % of the drug concentration (simulation box BW3) has higher intensity than the peaks for other concentrations (simulation boxes BW1 and BW2), which indicates that due to hydrogen bond formation in the simulation box BW3, more oxygen atom (Ow) in water molecules are located at a distance of 2.8 Å from oxygen atom of gemcitabine (Og). In Fig. 10b, the calculated gNg−Ow(r) for atom pairs Ng and Ow is presented, and it demonstrates that the highest peak belongs to the simulation box BW3. This means that more hydrogen bonds are formed between gemcitabine and water molecules in this simulation box. However, the RDF for simulation box BW3 in Fig. 10b indicates that by increasing the drug concentration, no more

Fig. 10 a the RDF, gOg − Ow(r), for oxygen atom of hydroxyl group of gemcitabine and oxygen atom of water. b the RDF, gNg − Ow(r), for nitrogen atom of amino group of gemcitabine and oxygen atom of water

165 Page 10 of 14

J Mol Model (2015) 21: 165 Table 6 First maximum and minimum in g Nc − Og (r) and the coordination number Simulation box

B1

B2

B3

BW1

BW2

BW3

First maximum (Å) First minimum (Å) < nNc-Og>

3.15 4.05 0.08

2.8 4.35 0.39

2.8 3.75 0.39

4.95 6.15 0.17

2.8 4.05 0.2

2.8 3.75 0.16

that, 40 % of drug concentration is the most appropriate concentration for this drug delivery system.

Vectorial autocorrelation function (VACF) The vectorial autocorrelation function (VACF) is a gage of angle by which the orientation of molecules changes within

Fig. 11 a the RDF, gOc − Ow(r), for oxygen atom of hydroxyl group of chitosan and oxygen atom of water. b the RDF, gNc − Ow(r), for nitrogen atom of amino groups of chitosan and oxygen atom of water

from the origin in g(r) plot and 〈nAB〉 is the average coordination number of B atoms around the A atom. In this study, andwere calculated and the results are listed in Tables 5 and 6. These Tables show that for 40 % drug concentration (the simulation boxes B2, BW2), < nOc-Og >andhave higher values than those for the simulation boxes B1and BW1. The values ofand corresponding respectively to 40 and 80 % drug concentration indicate that increasing concentration of the drug (from 40 to 80 %) has no effect on the coordination numbers and this confirms our previous results obtained for loading efficiency (LE) and radial distribution function (RDF)

Table 5 First maximum and minimum in g Oc − Og (r) and the coordination number Simulation box

B1

B2

B3

BW1

BW2

BW3

First maximum (Å) First minimum (Å)

2.6 3.45 0.05

2.6 3.45 0.24

2.6 3.15 0.2

3.15 3.45 0.003

2.8 3.15 0.08

2.8 3.15 0.082

Fig. 12 VACF for chitosan chain in the simulation boxes containing only chitosan and drug (Gemcitabine) in absence of water molecule for three concentration of the drug: a 10z, b 40 %, c80 %

J Mol Model (2015) 21: 165

Page 11 of 14 165

the chitosan chains with 10 % drug concentration are more flexible. However, comparing these results with those of other simulation boxes it is seen that by increasing the drug concentration to 40 and 80 % then to no significant effect occurs on the chitosan chain flexibility. The gyration radius The gyration radius is defined as [44]: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1ffi u0 X 2 u ðjri jÞ mi uB C uB i C Rg ¼ u@ X A t mi

ð8Þ

i

Fig. 13 VACF for chitosan chain in the simulation boxes in presence water molecule for three concentration of the gemcitabine drug: a 10 %, b 40 %, c 80 %

a given time and is defined as [43]: mðt Þ ¼< uðt 0 Þ:uðt 0 þ t Þ >

ð7Þ

where u(t0)) and u(t0 +t) are respectively the vectors characterizing the orientation of the polymer backbone or side chain at the beginning of the simulation run, t0 and t0 +t. In this work we used the vector indicating the end-to-end distance of chitosan chain to calculate VACF. For calculating VACF, the analytical Module of material Studio Software was utilized and the trajectories were obtained for 600 ps dynamics of each polymer chain in the studied simulation boxes. The results of VACF analysis are presented in Figs. 12 and 13 for the three chitosan chains. As seen from Fig. 12 a, b, and c which are respectively for the simulation boxes B1, B2, and B3, the range of variations of VACF are respectively 0.141, 0.009, and 0.041. Also as seen from Fig. 13 a, b, and c which are respectively for the simulation boxes BW1, BW2, and BW3, the range of variations of VACF are respectively 0.07, 0.02, and 0.02 which indicate that in both simulation boxes B1 and BW1, the range of variations are higher which means that,

Fig. 14 a the RDF, gOg − Hc(r), variation with simulation time for oxygen atom of hydroxyl group of gemcitabine and hydrogen atom of hydroxyl group for chitosan. b the RDF, gOg − Ow(r), variation with simulation time for oxygen atom of hydroxyl group of gemcitabine and oxygen atom of water. c the RDF, gOc − Ow(r), variation with simulation time for oxygen atom of hydroxyl group of chitosan and oxygen atom of water (for 40 % of the drug concentration in the simulation box BW2)

165 Page 12 of 14 Table 7 The radius of gyration Rg (Å) chitosan polymer chain for the studied simulation boxes

J Mol Model (2015) 21: 165

Simulation box

B1

B2

B3

BW1

BW2

BW3

Radius of gyration (Å)

11.2±1.9

12.8±1.1

12.3±2

10.9±0.5

11.6±1.9

11.9±0.8

where mi is the mass of atom i and ri is its vector position from the center of mass of molecule. The Rg for each chitosan chain is calculated and is represented in Table 7 for the studied simulation boxes. As seen from the result in this table, the simulation boxes B1 and BW1 have smaller Rg, but with increasing, the concentration of the drug in the other simulation boxes the Rg values increase which means that the rigidity of polymer chain increases as a result of drug presence in higher concentrations. However this increase is not significant for higher concentration of the drug (40 to 80 %). These results confirm our previous conclusion that, up to 40 % drug concentration more drug molecules take location in

vicinity of polymer chain to from hydrogen bonds. However, at higher concentrations, since all the active sites on the polymer chain have already been occupied, by the drug molecules, the added drug molecules can not interact with the polymer chain and therefore increasing the drug concentration has no significant effect on the flexibility of chitosan polymer chain. RDF variation with simulation time The RDF (gHc−Og(r)) variation for 40 % drug concentration versus simulation time has been plotted in Fig. 14a, which indicates the intermolecular interactions between oxygen atoms of gemcitabine (Og) and hydrogen atoms of hydroxyl groups of chitosan (Hc) changes with the simulation time. The results show that the first peak is located at about 1.8 Å and its intensity rises with the simulation time. Therefore, the probability of finding an Og atom at a distance of 1.8 Å from a Hc atom increases as the simulation proceeds. This indicates that more gemcitabine molecules tend to interact with the chitosan chains. In contrast, as shown in Fig. 14b for gOg−Ow(r) and in Fig. 14c, for gOc−Ow(r)), the first peaks are located at 2.8 Å and their intensities decrease with time. This indicates that, the probability of Og atoms being at distance of 2.8 Å from Ow atom, and the probability that the Oc atom being at distance of 2.8 Å from Ow atom decreases with time. The above mentioned results for loading of drug molecules on chitosan chains can be seen in the snapshot presented in Fig. 15.

Conclusions

Fig. 15 Configuration of simulation box with 40 % of the drug concentration in presence of water. a before molecular dynamic simulation. (b) after molecular dynamic simulation (drug molecules in red, polymer in yellow and water molecules in purple colors)

Biopolymers are known as the novel carriers for the drug delivery. In the present research the effects of drug (gemcitabine) concentration on the drug loading on the polymer (chitosan) in presence and in absence of water molecules were studied. The molecules were constructed through computer aided molecular modeling methods which were followed by MD simulations. In order to characterize the best drug concentration for efficient drug loading on the polymer, various parameters such as; loading efficiency, RDF and coordination number were evaluated. The molecular simulation results showed that the maximum loading of the drug occurred at 40 % of the drug (gemcitabine) concentration and the strongest interaction between chitosan and gemcitabine was observed for 40 % of drug concentration. In absence of water, the results indicated that more drug molecules were located at shorter distance

J Mol Model (2015) 21: 165

from polymer chains compared with the systems where the water molecules were present. It was also seen that in absence of water molecules the drug, gemcitabine molecules, tend to locate themselves in close vicinity of the chitosan chains. In contrast in presence of water molecule the drug gemcitabine molecules were located at farther distant from chitosan chain and instead water molecules take up position at close vicinity of chitosan chain.

Page 13 of 14 165 18.

19.

20.

21.

References

22. 23.

1. 2. 3. 4.

5.

6.

7.

8.

9.

10. 11.

12. 13.

14.

15.

16.

17.

Haley B, Frenkel E (2008) Nanoparticles for drug delivery in cancer treatment. Urol Oncol Semin Orig Invest 26:57–64 Faraji AH, Wipf P (2009) Nanoparticles in cellular drug delivery. Bioorg Med Chem 17:2950–2962 Torchilin VP (2007) Micellar nanocarriers: pharmaceutical perspectives. Pharm Res 24:1–16 Rapoport N, Pitt WG, Sun H, Nelson JL (2003) Drug delivery in polymeric micelles: from in vitro to in vivo. J Control Release 91: 85–95 Nishiyama N, Kataoka K (2006) Current state, achievements, and future prospects of polymeric micelles as nanocarriers for drug and gene delivery. Pharmacol Ther 112:630–648 Lin J, Zhu J, Chen T, Lin S, Cai C, Zhang L, Zhuang Y, Wang XS (2009) Drug releasing behavior of hybrid micelles containing polypeptide triblock copolymer. Biomaterials 30:108–117 Tang YF, Du YM, Hu XW, Shi XW, Kennedy JF (2007) Rheological characterisation of a novel thermosensitive chitosan/ poly (Vinyl Alcohol) blend hydrogel. Carbohydr Polym 67:491– 499 Ruel-Gariépy E, Leroux JC (2004) In situ-forming hydrogels—review of temperature-sensitive systems. Eur J Pharm Biopharm 58: 409–426 Kim SJ, Park SJ, Kim SI (2003) Swelling behavior of interpenetrating polymer network hydrogels composed of poly (Vinyl Alcohol) and chitosan. React Funct Polym 55:53–59 Hamidi M, Azadi A, Rafiei P (2008) Hydrogel nanoparticles in drug delivery. Adv Drug Deliv Rev 60:1638–1649 Gulrez SK, Al-Assaf S, Phillips GO (2011) Hydrogels: methods of preparation, characterisation and applications. Progress in molecular and environmental bioengineering-from analysis and modeling to technology applications. Prof. Angelo Carpi; In Tech, pp 117– 150 Ahmed EM (2015) Hydrogel: preparation, characterization, and applications. J Adv Res 6:105–121 Kavyani S, Amjad-Iranagh S, Modarress H (2014) Aqueous poly (Amidoamine) dendrimer G3 and G4 generations with several interior cores at Phs 5 and 7: a molecular dynamics simulation study. J Phys Chem B 118:3257–3266 Riva R, Ragelle H, Des Rieux A, Duhem N, Jérôme C, Préat V (2011) Chitosan and chitosan derivatives in drug delivery and tissue engineering. Adv Polym Sci 244:19–44 Park JH, Saravanakumar G, Kim K, Kwon IC (2010) Targeted delivery of Low molecular drugs using chitosan and its derivatives. Adv Drug Deliv Rev 62:28–41 Dash M, Chiellini F, Ottenbrite R, Chiellini E (2011) Chitosan—a versatile semi-synthetic polymer in biomedical applications. Prog Polym Sci 36:981–1014 Ravi Kumar MN (2000) A review of chitin and chitosan applications. React Funct Polym 46:1–27

24.

25.

26.

27.

28.

29.

30.

31.

32.

33.

34.

35.

36.

37. 38.

Jayakumar R, Prabaharan M, Sudheesh Kumar P, Nair S, Tamura H (2011) Biomaterials based on chitin and chitosan in wound dressing applications. Biotechnol Adv 29:322–337 Chiu YL, Ho YC, Chen YM, Peng SF, Ke CJ, Chen KJ, Mi FL, Sung HW (2010) The characteristics, cellular uptake and intracellular trafficking of nanoparticles made of hydrophobically-modified chitosan. J Control Release 146:152–159 Mitra S, Gaur U, Ghosh P, Maitra A (2001) Tumour targeted delivery of encapsulated dextran–doxorubicin conjugate using chitosan nanoparticles as carrier. J Control Release 74:317–323 Gu C, Gu H, Lang M (2013) Molecular simulation to predict miscibility and phase separation behavior of chitosan/poly (ECaprolactone) binary blends. Macromol Theory Simul 22:377–384 Vinsova J, Vavrikova E (2008) Recent advances in drugs and prodrugs design of chitosan. Curr Pharm Des 14:1311–1326 Kim IY, Seo SJ, Moon HS, Yoo MK, Park IY, Kim BC, Cho CS (2008) Chitosan and its derivatives for tissue engineering applications. Biotechnol Adv 26:1–21 Rungnim C, Rungrotmongkol T, Hannongbua S, Okumura H (2013) Replica exchange molecular dynamics simulation of chitosan for drug delivery system based on carbon nanotube. J Mol Graph Model 39:183–192 Arsawang U, Saengsawang O, Rungrotmongkol T, Sornmee P, Wittayanarakul K, Remsungnen T, Hannongbua S (2011) How Do carbon nanotubes serve as carriers for gemcitabine transport in a drug delivery system? J Mol Graph Model 29:591–596 Arias JL, Reddy LH, Couvreur P (2012) Fe3O4/chitosan nanocomposite for magnetic drug targeting to cancer. J Mater Chem 22: 7622–7632 Subashini M, Devarajan PV, Sonavane GS, Doble M (2011) Molecular dynamics simulation of drug uptake by polymer. J Mol Model 17:1141–1147 Wei XH, Niu YP, Xu YY, Du YZ, Hu FQ, Yuan H (2010) Salicylic acid-grafted chitosan oligosaccharide nanoparticle for paclitaxel delivery. J Bioact Compat Polym 25:319–335 Macháčková M, Tokarský J, Čapková P (2013) A simple molecular modeling method for the characterization of polymeric drug carriers. Eur J Pharm Sci 48:316–322 Wang XY, Zhang L, Wei XH, Wang Q (2013) Molecular dynamics of paclitaxel encapsulated by salicylic acid-grafted chitosan oligosaccharide aggregates. Biomaterials 34:1843–1851 Suknuntha K, Tantishaiyakul V, Espidel Y, Cosgrove T (2008) Molecular modeling simulation and experimental measurements to characterize chitosan and poly (Vinyl Pyrrolidone) blend interactions. J Polym Sci B Polym Phys 46:1258–1264 Prathab B, Aminabhavi TM (2007) Atomistic simulations to compute surface properties of poly (N-Vinyl-2-Pyrrolidone) (PVP) and blends of PVP/chitosan. Langmuir 23:5439–5444 Jawalkar SS, Raju KV, Halligudi SB, Sairam M, Aminabhavi TM (2007) Molecular modeling simulations to predict compatibility of poly (Vinyl Alcohol) and chitosan blends: a comparison with experiments. J Phys Chem B 111:2431–2439 Sun H (1998) Compass: an ab initio force-field optimized for condensed-phase applications overview with details on alkane and benzene compounds. J Phys Chem B 102:7338–7364 Rigby D, Sun H, Eichinger B (1997) Computer simulations of poly (Ethylene Oxide): force field, Pvt diagram and cyclization behaviour. Polym Int 44:311–330 Qiang L, Li Z, Zhao T, Zhong S, Wang H, Cui X (2013) Atomicscale interactions of the interface between chitosan and Fe3 O4. Colloids Surf A Physicochem Eng Asp 419:125–132 Hestenes MR, Stiefel E (1952) Methods of conjugate gradients for solving linear systems. J Res Natl Bur Stand 49:409–436 Fletcher R, Powell MJ (1963) A rapidly convergent descent method for minimization. Comput J 6:163–168

165 Page 14 of 14 39. 40.

41. 42.

Andersen HC (1980) Molecular dynamics simulations at constant pressure and/or temperature. J Chem Phys 72:2384–2393 Berendsen HJ, Postma JPM, van Gunsteren WF, DiNola A, Haak J (1984) Molecular dynamics with coupling to an external bath. J Chem Phys 81:3684–3690 Mc Quarrie DA (1973) Statistical mechanics. Hanper & Row, New York, p 263 Mahajan CV, Ganesan V (2010) Atomistic simulations of structure of solvated sulfonated poly (Ether Ether Ketone) membranes and their comparisons to nafion: I. Nanophase segregation and hydrophilic domains. J Phys Chem B 114:8357–8366

J Mol Model (2015) 21: 165 43.

Hu N, Chen R, Hsu A (2006) Molecular simulation of the glass transition and proton conductivity of 2,2′-benzidinedisulfonic acid and 4,4 ′ -diaminodiphenylether 2,2 ′ -disulfonic acid-based copolyimides as polyelectrolytes for fuel cell applications. Polym Int 55:872–882 44. Golzar K, Amjad-Iranagh S, Amani M, Modarress H (2014) Molecular simulation study of penetrant gas transport properties into the pure and nanosized silica particles filled polysulfone membranes. J Membr Sci 451:117–134