Chitosan-covered liposomes as a promising drug transporter: nanoscale investigations

Abstract

Liposomes are small artificial vesicles spherical shaped of 50–1000 nm in diameter. They are created from natural non-toxic phospholipids membranes. Externally, they are decorated with biocompatible polymers. Chitosan, a natural polymer, demonstrates exceptional advantages in drug delivery, in particular, as liposome cover. In this paper, Molecular Dynamics simulations (MD) are performed in the coupled NPT-NPH and NVT-NVE statistical ensembles to study the static and dynamic properties of DPPC membrane-bilayer with grafted cationic chitosan chains, with added Cl⁻ anions to neutralize the environment, using the Martini coarse-grained force-field. From the NPT-NPH MD simulations we found a chitosan layer L_DM ranging from 3.2 to 6.6 nm for graft chains of a degree of polymerization n_p = 45 and different grafting molar fractions X_p = 0.005, X_p = 0.014 and X_p = 0.1. Also, the chitosan chains showed three essential grafting regimes: mushroom, critic, and brush depending on X_p. The DPPC bilayer thickness D_B and the area per lipid A_l increased proportionally to X_p. From the NVT-NVE MD simulations, the analysis of the radial distribution function showed that the increase of X_p gives a more close-packed and rigid liposome. The analysis of the mean square displacement revealed that the diffusion of lipids is anomalous. In contrast, the diffusion of chitosan chains showed a normal diffusion, just after 100 ps. The diffusion regime of ions is found to be normal and independent of time. For the three identified regimes, the chitosan showed a tendency to adhere to the membrane surface and therefore affect the properties of the liposomal membrane.

---

1 Introduction

Recently, liposomes have aroused scientific interest in many disciplines encompassing theoretical physics, biophysical chemistry, and colloidal science. They demonstrate an excellent ability to increase the localization of drugs in diseased cells.¹ Thus, liposomes represent promising systems for drug delivery to target the drug to the site of action or improve the circulation time of the drug in vivo, aiming at side effect reduction as well as the optimization of treatments.^2,3 Also, liposomes have been known for their huge potential as drug carriers for a variety of compounds like low molecular weight drugs, therapeutic proteins and peptides, vaccines and diagnostic agents, toxins, enzymes, antigens and antibodies, and nucleotides.^4–13 The liposomes have been introduced in the market as efficient spherical nanovesicles,^6–14 they are formed by a self-closed phospholipidic bilayer carrying the drug in its aqueous core suspended in water. During the formation of liposomes, the active drugs are merged either in the aqueous solution or in the phospholipids themselves. Thus, liposomes allowing the delivery of hydrophobic and hydrophilic agents. In vivo pharmacokinetics and bio-distribution can be controlled by the size and lipid composition of the liposomes. When using liposomes as drug carriers, they can release the contented drugs into the cytoplasm of cells that have internalized them by endocytosis. Therefore, they promote the intracellular penetration of active products. The liposomes surfaces are made by a lipid bilayer. Thus, the lipids nature can increase the contact time at the epithelial surface. For instance, transcorneal drug transport is limited by the intrinsic permeation characteristics of the corneal epithelium.¹⁵

Externally, liposomes are covered with polymers, which making them stable in the intravenous environment. Grafted polymers prolong the blood circulation time of the liposomes.^16–18 Thus, the accumulation of drugs is significantly improved within diseased sites. From the literature, liposomes coated with chitosan have been used as a mucoadhesive delivery system because their positively charged surface promotes adhesion to negatively charged cell membranes.^19–21 Besides, the adhesive ability plays a vital role in drug delivery. Mainly for prolonging retention in the gastrointestinal tract and promoting penetration into the mucus layer.²²

Nowadays, chitosan, a natural polymer that receives increasing attention as a new material,²³ may represent an alternative to neutral polymers, such as polyethylene glycol PEG and polyethylene oxide PEO,^6,7 for liposomes covering. Chitosan does not only possess advantageous non-toxicity, low allergenicity, biocompatibility and biodegradability^23–26 but also contains free charges in low pH water solution allowing interactions with the polar heads of the phospholipidic bilayers in the liposomes.²⁷ Regarding the applicability of liposomes containing chitosan, several studies have highlighted the use of chitosan as an absorption-enhancing agent and have mentioned his inherent properties of mucoadhesion and penetration.^28–30 Moreover, because of its bioadhesive properties, chitosan has also been worth in novel bioadhesive drug delivery systems for the treatment of wounds and burns as well as for local administration of drugs in treatments of diseases in the eyes, nasal, mucosa, or tumors.^31–34 Besides, chitosan is reported to have other biological and biomedical properties, such as antitumor,³⁵ antimicrobial,³⁶ and antioxidant activities.³⁷ Chitosan has been used in drug delivery systems in different forms, like tablets, microspheres, micelles, vaccines, nucleic acids, hydrogels, nanoparticles, and conjugates. Chitosan can be used in drug delivery systems in both implantable as well as injectable forms through oral, nasal and ocular routes. Besides, it facilitates transmucosal absorption which is important in nasal and oral delivery of some polar drugs like peptides along with protein vaccines for their administration.³⁸ It is commonly used as an excipient in tablet formulation for oral medication. Chitosan microspheres have been extensively investigated for the controlled release of drugs and vaccines through oral and nasal delivery. They were prepared by complexation between the cationic chitosan and anionic compounds. For example, tripolyphosphate or alginates.^39,40 The chitosan hydrogels are three-dimensionally structured hydrophilic polyelectrolyte which can absorb and hold up to thousands of times more fluids than their dry weights and use in drug delivery.⁴¹ Thus, the drugs can be loaded in the hydrogels by diffusion, entrapment, and tethering. After injection of the loaded hydrogel into the body, the drug diffuse into the other tissues.⁴² The recent development of in situ forming depots using chitosan-based hydrogel has attracted much attention as a new method for controlled drug release.⁴³

Chitosan-coated liposomes can be investigated using different tools. In this context, Ketzasmin et al. reported a Dissipative Particle Dynamics (DPD) mesoscopic simulations of nanoliposomes. The nanoliposomes are constituted by lecithin and coated with a shell of chitosan.⁴⁴ The DPD simulation showed that chitosan is deposited well on the surface of the nanoliposome, as has been reported in some experimental works.^45,46 Mengoni et al. performed a dynamic light scattering with non-invasive backscattering (DLS-NIBS) experiments on a chitosan-coated liposomes system. Their results confirm the presence of the chitosan coating, which is remarked in the increase of the average size of liposomes and in the inversion of the zeta potential from negative value, for the uncoated liposomes, to a positive value for the coated ones. In addition, their results show that chitosan-coating improves the stability of liposomes. Indeed, in contrast to uncoated microemulsions,⁴⁷ uncoated liposomes are thermodynamically unstable. Thus, they are rarely used without coating in drug delivery since they aggregate easily during storage, degradation, and delivery of drugs upon administration.⁴² In another work, Mady and Darwish investigated the chitosan-coated DPPC liposomes physicochemical properties, using drug release rate, transmission electron microscopy (TEM), zeta potential and turbidity measurement. The electron microscopy and the zeta potential confirmed the coating of DPPC liposomes with a chitosan layer. In addition, the authors concluded that chitosan increases liposome stability during drug release.⁴⁸ M. Hasan et al. applied various techniques to study the interactions between phospholipid, chitosan, and curcumin showing that liposome coating has interesting advantages.⁴⁹ After liposome coating, a flow behavior change from Newtonian to shear-thinning is demonstrated. The hysteresis loop area decreased significantly and therefore liposomal dispersion stability increased significantly with the addition of the chitosan layer. Also, improvement of mechanical stability is demonstrated. In addition, the liposome membrane structure was not affected by the chitosan layer.

The chitosan has interesting applications in drug delivery. Thus, it is necessary to better understand the microscopic properties of this promising polyelectrolyte and their dependence on the grafting molar fraction. As schematically presented in the Fig. 1, chitosan can be found in different conformation: mushroom, critic or brush as observed for neutral polymers.⁵⁰ Thus, Molecular Dynamics simulations (MD) have a great advantage in the analysis of polymer conformation. Simultaneously, molecular dynamics is a powerful tool for studying the interaction between particles under different microscopic conditions. The liposome diameter range is 50–1000 nm. To mimic two adjacent liposomes the simulated system consists of two opposed lipid bilayers, where the opposed monolayers containing incorporated DPPC-chitosan lipo-polyelectrolyte. The confined area between the bilayers is filled with water molecules and Cl⁻ counterions. This system is supposed to mimic portions of two adjacent liposomes. A schematic representation is presented in Fig. 2.

Fig. 1 Schematic representation of liposomes covered with polymers: the polymers form a shield that protects the liposomes and play a lubrification role in the intravenous environment. The polymers can show different conformation depending on the grafting molar fraction.

Fig. 2 Coarse-grained representation of the simulated system: (a) DPPC molecule, (b) chitosan chain and (c) the initial distribution, two bilayer membranes opposite to each other with grafted chitosan chains. This initial configuration was generated using the Moltemplate package (https://moltemplate.org/), and visualized using the Ovito vizualisation tool (https://ovito.org).

In this paper, we present and analyze the results obtained from NPT-NPH and NVT-NVE molecular dynamics simulations. The structure and dynamics of the system components are found depending on the grafting molar fraction. The chitosan layer increases and the chitosan dynamics decrease proportionally to the grafting molar fraction X_p. As interesting results, the DPPC area per lipid and the DPPC thickness increase with increasing X_p. In contrast, the DPPC dynamics showed a slight decrease with increasing X_p. Thus, chitosan incorporation can affect slightly liposome properties. Chitosan chains and ions are adhered to the liposome surface. All these features indicate that chitosan can be used as a promising covering for liposomes designed for drug delivery.

2 Simulation method

2.1 Interaction potentials

We used the Martini force field that reproduces the free partitioning energies between the polar and nonpolar phases of a large number of chemical compounds, to describe the interaction potentials of DPPC lipids, water, anions and chitosan.^51–54 In the Martini coarse-grained model, approximately four heavy atoms have been grouped into a single coarse-grain bead, whose chitosan repetitive unit is composed of three sub-monomers (R1–R2–R3). Such units are bonded repetitively to construct chitosan chains with the desired polymerization degree. As described in detail in the original work,⁵⁴ the R1 bead comprises atoms of the hydroxyl side chain, the bead R2 comprises most of the atoms of the polymer backbone and the bead R3 maps the hydroxyl and amino groups opposite to the side chain R1. The three beads have different hydrophilicity, and the R3 bead can become protonated and is then well presented as R3⁺. The bonded interactions are parameterized by fitting the distributions of corresponding bonds, angles and dihedrals to the all-atom MD simulation results. Also, the coarse-grained bond and angle distributions match well with corresponding distributions in all-atom MD simulation. The non-bonded interactions are parameterized using the MARTINI CG force field of carbohydrates.

In the coarse-grained model, all particle pairs i and j at distance r interact via Lennard-Jones V_LJ and coulombic V_C interactions:

with ε_ij the potential depth, and σ_ij the bead radius.with, , where ε₀ the vacuum permittivity, and ε = 15 the medium dielectric-constant. For the DPPC lipids and chitosan chains, the bond interactions between adjacent pairs of beads are described by a harmonic potential V_bond(r), and a cosine squared harmonic potential V_angle1(θ) is used for angle interactions between triplets of DPPC beads:

V_bond = K_b(r − r₀)² (3)

V_angle1 = K_a(cos(θ) − cos(θ₀))² (4)

For the chitosan polyelectrolyte, the angle interaction is computed using a harmonic potential V_angle2. Taking account the chitosan chains flexibility, a Fourier dihedral interaction potential V_dihedral is used to compute the interaction between quadruplets of bonded beads.

V_angle2 = K_a2(θ − θ₀)² (5)

The parameters for bonded and nonbonded interactions were taken from the recently available study for DPPC, anions, and water from ref. 51–53, and for the chitosan polyelectrolyte from ref. 54. The nonbonded interaction between the chitosan polyelectrolyte monomers, water beads, Cl⁻ anions, and DPPC beads are calculated using the Lorentz–Berthelot mixing rules,^55,56 where and .

2.2 Simulation details

The simulated system consists of two membranes bilayer, with a grafted chitosan model, where each monolayer consists of 4624 DPPC lipid molecules with an area per lipid A_l = 0.49 nm². Different grafting molar fractions are tested (0.005, 0.014, and 0.1), with a degree of polymerization n_p = 45. These parameters are used in previous experimental and simulation works for the neutral polymer polyethylene glycol (PEG).⁵⁰ Thus, if chitosan has the same physical properties as PEG in term of the length of the polymer layer and the effect on the liposome membrane, it will be preferable that it be used as an alternative to PEG to cover liposomes for drug delivery applications because of its pharmaceutical advantages.

The confined area between the membranes bilayer is hydrated by 100 000 explicit water beads. The neutral charge was realized by adding Cl⁻ anions to compensate the positive charge of chitosan sub-monomer R3⁺. Molecular dynamics simulations were performed using the LAMMPS software package.⁵⁷ Initially, the two DPPC membrane-bilayers with incorporated DPPC-chitosan lipopolyelectrolyte, opposed to each other, are generated in a cubic simulation box of lengths (L_x = L_y = L_z = 50 nm). Then the confined volume between the two membrane-bilayers is filt with water and Cl⁻ ions. This preformed initial configuration was optimized using an energy minimization algorithm, by iteratively adjusting particle coordinates, for a sufficiently long time (10⁶ time-steps), a finite unitless stopping tolerance for energy (E_stop = 10⁻⁷) that represents the energy change between successive iterations divided by the energy magnitude, and a stopping tolerance for force (F_stop = 10⁻⁸ kcal Å⁻¹. The iterations are stopped when the energy change between outer iterations is less than E_stop, or the length of the global force vector is less than the F_stop. The stopping time, at which the equilibrium is achieved, is found of the order of (10⁵) time-steps for all simulated cases. Thus, the pre-formed initial distribution of the particles in the simulation box is corrected to be in local potential energy minimum. By using this minimization algorithm, the divergence is avoided and the thermodynamics equilibrium, in the production ensembles NPT-NPH and NVT-NVE, is reached quickly.

Then, molecular dynamics simulations are performed in NPT-NPH and NVT-NVE ensembles. Where N is the number of coarse-grained beads, V is the volume of the simulated system, T is the temperature, H is the enthalpy, E is the internal energy and P is the external pressure. In the NPT-NPH simulation, which updates the position and velocity of the system particles and allows the change in volume using time integration on Nosé–Hoover style non-Hamiltonian equations of motion of the system particles of molar mass m_i in Dalton, without energy divergence. We calculated the equilibrium distance between the two membranes, which correspond to the length of the chitosan brush bilayer, and the thicknesses of the membranes, for different molar fractions of lipopolyelectrolyte at ambient conditions of temperature and pressure (T = 300 K, P = 1 bar), by analyzing the density profiles. For all simulations, periodic boundary conditions are used.⁵⁸ The proposed initial distribution of the simulated system (Fig. 2) was equilibrated under the NPT-NPH conditions, by coupling the NPH barostat to the Langevin thermostat,^57,59 with the equilibrium temperature, which is proportional the particles velocities, T = 300 K and the damping parameter (damp = 100). This simulation updates the positions and the velocities of the particles for each time step. Also, it modifies the volume to reach the imposed equilibrium conditions on the external pressure and the temperature. From 10⁵ to 10⁶ time-step (δt = 10 fs), the area between the opposed membranes bilayer remains practically constant and then the layer of the polymer brush, along the z-direction, is calculated from the density profiles.

The second simulation is a Brownian dynamics (BD) performed on the NVT-NVE ensemble by coupling the Langevin thermostat which models an interaction with an implicit background solvent with the equilibrium temperature, which is proportional the particles velocities, T = 300 K and the damping parameter (damp = 100), to the NVE integration scheme which performs a constant NVE simulation to update the position and the velocity of coarse-grained beads.⁶⁰ In this simulation, we investigate the micro-distribution (through the radial distribution function g(r)) and dynamics (through the mean square displacement and the diffusion coefficients) of DPPC lipids, chitosan molecules, and Cl⁻ anions. The simulations were run for 10⁶ time-step, with a time step δt = 10 fs.

The visualizations of the NPT-NPH simulations, the snapshots shown in Fig. 3–5 correspond to the final step are also the most probable conformations, as we notice a thermodynamics equilibrium for all simulated cases after 10⁵ time-step, i.e., one nanosecond.

Fig. 3 Snapshot of the NPT-NPH simulation at T = 300 K and P = 1 atm of the polyelectrolyte conformation in the mushroom regime using the Ovito visualization tool.⁶¹ The model membranes composed of DPPC lipid bilayer are decorated with a molar fraction X_p = 0.005 of DPPC-chitosan lipopolyelectrolyte with the same degree of polymerization n_p = 45. The whole is hydrated by water molecules. After one nanosecond (10⁵) time-step δt = 10 fs, the simulated system is equilibrated and the membrane–membrane distance remains constant. The chains do not interact and behave almost like isolated grafted chains.

Fig. 4 Snapshot of the NPT-NPH simulation at T = 300 K and P = 1 atm of the polyelectrolyte conformation in the critic regime. The model membranes composed of DPPC lipid bilayer are decorated with a molar fraction X_p = 0.014 of DPPC-chitosan lipopolyelectrolyte with the same degree of polymerization n_p = 45 and the whole is hydrated by water molecules. The transition between mushroom and brush regimes take place when the chitosan chains begin to overlap.

Fig. 5 Snapshot of the NPT-NPH simulation at T = 300 K and P = 1 atm of polyelectrolyte conformation in the brush regime. The model membranes composed of DPPC lipid bilayer are decorated with a molar fraction X_p = 0.1 of DPPC-chitosan lipopolyelectrolyte with the same degree of polymerization n_p = 45 and the whole is hydrated by water molecules. The chitosan chains are extended from the grafting surface forming a brush-like conformation.

More pieces of information on the energy minimization algorithm and the different statistical ensembles, used in this work, are available in detail on the online version of the LAMMPS package (http://lammps.sandia.gov). All interaction potentials parameters, needed for the MD simulation, are listed in detail in Tables 1–4. More detailed information about the MARTINI coarse-grained representation and parametrization of the coarse-grained interaction potentials of lipids, polymers (monomers), water and ions, can be found on the MARTINI home page (http://cgmartini.nl/).

Table 1 Non-bonded interaction potential parameters

Pair type	σ (Å)	ε (kcal mol^−1)
Q0–Q0, C1–C1	4.7	0.836521
Qa–Qa, MW–MW	4.7	1.195030
Na–Na, Q0–Na, Qa–Na	4.7	0.956024
Q0–Qa	4.7	1.075527
Q0–C1, Qa–C1	6.2	0.478012
MW–C1	4.7	0.478012
Na–C1	4.7	0.645316
Q0–MW, Qa–MW	4.7	1.338434
Na–MW	4.7	0.956024
R1–R1, R2–R2, R1–R2	4.7	1.195030
R1–R3^+, R2–R3^+, R3^+–R3^+	4.7	1.338434
R1–MW, R2–MW, R1–Qa, R2–Qa	4.7	1.195030
R1–C1, R2–C1, R1–Q0, R2–Q0	4.7	0.999834
R1–Na, R2–Na	4.7	1.068867
R3^+–MW, R3^+–Qa	4.7	1.264701
R3^+–C1, R3^+–Q0	4.7	1.058124
R3^+–Na	4.7	1.131183
Cl^−–Cl^−	4.83	0.012785
Cl^−–C1, Cl^−–Q0	4.765	0.103416
Cl^−–MW, Cl^−–Qa	4.765	0.123606
Cl^−–Na	4.765	0.110556
Cl^−–R2, Cl^−–R2	4.765	0.123606
Cl^−–R3^+	4.765	0.130812

---

Table 2 Bonded interaction potential parameters

Bond type	r0 (Å)	Kb (kcal mol^−1)
Q0–Qa	4.7	2.98
Qa–Na	4.7	2.98
Na–Na	3.7	2.98
Na–C1	4.7	2.98
C1–C1	4.7	2.98
R2–R2	5.212	7.17
R2–R3	1.934	7.17
R1–R1	2.494	7.17

---

Table 3 Angle interaction potential parameters

Angle type	θ0 (°)	Ka1 (kcal mol^−1)
Qa–Na–Na	120	5.97
Na–C1–C1	180	5.97
C1–C1–C1	180	5.97

---

Angle type	θ0 (°)	Ka2 (kcal mol^−1)
R10–R20–R30	123.76	107.553
R20–R21–R22	163.78	83.652
R20–R21–R11	70.94	107.553
R10–R20–R21	90.88	23.9006
R30–R20–R21	135.94	71.7017
R31–R21–R20	58.01	191.205

---

Table 4 Dihedral interaction potential parameters

Dihedral type	ϕ0 (°)	Ki (kcal mol^−1)	ni
R10–BR0–R21–R11	175.62	1.91205	3
R10–R20–R21–R31	−28.59	2.39006	3
R30–R20–R21–R11	36.57	4.3021	3

---

3 Results and discussions

3.1 Length of the chitosan layer

Chitosan chains have n_p Kuhn segments and carry n_p elementary charges e, where the so-called Kuhn segment is composed of three attached beads (R1–R2–R3⁺). For each chain, the same number n_p of mobile monovalent Cl⁻ counterions distributed in the solution above the membranes. In a recent report,⁵⁰ we studied the grafted PEG conformations using MD simulations. We showed that, when the lipo-polymers are incorporated into a bilayer membrane, the PEG neutral polymers exhibit three conformations, namely mushroom, critic and brush regimes, for a fixed degree of polymerization n_p = 45. These conformations are found depending on the lipo-polymer grafting molar fraction (X_p). Thus, the chitosan chains are expected to show the same conformations.

3.1.1 Length of the chitosan layer for low grafting molar-fraction. At high salt concentration and low grafting molar fraction, the conformation of chitosan chains is independent of the grafting molar fraction and salt concentration, and each chain is isolated from its neighbors. Therefore, there are no collective stretching effects. For simplicity, we consider the chitosan chains as linear homopolymer chains in a homogeneous environment and disregard all monomer–monomer interactions, including the hard-core ones that prevent the chains from overlapping and the coulombic interaction between charged monomers. For high molecular weights, all such polymers fall into the same universality class of non-avoiding random walks. Thus, they exhibit a degree of common behavior. For instance, their length scales with their molecular weight with a universal exponent, ν = 1/2. This behavior is common for natural and artificial polymers. For maximum simplicity, we select the freely-jointed chain, where each monomer has a fixed length a, and the joint between sequential monomers is completely flexible. In this regime, the polymer chains have a random distribution.where r_i the modulus of the vector r_i. The factor in the denominator ensures that the distribution is properly normalized such that:

The size of the particular conformation can be characterized by the length of its end-to-end vector:

by averaging the length of this vector over all possible conformations, R₀, can then be used as a measure for the length of the polymer. Since the simple average, 〈|R|〉, is mathematically difficult to evaluate, then it is common practice to use the alternative root-mean-square average:

3.1.2 Length of the chitosan layer for the critical grafting molar-fraction. The critic regime is an intermediate regime between the mushroom and brush regimes. In this regime, the chitosan chains are not totally stretched but the inter-chains interaction is significant. Thus, the chains begin to overlap. For chitosan chains with a fixed polymerization degree n_p, the critical regime depends only on the grafting molar fraction X_p. Thus, the transition mushroom-brush occurs at such as:where A_l is the area of the membrane by lipids, as we will calculate in Section. 3.4, is found of the order of (A_l = 0.62 nm²) for the critical grafting molar fraction . For a graft chitosan with (n_p = 45), the transition between the mushroom and brush conformations must take place at for fluid phase membranes. The effect of the chitosan coating on the area per lipid will be discussed in detail in Section 3.4.

3.1.3 Length of the chitosan layer for high grafting molar-fraction. The layer of the chitosan brush in contact with a salt solution depends on the Bjerrum length l_B defined as the distance at which two elementary charges e interact with thermal energy k_BT = e²/(εl_B) in a solvent with dielectric constant ε, and, of course as well as for neutral brush, to the molar fraction of graft polyelectrolyte X_p, as presented firstly by Milner^62,63 and proved experimentally by Kenworthy et al.⁶⁴ We note that in this work we simulate a neutral system (n(R3⁺) = n(Cl⁻)), so the dependence of the polyelectrolyte layer on the salt concentration per chain is not discussed in this paper. This dependence is discussed in detail by Zhulina and Rubinstein.⁶⁵ For very high molar fraction, 1 chains create an electric field which stretches each chain beyond the size of a free chitosan R₀ (eqn (10)). The size of the chain in this regime, called the Pincus regime,⁶⁵ is determined by the balance of the electrostatic energy per chain F_elecwhere f is the degree of chain ionization (fraction of charged Kuhn segments), and F_elast is the elastic free energy per chainwith ν = 3/5 in good solvent conditions, giving a brush height

While the chitosan is highly charged (f = 1) with condensed Cl⁻ counterions, there is no Pincus regime and X^C_p > X^osm_p, with X^osm_p the osmotic molar fraction. For osmotic brushes, the grafting molar fraction is large enough that the counterions become highly localized, thus determining the properties of the brush. Therefore, in this case, the chitosan layer is determined by the balance of the osmotic pressure of the counterions and the elastic energy of the chains, which gives the following formula.⁶⁵

R ∼ an_p (15)

On the other hand, the energy of short-range excluded volume interactions between monomers becomes higher than the counterion osmotic contribution to the brush free energy. In this case, Ekaterina and Rubinstein demonstrate that PE brushes are in a quasi-neutral brush regime and have the same properties as neutral brushes.⁶⁵ Then the chitosan layer can be determined through the mean field theory. In the brush regime, the chitosan chains are extended and the problem is essentially unidimensional. The chains are confined in the normal direction to the grafting surface. In the framework of the mean-field theory, for a grafted polymer chains, the one-dimensional version of the average field expression for free energy is:⁶⁶

Here, for the polyelectrolyte brush, ν = a³ the excluded volume per repetitive unit (R1–R2–R3), X_p is the molar fraction of the graft polyelectrolyte and A_l is the area per lipid. Minimizing the free energy of eqn (16) with respect to R (the end-to-end distance of the chain), i.e., gives the following expression for the equilibrium length of the polymer brush:

3.2 Chitosan monomers and Cl⁻ counterions density distribution

Fig. 6a and b depict the density profile of grafted chitosan chains and Cl⁻ counterions, respectively, using a degree of polymerization n_p = 45, for three different molar fraction. It has been observed that the distance between two parts of adjacent liposomes, which are supposed as flat surfaces, increases with the increase of the grafting molar fraction. For bilayer of graft chitosans, in the ambient conditions of temperature and pressure, the density profile of the chitosan monomers present as symmetric curves centered in the middle of the z-direction of the simulation box (L_z = 0), giving that there is no remarkable interpenetration between the two polyelectrolytes layers. Thus, the chitosan bilayer is considered as polyelectrolyte brush in weak interpretation regime.⁶⁶ Then, the length of the chitosan layer can be estimated as L_DM = d/2. According to Fig. 6a, the monomer density decreases away from the grafting surface. Noted that, larger X_p implies a larger number of monomers. Therefore, an increase in X_p implies that we find a larger value of the monomer concentration at a given distance from the wall, regardless of the extent of lowering of the monomer concentration away from the grafting surface. Fig. 6a demonstrating the existence of finite monomer concentration near the end of the chitosan chains. Thus, chitosan adheres to biological membranes. Fig. 6b shows the counterion concentration distribution for different X_p. The counterion concentration distribution almost obeys the monomer concentration distribution for the mushroom and critic regimes. In contrast, for the brush regime, the counterions are pushed away from the grafting surface by sterical effects. Therefore, the counterion concentration decreases near the end of the chains as it is observed for the monomer concentration distribution. For the brush regime, we witness a more uniform counterion density distribution, with a small increase at the end of the chitosan chains, such an increase is observed also in the monomer's density. Thus, the counterions remain tightly bound to the brushes. Obviously, there is a quantitative variation between the monomer and counterion density profile, particularly near the end of chains, at that location the monomer density decreases steeply while the corresponding counterion density shows a much flatter decrease.

Fig. 6 The density profile in the z-direction of the confined chitosan chains (a) and Cl⁻ counterions (b).

In Table 5 we present the chitosan layers obtained by the present MD simulations at the NPT-NPH statistical ensemble under the ambient condition of pressure and temperature, for different molar fractions of chitosan (0.005, 0.014 and 0.1) using a lipopolyelectrolyte composed of a DPPC molecule bound to a chitosan polyelectrolyte chain of degree of polymerization n_p = 45, for these parameters the transition mushroom-brush lipopolyelectrolyte molar fraction is X_p = 0.014, from Fig. 3–5, it is clear that for X_p = 0.005 the lipopolyelectrolyte show a mushroom regime, for X_p = 0.014 show a critical regime and for X_p = 0.1 show a brush regime. From molecular dynamics simulation, we find that the thickness of the chitosan layer increases with the increase of the grafting molar fraction. The chitosan layer for n_p = 45 and X_p = 0.005, 0.014 and 0.1 are found around 3.2 nm, 3.5 nm and 6.6 nm, respectively. Previous experimental results for the PEG layer with the same parameters, using the X-ray diffraction method, give a polymer layer in the range of 3.75–6.6 nm.^64,67 The MD results for chitosan polyelectrolyte are in good agreement with the experimental findings for PEG neutral polymer. Thus, the incorporation of DPPC-chitosan lipopolyelectrolyte prevents the adhesion liposome–liposome even with a small molar fraction of grafted polyelectrolyte, and they form a layer which protects the liposome from adhesion phenomena and from supposed system immune attack. Then, the chitosan can be used as an alternative for liposome covering.

Table 5 Thickness of the polyelectrolyte layer obtained by MD simulations L_DM for DPPC–chitosan_np lipopolyelectrolyte grafted in DPPC bilayer membranes, and from theoretical prediction L_Th

Xp	np	Regime	LDM (nm)	LTh (nm)
a Polymer layer from eqn (10), with a = 0.47 nm.b Polymer layer from eqn (17), with a = 0.47 nm and Al = 0.7 nm^2.
0.005	45	Mushroom	3.2 ± 0.2	3.15^a
0.014	45	Critic	3.5 ± 0.2	3.47^b
0.1	45	Brush	6.6 ± 0.2	6.68^b

---

3.3 DPPC bead density distribution

The bead density profile along the normal of the membranes-bilayer is analyzed to investigate the effect of chitosan grafting density on the DPPC bilayer-membrane thickness D_B and the hydrocarbon thickness D_HH. The D_B and D_HH were calculated for the three grafting regimes, at the same conditions of external-pressure and temperature (T = 300 K and P = 1 bar). Fig. 7a, b and c represent the density profiles of all bead types (C, N_a, Q_a, and Q₀) for the mushroom, critic, and brush regimes, respectively. The beads density profiles across the membrane show a slight extent and a slight decrease in the monolayer in which the lipo-polyelectrolyte is incorporated for comparison to the liposome intra-monolayer. This deviation increase proportionally to the grafting molar fraction X_p. Previous MD simulation study on DPPC bilayer-membranes, without polymer coating, revealed D_B = 3.85 nm.⁶⁸

Fig. 7 Reduced density profile of DPPC beads: (a) reduced density in the mushroom regime, (b) reduced density in the critic regime, (c) reduced density in brush regime.

In the present study, the D_B, given by the distance between the maxima of the density profile of Q0 beads, amounts to 4.1, 4.16 and 4.2 nm for mushroom, critic, and brush regime, respectively. Thus, the incorporation of lipo-polyelectrolyte increases the membrane thickness from 0.15 to 0.35 nm depending on X_p. For the hydrocarbon, the bilayer beads distribution is a Gaussian like, since the two density peaks have overlapping, so they are unresolvable and they look like one peak. A closely related quantity is the half-width at half maximum (FWHM). Thus, the FWHM is used to define hydrocarbon thickness. According to this concept, the hydrocarbon thickness is found 2.75 nm, 2.79 nm and 3.05 nm for mushroom, critic and brush regimes, respectively. From the literature, Norbert et al. have performed an experimental study combining the analysis of X-ray and neutron ULV data to measure DPPC bilayer D_B and DPPC hydrocarbon D_HH thicknesses.⁶⁹ They found D_B around 3.9 nm and D_HH around 2.85 nm. In comparison to their works, our results are in good agreement with the experimental findings, in terms of D_B and D_HH.

3.4 Area per DPPC-lipid in DPPC-liposome with incorporated DPPC-chitosan lipo-polyelectrolyte

The simulations were performed at the ambient conditions of temperature and pressure (T = 300 K, P = 1 bar). For these conditions, the membrane bilayer is found in the lamellar phase giving an area per lipid A_l around 0.6 nm², 0.62 nm², and 0.66 nm² for mushroom, critic, and brush regimes, respectively. This change of the area per lipid is due to the steric interactions between the chitosan chains, which are collected in the brush regime. Thus, giving rise to the lateral pressure which will tend to stretch the lipid membrane. Among other things, this will change the membrane permeability properties, by chain fusion. Experimental proof of such lateral expansion, is provided by the increase of the movement freedom of lipid chains with the incorporation of lipopolymer into bilayer membranes.⁷⁰

The extent of the resulting membrane expansion is determined by Marsh^71,72 using a precise theoretical treatment using a virial equation of state as:

with, , such as, A_l the area per lipid in the presence of grafted chitosan brush, A_l,0 the area per DPPC lipid in the absence of the chitosan brush, X_p the grafting molar fraction and m_F a critical exponent, which can take two values according to the used theoretical approach, for the scaling theory or for the mean-field theory.^71,72 From the present MD simulation, the membrane expansion is found around 0, 0.03 and 0.1 for mushroom, critic and brush regimes, respectively. The obtained values cover the theoretical predictions (Fig. 4 of ref. 71).

3.5 Radial distribution function (RDF)

The radial distribution function (RDF) shows the probability of finding a particle distributed around another central particle. It provides information on the local structure of coarse-grained beads in a group defined among the various beads of the system.where, ρ = N/V is the density of the system, N the number of particles and V the volume of the simulation box.

Fig. 8 shows the radial distribution function (RDF), which gives the probability of finding a coarse-grained bead at a distance r from another bead of the same considered group of molecules. For polyelectrolyte beads group (Fig. 8a), in the brush regime, when the chains are more stretched, the polyelectrolyte RDF has one narrow peak, at r = 5.2 Å, corresponding to an ordered monomers distribution, which represents the inter-chain interaction. In the mushroom regime, when the chains are congested, the polyelectrolyte RDF has two narrow peaks, the first peak, appearing at r = 5.2 Å, represent the intra-chains interaction and a second large peak centered at r = 9 Å due to the strong correlation between the chitosan monomers. Inspection of the radial distribution function confirms that, for brush, critic and mushroom regimes, the membrane bilayer is in a fluid phase, the latter is characterized by an ordered distribution at a short distance and a disordered distribution at long distance. The two main peaks, located at r = 5.2 Å and 9 Å, respectively, of the RDF between lipids, increase proportionally to the graft polyelectrolyte molar fraction, i.e., the membrane lipids are more congested for the brush regime (Fig. 8b). The incorporation of lipopolymer affects the liposome structure positively, the stability of DPPC lipids bilayer increases with the increase of the grafting molar fraction. Fig. 8c shows the RDF of the charged group in the system (polyelectrolyte + ions), two main peaks are observed at r ≈ 5 Å and r ≈ 9 Å, respectively. The RDF main peaks are higher for the low PE incorporation, in the mushroom and critic regimes, because the ions come close to the PE, while in the brush regime the ions are pushed away because of the sterical effect. Fig. 8d shows the RDF of Cl⁻ ions, the RDF present two main peaks located at r ≈ 5 Å and r ≈ 8.5 Å, respectively. The first peak is higher for the brush regime, in contrast, the second peak is higher for the mushroom and critic regimes, in the mushroom and critic regimes the ions are distributed in the whole volume. However, in the brush regime, the ions can form confinements, near the end of the chains, because of the PE-chains steric effect.

Fig. 8 RDF of membrane grafted chitosan components: (a) RDF of the group of chitosan molecules, (b) RDF of the group of lipid molecules, (c) RDF of the group of (chitosan molecules + Ions), (d) RDF of the group of ions.

3.6 Dynamics properties

3.6.1 Theoretical background of diffusion. Liposomes are constructed from lipid bilayers decorated with polymers and dispersed in a solvent. The membrane of considered liposome is formed by DPPC phospholipids. The dynamic behavior of these lipid bilayer systems has been studied in-depth in terms of molecular dynamics simulations using both all-atoms and coarse-grained approaches and using different experimental tools.^73–79 Within the framework of MD simulations, the diffusion is studied by analyzing the mean-square-displacement expressed as follows:

In the expression above, r_i(t) represents the temporal position of a random walker (i), which can be a DPPC bead, a chitosan monomer, a water bead, or an ion. Here, t indicates the time and t₀ is an initial time at which the random walker begins to move. We recall that, for the standard diffusion of simple liquid (water, ions and so on), the particles move ballistically at short times. Thus, the mean-square-displacement 〈Δr²(t)〉 ∼ t², which is followed by a crossover to Fickian diffusion, characterized by 〈Δr²(t)〉 ∼ t for long times. On the other hand, for molecular systems such as polymers and lipids, a caging effect causes a subdiffusion regime that intermediates the ballistic and the diffusive regimes. In this case, the particles are trapped by their neighbors. Thus, this regime is characterized by 〈Δr²(t)〉 ∼ t^α, with 0 ≤ α ≤ 1. The value of α in the present of the cage effect remains astounding. In fact, there are some theoretical approaches based on some complex memories functions developed to discuss the subdiffusion laws. From the literature, Zwanzig in a series of works developed a theoretical approach based on a generalized Langevin equation (GLE),^80,81 this approach can be adopted to discuss the subdiffusion observed in polyelectrolyte and lipids bilayer. From a mathematical viewpoint, this is simply an extension of the standard Langevin equation for simple liquid, where the friction is assumed to be determined by the instantaneous velocity of the particle. It is noted that the difficulty of understanding the subdiffusion phenomena is to mathematically explicit the cage effect that, physically speaking, depends on various physical parameters, namely, the temperature, the particle density, the particles arrangement nature and so on. In this context, the GLE is expressed as:

where v is the velocity of a moving particle, is the relaxation rate, where ζ is the friction coefficient and m is the mass of the tracer, κ is the memory-function that expresses the friction retardation, and F_s(t) is a random force felt by the moving particle due to its collisions with other particles. Thus, the random force verify the equation:

〈F_s(t)·F_s(0)〉 = 6mk_BT[γδ(t) + κ(t)], t > 0 (22)

For this considerations VACF solves the following differential equation:

To resolve the GLE equation in the coarse-grained representation of DPPC lipids and chitosan chains, we adopt a recent memory function, proposed by Flenner,⁸² developed to study the dynamic aspect of lipid atoms. Flenner's approach is based on the Zwanzig–Mori projection method for modeling the 〈Δr²(t)〉. The starting point is to consider the equation of motion for density autocorrelation function Φ_s(q,t) = 〈n(−q,0)n(q,t)〉 of a selected particle at wave vector q,⁸³

where , k_B stand for the Boltzman's constant, T denote the temperature and m is the coarse-grained bead mass.

The equation of motion for the MSD in two dimension can be obtained from

which giveswhere .

Flenner et al. propose a memory function, with explicit crossover times scales and which reproduce the three diffusion regimes, namely, ballistic, normal, and subdiffusive⁸²

where δ(t) is the Dirac delta function, B is a dimensionless parameter, τ₁ is the characteristic time for the crossover from the subdiffusion to the normal diffusion, τ₂ is the onset time of the subdiffusion regime, and τ₃ is the characteristic time for the crossover from the ballistic to the subdiffusive regime.

Analysis of this memory function revealed that, for B = 0, the Brownian diffusion laws, expected for ions, are recovered and the 〈Δr²(t)〉 behaves like:

If (t/τ₂)^α ≪ 1

Generally, the long-term mean-square-displacement 〈Δr²(t)〉 shows a normal diffusion with α = 1 common for all constituent particles, in good agreement with the Langevin theory. But the peculiarity of congested systems (chitosan chains and DPPC lipids) is that the mean-square-displacement exhibits a sub-diffusion regime, with 0 < α < 1, extended over a large time interval. This is explained by the fact that the coarse-grained beads, instead of performing a random walk, its movement is blocked by its neighbors and move with a concerted manner, in particular for DPPC lipids at in vivo and in vitro temperatures, which are generally below their melting temperature found around 314 K.

The subdiffusion phenomenon is represented mathematically in the memory function by the term of the power law. The analytical determination of the diffusion exponent α seems to be complex. However, a numerical computation, taking into account the term of power law in the memory function,⁸² revealed that the dynamics pass slowly from the ballistic regime to the subdiffusive regime with α < 1. The region of subdiffusion separating the ballistic and the normal regimes is due to the cage effect. Thus, in this region the mean-square-displacement, 〈Δr²(t)〉, behaves like:

〈Δr²(t)〉 ∼ t^α, (α < 0),  t > τ₂ (30)

The dynamics of the chitosan chains, in the mushroom and critic regimes, can be also discussed using the Rouse model,⁸⁴ which is a simple model based on the entropic spring force acting between monomers. Thus, the equation of motion of the chitosan-beads is described by the following Langevin equations:

where, (R₁, R₂, …, R_N ≡ R_n) represent the beads positions, k = 3k_BT/σ_R2² depicts the entropic spring force and f_n a random force representing the random collisions between monomers, ξ is a friction coefficient related to the diffusion coefficient, at large time, by D_α=1 = k_BT/ξ.

The solution of the Langevin equation gives the following behavior of the mean-square-displacement as a function of time.

with τ is a crossover time, d the space dimension d, and D_α the normal diffusion coefficient.

3.6.2 NVT-NVE MD simulation results. The dynamics laws of the system components (ions, chitosan, and lipids) are identified at the three studied chitosan conformations, using an NVT-NVE MD simulation. Fig. 9, 10 and 11 show the mean-square-displacement as a function of time in the linear–linear scale and in the log–log scale, for chitosan chains, DPPC lipids, and Cl⁻ ions, respectively. At a fixed temperature T = 300 K. The diffusion coefficients, that are the slopes of the best linear log–log plot of the MSD as a function of time, are given in Tables 6–8. In Fig. 9b, we depict the log–log plot of MSD against time for chitosan, for three values of the grafted chitosan molar fraction correspond to the mushroom, critic and brush regimes. We first remark that after a short initial regime of subdiffusive motion (t = 100 ps), the chitosan dynamics reach the normal diffusion, after this time MSDs are straight lines with a slope α ≃ 1, but decrease as the polymer pass from the mushroom to the brush regime. We note that α < 0.5 for the short time, this deviation from the Rouse model is attributed to the cage effect provoked by the high density. In a recent experimental study, the diffusion coefficient D_c of the polyelectrolyte inside a gel matrix was measured using dynamic light scattering and was found to be in the order of 10⁻⁸ cm² s⁻¹.⁸⁵ The MD results cover the experimental findings and show a dependence of the diffusion coefficients to the grafting densities of chitosan chains.

Fig. 9 Mean square displacement of confined chitosan chains: (a) 3 dimensions MSD as a function of time, (b) log–log of MSD as a function of time.

Table 6 Diffusion coefficient D_α of chitosan polyelectrolyte obtained by MD simulations in the NVT-NVE statistical ensemble

Regime	t/fs	α	(Dα) (cm^2 s^−α)
Brush	[0,10^5.3]	0.22	0.09 × 10^−8
	[10^5.3,10^7]	0.71	2.94 × 10^−8
Critic	[0,10^5.32]	0.45	0.10 × 10^−8
	[10^5.32,10^7]	0.80	8.12 × 10^−8
Mushroom	[0,10^5.31]	0.45	0.10 × 10^−7
	[10^5.31,10^7]	1	1.43 × 10^−7

---

Table 7 Diffusion coefficient D_s of DPPC lipids obtained by MD simulations in the NVT-NVE statistical ensemble

Regime	t/fs	α	(Dα) (cm^2 s^−α)
Brush	[0,10^5.2]	0.44	0.50 × 10^−7
	[10^5.2,10^7]	0.6	1.77 × 10^−7
	[0,10^5.1]	0.44	0.62 × 10^−7
Critic	[10^5.1,10^6.7]	0.6	1.25 × 10^−7
	[10^6.7,10^7]	1	3.61 × 10^−7
	[0,10^5.1]	0.44	0.66 × 10^−7
Mushroom	[10^5.1,10^6.5]	0.6	2.46 × 10^−7
	[10^6.5,10^7]	1	4.57 × 10^−7

---

Table 8 Diffusion coefficient D_α of confined ions obtained by MD simulations in the NVT-NVE statistical ensemble

Regime	t/fs	α	(Ds) (cm^2 s^−α)
Brush	[0,10^7]	1	1.52 × 10^−6
Critic	[0,10^7]	1	1.64 × 10^−6
Mushroom	[0,10^7]	1	2.47 × 10^−6

---

In Fig. 10b we depict the log–log plot of MSD against time for DPPC lipids. The dynamics behavior of DPPC lipids can be affected by the presence of incorporated DPPC-chitosan lipo-polyelectrolyte with different grafting molar fraction. For a pure DPPC lipid bilayer with different molar fractions of incorporated DPPC-PEG lipo-polymer, we observed anomalous diffusion of lipids with a scaling exponent α < 1, depending on time. Also, we found that the calculated diffusion coefficient of lipids, at large time t ≥ 10^6.6 fs, is of the order of 10⁻⁸–10⁻⁷ cm² s⁻¹ and decreases proportionally to the DPPC-PEG molar fraction. This dependence on the molar fraction is proved experimentally by M. L. Wagner and L. K. Tamm using fluorescence recovery after photobleaching (FRAP).⁸⁶ In their work, the authors found that high lateral lipid diffusion is observed, in supported lipid bilayers on a polyethyleneglycol (PEG) cushion, when the molar fraction of graft PEG kept slightly below the transition to the brush regime. In this work, we found that the scaling exponent for DPPC lipid follows the same instance, α < 1. But, we note that α found for DPPC decorated with chitosan polyelectrolyte is slightly lower for comparison to that found for DPPC membranes decorated with PEG neutral polymer. This can be explained by the fact that the PEG polymer chains are more flexible compared to those of chitosan. On the first hand, the calculated diffusion coefficient of lipids, in the present study, is also found in good agreement with the previous MD simulations of DPPC membranes with incorporated neutral polymers. On the other hand, the obtained values of D_α cover the experimental findings.^86–88

Fig. 10 Mean square displacement of membranes: (a) 2 dimensions MSD as a function of time for DPPC molecules (the component of MSD in the z-direction is neglected), (b) log–log of MSD as a function of time for DPPC molecules.

In Fig. 11b, we depict the log–log plot of MSD against time for ions, for three values of the grafted chitosan molar fraction correspond to the mushroom, critic and brush regimes. We first remark that the ions dynamics follow a normal diffusion, the MSDs curves are straight lines of the same diffusion exponent α = 1, but decrease as the polymer pass from the mushroom to the bush regime. The calculated D_c for Cl⁻ ions in the chitosan branched polyelectrolyte, confined between two adjacent liposomes, is of the order of 10⁻⁶ cm² s⁻¹, these values are in good agreement with the experimental measurements of diffusion coefficients in aqueous solutions,⁸⁹ and are in accord with the calculated values from recently published works using all-atom MD simulations.⁹⁰

Fig. 11 Mean square displacement of confined ions: (a) 3 dimensions MSD as a function of time, (b) log–log of MSD as a function of time.

We note that the ballistic regime is not observed for all constituents, this can be explained by the fact that the simulated system is previously equilibrated in the NPT-NPH MD simulation.

4 Conclusion

The aim of this study is to examine the structural and dynamic properties of DPPC membranes covered by the chitosan chains using a coarse-grained MD simulation method, with particular emphasis on the application of this system in drug delivery. Chitosan chain conformations are characterized in the NPT-NPH conditions, as a function of the grafting molar fraction X_p, through calculating the length of the chitosan layer. On one hand, the obtained results are found to be consistent with the theoretical investigations, which predicts three different conformations of chitosan chains depending on X_p, namely, mushroom, critic and brush regimes. On the other hand, the calculated lengths of chitosan layers are found in agreement with the available studies on neutral PEG polymer, from MD simulation and X-ray experiment. The dynamics properties have been investigated in the NVT-NVE conditions, the dynamics of the membrane slowing during the transition from the mushroom regime to the brush regime. The chitosan chains move freely following a normal diffusion and increase the dynamic flexibility of the liposomes. The ions follow a normal diffusion for all the conformation regimes independent of diffusion time. The calculated diffusion coefficients for all constituents namely, DPPC lipids, chitosan chains, and ions, are found in good agreement with the experimental ones, available in the literature. Regarding the application of chitosan covering liposomes in drug delivery, the present findings confirm the role of coated chitosan in the protection of liposomes against the adhesion phenomena expected between liposomes–liposomes and liposomes–target cells. Also, these findings provide additional information about the role of chitosan in the stabilization and the lubrication of liposomes. Thus, these results will provide a quantitative basis for designing liposomes covered with chitosan for drug delivery applications. Finally, we underline that grafting chitosan to the liposomes can also affect other properties of its membranes bilayer, for example, the melting temperature and the elastic properties.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors would like to express their sincere gratitude to Professor Mabrouk Benhamou, from Moulay Ismail University, for several helpful discussions.

Notes and references

1. H. Daraee, et al., Artif. Cells, Nanomed., Biotechnol., 2016, 44, 381–391.
2. C. Ross, et al., Int. J. Nanomed., 2018, 13, 8507–8522.
3. M. Houshmand, et al., Nanomaterials, 2020, 10, 276.
4. F. N. Kök et al., Biomimetic Lipid Membranes: Fundamentals Applications and Commercialization, Springer, 2019.
5. C. I. Nkanga et al., General Perception of Liposomes: Formation, Manufacturing and Applications, Liposomes-Advances and Perspectives, IntechOpen, 2019.
6. D. C. Drummond, et al., Pharmacol. Rev., 1999, 51, 691–744.
7. Y. Okamoto, et al., J. Drug Delivery Sci. Technol., 2020, 55, 101381.
8. N. S. Templeton and D. D. Lasic, Mol. Biotechnol., 1999, 11, 175–180.
9. G. Gregoriadis, et al., Vaccine, 2002, 20, B1–B9.
10. E. Fattal, et al., Adv. Drug Delivery Rev., 2004, 56, 931–946.
11. S. M. B. Souza, et al., Colloids Surf., B, 2004, 36, 13–17.
12. M. S. Martina, et al., J. Am. Chem. Soc., 2005, 127, 10676–10685.
13. F. Frézard, et al., Quim. Nova, 2005, 28, 511–518.
14. G. J. Charrois and T. M. Allen, Biochim. Biophys. Acta, Biomembr., 2003, 1609, 102–108.
15. U. Pleyer, J. Grammer, P. Kosmidis and D. G. Ruckert, Surv. Ophthalmol., 1995, 39, S3–S16.
16. Y. Xia, et al., Biomaterials, 2016, 79, 56–68.
17. S. Deodhar and A. K. Dash, Ther. Delivery, 2018, 9, 857–872.
18. W. Yan, et al., Nanomedicine, 2020, 15, 303–318.
19. M. L. Gonzalez Rodriguez, L. B. Barros, J. Palma, P. L. Gonzalez Rodriguez and A. M. Rabasco, Int. J. Pharm., 2007, 337, 336–345.
20. S. S. Feng, G. Ruan and Q. T. Li, Biomaterials, 2004, 25(21), 5181–5189.
21. J. Guo, Q. Ping, G. Jiang, L. Huang and Y. Tong, Int. J. Pharm., 2003, 260(2), 167–173.
22. Z. H. Wu, Q. N. Ping, Y. Wei and J. M. Lai, Acta Pharmacol. Sin., 2004, 25(7), 966–972.
23. K. Tai, et al., Food Hydrocolloids, 2020, 99, 105355.
24. A. A. Barba, et al., Coatings, 2019, 9, 491.
25. M. Jo, et al., Food Hydrocolloids, 2019, 94, 603–612.
26. M. Tekade et al., Biomaterials and Bionanotechnology, Academic Press, 2019, pp. 459–493.
27. M. Qiang, et al., Molecules, 2020, 25, 610.
28. M. Gao, et al., J. Pharm. Sci., 2019, 108, 3623–3629.
29. A. Haasbroek, et al., Pharmaceutics, 2019, 11, 36.
30. A. Kumar, et al., Int. J. Biol. Macromol., 2016, 91, 615–622.
31. W. Chen, et al., Int. J. Pharm., 2020, 119115.
32. N. Morin-Crini, et al., Environ. Chem. Lett., 2019, 1–26.
33. V. Nikolić et al., Characterization and Biology of Nanomaterials for Drug Delivery, Elsevier, 2019, pp. 587–625.
34. P. Wu, et al., J. Biomed. Nanotechnol., 2009, 5, 557–564.
35. M. Z. Karagozlu and S. K. Kim, Adv. Food Nutr. Res., Academic Press, 2014, vol. 72, pp. 215–225.
36. A. Martins, et al., Int. J. Mol. Sci., 2014, 15, 20800–20832.
37. D. H. Ngo and S. K. Kim, Adv. Food Nutr. Res., 2014, 73, 15–31.
38. B. Homayun, et al., Pharmaceutics, 2019, 11, 129.
39. H. Mndlovu, et al., Carbohydr. Polym., 2019, 222, 114988.
40. S. Monti, et al., Phys. Chem. Chem. Phys., 2019, 21, 13099–13108.
41. B. Cheng, B. Pei, Z. Wang and Q. Hu, RSC Adv., 2017, 7(67), 42036–42046.
42. T. Mengoni, et al., Pharmaceutics, 2017, 9, 56.
43. S. Supper, et al., Eur. J. Pharm. Biopharm., 2014, 88, 361–373.
44. K. A. Terrón-Mejía, et al., Nanomaterials, 2018, 8, 425.
45. S. Ilk, et al., Artif. Cells, Nanomed., Biotechnol., 2017, 45, 907–916.
46. C. Laye, et al., J. Food Sci., 2008, 73, N7–N15.
47. M. Lemaalem, R. Ahfir, A. Derouiche and M. Filali, RSC Adv., 2020, 10(59), 36155–36163.
48. M. M. Mady and M. M. Darwish, J. Adv. Res., 2010, 1, 187–191.
49. M. Hasan, et al., RSC Adv., 2016, 6, 45290–45304.
50. M. Lemaalem, et al., RSC Adv., 2020, 10, 3745–3755.
51. S. J. Marrink, et al., J. Phys. Chem. B, 2007, 111, 7812–7824.
52. C. A. López, et al., J. Chem. Theory Comput., 2009, 5, 3195–3210.
53. X. Periole and S. J. Marrink, Biomolecular Simulations, Humana Press, Totowa, NJ, 2013, pp. 533–565.
54. H. Xu and S. Matysiak, Chem. Commun., 2017, 53, 7373–7376.
55. H. A. Lorentz, Ann. Phys., 1881, 248, 127–136.
56. D. Berthelot “Sur le mlange des gaz”, Comptes rendus hebdomadaires des sances de l'Acadmie des Sciences, 1898, vol. 126, pp. 1703-1855.
57. S. Plimpton, et al., Sandia Natl. Lab., 2007, 18, 43.
58. M. P. Allen, Introduction to molecular dynamics simulation, Computational soft matter: from synthetic polymers to proteins, 2004, vol. 23, pp. 1–28.
59. R. L. Davidchack, et al., J. Chem. Phys., 2015, 142, 144114.
60. T. F. Miller III, et al., J. Chem. Phys., 2002, 116, 8649–8659.
61. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012.
62. S. T. Milner, et al., Europhys. Lett., 1988, 5, 413.
63. S. T. Milner, et al., Macromolecules, 1988, 21, 2610–2619.
64. A. K. Kenworthy, et al., Biophys. J., 1995, 68, 1921–1936.
65. E. B. Zhulina and M. Rubinstein, Macromolecules, 2014, 47, 5825.
66. P. R. Desai, et al., Phys. Rev. E, 2018, 97, 032503.
67. Y. Schilt, et al., Biochim. Biophys. Acta, Gen. Subj., 2016, 1860, 108–119.
68. S. J. Marrink, et al., J. Phys. Chem. B, 2004, 108, 750–760.
69. K. Norbert, et al., Gen. Physiol. Biophys., 2009, 28(2), 117–125.
70. E. A. Disalvo and S. A. Simon, Permeability and stability of lipid bilayers, CRC Press, 2017.
71. D. Marsh, et al., Biochim. Biophys. Acta, Biomembr., 2003, 1615, 33–59.
72. D. Marsh, Biophys. J., 2001, 81, 2154–2162.
73. N. Hadrioui, et al., RSC Adv., 2020, 10, 8568–8579.
74. N. Kučerka, et al., Biochim. Biophys. Acta, Biomembr., 2011, 1808, 2761–2771.
75. P. Schwille, et al., Cytometry, 1999, 36, 176–182.
76. N. Kahya, et al., J. Biol. Chem., 2003, 278, 28109–28115.
77. K. Bacia, et al., Biophys. J., 2004, 87, 1034–1043.
78. J. Yguerabide, et al., Biophys. J., 1982, 40, 69–75.
79. E. A. J. Reits and J. J. Neefjes, Nat. Cell Biol., 2001, 3, E145–E147.
80. R. Kubo et al., Statistical Physics II: Nonequilibrium Statistical Mechanics, Springer Science and Business Media, 2012, vol. 31.
81. R. Zwanzig, Nonequilibrium Statistical Mechanics, Oxford University Press, 2001.
82. E. Flenner, et al., Phys. Rev. E, 2009, 79, 011907.
83. J. HansenI. McDonaldTheory of Simple LiquidsElsevier, Amsterdam, 2006.
84. C. Bennemann, et al., Comput. Theor. Polym. Sci., 1999, 9, 217–226.
85. A. Rahalkar and M. Muthukumar, Macromolecules, 2017, 20, 8158–8168.
86. M. L. Wagner and L. K. Tamm, Biophys. J., 2000, 79, 1400–1414.
87. T. Wang, et al., Langmuir, 2010, 26, 11157–11164.
88. M. A. Deverall, et al., Biophys. J., 2005, 88, 1875–1886.
89. E. R. Nightingale, J. Phys. Chem., 1958, 63, 1381–1387.
90. Y. Ruan, et al., Langmuir, 2016, 32, 13778–13786.

---