Computational Study of the Interaction of a PEGylated Hyperbranched Polymer / Doxorubicin Complex with a Bilipid Membrane

Fully atomistic molecular dynamics simulations are employed to study in detail the interactions between a complex comprised by a PEGylated hyperbranched polyester (HBP) and doxorubicin molecules, with a model dipalmitoylphosphatidylglycerol membrane in an aqueous environment. The effects of the presence of the lipid membrane in the drug molecules’ spatial arrangement were examined in detail and the nature of their interaction with the latter were discussed and quantified where possible. It was found that a partial migration of the drug molecules towards the membrane’s surface takes place, driven either by hydrogen-bonding (for the protonated drugs) or by hydrophobic interactions (for the neutral drug molecules). The clustering behavior of the drug molecules appeared to be enhanced in the presence of the membrane, while the development of a charge excess close to the surface of the hyperbranched polymer and of the lipid membrane was observed. The uneven charge distribution created an effective overcharging of the HBP/drug complex and the membrane/drug surface. The translational motion of the drug molecules was found to be strongly affected by the presence of the membrane. The extent of the observed changes depended on the charge of the drug molecule. The build-up of the observed charge excesses close to the surface of the polymeric host and the membrane, together with the changes in the diffusional behavior of the drug molecules are of particular interest. Both phenomena could be important at the latest stages of the liposomal disruption and the release of the drug cargo in formulations based on relevant liposomal carriers.


Introduction
Vesicles based on lipids are recognized as promising non-viral vectors for drug or gene delivery purposes [1][2][3][4][5].In particular, phospholipid liposomes are considered as suitable vehicles for pharmaceutical compounds due to their biodegradable and biocompatible nature, which is combined with low toxicity levels [6].In such liposome/drug formulations, the nature of the interactions between the guest molecules and the lipid vesicle is of paramount significance since it is intimately related to the ability of the carrier to minimize the loss of the bioactive cargo during the transfer to the desired destination, as well as to the mechanism of its liposomal escape at the target site [7,8].
Molecular simulations have started playing an increasingly important role in the description of the interactions of bioactive molecules with lipid layers [9][10][11][12], due to their ability to provide information at the molecular or even at the atomic level.Among other model phospholipid layers, dipalmitoylphosphatidylglycerol (DPPG)-based systems have been recently studied in drug/liposome formulations [13].DPPG membranes (formed by a double layer of lipids, i.e., a bilipid layer) carry a Fluids 2019, 4, 17 2 of 15 negative charge per lipid at physiological pH conditions [14,15] which renders them good candidates when complexation with cationic bioactive compounds is desired [13,14,16].This attribute of DPPG lipids can be exploited in more complex drug/liposome formulations, where optimization of the drug loading conditions is promoted by the presence of multifunctional molecules such as hyperbranched polymers [4, 17,18].
In this work we examine by means of fully atomistic molecular dynamics simulations, the interactions of a complex comprised by a functionalized (PEGylated) hyperbranched polymer (HBP) and an anticancer drug, doxorubicin, with a model bilipid DPPG membrane in an aqueous environment.Structural and dynamic characteristics of the hyperbranched-polymer/doxorubicin aqueous solution have been investigated in our previous work [19].We henceforth aim at elucidating the effects of the presence of the DPPG bilayer to the drug spatial distribution with respect to the polymer and the lipid membrane, the nature of the drug/DPPG interactions, the effective charge distributions related to possible overcharging phenomena of the polymer and of the bilipid DPPG layer, and the changes in the translational behavior of the drug molecules imparted by the presence of the bilipid membrane.

Materials and Methods
Initial configuration of the DPPG bilipid membrane was taken from a pre-equilibrated system of a several hundred ns long simulation, from a publicly available repository [20].Pre-equilibration of the hydrated DPPG membranes were performed by Jämbeck and collaborators [21][22][23] in explicit TIP3P [24] water solution.Figure 1 illustrates a single DPPG lipid, a doxorubicin molecule, and the structure of the PEGylated hyperbranched polymer.
Fluids 2019, 4 FOR PEER REVIEW 2 drug/liposome formulations [13].DPPG membranes (formed by a double layer of lipids, i.e., a bilipid layer) carry a negative charge per lipid at physiological pH conditions [14,15] which renders them good candidates when complexation with cationic bioactive compounds is desired [13,14,16].This attribute of DPPG lipids can be exploited in more complex drug/liposome formulations, where optimization of the drug loading conditions is promoted by the presence of multifunctional molecules such as hyperbranched polymers [4, 17,18].
In this work we examine by means of fully atomistic molecular dynamics simulations, the interactions of a complex comprised by a functionalized (PEGylated) hyperbranched polymer (HBP) and an anticancer drug, doxorubicin, with a model bilipid DPPG membrane in an aqueous environment.Structural and dynamic characteristics of the hyperbranched-polymer/doxorubicin aqueous solution have been investigated in our previous work [19].We henceforth aim at elucidating the effects of the presence of the DPPG bilayer to the drug spatial distribution with respect to the polymer and the lipid membrane, the nature of the drug/DPPG interactions, the effective charge distributions related to possible overcharging phenomena of the polymer and of the bilipid DPPG layer, and the changes in the translational behavior of the drug molecules imparted by the presence of the bilipid membrane.

Materials and Methods
Initial configuration of the DPPG bilipid membrane was taken from a pre-equilibrated system of a several hundred ns long simulation, from a publicly available repository [20].Pre-equilibration of the hydrated DPPG membranes were performed by Jämbeck and collaborators [21][22][23] in explicit TIP3P [24] water solution.Figure 1 illustrates a single DPPG lipid, a doxorubicin molecule, and the structure of the PEGylated hyperbranched polymer.Figure 2 depicts the equilibrated DPPG membrane comprised by 128 lipids.As a next step, the bilipid membrane was combined with the HBP/Doxorubicin complex that was taken from our previous study [19] with the aid of the VEGA ZZ software [25].The HBP/Doxorubicin complex was initially placed on top of one side of the DPPG bilayer with a separation of approximately 65 Å between the center of mass of the HBP and that of the membrane (no constraints have been imposed to any of the constituents during the equilibration or the production runs).For comparison purposes, we have kept exactly the same number of neutral and protonated drug molecules and the same ionic strength of 165 mM as in Reference [19].Figure 3 depicts the initial configuration of the bilipid/hyperbranched/drugs/counterions system, with (a) and without (b) the water molecules.The number of the different components comprising the examined system are presented in Table 1. Figure 2 depicts the equilibrated DPPG membrane comprised by 128 lipids.As a next step, the bilipid membrane was combined with the HBP/Doxorubicin complex that was taken from our previous study [19] with the aid of the VEGA ZZ software [25].The HBP/Doxorubicin complex was initially placed on top of one side of the DPPG bilayer with a separation of approximately 65 Å between the center of mass of the HBP and that of the membrane (no constraints have been imposed to any of the constituents during the equilibration or the production runs).For comparison purposes, we have kept exactly the same number of neutral and protonated drug molecules and the same ionic strength of 165 mM as in Reference [19].Figure 3 depicts the initial configuration of the bilipid/hyperbranched/drugs/counterions system, with (a) and without (b) the water molecules.The number of the different components comprising the examined system are presented in Table 1.The forcefield parameters for the hyperbranched polymer and the drug molecules were taken from the generalized Amber forcefield (GAFF) [26] and for the water molecules from the TIP3P model [24] as described in Reference [19].For the bilipid membrane, the GAFF lipid set of parameters was utilized, which was developed specifically for the description of lipids [27].All simulations were performed by the aid of the NAMD simulation package [28] with a timestep of 1 fs and a saving frequency of 2 ps, using periodic boundary conditions.Temperature control was performed using the Langevin method (with a damping coefficient of 5 ps −1 ) [29].Electrostatic interactions were calculated via the particle mesh Ewald (PME) scheme [30], while the distance cutoff for the van der Waals interactions was set at 12 Å.The percentage of the protonated doxorubicin molecules was calculated through the Henderson−Hasselbalch relation [31] taking into account the pKa of the drug's ionizable group at 37 °C and considering physiological pH conditions.An appropriate amount of chlorine and sodium ions was added for the overall neutralization of the system and for accomplishing an ionic strength of 165 mM [32].More detailed information regarding the construction and the equilibration of the hyperbranched polymer/drug complex can be found in Reference [19].Following the construction of the system, energy minimization with steepest descent and conjugate gradient methods was applied, prior to 40 ns long MD equilibration runs in the canonical (NVT) ensemble at T = 310 K.This procedure ensured stabilization of energetic (total energy), thermodynamic (pressure), and conformational properties (i.e., average size of the PEGylated hyperbranched polymer) of the polymeric component of the system.The equilibrium area per lipid was approximately 75 Å 2 and the The forcefield parameters for the hyperbranched polymer and the drug molecules were taken from the generalized Amber forcefield (GAFF) [26] and for the water molecules from the TIP3P model [24] as described in Reference [19].For the bilipid membrane, the GAFF lipid set of parameters was utilized, which was developed specifically for the description of lipids [27].All simulations were performed by the aid of the NAMD simulation package [28] with a timestep of 1 fs and a saving frequency of 2 ps, using periodic boundary conditions.Temperature control was performed using the Langevin method (with a damping coefficient of 5 ps −1 ) [29].Electrostatic interactions were calculated via the particle mesh Ewald (PME) scheme [30], while the distance cutoff for the van der Waals interactions was set at 12 Å.The percentage of the protonated doxorubicin molecules was calculated through the Henderson−Hasselbalch relation [31] taking into account the pK a of the drug's ionizable group at 37 • C and considering physiological pH conditions.An appropriate amount of chlorine and sodium ions was added for the overall neutralization of the system and for accomplishing an ionic strength of 165 mM [32].
More detailed information regarding the construction and the equilibration of the hyperbranched polymer/drug complex can be found in Reference [19].Following the construction of the system, energy minimization with steepest descent and conjugate gradient methods was applied, prior to 40 ns long MD equilibration runs in the canonical (NVT) ensemble at T = 310 K.This procedure ensured stabilization of energetic (total energy), thermodynamic (pressure), and conformational properties (i.e., average size of the PEGylated hyperbranched polymer) of the polymeric component of the system.The equilibrium area per lipid was approximately 75 Å 2 and the parallel orientation of the lipids was found to persist throughout the simulation.Production runs were also performed in the NVT ensemble and at the same (T = 310 K) temperature.

Results
The analysis performed in the simulation data aimed at investigating the spatial arrangement of the drug molecules with respect to the hyperbranched carrier and the negatively charged DPPG membrane, the degree of their association with the latter, and the driving forces responsible for the observed association.

Spatial Distribution of Drug Molecules
To examine the spatial arrangement of the drug molecules with respect to the hyperbranched nanoparticle and the bilipid membrane, we have calculated relevant average distances between centers of mass as well as relevant mass distributions.Figure 4 shows the average distance between the center of mass of the DPPG membrane and the center of mass of the other moieties.
If we take into account that the boundary of the HBP periphery was found to be at approximately 20 Å distance from its center of mass, it appears that drug molecules remain on average close to the HBP's surface.However, the protonated doxorubicin molecules together with the Na + counterions appear to be somewhat closer to the membrane.The closer proximity of the positively charged moieties should be attributed to the favorable Coulombic interactions with the membrane.Lack of such favorable electrostatic interactions between the DPPG membrane and the HBP, as well as the repulsive interactions with the Cl − counterions, can account for the larger distances between those moieties and the bilipid membrane.
To elaborate more on the spatial arrangement of the drug molecules with respect to the hyperbranched nanoparticle, we have calculated their mass distribution as a function of their distance from the center of mass of the HBP, as illustrated in Figure 5.If we take into account that the boundary of the HBP periphery was found to be at approximately 20 Å distance from its center of mass, it appears that drug molecules remain on average close to the HBP's surface.However, the protonated doxorubicin molecules together with the Na + counterions appear to be somewhat closer to the membrane.The closer proximity of the positively charged moieties should be attributed to the favorable Coulombic interactions with the membrane.Lack of such favorable electrostatic interactions between the DPPG membrane and the HBP, as well as the repulsive interactions with the Cl − counterions, can account for the larger distances between those moieties and the bilipid membrane.
To elaborate more on the spatial arrangement of the drug molecules with respect to the hyperbranched nanoparticle, we have calculated their mass distribution as a function of their distance from the center of mass of the HBP, as illustrated in Figure 5.The main feature characterizing the mass distributions of the drug molecules is their bimodal behavior.No such bimodal behavior for the protonated doxorubicin was noted in the absence of the membrane [19].The second peak of the drugs' distribution located beyond the HBP boundary exhibits a partial overlap with the membrane distribution.The degree of this overlap is considerably higher between the distribution describing the protonated drug molecules and that of the membrane.Quantification of the degree of penetration of the drug molecules within the HBP's structure and within the bilipid membrane can be made by calculating the area of the profiles within the boundaries of the HBP and the membrane respectively.Tables 2 and 3 list the so-calculated degrees of drug penetration.
Table 2. Average number and degree of penetration of drug molecules within the HBP's structure.The main feature characterizing the mass distributions of the drug molecules is their bimodal behavior.No such bimodal behavior for the protonated doxorubicin was noted in the absence of the membrane [19].The second peak of the drugs' distribution located beyond the HBP boundary exhibits a partial overlap with the membrane distribution.The degree of this overlap is considerably higher between the distribution describing the protonated drug molecules and that of the membrane.Quantification of the degree of penetration of the drug molecules within the HBP's structure and within the bilipid membrane can be made by calculating the area of the profiles within the boundaries of the HBP and the membrane respectively.Tables 2 and 3 list the so-calculated degrees of drug penetration.The percentages appearing in Table 3 indicate that a significant amount of drug molecules have migrated from the HBP/drug complex towards the DPPG surface.This percentage is considerably higher for the protonated drug molecules, apparently due to the attractive Coulombic interactions with the negatively charged lipids.It must be noted here, that the percentages appearing in Tables 2 and 3 reflect equilibrium conditions, since the timescale of the simulation is orders of magnitude longer compared to the timescale necessary for the drug molecules to reach the diffusive regime at conditions similar to the ones examined here [19], while the average distances between the different molecules remain stable (see Figure 4).

Charge Distributions
In our effort on one hand to discuss further the driving forces which are responsible for the mass distributions described before, and on the other hand to explore overcharging phenomena which are known to play a significant role in drug or gene-delivery formulations [33][34][35] and in the conformational properties of the hyperbranched molecules [36], we calculated component and overall charge distributions close to the HBP/drug complex and close to the membrane's surface.Figure 6 portrays the overall charge profile (arising from all components in the system) when moving radially outward with respect to the center of mass of the HBP (all charges henceforth are presented in units of an electron charge).It is shown that the overall charge profile exhibits a strong fluctuation around zero at distances extending from the center of mass of the PEGylated hyperbranched polyester to a region close to its periphery.At the outermost region, however, a net positive charge develops.As it was found, this charge arises from the positively charged drug molecules which are concentrated close to the HBP's outer surface.Calculation of the overall charge of the complex rendered a net charge of approximately 8.5e (integration from the polymer's center of mass to its boundary at ~22.5 Å).
Fluids 2019, 4 FOR PEER REVIEW 8 conditions similar to the ones examined here [19], while the average distances between the different molecules remain stable (see Figure 4).

Charge Distributions
In our effort on one hand to discuss further the driving forces which are responsible for the mass distributions described before, and on the other hand to explore overcharging phenomena which are known to play a significant role in drug or gene-delivery formulations [33][34][35] and in the conformational properties of the hyperbranched molecules [36], we calculated component and overall charge distributions close to the HBP/drug complex and close to the membrane's surface.Figure 6 portrays the overall charge profile (arising from all components in the system) when moving radially outward with respect to the center of mass of the HBP (all charges henceforth are presented in units of an electron charge).It is shown that the overall charge profile exhibits a strong fluctuation around zero at distances extending from the center of mass of the PEGylated hyperbranched polyester to a region close to its periphery.At the outermost region, however, a net positive charge develops.As it was found, this charge arises from the positively charged drug molecules which are concentrated close to the HBP's outer surface.Calculation of the overall charge of the complex rendered a net charge of approximately 8.5e (integration from the polymer's center of mass to its boundary at ~22.5 Å).To explore the charge profiles with respect to the bilipid's surface, we considered a reference frame with its origin in the center of mass of the bilayer and with the directions of its axes defined by the principal axes of inertia of the membrane.The profiles, shown in figure 7, were calculated along the direction normal to the plane defined by the directions of two of the principal axes of inertia  To explore the charge profiles with respect to the bilipid's surface, we considered a reference frame with its origin in the center of mass of the bilayer and with the directions of its axes defined by the principal axes of inertia of the membrane.The profiles, shown in Figure 7, were calculated along the direction normal to the plane defined by the directions of two of the principal axes of inertia which lie parallel to the membrane's surface (here denoted as z). .Symmetrized charge distributions of all the moieties in the examined system, as a function of the distance from the plane, which is defined from the center of mass of the lipid membrane (considered as an individual molecule) and the two principal inertia axes, which are parallel to the lipid layer."Q total" denotes the sum of all the component charge profiles.
Focusing on the overall sum of the distributions, it appears that there exists a net negative charge close to the membrane's surface.The membrane boundary is located close to 30 Å as this has been determined independently by the corresponding mass distribution (not shown here).As anticipated, the net charge arises by the negatively charged lipids.The concentration of positive charges arising from the Na + and the DOX+ drug molecules does not suffice to counterbalance the lipids' negative charge.Estimation of the overall charge of the membrane complex results to a value of approximately −100 (integration from the membrane's center of mass to its boundary in both directions), which is appreciably reduced (by about 22%) compared to the entire membrane charge (−128).

Effects of the Presence of the Membrane in the Clustering Behavior of Doxorubicin
As has been demonstrated in past studies (see Reference [19] and references therein) Doxorubicin tends to self-associate in aqueous solutions forming oligomeric clusters.This behavior is attributed to hydrophobic forces and to the π−π electron interaction between the planar aromatic parts of the drug molecules.The formation of such drug aggregates and the structural characteristics of these clusters may be of paramount importance regarding their pharmacological efficiency, since they can affect their transport properties, their permeation through a membrane, and ultimately the final drug biodistribution.To monitor the changes imparted in the associative behavior of doxorubicin molecules in the presence of the polymeric host and the bilipid membrane, the relative arrangement of the drug molecules was examined through radial distribution functions arising from their centers of mass.
Figure 8 shows the radial distribution functions of the center of mass of the drug molecules (neutral and protonated) with and without the presence of the lipid membrane (the behavior without the presence of the membrane was taken from Reference [19]).The distributions without the membrane are almost bimodal, characterized by a sharp peak close to an intermolecular distance of 5 Å and a second maximum close to 7-8 Å.The short-distance peak arises from the closest neighbors of a drug molecule and is consistent with the formation of dimers, while the peak corresponding to longer distances arises from the presence of neighboring clusters.It must be noted that in the case of the protonated drugs, the electrostatic repulsions due to the presence of the charged amine groups affect the spatial arrangement of the molecules (i.e., the relative intensity and the number of the Symmetrized charge distributions of all the moieties in the examined system, as a function of the distance from the plane, which is defined from the center of mass of the lipid membrane (considered as an individual molecule) and the two principal inertia axes, which are parallel to the lipid layer."Q total" denotes the sum of all the component charge profiles.
Focusing on the overall sum of the distributions, it appears that there exists a net negative charge close to the membrane's surface.The membrane boundary is located close to 30 Å as this has been determined independently by the corresponding mass distribution (not shown here).As anticipated, the net charge arises by the negatively charged lipids.The concentration of positive charges arising from the Na + and the DOX+ drug molecules does not suffice to counterbalance the lipids' negative charge.Estimation of the overall charge of the membrane complex results to a value of approximately −100 (integration from the membrane's center of mass to its boundary in both directions), which is appreciably reduced (by about 22%) compared to the entire membrane charge (−128).

Effects of the Presence of the Membrane in the Clustering Behavior of Doxorubicin
As has been demonstrated in past studies (see Reference [19] and references therein) Doxorubicin tends to self-associate in aqueous solutions forming oligomeric clusters.This behavior is attributed to hydrophobic forces and to the π−π electron interaction between the planar aromatic parts of the drug molecules.The formation of such drug aggregates and the structural characteristics of these clusters may be of paramount importance regarding their pharmacological efficiency, since they can affect their transport properties, their permeation through a membrane, and ultimately the final drug biodistribution.To monitor the changes imparted in the associative behavior of doxorubicin molecules in the presence of the polymeric host and the bilipid membrane, the relative arrangement of the drug molecules was examined through radial distribution functions arising from their centers of mass.
Figure 8 shows the radial distribution functions of the center of mass of the drug molecules (neutral and protonated) with and without the presence of the lipid membrane (the behavior without the presence of the membrane was taken from Reference [19]).The distributions without the membrane are almost bimodal, characterized by a sharp peak close to an intermolecular distance of 5 Å and a second maximum close to 7-8 Å.The short-distance peak arises from the closest neighbors of a drug molecule and is consistent with the formation of dimers, while the peak corresponding to longer distances arises from the presence of neighboring clusters.It must be noted that in the case of the protonated drugs, the electrostatic repulsions due to the presence of the charged amine groups affect the spatial arrangement of the molecules (i.e., the relative intensity and the number of the peaks).
enhancement of their tendency to cluster close to the lipid surface.This behavior could also be related to the effective higher concentration of the drug molecules in the system with the membrane.As far as the neutral drugs are concerned, a similar behavior can be observed.It is noticed that in the presence of the DPPG bilayer, the third peak arising from neighboring clusters exhibits a higher intensity compared to the analogous peak without the presence of the membrane.In other words, the limited diffusion of the drug clusters due to the geometric restriction imposed by the membrane, combined with hydrophobic forces, results in a higher degree of the drug molecules' localization.In general, therefore, for both cases (protonated and neutral drugs), the presence of the bilipid membrane appears to reinforce the drug clustering behavior.

Hydrogen Bonding
In the absence of the DPPG membrane, it was found that apart from the geometric entrapment of the drug molecules within the polymeric structure, another mechanism, that of hydrogen bonding, was also responsible for the association of doxorubicin with the polymeric host [19].This route for complexation can play a key role in drug delivery processes [37].In the presence of the bilipid In the case of the protonated drugs in the presence of the membrane, it is observed that on average the separation between the formed clusters decreases (see the distance between the first and the second peak with and without the presence of the membrane).This can be associated with the electrostatic attraction of the protonated drug molecules by the membrane, which results in an enhancement of their tendency to cluster close to the lipid surface.This behavior could also be related to the effective higher concentration of the drug molecules in the system with the membrane.As far as the neutral drugs are concerned, a similar behavior can be observed.It is noticed that in the presence of the DPPG bilayer, the third peak arising from neighboring clusters exhibits a higher intensity compared to the analogous peak without the presence of the membrane.In other words, the limited diffusion of the drug clusters due to the geometric restriction imposed by the membrane, combined with hydrophobic forces, results in a higher degree of the drug molecules' localization.In general, therefore, for both cases (protonated and neutral drugs), the presence of the bilipid membrane appears to reinforce the drug clustering behavior.

Hydrogen Bonding
In the absence of the DPPG membrane, it was found that apart from the geometric entrapment of the drug molecules within the polymeric structure, another mechanism, that of hydrogen bonding, was also responsible for the association of doxorubicin with the polymeric host [19].This route for complexation can play a key role in drug delivery processes [37].In the presence of the bilipid membrane, additional hydrogen-bonding-capable groups are introduced in the system, creating thus antagonistic to the HBP centers for the association of drug molecules.Recalling the image discussed in Section 3.1 (Figure 5 and Table 3), it is of interest to examine the role of hydrogen bonding in the membrane/drug association.To this end, we have examined pair correlation functions of hydrogen-bonding-capable atoms.The criteria we followed for the identification of a hydrogen bond were based on the hydrogen-acceptor distance in combination with the angle formed by the donor-hydrogen-acceptor triplet.The maximum distance considered, was the extent of the corresponding peak in the pertinent pair correlation function, while the minimum donor-hydrogen-acceptor angle was taken to be 120 • [38,39].
Such pair correlation functions between atoms belonging to the protonated doxorubicin molecules and to the lipids of the membranes are shown in Figure 9.In all cases, the existence of a sharp peak close to 2.5 Å implies the presence of hydrogen bonds between the examined atomic pairs.The intensity of the hydrogen-bonding peak is indicative to the abundance of each kind of hydrogen bond.Therefore, the pair between the amine hydrogen of the DOX+ and the phosphorus atom of the lipid, appears to be the more abundant.Analogous pair correlation functions between the neutral drug molecules and the lipid are shown in Figure 10.membrane, additional hydrogen-bonding-capable groups are introduced in the system, creating thus antagonistic to the HBP centers for the association of drug molecules.Recalling the image discussed in Section 3.1 (Figure 5 and Table 3), it is of interest to examine the role of hydrogen bonding in the membrane/drug association.To this end, we have examined pair correlation functions of hydrogenbonding-capable atoms.The criteria we followed for the identification of a hydrogen bond were based on the hydrogen-acceptor distance in combination with the angle formed by the donorhydrogen-acceptor triplet.The maximum distance considered, was the extent of the corresponding peak in the pertinent pair correlation function, while the minimum donor-hydrogen-acceptor angle was taken to be 120° [38,39].Such pair correlation functions between atoms belonging to the protonated doxorubicin molecules and to the lipids of the membranes are shown in Figure 9.In all cases, the existence of a sharp peak close to 2.5 Å implies the presence of hydrogen bonds between the examined atomic pairs.The intensity of the hydrogen-bonding peak is indicative to the abundance of each kind of hydrogen bond.Therefore, the pair between the amine hydrogen of the DOX+ and the phosphorus atom of the lipid, appears to be the more abundant.Analogous pair correlation functions between the neutral drug molecules and the lipid are shown in Figure 10.
To quantify the tendency for hydrogen-bond formation between the examined pairs, we have calculated the average number of each kind of hydrogen bond, per timestep, as presented in Table 4.    Inspection of the values in Table 4 verifies that hydrogen bonding between DOX+ amine hydrogens (HN) and lipid phosphorus atoms (P5) is the most frequent.A lower degree of hydrogen bonding which appears to be marginally statistically significant, is that of the DOX+ hydroxyl hydrogen (HO) and lipid hydroxyl oxygen (OH), while an even lower degree of hydrogen bonding is observed between the DOX+ amine hydrogen (HN) and the lipid's hydroxyl oxygen (OH).In conjunction with the results from Section 3.1 it can be concluded that one of the mechanisms which is responsible for the association between the protonated drug molecules and the membrane's surface is the formation of hydrogen bonds (another mechanism is related to the electrostatic interactions, as discussed in Section 3.2).As far as the neutral drugs are concerned, their association with the lipid membrane discussed in Section 3.1 does not appear to be driven by hydrogen-bond formation (the average number of hydrogen bonds is not statistically significant).A more probable driving force for the latter is therefore the hydrophobic interactions of the neutral drug molecules with water, which can result in a more favorable interaction with the membrane's surface.To quantify the tendency for hydrogen-bond formation between the examined pairs, we have calculated the average number of each kind of hydrogen bond, per timestep, as presented in Table 4. Inspection of the values in Table 4 verifies that hydrogen bonding between DOX+ amine hydrogens (HN) and lipid phosphorus atoms (P5) is the most frequent.A lower degree of hydrogen bonding which appears to be marginally statistically significant, is that of the DOX+ hydroxyl hydrogen (HO) and lipid hydroxyl oxygen (OH), while an even lower degree of hydrogen bonding is observed between the DOX+ amine hydrogen (HN) and the lipid's hydroxyl oxygen (OH).In conjunction with the results from Section 3.1 it can be concluded that one of the mechanisms which is responsible for the association between the protonated drug molecules and the membrane's surface is the formation of hydrogen bonds (another mechanism is related to the electrostatic interactions, as discussed in Section 3.2).As far as the neutral drugs are concerned, their association with the lipid membrane discussed in Section 3.1 does not appear to be driven by hydrogen-bond formation (the average number of hydrogen bonds is not statistically significant).A more probable driving force for the latter is therefore the hydrophobic interactions of the neutral drug molecules with water, which can result in a more favorable interaction with the membrane's surface.

Diffusional Behavior
The effects of the presence of the DPPG membrane in the clustering characteristics and in the associative behavior of the drug molecules are expected to be reflected in their transport properties as well.To obtain information regarding the diffusional motion of the different molecular components of the examined system, we monitored the mean squared displacement (MSD) of their center of mass as depicted in Figure 11.For comparison purposes we have also included the diffusional behavior of the drug molecules in the membrane-free case [19].

Diffusional Behavior
The effects of the presence of the DPPG membrane in the clustering characteristics and in the associative behavior of the drug molecules are expected to be reflected in their transport properties as well.To obtain information regarding the diffusional motion of the different molecular components of the examined system, we monitored the mean squared displacement (MSD) of their center of mass as depicted in Figure 11.For comparison purposes we have also included the diffusional behavior of the drug molecules in the membrane-free case [19].Focusing on the behavior of doxorubicin molecules, it appears that the presence of the DPPG membrane drastically affects the drug diffusional motion.Moreover, it appears to do so to a different degree depending on whether the drug molecules are neutral or protonated.The diffusional motion of the neutral drug at all the examined timescales is slower compared to the analogous behavior in the membrane-free system, remaining very close to that of the HBP.This observation can be rationalized on one hand by the fact that the majority of the neutral drug molecules are located within the HBP's structure (see Figure 5 and Table 2) and on the other hand by the more restrictive environment (i.e., higher effective concentration compared to the membrane-free case, enhanced clustering behavior, partial association with the bilipid membrane) experienced by these molecules in the presence of the bilipid layer.
In contrast, the protonated drug molecules diffuse much faster compared to their behavior in the absence of the membrane as far as it concerns the short timescales.This trend changes at larger timescales where the transport of the center of mass of the latter becomes slower than the membranefree case.The origin for the more complex behavior exhibited by the protonated drug should be related to the electrostatic interactions between the cationic doxorubicin molecules and the anionic part of the DPPG lipids.At short timescales, Coulombic interactions between the two oppositely charged moieties provide the driving force for an enhanced diffusive behavior.At larger timescales, where the effects of the more constrictive environment become more important, all those factors affecting the neutral drug molecules start affecting the clusters of the protonated drug molecules as well, resulting in the slowing down of their diffusional motion.Focusing on the behavior of doxorubicin molecules, it appears that the presence of the DPPG membrane drastically affects the drug diffusional motion.Moreover, it appears to do so to a different degree depending on whether the drug molecules are neutral or protonated.The diffusional motion of the neutral drug at all the examined timescales is slower compared to the analogous behavior in the membrane-free system, remaining very close to that of the HBP.This observation can be rationalized on one hand by the fact that the majority of the neutral drug molecules are located within the HBP's structure (see Figure 5 and Table 2) and on the other hand by the more restrictive environment (i.e., higher effective concentration compared to the membrane-free case, enhanced clustering behavior, partial association with the bilipid membrane) experienced by these molecules in the presence of the bilipid layer.
In contrast, the protonated drug molecules diffuse much faster compared to their behavior in the absence of the membrane as far as it concerns the short timescales.This trend changes at larger timescales where the transport of the center of mass of the latter becomes slower than the membrane-free case.The origin for the more complex behavior exhibited by the protonated drug should be related to the electrostatic interactions between the cationic doxorubicin molecules and the anionic part of the DPPG lipids.At short timescales, Coulombic interactions between the two oppositely charged moieties provide the driving force for an enhanced diffusive behavior.At larger timescales, where the effects of the more constrictive environment become more important, all those

Conclusions
As was demonstrated by the analysis presented, the presence of the DPPG membrane imparts significant changes regarding the spatial distribution and the clustering behavior of the drug molecules.Almost 1/3 of the protonated and 1/5 of the neutral doxorubicin molecules appear to migrate from the drug/HBP complex towards the membrane's surface.Localization of drug molecules close to the membrane's surface enhances their tendency to form clusters, particularly for the positively charged drug molecules.The latter form persistent hydrogen bonds mainly with the phosphorus atoms of the hydrophilic part of the lipids.The affinity of the neutral doxorubicin molecules to the lipid bilayer appears to be driven by their hydrophobic nature.The uneven distribution of the DOX+ molecules and that of the counterions' between the HBP and the bilipid membrane results in the development of positive and negative surface charge excesses, respectively.This may eventually lead to the coalescence of the HBP/drug complex with the membrane and thus towards a mechanism for a liposomal escape in formulations in which such HBP/doxorubicin complexes are encapsulated in a DPPG-based liposomal carrier.
The timescales associated with the migration of the neutral and the protonated doxorubicin molecules from the HBP/drug complex to the membrane's surface, and thus the build-up of the final charge distribution, relates to the respective diffusional behavior.It appears that the presence of the DPPG membrane drastically affects the translational motion of doxorubicin with respect to the membrane-free case.The extent of the observed changes depends on the charge carried by the drug molecules.The neutral drug diffuses slower at all the examined timescales while the protonated drug diffuses faster at short timescales and slower at longer timescales, compared to the behavior observed in the absence of the DPPG bilayer.
The elementary mechanisms characterizing the interactions of doxorubicin molecules with a hyperbranched polymeric host in the presence of an anionic bilipid membrane may serve as a basis for a better understanding of the microscopic picture in other liposomal-based formulations for drug-delivery purposes, where similar interactions are expected to be present.

Figure 2 .
Figure 2. The pre-equilibrated model of the DPPG membrane that was used in our study.

Figure 2 .Figure 3 .
Figure 2. The pre-equilibrated model of the DPPG membrane that was used in our study.Fluids 2019, 4 FOR PEER REVIEW 5

Figure 3 .
Figure 3. (a) The constructed system including water molecules (shown in red-white).(b) The same system without the water molecules.Blue beads represent the chlorine ions and orange beads the sodium ions.Each drug molecule is shown in different color and in rod representation, while the PEGylated hyperbranched molecule is shown in red color.In the lower part of the snapshot, the lipids of the membrane are shown in wire representation.

FluidsFigure 4 .
Figure 4. Average distance as a function of time between the center of mass of the lipid membrane (considered as an individual molecule including all the lipids) and those of the hyperbranched polymer, the counterions (Cl − , Na + ), and the neutral (DOX) and the protonated (DOX+) drug molecules.

Figure 4 .0Figure 5 .
Figure 4. Average distance as a function of time between the center of mass of the lipid membrane (considered as an individual molecule including all the lipids) and those of the hyperbranched polymer, the counterions (Cl − , Na + ), and the neutral (DOX) and the protonated (DOX+) drug molecules.Fluids 2019, 4 FOR PEER REVIEW 7

Figure 5 .
Figure 5. Weight distribution of the drug molecules and of the other moieties with respect to the center of mass of the hyperbranched polymer (a spherical symmetry of the hyperbranched polyester (HBP) molecule has been considered).The notation in the legend is as in Figure 4, with LIP denoting the DPPG membrane.The vertical dashed line denotes the location of the radius of gyration of the polymer.

Figure 6 .
Figure 6.Overall charge distribution with respect to the HBP center of mass arising from all the components in the system.The vertical dashed line denotes the surface boundary of the hyperbranched nanoparticle.

Figure 6 .
Figure 6.Overall charge distribution with respect to the HBP center of mass arising from all the components in the system.The vertical dashed line denotes the surface boundary of the hyperbranched nanoparticle.

FluidsFigure 7
Figure7.Symmetrized charge distributions of all the moieties in the examined system, as a function of the distance from the plane, which is defined from the center of mass of the lipid membrane (considered as an individual molecule) and the two principal inertia axes, which are parallel to the lipid layer."Q total" denotes the sum of all the component charge profiles.

Figure 7 .
Figure7.Symmetrized charge distributions of all the moieties in the examined system, as a function of the distance from the plane, which is defined from the center of mass of the lipid membrane (considered as an individual molecule) and the two principal inertia axes, which are parallel to the lipid layer."Q total" denotes the sum of all the component charge profiles.

Figure 8 .
Figure 8. Radial distribution functions of the center of mass of the doxorubicin molecules in the presence (solid lines) and in the absence (dashed lines) of the bilipid membrane, for the (a) neutral and the (b) protonated drug molecules.

Figure 8 .
Figure 8. Radial distribution functions of the center of mass of the doxorubicin molecules in the presence (solid lines) and in the absence (dashed lines) of the bilipid membrane, for the (a) neutral and the (b) protonated drug molecules.

Figure 9 .
Figure 9. Pair correlation functions among atoms belonging to the protonated drug molecules (DOX+) and the DPPG lipids (LIP).(a) Pairs between amine hydrogens of DOX+ (HN) and hydroxyl oxygens of the lipid (OH).(b) Pairs between hydroxyl hydrogens (HO) of DOX+ and hydroxyl oxygens (OH) of the lipid.(c) Pairs between amine hydrogens of DOX+ (HN) and phosphorus (P5) of the lipid.

Figure 9 .
Figure 9. Pair correlation functions among atoms belonging to the protonated drug molecules (DOX+) and the DPPG lipids (LIP).(a) Pairs between amine hydrogens of DOX+ (HN) and hydroxyl oxygens of the lipid (OH).(b) Pairs between hydroxyl hydrogens (HO) of DOX+ and hydroxyl oxygens (OH) of the lipid.(c) Pairs between amine hydrogens of DOX+ (HN) and phosphorus (P5) of the lipid.

Figure 10 .
Figure 10.Pair correlation functions among atoms belonging to the neutral drug molecules (DOX) and the DPPG lipids (LIP).(a) Pairs between amine hydrogens of DOX (HN) and hydroxyl oxygens of the lipid (OH).(b) Pairs between hydroxyl hydrogens (HO) of DOX and hydroxyl oxygens (OH) of the lipid.(c) Pairs between amine hydrogens of DOX (HN) and phosphorus (P5) of the lipid.

Figure 10 .
Figure 10.Pair correlation functions among atoms belonging to the neutral drug molecules (DOX) and the DPPG lipids (LIP).(a) Pairs between amine hydrogens of DOX (HN) and hydroxyl oxygens of the lipid (OH).(b) Pairs between hydroxyl hydrogens (HO) of DOX and hydroxyl oxygens (OH) of the lipid.(c) Pairs between amine hydrogens of DOX (HN) and phosphorus (P5) of the lipid.

Figure 11 .
Figure11.Mean square displacement arising from the center of mass of the drug molecules, the hyperbranched polyester, and the bilipid membrane (considered as a single moiety including all the lipids).The behavior of the drug molecules in the membrane-free case is also shown for comparison purposes.

Figure 11 .
Figure11.Mean square displacement arising from the center of mass of the drug molecules, the hyperbranched polyester, and the bilipid membrane (considered as a single moiety including all the lipids).The behavior of the drug molecules in the membrane-free case is also shown for comparison purposes.

Fluids
neutral drug molecules start affecting the clusters of the protonated drug molecules as well, resulting in the slowing down of their diffusional motion.

Table 1 .
Number of molecules/ions of each type included in the simulated system.

Table 1 .
Number of molecules/ions of each type included in the simulated system.

Table 2 .
Average number and degree of penetration of drug molecules within the HBP's structure.

Table 3 .
Average number and degree of penetration of drug molecules within the membrane's structure.

Table 4 .
Average number of hydrogen bonds per timestep for the examined pairs.

Table 4 .
Average number of hydrogen bonds per timestep for the examined pairs.