1 Introduction

Dendrimers are dendritic polymers consisting of core, branch,and terminal groups where each branch layer is called ‘‘generation (G)’’ ^{[1, 2, 3]}. The architecture of dendrimers with hollow interior spaces allows their applications as potential nanoplatforms for drug delivery ^{[4, 5]}. Dendrimers help to improve solubility and pharmaceutical kinetics. Furthermore,dendrimer nanoplatform can be decorated with targeted functional groups and/or imaging agents for targeted therapy and molecular diagnosis ^{[6, 7, 8, 9]}.

Poly (amidoamine) (PAMAM) dendrimers are the first commercialized dendrimer family synthesized by Esfand and Tomalia ^[10]. PAMAM consists of a series of highly dendritic macromolecules with typical molecular weights 657-935 kg/mol (G0-G10). PAMAM consists of three basic parts with various functions. Firstly,the internal core leads to the differences in symmetry and has the electrostatic interactions with the dendritic repeated layers ^[11]. Secondly,the middle repetitive units can provide hollow spaces to encapsulate drug molecules. The middle parts are distinguished by the same concentric circles determined by the repetitive units as ‘‘generation’’ ^[12] (see Fig. 1),e.g.,G0, G1,G2,...,Gn are defined from the central core to the peripheral surface. The diameters of PAMAM dendrimers have approximately linear correlation with the generation of the dendrimers ^{[13, 14]}. Thirdly,the external surface terminated with amine groups can compound with various guest molecules,such as drugs,genetic materials,fluorescent molecules,and HBVc (FA,etc.) ^{[15, 16, 17, 18]}. The primary amino groups on the surfaces of PAMAM can be modified by different functional groups to be used as a multifunctional vector molecule in drug delivery systems ^{[19, 20, 21, 22]}. However, the charged surface groups may destroy the structures of cell membranes and induce the hemolysis of human red blood cells. Wang et al. ^[23] modified the surface of dendrimers with poly (ethylene glycol) (PEG) to reduce the cytotoxicity of dendrimers significantly. The PEGylation of the denrimer can also improve the aqueous solubility and circulation time in blood system of human body ^{[24, 25, 26]}.

Many anticancer drugs (such as CE6,DOX,MTX,and SN38) are limited in clinical use due to their poor solubility and biocompatibility ^{[27, 28]}. It is shown that the dendrimer-based drug-delivery systems can improve drug’s solubility and biocompatibility,realizing the targeted delivery and the controlled release of drugs ^{[13, 29]}. The influence of pH environment of human body can release the drugs carried in dendrimers in a controlled manner ^{[30, 31]}.

In dendrimer-based drug delivery systems,drug molecules can either bind covalently with dendrimers to form dendrimer-drug conjugates or incorporate non-covalently into dendrimers to form complexes ^[32]. In drug@dendrimer complex,drug molecules can be released more readily than conjugated system. On the other hand,it is not desirable if dendrimers release drugs too soon before reaching targeted cells ^[33]. Therefore,it is critical to understand the binding mechanisms of drug@dendrimer complex including binding strength and binding sites as well as the effects of surface functionalization. This work focuses on how binding strength and conformation vary with the number of loaded drug molecules.

In order to understand the mechanism of the interaction between the drug molecules and dendrimers at an atomic level,we carried out all-atom classical molecular dynamics (MD) simulations to study the dendrimers with and without the PEGylation. We load four types of anticancer drugs, including CE6,DOX,MTX,and SN38,respectively, gradually to saturated doses. In this work we address the issues of the binding strength and locations of the drug molecules as well as the sizes of complexes via the analyses of the binding energy (E_b) of drug molecules incorporated into dendrimers,the distribution of drug molecules, and the radius of gyration (R_g). 2 Computation methods 2.1 Computation models and simulation procedures

To study the influences of PEGylation of dendrimers as well as the types and numbers of drug molecules,we built full atomic models of non-covalently bound drug-dendrimer complexes,denoted as DN@PAMAM (or DN@GN) and DN@PAMAM-PEG (or DN@GN-PEG),respectively, for the dendrimers of generations G0-G3 with and without surface PEGylation. In the case of PEGylation,each terminal group binds covalently with a single PEG chain with atomic weight*2 000. DN represents drug molecules among one of the four types: CE6,DOX,MTX,and SN38, where N(0,1,2,4,...) is the number of drug molecules that compound with PAMAM. The drug@dendrimer complexes were first relaxed by geometry optimization (GO), followed by MD at the NVT ensemble at 298 K in vacuum. To examine the solvent effects,some systems,e.g., D₃₆@G3G3 and D₄₂@G3G3-PEG,were simulated in water solutions using MD at the NPT and NVT ensembles after geometry optimization.

Four types of anticancer drug molecules were studied in this work: Chlorin e6 (CE₆,C₃₄H₃₆N₄O₆,M=596),Doxorubicin (DOX,C₂₇H₂₉NO₁₁,M=543),Methotrexate (MTX,C₂₀H₂₂N₈O₅,M=454) and 7-ethyl-10-hydroxycampothecin (SN38,C₂₂₆H₂₀N₂O₅,M=392). The dendrimers G0(C₂₂H₂₈N₁₀O₄,M=516),G1(C₆₂H₁₂₈N₂₆O₁₂ , M=1 428),G2(C₁₄₂H₂₈₈N₅₆O₂₈,M=3 252) and G3(C₃₀₂H₆₀₈N₁₂₂O₆₀,M=6 900) were built using dendrimer module in materials studio (MS) software package. The charges of the dendrimers and drugs were assigned by Gasteiger method. And they were electrically neutral as a whole. The PEG (C90H182O45,M=1 982) chain was covalently bonded with the primary amine groups on the surface of dendrimer to build G0-PEG,G1-PEG,G2-PEG and G3-PEG. Then we added drug molecules to the hollow spaces of PAMAM to form the DN@PAMAM complex. To study the influence of solvation,we built the periodic cubic cell structure with a number of drug molecules carried in various generations of dendrimer where the distances between the adjacent complexes are 1 nm. H2O molecules were added by amorphous cell module of MS. The density of H2O was 1.0 g/cm³ . The CE616@G2-PEG,DOX16@G2-PEG,CE642@G3-PEG and DOX42@G2-PEG complex systems have 7 870,8 704,15 097 and 14 385 H2O molecules,respectively.

In this work we used the classical MD simulation software (forcite module implemented in MS ^[34]) to execute the GO and MD calculations. We used the Dreiding force field in both GO and MD computations ^[35]. The convergence tolerance of the energy and force in GO is 0.836 J/mol and 41.8 J/(mol·nm),respectively. The NPT MD was executed at 298 K and 1.013×10⁵ Pa,followed by NVT MD runs at 298 K for 1 ns. The other details of the MD simulation conditions are set as follows. (1) Dynamics parameters: in the NPT and NVT ensemble we use the thermostat (nose) to keep the target temperature equal to 298 K. The time step is 1.0 fs and the total simulation time is 1 ns. (2) Energy parameters: we use Dreiding force field to calculate the electrostatic,van der Waals and Hydrogen bond force. For these three long-range force terms,the cut-off of the distances are 1.55 nm,1.55 nm and 0.45 nm,spline widths are 0.1 nm,0.1 nm and 0.05 nm,respectively. The buffer widths are 0.05 nm and truncation methods are cubic spline. After equilibrium is reached,the last 0.5 ns in the MD trajectories was used for statistical analyses. 2.2 Definitions of binding energy and its associated components

To evaluate binding strength quantitatively,we calculated binding energy and its associated components including interaction energy and deformation energy. The relationships among these energies can be written

where E_b-PD is the instantaneous binding energy of drug molecules and PAMAM of the last snapshot of NVT MD; E_i-PD is the instantaneous interaction energy of PAMAM and drug molecules;E_d-p is the instantaneous deformation energy of PAMAM;E_di-D is the instantaneous deformationand-interaction energy of drug molecules. These energy components are defined in Eqs. (2)-(6). 2.2.1 Binding energy ofDN@PAMAMcomplex

where E_b is the average binding energy of DN@PAMAM complex;E_P-D is the total energy of DN@PAMAM complex averaged over equilibrium MD trajectories;E_D and E_p are the total energies of an isolated drug molecule and PAMAM averaged over equilibrium MD trajectories, respectively; Nis the number of drug molecules. 2.2.2 Instantaneous binding energy ofDN@PAMAM complex

where E_1-PD is the single-point energy of DN@PAMAM complex based on a single snapshot structure. 2.2.3 Instantaneous interaction energy betweenPAMAM and drug molecules inDN@PAMAMcomplex

whereE_1-D is the single-point energy of all drug molecules as a whole;E_1-P is the single-point energy of PAMAM. 2.2.4 Deformation energy ofPAMAMin DN@PAMAM complex

2.2.5 Deformation-and-interaction energy of drug molecules inDN@PAMAMcomplex

The binding energy (E_b-PD) defined in Eq. (3) is obtained based on the energy difference between thermodynamic equilibrium structures and often involves the deformation energies,while the instantaneous interaction energy (E_i-PD ) defined in Eq. (4) is calculated as the difference of singlepoint energies of the frozen snapshot structures so that pure interaction of drug-PAMAM can be evaluated alone. Therefore,both the binding energy and interaction energy are used in this work to understand the drug-dendrimer bindings from different points of views. The deformation energy (E_d-p) defined in Eq. (5) measures the energy change of PAMAM structures caused by deformation with respective to the corresponding equilibrium MD structures. The deformation-and-interaction energy of drug molecules (E_di-D ) defined in Eq. (6) evaluates the energy changes caused by both the deformation of drugs and the interactions among drug molecules. 2.3 Radius of gyration

R_g represents the size of DN@PAMAM complex,defined as follows

where M is the mass of the whole complex,m_i the mass of an atom i,r_i the position of an atom i,and R the center of mass of the complex ^[36].<>refers to average value in a NVT ensemble (T=298 K). 2.4 Diffusion coefficient of DN@PAMAM complex in water

where D is the total diffusion coefficient of the drug-dendrimer complex;r_i (t) and r_i(0) stand for the position at time t and 0,respectively. The coordinates of the last 500 ps in MD trajectories were used when the systems reached equilibrium in NVT simulations. 3 Results and discussions 3.1 Binding energy

Binding energy (E_b),defined in Eq. (2),measures the binding strength of dendrimer-drug complex including the interaction energies of drug-dendrimer and drug-drug (if any) as well as the deformation energies of drug and dendrimers. SmallerEbwith larger absolute value means stronger binding. When a single drug molecule is loaded into PAMAM or PAMAM-PEG,theEbscatters at various binding positions (see Figs. 2 and 3). The average E_b of the studied drug molecules decreases as DOX ((-16.23±3.09)×10⁴ J/mol)≈CE6 ((-15.94±3.17)×10⁴ J/mol)[MTX ((-13.49±2.16)×10⁴ J/mol)≈SN38 ((-13.44±3.16) ×10⁴ J/mol) in the D₁@G3 complexes. The negativeEb indicates that it is energetically favorable to incorporate these drug molecules into PAMAM. The standard deviations ofEbare*20% around their averaged value indicating the large site dependence of binding energy. Considering both interaction energy and deformation energy (confirmed later in Fig. 5 when N is small),we find that the interaction energy between PAMAM and drug is larger than the deformation energy of PAMAM leading to overall attractive PAMAMdrug interactions. This indicates that PAMAM-drug interactions are stronger than the intramolecular interactions among PAMAM branches.

In the case of D₁@G3-PEG (see Fig. 3),the averageEbof the studied drug molecules decreases as CE6 ((6.28±23.91) ×10⁴ J/mol) >DOX ((8.60±21.30)×10⁴ J/mol)> (MTX (11.06±17.76)×10⁴ J/mol)[SN38 ((11.83±23.30) ×10⁴ J/mol). The positive average E_b indicates that it is energetically unfavorable to incorporate drug molecules into PAMAM-PEG. Further analyses on deformation energies show thatEbbecomes positive mainly due to large deformation energy of PAMAM-PEG upon drug insertion. The drug molecules can bind with either PAMAM or PEG in D₁@G3-PEG. The earlier results show thatEbof drug-PAMAM is negative. Therefore,we conclude that positiveEbof drugPEG dominates and drug-PEG interactions are weaker than the intramolecular interactions among PEG chains.

When multiple drug molecules are encapsulated into PAMAM,we calculatedEbas functions of the number of drug molecules (N) shown in Fig. 4.The E_b fluctuates whenNis small as they scatter depending on various possible binding sites. However,E_b tends to approach towards constant values between (10.87-15.05)×10⁴ J/mol asNincreases beyond certain values,e.g.,~20 in DN@G3 or~30 in DN@G3-PEG. Both systems have similar trendsEb (CE6)≈Eb (DOX)>Eb(MTX)>Eb(SN38) when N is large.

On the basis of analyses above,the drug-PAMAM binding is stronger than the intramolecular interactions of PAMAM,but the drug-PEG binding is weaker than the intramolecular interactions among PEG chains. Although drugs bind with PAMAM energetically more favorable, when the number of drug molecules increase,drug should bind with PEG chains that provide more hollow spaces. Even though the drug-PEG binding energy is thermodynamically unfavorable,the complexes remain stable in the MD simulations,suggesting that the complexes are kinetically metastable once they form. 3.2 Instantaneous binding energy,interaction energy,and deformation energy of DN@PAMAM complex

In order to better understand the relationship among the four types of energies (see Eq. (1)),Eb-PD,Ei-PD,Ed-P,and Edi-D,we calculated and plotted theses energies for the DN@PAMAM in Fig. 5 and DN@PAMAM-PEG complexes in Fig. 6. The studied systems correspond independently to different numbers and types of drug molecules,so the trends obtained from these results are statistically meaningful.

As the number of drug moleculesNincreases,the instantaneous binding energies approach to constant values quickly after initial fluctuations at smallN(\12). We find that the trends of deformation energies of dendrimers coincide well with those of instantaneous binding energies,indicating that the conformations of dendrimers within the complexes are critical to determine the binding strength of drug-dendrimer. As the number of loaded drug moleculesNincreases,the deformation energies of drug molecules decrease and the interaction energies become weaker until both level off. The trends discussed above are qualitatively similar in both DN@G3 and DN@G3-PEG complexes. Quantitatively,the interaction energies of drug-(G3-PEG) are stronger than those of drug-G3,but the deformation energies of G3-PEG are larger than those of G3. 3.3 R_g

To examine the size changes,the radii of gyration of the drug molecules,dendrimers,and DN@G3 complexes are calculated and averaged over the equilibrium MD trajectories at 298 K (see Table 1). The R_g of isolated drug molecules (D) studied in this work is R_g (CE6)= 0.465 nm,R_g (DOX)=0.479 nm,R_g (MTX)=0.525 nm, and R_g (SN38)=0.434 nm,respectively. The R_g values of these drug molecules change littlewhen they are incorporated into PAMAM regardless of the number of molecules and PEGylation except that the largest floppy MTX molecules expand slightly 3%-8%. The whole D₃₆@G3G3 complexes expand 34%-48% relative to the pristine G3. The R_g of D₃₆@G3G3 is slightly larger or smaller than that of G3 embedded within the complex [G3(D₃₆@G3G3)] indicating that some drug molecules reside on the periphery of PAMAM. These results show that pristine PAMAM expands significantly when drug molecules are inserted while the drug molecules themselves are relatively rigid withoutlarge deformation.

When PAMAM is PEGylated,the R_g of D₃₆@G3G3-PEG increases only 2%-5% after 36 drug molecules are inserted. These results imply that PEG chains have more empty interior spaces to accommodate drugs than pristine PAMAM. The R_g values of D₃₆@G3G3-PEG are always smaller than those of G3-PEG within D₃₆@G3G3-PEG complexes,indicating that the drug molecules are encapsulated deeply inside G3-PEG.

When a single molecule is inserted,most G3 and G3-PEG shrink slightly (<2%). When many molecules are inserted (e.g.,N=36),the R_g of G3 increases 37%-46%, while that of G3-PEG increases only 8%-13%. These results indicate that small number of drug molecules shrink PAMAM as glue while many drug molecules expand PAMAM.

R_g values as functions of the number of drug molecules are plotted in Fig. 7 for DN@G3 and DN@G3-PEG. The R_g of DN@G3 is linearly proportional to the number of the drug molecules. The R_g of DN@G3-PEG fluctuates and has weaker dependence on the number of drugs. When the drugs are encapsulated into the G3,the R_g increases 34.3%-48.4% specifically from 1.151 nm (N=0) to 1.546-1.708 nm (N=36). The R_g of G3-PEG complexes increase 1.6%-5.0% specifically from 2.626 nm (N=0) to 2.668-2.756 nm (N=36). The R_g of G2(G2-PEG) is 0.896 nm (1.996 nm). The other results of G2(G2-PEG) system are similar to G3(G3-PEG) system.

3.4 Distributions of the distances between drugs and dendrimers

It is crucial to understand the binding sites of drug molecules embedded in dendrimers. The probability and density distributions of the distances between the centers of masses (COM) of the dendrimers and drug molecules are calculated based on the averaged equilibrium MD trajectories of D₃₆@G3G3 and D₄₂@G3G3-PEG complexes (see Figs.8 and 9). The results show that most drug molecules stay 1.5 nm away from the COM of G3. These distances are close to the R_g of complex 1.6 nm,indicating that the positions of drug molecules determine the R_g of the complexes and many of drug molecules expose on the surface of D₃₆@G3G3. On the other hand,most drug molecules are located 2.0 nm away from the COM of G3-PEG. These distances are smaller than theR_g of complexes 2.7 nm, implying that drug molecules are encompassed by the PEG chains and deeply embedded. Indeed,visual inspection of MD snapshots shown in Figs.10 and 11 supports the results of R_g. The information about the binding sites suggests that PEGlyation not only can reduce the cytotoxicity and improve the aqueous solubility of the dendrimers but also can directly interact with drug molecules. Therefore,designing functional groups of PEG chains is critical to develop efficient drug-delivery system.

3.5 Diffusion coefficient of DN@PAMAM complex in water

To examine the diffusion behaviors of dendrimer-based drug-delivery system (DDS),we calculated the diffusion coefficients Dof DN@G3 in water (see Table 2). The calculations are converged as shown in Figs. S1 and S2 in Supporting Information. The diffusion coefficients of DN@G3 and DN@G3-PEG complexes are 5.0×10^-4 -8.0×10^-4 ,which are smaller than the diffusion coefficient of water 2.0×10^-2 . Therefore,the dendrimer-based DDS diffuse slowly in solvents,which elongates the retention time of DDS in human body.

4 Conclusions

This work carried out molecular dynamics study on the dendrimer-based drug delivery systems focusing on the binding strength and conformation of DDS as functions of the numbers and types of drug molecules. The studied DDSs are PAMAM dendrimers loaded with the anticancer drug molecules DOX,MTX,SN38,and CE6,respectively. The PEGylated PAMAM was also studied to examine the effects of functionalization. On the basis of analyses of binding energies and their associated components including interaction energies and deformation energies,we find that the drugs bind with PAMAM more strongly than the intramolecular bindings of PAMAM,while the drugs bind with PEG chains more weakly than the intramolecular bindings among PEG chains. The deformed conformations of dendrimers embedded in the complexes are crucial to determine the binding strength of the drug-dendrimer complex. These results indicate that weakening the intramolecular interactions among functionalization groups of PAMAM can facilitate the incorporation of drug molecules.

Loading multiple drug molecules enlarges PAMAM dendrimers significantly but has much weaker effects on the PEGylated PAMAM. The drug molecules themselves remain similar sizes even after loaded into dendrimers. The PAMAM-PEG can accommodate more drug molecules and encompass drugs more deeply and strongly than pristine PAMAM. Our studies show that the drug molecules directly bind with PEG chains suggesting that designing the functionalization groups of PAMAM is critical to improve its drug bindings. The information extracted from this work provides helpful theoretical evidence to develop more effective dendrimer-based drug delivery systems.

Acknowledgments The authors are grateful for financial supports from ‘‘Shanghai Pujiang Talent’’ program (Grant No. 12PJ1406500), ‘‘Shanghai High-tech Area of Innovative Science and Technology (Grant No. 14521100602)’’,STCSM; ‘‘Key Program of Innovative Scientific Research’’ (Grant No. 14ZZ130) and ‘‘Key Laboratory of Advanced Metal-based Electrical Power Materials’’,the Education Commission of Shanghai Municipality; State Key Laboratory of Heavy Oil Processing,China University of Petroleum (Grant No. SKLOP201402001); National Natural Science Foundation of China (Grant Nos. 51202137,61240054,and 11274222). Computations were carried out at Hujiang HPC facilities at USST,Shanghai Supercomputer Center,and National Supercomputing Center in Shenzhen,P.R. China.