Hydrogen Bonding Regulates the Rigidity of Liposome‐Encapsulated Chlorin Photosensitizers

Abstract Liposomal formulations facilitate the administration of hydrophobic drugs, avoiding precipitation and aggregation phenomena when injected in polar solvents. The integration of the photosensitizer into the liposome may alter the fluidity of the system and, thus, modify the delivery process of the drug. Such a change has been observed for the liposomal formulation of Temoporfin, which is one of the most potent chlorin photosensitizers employed in photodynamic therapy. Here, all‐atom molecular dynamics simulations have been performed to identify the nature of the intermolecular interactions that might be responsible of the different lipids freedom of motion when the drug is introduced in the bilayer. It is found that Temoporfin participates as a hydrogen donor in strong hydrogen‐bonding interactions with the polar groups of the phospholipids. The theoretical analysis suggests that the rigidity of drug/liposome complexes can be modulated by considering the different hydrogen‐bond ability of the photosensitizer and the carrier material.


Introduction
Photodynamic therapy is at reatment that exploits the interaction of light with ap hotosensitive drug in the presence of molecular oxygen. [1][2][3] Such an interaction produces cytotoxic reactive oxygen speciest hat may induce the deatho ft he targeted tumor cells. [4,5] The advantage of photodynamic therapyl ies in its potential to achieve high selectivity towards the diseased tissues through the following therapeutic mechanism:t he drug is administered to the patient and, after as pecific interval, in which the drug accumulates in the tumor tissues,the interested region is irradiated with as ource of light with appropriate wavelength. [2,3,5] Chlorin derivatives are second-generation photosensitizers widely employedi na pproved clinical treatments and also in clinicalt rials. [1,6,7] These compounds are porphyrin-based species, [8] in which one of the double bonds has been reduced to enhancethe light absorption at lower energies, [9] in the so-called therapeuticw indow. [4] The best performance of second-generation photosensitizers over first-generation photosensitizers is well illustrated when comparing Te -moporfin (mTHPC) to Photofrin,w ith the former being 100 times more photoactive than the latter. [10,11] However,p orphyrin derivatives present severald isadvantages, especiallya tt he time of administration into patients. In polar media, likew ater or blood, porphyrins tend to form aggregates, [11,12] which are less photoactive [13] and considerably diminish the efficacy of the treatment. [14] This is the case for mTHPC( see Figure 1a), currently employed in photodynamic therapy to treat certain types of cancera nd precancerous conditions. [11,15,16] This drug is nowadays administrated as Foscan (biolitecp harma Ltd.) within an ethanol/propylene glycol mixture and it has shown aggregation and precipitation phenomena at the moment of injection, when the photosensitizer enters the blood stream. Liposomes have been identified as valuable nanocarriers of hydrophobic drugs, favoring the accumulationo ft he photosensitizer in the tumor tissues and solving the solubility issue in the administration phase of the treatment. [17] As shown in Figure 1b,liposomes consist of an assembly with ah ydrophilic inner cavity formed by the phosphate groups of the lipids, which contains water inside, ah ydrophobic inner part of the double layer formed by the alkyl chains, and an outer polar spherical surface also formed by the phosphate groups.
The efficacy of liposomal formulations of chlorin derivatives has been demonstrated for mTHPC, in which the drug has been embedded in the non-polar region of the liposome, as it is the case of the Foslip formulation. First, thesef ormulations solve the solubility problem, providing an easier administration of the drug to the patient,a st he final product can be dispersedi nw ater and administrated orally,i ntravenously,o r through skin permeation. [18,19] Second, Foslip and other liposome formulations show ah ighere fficiency with respectt o Foscan,p resentingalarger selectivity for tumor tissues with respect to healthy ones and provoking faster dissociation of Liposomal formulations facilitatet he administration of hydrophobic drugs, avoiding precipitation and aggregation phenomena when injected in polar solvents. The integration of the photosensitizer into the liposome may alter the fluidity of the system and, thus, modify the delivery process of the drug. Such ac hange has been observed for the liposomal formulation of Te moporfin, which is one of the most potent chlorin photosensitizers employedi np hotodynamic therapy.H ere, allatom molecular dynamics simulations have been performed to identifyt he nature of the intermoleculari nteractions that might be responsible of the different lipids freedom of motion when the drug is introduced in the bilayer.I ti sf ound that Temoporfin participates as ah ydrogen donor in strong hydrogen-bonding interactions with the polar groups of the phospholipids. The theoretical analysis suggestst hat the rigidity of drug/liposome complexes can be modulated by considering the different hydrogen-bond ability of the photosensitizer and the carrier material. the tumor cells, which resultsi nm ore effective growth inhibition. [17,[20][21][22][23] More importantly,l iposomal formulations of Temoporfin do not manifest any dark toxicity. [21,24] Finally,t he liposomal formulation ensures the presence of the drug in its monomeric form, which is the most photoactive species. [25,26] The factors that most influence the success of al iposomal formulation are the loading capacity of the carrier andt he release of the drug to the targeted cells, that is, the actual delivery process. Both factorsa re regulated by the interactions between the vehicle material and the drug. To provide an optimal formulation,s uch interactions should be of medium strength so that the drug is safely hosted in the carrier, but its releasei sn ot hindered. Intermolecular interactions, therefore, play av ery important role in the transport andd elivery of the drug, [27] and the elucidation of the nature of such interactions is relevant to designn anocarrierswith optimal properties. Different molecular compositions for possible mTHPCl iposomal formulationsh ave been investigated by differentials canning calorimetrya nd electron spin resonances pectroscopy. [28] The latter has shown that the photosensitizeri sl ocated in the hydrophobic part of the double layer and that its presence hinders the motiono ft he phospholipids, giving more rigid vesicles. [28] Thus, the drug-lipid interactionsl ower the molecular motion of the phospholipids. However,t he nature of these interactions, responsible for the high loading capacity of the carriera nd its different fluidity is unknown.
Coarse-grained methods are typically used to investigate the structuralb ehaviour of loaded and pure liposomes. [29,30] However,t hese methods are not able to provide atomistic details, such as interatomic interactions responsible for the final structural properties of the system. Therefore, in this work, we have performed all-atom classical molecular dynamics simulations to unveil the nature of the lipid-drug interactions that might be relevant in the liposomal formulation Foslip as am odel for chlorin-based photosensitizers. We show that hydrogen bonding between the photosensitizer and the phosphate heads of the lipids is the main interaction responsible for the large loading capacity of the carrier and the different mobility shownb y the phospholipids in calorimetric experimentsa nd electron spin resonance measurements. [28,31] 2. Results and Discussion

Convergencet oaSteady State
The first step in any molecular dynamics simulation of al arge system is to obtain an equilibrated structure. However,t his is very hard to achieve in liposomes of small size, whicha re metastable systems that converge very slowly to as teady state. [32] Liposomes with diameters smallert han 20 nm have a very high propensity to fuse with other liposomes, forming larger vesicless uch that the stress imposed by the highc urvature of small spherical shapes, which constrains the lipids to be very tightly packed, is released. [33] Therefore, the relatively small-sized liposomem odelled in this work (diametero f 12.6 nm;s ee Figure 2a)u ndergoes very strong curvature effects, which influence the packing of the phospholipids and make equilibration very difficult.
Commonly,c oarse-grained methodsa re employed to study lipid vesicles, allowing larger sizes (20-34 nm of diameter) and longer simulation times, [32,34] at the price of losing atomistic insight. Even in such cases,t he definition of at rue steady state is extremelyd ifficult and requires simulation times that can reach hundreds of microseconds. [34] Achieving atrue equilibrated structure, in fact, would requiret he simulationo fp henomena like the flip-flop of phospholipids between the inner and outer shell of the double layer as well as water exchange between the inner cavity of the liposomea nd the surrounding solvent. [32] In this work, we have monitored two parameters to analyze the structural quasi-equilibration of the system:t he area per lipid A L of the two layers of the membrane and the hydration number of the polar heads of the phospholipids. These two properties converge very slowly to as teady value in flat lipid membranes and are, therefore, employed as criteria to define the equilibrium of the system. [35] The area per lipid A L hasb een calculated considering the two spherical surfaces formed by the phosphorus atoms of the polar heads of the lipid chains, as shown in Figure 2a.T he radius of such spheres was calculated at intervals of 0.5 ns alongt he heatings imulation in the canonical ensemble (constant number of atoms, volume and temperature [NVT]) and at intervals of 10 ns along the simulation in the isothermal-isobarice nsemble( constant number of atoms, pressure and temperature [NPT]). The average radius (r) obtained for the inner and outer leaflets was employed to calculate the area of the two surfaces defined by the phosphate groups of the lipids (Figure 2a). Then, the area per lipid is given by the area of the surfaced ivided by the number of lipid moleculeso ft he leaflet (N L )[ Eq. (1)]: As seen in Figure 2b, A L presents strong oscillations along the NVT heating step, especially for the outer shell of lipid molecules. Instead, in the simulation at constant pressure (NPT), A L reaches ac onstant value almosti mmediately.T he area per lipid of the inner leaflet oscillates around 0.4 nm 2 with as mall deviation in the last 20 ns of the simulation, where it slightly increases up to 0.49 nm 2 .T his shows that af ull equilibrated structure is hard to achievef or relatively small liposomes,a s explained above. The A L value for the outers hell molecule is 1.02 nm 2 and remains stable throughout the whole NPTsimulation. The observed inner-leaflet structural fluctuationsi saconsequence of the curvature of the double layer,w hich influence the stabilityo ft he membrane inner leaflet to al arger extent. [32] The experimental A L for flat phospholipids double layers is found in ar ange of 0.55-0.72 nm 2 . [35] However,s uch values are only valid for very large vesiclesw here the surface locally approximate af lat bilayer.T he smalls ize of our theoretical model induces ah ighc urvature that constrains the inner shell lipids to be very tightly packed, whereas the ones in the outer layer are more loosely arranged. As there is no free lunch, such geometricala lterations introduce an error in the simulated A L .P revious molecular dynamics simulations, [29,30] employing a coarse-grained model tod escribe al iposome, reported A L valuesf or the inner and outer shells of 0.5 and 0.8 nm 2 ,respectively,t hat is, in better agreement with the experimental estimation, but still affected by the spherical shapeo ft he bilayer.
As intermoleculari nteractions are relevant for the structural properties of the bilayer, [36] we have selected the hydration number as the second parameter to monitor the equilibration process, whichd irectly correlatesw ith A L .T he hydration number represents the total number of water molecules per polar lipid head present in the first solvation shell of the liposome. This has been calculated by integrating the first peak of the radial distribution function of the water molecules and the negativelyc harged phosphate oxygen atom of the DPPC lipids. [35] As seen in Figure 2c,t he hydration number shows important oscillations during the heating simulation, but it convergest oavalue of 1.52 water molecules after 20 ns in the NPT simulation. We conclude, therefore, that the drug/liposome assembly is relativelyw ell equilibrated after 20 ns of dynamics in the NPT ensemble, attending to the time evolution of the area per lipid and the hydration number,w hich were recommended convergencec riteria. [35] Accordingly,t he following analyses are performed using framesf rom the dynamics between 20 and 100 ns.

Diffusion Coefficients
As shown by electron spin resonance spectroscopy experiments, [28] the integrationo fm THPC into the liposomei nduces al owering in the molecular freedom of motion of the lipids. This change in the mobility is ac onsequence of the different chemicale nvironmente xperienced by the lipid chains when they are in the presence of the drug versus when they are located in the pure liposome.T he appearance of new intermolecular interactions might induce an increaseo rd ecrease of the lipid motion, depending on the strength of the drug-lipid interactions. The lowering of the motion,e xperimentally observed, [28] indicates that the interactions in the loaded liposome are stronger than in the pure one.
To theoretically investigate the effect of the drug on the mobility of the lipids, we have calculated the root-mean-square displacement( RMSD) of the lipids along the 100 ns of the simulation in the NPTe nsemble by considering different regions of the liposome, definedb yt he separation between the lipids and ad rug molecule (see Figure 3a). We expect that the lipids in the immediate surroundings of the drug present less fluidity because of possible strong drug-lipid intermoleculari nteractions. First, the RMSD of the mTHPC molecules (without lipids) has been computed. Then, we have extended the analysist o the surrounding lipid molecules, gradually increasingt he distance from the drug molecules between2and 10 ,a ss ketched in Figure 3a.T he RMSD of the lipids located in the different regions,p lotted in Figure 3b,c learly evidences that the lipid molecules in the immediate surroundingso ft he drug molecules show less freedom of motion. In other words, the farther the residues are from the drug, the larger is their freedom of motion and, therefore, their displacementa long the simulation. This behaviori ndicatest hat the lipid molecules interact stronger with the drug than with other lipid molecules, and that their motion is reduced with respectt ot he situation where they are surrounded by only other lipid molecules. A more quantitative analysisc an be performed by calculating the diffusionc oefficient D from the mean-squared isplacement (MSD),a long the simulation time t,u sing the so-called Einstein relation[Eq. (2)]: [37,38] MSD ¼ 6Dt ð2Þ The MSD of the last 80 ns of the simulation is linearly fitted and the slope provides the diffusion constant D of the lipid residues located at different distances from the drug. Figure 3c shows that the variation of the diffusion coefficient with the separation r from the drug residues can be interpolated by the increasing exponentialf unction D(r) = Ae B(r) + D 0 ,w here A, B, and D 0 are fitting parameters with values of À0.008 2 ps À1 , À0.289 À1 ,a nd 0.01 2 ps À1 .T his indicates that the diffusion coefficient of the lipid chains increases with the separation between the lipids and the drug, and it reaches the value of 0.01 2 ps À1 at larger distances, where drug-lipid interactions are less significant or absenta nd are replaced by lipid-lipid interactions. In addition, the fast increase of D with the distance r from the drug residues (D reaches half of its maximum value already at ad rug-lipid separation of only 2.5 )c onfirmst hat the interactions controlling the system are of short-range character.
After confirming that the liposomep resents larger rigidity due to the intermolecular interactions between the drug and the phospholipids, the next step was to identify the nature of such interactions. To this aim, first, we investigated which fragments of the lipid chains,n amely heads or tails, play am ore significant role in the interactions with the drug molecules. Thus, the RMSD in the NPT ensemble is computed for the polar heads and the nonpolar tails of the phospholipids. In addition, this analysis ( Figure 4) is performed by considering two Figure 3. a) Schematic representation of the drug Temoporfin and the surrounding residues used in the RMSD analysis. The molecules represented with bolder lines are found in asphere of 2 radius centered at the drug. Thinnerl ines represent instead residues inside as phere of 10 radius. b) RMSD of the residues containedint he different size spheres surrounding the drug residuesd uring the NPT simulation.c)Diffusion coefficient reported against the separationfrom the drug molecules. The color scheme is the same as in (b).  Figure 3a). In the case of the residues in the immediate surroundings of the drug molecules (cutoff of 2 ,F igure 4a), the RMSD of the drug and the polar heads almosto verlap and oscillates aroundaRMSD of 13 . The nonpolar tails, instead, show al arger mobility with an average RMSD value of 20 throughout the simulation. The smaller mobility of the polar heads in comparison with the nonpolart ails clearly points to an existing interaction between the polar heads and mTHPC. Consideringn ow the RMSD of the residues at farther distances from the drug (cutoff of 10 , Figure 4b), the coupling between the polar heads and drug molecules is lost, and the tail and head residues show approximately the same freedomo fm otion.C learly,t he drug-head interactions, which are evident at short distances, vanish at larger separations, where the motion of the lipids is less perturbed by the presence of the drug.
The large number of hydrogen donora nd acceptor atoms present in both the drug and polar heads of the phospholipids (Figure 1a)a nd the short-range character of the drug/head interactions, as indicated by the rapid change of the diffusion coefficient of the lipids with the distance from the drug (Fig-ure 3c), suggest that hydrogen bonding is responsible for the rigidity in the proximity of the Te moporfin molecules.

Hydrogen-Bond Analysis
The presence of an etwork of hydrogen bonds between mTHPCm olecules and the polar heads of the lipids would explain the rigidity of the phospholipids in the direct proximity of the drug molecules. Hydrogen bonds are, in fact, interactions that are particularlyr elevant in ah ydrophobic environment, where water molecules cannot easily interfere. Hydrogen bonds are also somewhat directional interactions and can form rather ordered structures. Finally,t he presence of hydrogen bonds between the photosensitizer and the carrierm aterial would assure ag ood loading capacity of the vesicle. [27] The polar heads of the DPPC lipids are hydrogen acceptors, whereas the polar heads of the DPPG lipids and the drug moleculesc an act as both hydrogen donors and acceptors in the formation of hydrogen bonds. Therefore, we can classify the hydrogenb onds as drug-drug, head-drug, and head-head hydrogen bond pairs, see Scheme1.
According to the coupling of the polar head and the drug motion found in the analysis of the RMSD reported in The strength of ah ydrogen bond is intimately related with the geometrical parameters of the atoms involved in the interaction, specifically, with the hydrogen donor-hydrogen acceptor distance r D-A and with the angle q D-H-A formed by hydrogen donor,h ydrogen atom, andh ydrogen acceptor shown in Scheme 1d.T able 1s hows the total number of hydrogen bonds found along the full simulation, the number of hydrogen bondst hat last for more than half of the simulation time, the total number of fragment pairs that can potentially form Scheme1.Schematicrepresentation of the three categories considered for the hydrogen-bondanalysis:drug-drug pair forming ah ydrogen bond (a), drug-head pair involved in ah ydrogen bond (b), and two phospholipids connected through ah ydrogen bondbetween the polar head residues (c). The tails of the phospholipids are pictured in gray,whereas, for the residues involved in the bonds, the color scheme is red for oxygen atoms, blue for nitrogen,tan for phosphate,c yan for carbons, and white for hydrogens. d) A zoom on the atomsi nvolved in the hydrogenb ond(encircled in yellow in the previous panels), with the geometrical parametersshown.  Table 1, the number of long-lived head-drugh ydrogen bonds (27) is larger than the number of long-lived head-head hydrogen bonds (12), despite the number of head-drug pairs (113050) being smaller than the number of head-head pairs (116886). Althought he number of hydrogen bondsd epends on the number of speciesc onsidered in the box, the present numbers clearly indicatet hat the lipid chains are more prone to participate in hydrogen bondingw ith the drug than with other lipid chains. In addition, the head-drug hydrogen bondsa re not only more abundant buta lso stronger,a si ndicated by the smallest r D-A value (2.67 )a nd the largest q D-H-A value (1628).
Am ore quantitative picturei so btained by considering the distribution of the hydrogen-bond lifetimes and geometrical parameters obtained by the analysiso fa ll hydrogen bonds formed along the simulation. Figure 5a shows the number of drug-drug, head-drug, and head-head hydrogen bonds, grouped according to their lifetimes as ap ercentage of the total simulation time. Head-head hydrogen bonding dominates at short lifetimes (less than 10 %o fs imulation time), whereas they are surpassed by the drug-head hydrogen bonds at longer lifetimes, especially at lifetimes over 70 %o f the simulation time. The drug-head distribution in Figure 5a presents am inimum at al ife time of 50 %; such af eature is just ar esult of the way the distribution was built and has no physicalmeaning.
To comparet he strength of the interactions between the three different pairs of residues, the probability distribution of r D-A and q D-H-A is depicted in Figures 5b and 5c,r espectively. The r D-A distribution peaks at as horter distance for the drughead interaction (2.67 )t han those for the head-head (2.69 ) and the drug-drug (2.83 ), whereas the q D-H-A distribution peaks occur at larger angles for the drug-head interaction (1628 and 1648)t han for head-head (1618)a nd the drug-drug (1608)o nes.T herefore, both the hydrogen-bond lifetimes and the probability distributiono ft he geometrical parameters indicate that the hydrogen-bond strengths follow the order drughead > head-head > drug-drug. The strongerh ydrogen bonds found between Te moporfin and the polar headse xplain why the mobility of the phospholipids is smaller in close prox-imity to the drug. In addition, the hydrogen-bond behavior also explains the good loading capacity of the liposomes and the loss of mobility of severall iposomal formulations of mTHPC. [28] Finally,w ediscuss the reasonsb ehind the stronger interactions between the photosensitizer and the lipid heads than between different lipid heads. Ac loser inspection of the atoms involved in hydrogen bonding reveals that, in the drug-head interaction, mTHPCp lays the role of hydrogen donor from its phenolm oieties virtually all the time, whereas, in the headhead interaction, the hydrogen donor atoms are the Oa toms  . Owingt ot he electron-withdrawing effect of the aromatic ring of the phenolic substituent, mTHPCi sabetter hydrogen donor than the glycol moiety of the DPPG.I ng eneral, aromatica lcohols are better hydrogen donors than aliphatic ones. This is the reason why phenol acts as ah ydrogen donor and methanola sahydrogen acceptori nt he phenol-methanol dimer. [39] Based on our analysis, one can conclude that the rigidity of the liposomal formulation of mTHPC, closely related with the loading and release capacity of the liposome,c an be modulated considering the strength of the hydrogen bonds that the phenol moieties of the drug could form with phospholipids suitable for liposomal formulations. This consideration can be extended to different photosensitizers baring H-bonding groups in their substituents. In case of phenol moieties, furtherf unctionalization of the photosensitizer can be pursued. If the functionalization enhances the electron-withdrawing ability of the aromatic ring, the drug would be ab etter hydrogen donora nd the drug/liposome assembly would be more stable. Contrarily,i ft he aromatic ring loses its electron-withdrawinga bility upon functionalization, the drug is al ess effective hydrogend onor andt he drug/liposome assembly would be less stable. However,i ti si mportant to note that such af unctionalization should not modify other properties of the drug, for example, its hydrophobicity or the energy of the electronically excited states and, thus, the photosensitization efficiency and mechanism of the drug should not be altered.

Conclusions
In this work, we have performed all-atom molecular dynamics simulations to unravel the reasonsb ehind the smaller membrane fluidity experimentally found for several liposomeformulations of the photosensitizer Te moporfin, [28] ap rototypical chlorin photosensitizer.T he analysis of the motion of the different residues composing the liposome andt he calculation of diffusionc oefficients of the lipid chains in the immediate surroundings of the photosensitizer have shown that the interactions that confer rigidity to the drug/liposome assembly are of short-range character.T hese interactions are identified as relatively strong hydrogen bonds between the photosensitizer and the polar heads of the phospholipids. Owing to electron withdrawing effects, the drug is ap owerful hydrogen donor and, thus, the rigidity of the liposomei nduced by the drug into the phospholipids in its immediate surroundings derives from the tenacity of such interactions. Our theoretical analyses suggest that the stability of Temoporfin liposomal formulation could be also ascribed to the strengtho ft he hydrogen bond formed by the drug and the carrier material. Such considerations on the hydrogen-bonding capability of the drug and the phospholipids can be extended to other chlorinp hotosensitizers and, thus, can be relevant to design new drug formulations, as intermolecular interactions between the drug and the nanocarrier affect the loading capacity of the carrier, the drug releasep rocess, and the pharmacokinetic properties of the drug. [40] Computational Details The liposomei nitial structure was built with the Packmol software package, [41] with ac omposition of DPPC/DPPG phospholipids with ar atio of 9:1( w/w) and with am olar ratio of lipid/ mTHPCo f1 2:1. The inner leaflet of the liposomeb ilayer is formed by 118p hospholipid molecules confining4 72 molecules of water in the inner cavity.T he number of water molecules placed in the inner cavity was calculated so that the density of the solventi nt he inner cavity equals the one that surrounds the liposome. The second leaflet of the vesicle is assembled from 366 lipid molecules, resulting in al iposome with ad iameter of about 126 .T he liposome also contains 40 moleculeso fT emoporfin, placed in the hydrophobic region of the double layer,b ased on previous electron spin resonance spectroscopym easurements. [28] The resulting drug/liposome complex was solvated with 109 387 water molecules confined in a cubic box of 83 side. Na + counter ions were added to neutralize the DPPG negative charges and NaClw as added to reach as alt concentration of 0.15 M.
The energy of the system was first minimized in seven steps. In the first step, only 200 residues of the liposomew ere minimized,w hereas the other residues were spatially constrained by employing ah armonic force constant of 50 kcal mol À1 À2 . In each of the next fives teps, 200 more residues were allowed to move during minimization togetherw ith the ones already optimized in the previous steps.I nt he final minimization step, the whole liposomew as minimized withouta ny constrains. Each minimization calculation consists of 10 000 steps, where the first 5000 were performed by employing the steepestd escent algorithm whilet he last 5000 were performed by employing the conjugated gradient algorithm. The system was then slowly heated in the canonical ensemble (NVT ensemble) from 0t o3 00 Ki nf our different simulations of 1nse ach. In the first three heating simulations, the temperature was increasedb y1 00 Ki ntervals from 0t o3 00 K. Duringt he three heatings teps, the motion of the residues forming the liposome was constrained with decreasing harmonic force constants of 10, 5a nd 1kcal mol À1 À2 .T he last heating simulation was performed at af ixed temperature of 300 Ka nd with no constraints applied on the motion of the drug/liposome assembly.A fterwards, the system was run in the isothermal-isobaric ensemble( NPTe nsemble) to equilibrate the density of the system and the structure of the liposomef or 100 ns, and the snapshots were printed every 20 ps for analysis. As discusseda bove, after 20 ns, the system reached am etastable state and, therefore, the last 80 ns of the simulation were employedi nt he analyses.
For all of the simulations, periodic boundary conditions were applieda nd at ime step of 2f sw as employed.T he SHAKEa lgorithm [42] was used to freeze bonds involvingHa toms. The Langevint hermostat was used with ac ollision frequency gammao f1 .0 ps À1 while the pressure was kept fixed at 1atm with the Berendsen barostat and isotropic positions caling. For the non-bondedi nteractions calculation, ac utoff of 10.0 was employed together with the particle-meshE wald method to calculate the coulombic interactions. [43] The phospholipids were described by the Lipid17 [44] force field, available in the AmberTools17. [45] The intramoleculara nd Lennard-Jones parameterso ft he drug were taken from the generalA mberf orce field. [46] The restrained electrostatic potential( RESP) charges were computed by density functional theory, employing the B3LYP [47] functional andt he 6-31G* basis set [48] using the Gaus-sian09 [49] softwarep ackage. We employed the TIP3P model [50] for the water molecules and appropriate Amberp arameters to describe the NaCl ions. [51] All the simulations were performed with the GPU-based module PMEMD implemented in AMBER16. [45] The visualization of the simulationt rajectories and the RMSD calculations were performed with the visual molecular dynamics (VMD) program, [52] the hydrogen-bond analysis was performed with the CPPTRJ software [53] by employing ac utoffo f 3 for the D-A distance and 1358 for the D-H-A angle. By using the Einstein relation, Equation (2), the RMSD curves were linearly fitted to calculate the diffusion coefficient of the drug and lipid molecules.