Xiaohui Wang ,Jiachen Li ,Li Zhang ,*

1 School of Economicsand Management,Tianjin University of Technology and Education,Tianjin 300222,China

2 Department of Chemistry,Zhejiang Sci-Tech University,XiaSha Higher Education Zone,Hangzhou 310018,China

Keywords:Water diffusion coef fi cient Controlled release Molecular dynamic simulation

ABSTRACT Poly-lactic acid(PLA)is w idely used as a controlled drug release material and the diffusion property of water w ithin the polymer matrix is closely related to the drug release pro fi le.This paper studies the w ater diffusion in PLA by molecular dynamic simulations.Free volume analysis indicates that w ater molecules are expected to fi ll in the free volumes of the polymer matrix forming w ater clusters at low w ater content.Along w ith the increase of the water concentration,the polymer starts to sw ell and the density of the system starts to drop.Due to the high mobility of water within water cluster,the calculated diffusion coef fi cient dramatically increases along with the incensement of w ater content.Thus,we conclude that the diffusion of water is a self-accelerate process,w ith higher mobility of water in the case where more water exists.

1.Introduction

During the last two decades,bio-based materials have aroused extensive interest in controlled drug delivery system[1-7].In contrast w ith traditional pharmaceutical preparations,controlled release formulation could release activeagent w ith appropriate rateaswell asbalance betw een the minimum effective dose and toxic threshold dose for a sustained interval of time[8].Thus,along-term activity can be achieved w hich w as essential especially for chronic diseases.A w idely used controlled release material is poly-lactic acid(PLA)[9]and its co-polymers,such as poly(lactic-co-glycolic acid)[10,11],poly(lactic acid)-poly(ethylene oxide)-poly(lactic acid)[12],and poly(lactic acid)-poly(ethylene glycol)[13].The most useful property of PLA and its copolymer is biocompatibility w ith tissues and biodegradability at a given time,leading to non-toxic sub-products which can be released or metabolized by the organism[14,15].

The release speed of drugs in polymeric drug delivery system can be in fl uenced by many aspects of the controlled release material,including polymer degradation,mass erosion as w ell as the diffusion of drugs in the polymer matrix[16].The insight of the drug release mechanism from thepolymer matrix can help usunderstand therelease behavior of controlled release systems.Mathematical modeling has been employed to investigate the releasing mechanism of drug molecules in controlled release materials[13-19].

Meanwhile,molecular dynamics(MD)simulation aroused our interest since it is effective to explore the drug releasing mechanism,especially for the diffusion properties of molecules in the polymer matrix,at the molecular level.The diffusion properties of small molecules in variouspolymer matrixes have been w idely studied in the past decades[20-26].More recently studies gradually start to investigate diffusion properties of larger molecules(drugs),like property of cholesterol in poly(methacrylic acid)hydrogel[27,28]and tetracycline in poly(styrene-b-isobutylene-b-styrene)[29].

Until now,the“hopping mechanism”is the mostly investigated mechanism for gas molecule diffusion w ithin the polymer matrix.Generally speaking,the diffusion of small molecule in the polymer matrix w as achieved by its hopping betw een adjacent free volumes.Zhao et al.have revealed the diffusion properties of(model)drug molecules w ith different sizes,like benzene,aspirin or nifedipine in poly(lactic acid-co-ethylene glycol)[30,31].They proposed a competitive mechanism depending on the size of drug molecules:small molecules tend to follow the“hopping mechanism”w hile larger moleculesmainly diffuse along w ith the w riggling of the polymer matrix.

There is abundant amount of w ater in human body environment;w ater could penetrate into the polymer matrix and lead to polymer sw elling during the release process of drug molecules.In other w ords,the release speed of drugs would be eventually in fl uenced by the diffusion process of drugs in the polymer matrix as w ell as the polymer sw elling.The molecular level mechanism underlying the w ater diffusion procedure in the PLA polymer matrix is still limited,as w ell as the study for large drug molecule diffusion in PLA w ith different w ater contents.In this w ork,the diffusion behavior of w ater molecules in PLA w as investigated through the MD simulation method.The aim of this study is to uncover the factors that in fl uence the w ater diffusion in the PLA polymer matrix.

2.Computational Methods and Simulation Details

2.1.Simulation details

Water penetration into PLA is a complicated procedure w hich is mostly composed of the follow ing tw o steps:in the fi rst step,w ater molecules diffuse into the polymer matrix(controlled release material)from the outside environment;in the second step,w ater molecules diffuse inside the polymer matrix.In this study,w e only focused on the second step by simulating w ater diffusion w ith a series of fi xed water concentrations.

The simulation box,w ith periodic boundary conditions,contained three PLA polymer chains with a polymerization degree of 500.MD simulations and data analysis w ere carried out w ith GROMACS 4.6.1 package[32].We used the PLAFF3 force fi eld which is developed specifically for PLAby McAliley and Bruce[33].Acut-off distance of 1.25 nm w as employed for the van der Waals interaction.The Particle-mesh Ewald method[34]was used to calculate the long-ranged electrostatic interactions w ith the same cut-off distance of 1.25 nm.The TIP4P w ater model[35]w as employed to simulate the explicit w ater in the system.The cross interaction parameters w ere obtained by geometry averages,

MD simulations w ere performed at 298 K,and the simulation time step is set to be 1 fs.Firstly,the PLAchains,together w ith the H2Omolecules,w ere optimized locally with a steep descent method at low density(typically＞10 timeslarger unitcell than theequilibrium unitcell after Np T simulation).Then w e employed an isothermal-isobaric ensemble(Np T)for 3 ns w ith 101.3 k Pa system pressure in order to compressthesystem to itsequilibrium density.Subsequently,canonical ensemble(NVT)simulation w as performed for 50 ns at 310 K and the trajectories of the last 30 ns NVT simulation w ere used to calculate the statistical properties.In order to reduce the statistical error,all the calculations w ere repeated three times from different initial structures and velocities.The w eight percentage of water was set between 1.6 wt%and 50 wt%.

2.2.Diffusion coef fi cient

The diffusion coef fi cients of the w ater molecules in PLA w ere calculated from the slope of the mean square displacement(MSD)as follow s,

w here D is the diffusion coef fi cient,t represents the time,and r(t)is the position vector of the center-of-mass of the molecule at time t.The bracket denotes an ensemble average.Generally,the MSD w ith respect to time is linear except in the fi rst range of picoseconds.Therefore,in this w ork,only the linear part w as used in the fi tting.

2.3.Radial distribution function and coordination number

The microscopic structure can be characterized using the radial distribution function gAB(r)betw een two kinds of atoms A,B.It can be calculated from the following equation,

w here 〈ρB(r)〉is the atom density of species Bat a distance r around atom A,and 〈ρB〉localis the atom density of species Baveraged over all spheresaround atom Aw ith half of the simulation box length.Thecoordination number w as obtained by integration of radial distribution function from zero to a required atom-atom distance.

2.4.Free volume

The free volume of the polymer matrix can be estimated via the volume unoccupied by the atoms of the system.But in this study,small free volume in the polymer matrix has little effect on the diffusion of w aters and the method mentioned above cannot distinguish the size of the free volume.Another method w as used here to get a more accurate result.The simulated boxed w as divided by mesh points w ith about 0.05 nm along each direction.At each mesh point,a test ball with radius in range of 0.2-0.9 nm w as inserted to test w hether it overlaps w ith existing atoms.Thus,the free volume of each system can be calculated by the ratio of the mesh points which do not overlap with the polymer matrix.

3.Results and Discussion

3.1.Simulation results

3.1.1.Diffusion coef fi cient

The diffusion coef fi cients of w ater molecules in the polymer matrix w ere successfully obtained during the simulation from mean square displacement according to Eq.(3)(selected MSD pro fi les are show n in Supporting Information Fig.S1).As show n in Fig.1 and Table 1,the diffusion coef fi cientschange along with water content,but it isnot alinear curve,and this curve can be approximately divided into tw o phases.At low w ater content(＜5 w t%),the calculated diffusion coef fi cients are small and in a narrow range,betw een 0.03 and 0.04 × 10-9m2·s-1.The diffusion coef fi cients increase dramatically w hen the w ater content exceeds 5 w t%,and the curve increases nearly linearly w hen the w ater content is above 17 w t%.The maximum diffusion coef fi cient is 2.10 × 10-9m2·s-1,w hich is about tw o magnitudes larger than the value at the low est w ater content.How ever,it is still smaller than the diffusion coef fi cient of pure w ater(～4×10-9m2·s-1)[36],w hich indicates water mobility is still con fi ned owing to the existence of polymer even w hen the w ater content is larger than 50 w t%.

Fig.1.Simulated diffusion coef fi cient(10-9 m·s-1)of w ater in PLAw ith different w ater concentrations.

Table 1 Summary of calculated diffusion coef fi cient(D),density(ρ),coordination number(CN)and free volume of PLA w ith different contents of w ater

3.1.2.Density of the system

As show n in Fig.2,the calculated densities of PLA in different concentrations of w ater varied betw een the ranges of 1.09-1.19 g·cm-3.One interesting observation is that the density of the system slight increases at low w ater concentration,w hile dramatically decreases at high water concentration.The density at the lowest water concentration(1.6 w t%)w as 1.18 g·cm-3w hich w as very close to the calculated density of pure amorphous PLA,1.15 g·cm-3[25]and slightly low er than the experimental reported crystal PLA density 1.26 cm-1[37].The slightly low er simulated density might be due to the defect of the Np T compression procedure,w hich is unable to fi nd the global minimum of this system w ith higher density.The density increased slightly along w ith the w ater concentration rise and the maximum density(1.19 g·cm-3)is observed w hen w ater concentration increases to 6.26 w t%,w hich might be because of the fi lling of w ater molecules in the free volume of the PLA matrix.Afterw ards,the system density started to drop dramatically and decreased to 1.09 g·cm-3at 50 w t%w ater concentration,because of the low density of pure w ater.

Fig.2.Density of amorphous PLA with different w ater concentrations.

3.1.3.Radial distribution function and coordination number

The radial distribution function(RDF)is an important characteristic to illustrate the microscopic structure of liquids.The RDFs betw een oxygen atoms(OW)in w ater molecules are calculated to reveal the structure of water molecules around the polymer networks.The position and height of peaks in gOW-OW(r)curve for bulk w ater are in good agreement w ith Sarupria and Garde's w ork[38].From Fig.3,it can be found that the shape of gOW-OW(r)curves for bulk water and w ater in polymer netw ork is similar.The fi rst peak is at about 0.28 nm for all the systems.The second peak is located at around 0.45 nm.By comparing the peak height at 0.28 nm and 0.45 nm,w e observed that the relative peak height(@0.27 nm/@0.45 nm)is decreasing,indicating theordered rangeisincreasing along with the increaseof water content.

Fig.3.Selected RDFof w ater oxygen-w ater oxygen pairs in PLA with different w ater concentrations.

In order to further investigate the structure difference of w ater clusters in the polymer matrix,the accumulative coordination number of OW-OWpairsin the fi rst water shell wasfurther calculated,i.e.counting OW-OWpairsw hosedistanceissmaller than 0.37 nm,middleof the fi rst and second peaks in the RDF.At low w ater concentration(1.6 w t%),this value is only 2.3.But it dramatically increased to 3.7 at 9.1 w t%w ater concentration.Afterw ards,this parameter continues to increase although the increasing speed decreases(Fig.4).In the range of current simulation,it reaches to 5.6 at 50 w t%w ater concentration,still low er than the value in pure w ater(6.1,simulation results).

Fig.4.Calculated water oxygen-w ater oxygen coordination number in PLAwith different water concentrations.

3.1.4.Free volume

Fig.5.Percentage of free volume(test ball size is 0.2 nm)for PLA polymer matrix with different water concentrations.

The free volume of the system w as calculated according to the method described in Section 3.Two typesof free volume were calculated,either w ith w ater or w ithout w ater in the systems.Fig.5 show s the free volume in case contains w ater molecules,it could be observed the decreasing of the free volume w hen the w ater concentration is low.When the w ater concentration is higher than 6.3 w t%,the free volume of the system starts to fl uctuate around 4.2,w hich is expected to be the inaccessible area of w ater in PLA.

We further calculated the free volume without water to roughly estimate the(sphere)w ater cluster size in the system.When more w ater enters the system,large w ater clusters could be formed w hich is re fl ected by the existence of free volume in the case w ith large test balls.Detailed results are listed in Table 2.

3.2.Discussion

Previous studies[19,25,26]indicate small molecules tend to follow the“hopping mechanism”w hen they diffuse in the polymer matrix.Indeed,hopping phenomenon could be observed at the low water concentration case in our simulation.As show n in Fig.6a,w ater molecules tend to occupy the free volume of thepolymer matrix,and it can mostly fl uctuatein asmall freevolumemost of thetime,asindicated by thetwo black areas in the trajectory.The line w hich connects these tw o black areas indicates the water molecule hopping from one free volume to a neighbor free volume.Because the hopping chance is low,the w ater diffusion coef fi cients remain low in the similar range,c.a.0.04×10-9m2·s-1,w hen w ater concentration is low er than 5 w t%.

As discussed above,the diffusion coef fi cients increase dramatically w hen the w ater concentration is higher than 6.26 w t%,it suggests that the w ater diffusion process w ould accelerate.This may be caused by the relatively larger water cluster formation when more water exists in the PLA matrix,w hich might be due to the hydrophobic nature of the PLAchain.The polymer matrix starts to swell and polymer density starts to decrease ow ing to the relative strong interaction betw een w ater molecules and polymer matrix,this is consistent w ith the observation that the accessible free volumes are mostly occupied at higher w ater content(Fig.5).Meanw hile,analysisalso show sw ater cluster becomes larger w ith higher w ater concentration in the system(Table 2),w hich is also consistent w ith our coordination number analysis.As show n in Fig.6b,the trajectory of w ater molecules is more randomly distributed w ithin the PLA matrix,consistent w ith the larger free volumes at high w ater content(Table 2).At low w ater concentration,water cluster issmall with large percentage of“boundary”water molecules,w hich could not form a complete H-bond coordination netw ork with other water molecules as in bulk water.Thus,the average coordination number in this case is smaller than the value in bulk w ater.Forexample,in our low est w ater concentration,thisvalueisonly 2.3.When more water existsin PLA,larger water clustersform in the polymer matrix.Our calculation results fi nd thisvaluereachesto 5.6 at 50 w t%w ater concentration.How ever,due to the interface between w ater and polymer matrix,thisvalueisstill lower than 6.1 which istheresult from simulation of pure w ater.Increasing formed coordination bond means larger water clusters start to form in the polymer matrix.Due to the high mobility of w ater molecules w ithin w ater clusters,w ater diffusion becomes easier in the large w ater clusters and the diffusion coef fi cient starts to linearly increase after water concentration＞5 wt%.

Table 2 List of calculated free volume only counting polymer matrix(i.e.free volume+w ater volume)

Fig.6.Selected trajectory of one water molecule in(a)1.6 wt%and(b)50.0 wt%water content of PLA.The straight lines are auxiliary lines.

Based on our simulation,w e proposed a self-accelerate diffusion process.In the early stage,the w ater diffusion is a slow process,w ith w ater molecules mainly fl uctuating w ithin small free volume,w ith a small chance to hop betw een neighbor free volumes.Due to the hypothetic nature of PLA,w ater w ill form clusters even at very low w ater concentration.From the free volume analysis and coordination number discussed above,it could befound that thesizeof water clusters expands w ith the increase of the w ater content in the system.There larger water clusters force the PLA chains to move aw ay from each other and generate larger free volumes w ithin the PLA matrix.Meanw hile,larger free volumes serve as new channels for w ater diffusion w ithin the PLAmatrix.The diffusion of water within the PLA matrix is much faster and easier than in low w ater content.In other w ords,more water molecules facilitate w ater transportation within the PLA matrix.Thus,the w ater diffusion in the PLA matrix is a self-accelerate process.

4.Conclusions

In this study,w e successfully simulated the PLAw ith different w ater contents by molecular dynamic simulations,and calculated the diffusion coef fi cient of w ater moleculesin PLAw ith different w ater contents.We observed low water mobility,i.e.small diffusion coef fi cient,in the case w ith low w ater concentration.In the case of the highest w ater concentration,the calculated diffusion coef fi cient is about tw o magnitudes larger than the value with the lowest water amount.

Free volume analysis indicates that w ater molecules are expected to fi ll in the free volumes of the polymer matrix and it has less effect to expand the volume and low er the density of the system at low w ater content.Afterw ards,the polymer starts to sw ell and the density of the system starts to drop.In the meanwhile,the size of water cluster is increasing along w ith the incensement of w ater content,w hich is re fl ected by the increase of average w ater oxygen-w ater oxygen coordination number as w ell as free volume analysis in the case w ithout including w ater.Due to the high mobility of w ater w ithin w ater cluster,the calculated diffusion coef fi cient dramatically increases along w ith the incensement of water content.Thus,w e conclude that the diffusion of water isaself-accelerate process,with higher mobility of water in the case w here more w ater exists.

Nomenclature

D diffusion coef fi cie nt

ρ density

SupplementaryMaterial

Supplementary data to this article can befound online at https://doi.org/10.1016/j.cjche.2018.09.009.