Matrix interaction according to Eq. (2).
Liposomes are essential components in the development of functional materials for drug delivery; this is mainly due to its ability to self-associate spontaneously and form bilayer vesicles. In these potential applications, knowing the size of self-assembled liposomes is essential for optimal performance; however, this process still has many unanswered questions. Conventional experimental techniques to study self-assemblies of liposome nanoparticles still have a great challenge. Computational simulations emerge as an alternative to understand the role of thermodynamic properties responsible for the self-assembly, particularly when they are unreachable experimentally because of limited time and length resolutions. In this chapter, we present the advantages and disadvantages of dissipative particle dynamic method to explore the functioning of liposome self-assembly in the transport of drugs.
- dissipative particle dynamics
- liposomes self-assemblies
- drug delivery
- chitosan polymers
The scenario of the modern era of interdisciplinary research between the natural sciences, including mathematics and theoretical physics, biophysics, chemistry, colloid science, and biochemistry, among others, along with advances in computational sciences, makes possible the manipulation of matter on a molecular scale through nanotechnology and bionanotechnology. This interdisciplinary research is able to study nanosystems made of phospholipids can be self-assembled spontaneously. One of the most studied cases in recent years is the self-assembly of liposomes.
The liposomes, vesicular systems based on lipids, in the last decades have increased their potential of their applications in various scientific disciplines and extensively been used as carriers for numerous molecules in cosmetic and pharmaceutical industries and medical and agricultural fields due to their high capacity and efficiency in drug delivery [1, 2]. For purposes in drug delivery systems, the liposomes are suitable candidates, due to their multiple properties that make them unique in these processes, such as biocompatibility, biodegradability, low toxicity, ability to modify the pharmacokinetic profile of the drug loaded, etc. . In this way, liposomes have been successfully applied in the delivery of nucleic acid systems , in the intestinal lymphatic delivery , for enhanced oral delivery . Allen and Cullis documented that these advances have led to numerous clinical trials in such diverse areas as the delivery of anticancer, antifungal, and antibiotic drugs; the delivery of gene medicines; and the delivery of anesthetics and anti-inflammatory drugs . Despite the innumerable applications and advances of the drug delivery of liposomes, it presents enormous challenges still. For instance, He et al. argued that by modulating the compositions of the lipid bilayers and adding polymers or ligands, both the stability and permeability of liposomes can be greatly improved for oral drug delivery . In order to achieve direct, effective, and selective delivery of the liposome itself, a deeper understanding of the mechanisms for interaction between the water/ligands/polymers and the liposome is required . It has been reported in different sources that the mechanism of liposome formation is based on the unfavorable interactions occurring between amphiphilic compounds and water molecules, where the polar head groups of phospholipids are subjected to the aqueous phases of the inner, the hydrophobic hydrocarbon tails are associated with the bilayer, and spherical core shell structures are formed [10, 11, 12].
In addition to a large number of experiments that have been dedicated to the study of the mechanism of liposome formation, computational studies have made important contributions to the subject. These methodologies have been shown to be very effective in knowing in a very specific way the role played by molecular interactions in complex systems such as liposome self-assembly. Nowadays, it is increasingly common for research protocols to include a computational analysis prior to carrying out the experimental stage. Different techniques of computer modeling can facilitate quantitative understanding of experimental observations and secure fundamental interpretation of underlying phenomena . For instance, Thota and Jiang in an interesting review documented that with the use of simulation methods, it is possible to elucidate the mechanisms of drug loading/release, which are indispensable in the rational screening and designing new amphiphiles for high-efficacy drug delivery .
Most of the computational studies for drug delivery use molecular dynamics (MD) or Monte Carlo (MC) simulations to obtain structural and dynamic information. In the case of MD, for example, the study of interactions is carried out considering all atoms explicitly commonly. These interactions are evaluated by analytical functions that form the so-called force fields. The level of precision of a simulation with MD relies on the force field used. There are traditional force fields such as CHARMM , AMBER , and OPLS  used to study the self-assembly of complex systems such as liposomes. Due to the scheme of evaluation of the interactions, which is through both the intra- and intermolecular forces, the MD simulations are slow and computationally expensive when the number of atoms is very large. This limitation makes it almost impossible to reach the time scales (of the order of the microseconds and beyond) where the thermodynamic observables have greater physical significance. For this reason, this technique turns out to be little viable to study the systems planned in this work, although it is important to point out that diverse groups use it with very successful results.
Modeling technique alternatives to MD have been developed with the aim of overcoming the bottleneck when it is required to study systems composed of a large number of atoms, as well as the issue of length time scales. One of those techniques is the dissipative particle dynamics (DPD). The main philosophy of the DPD model is based on the so-called coarse-grained (CG) models. These models allow to diminish considerably the number of atoms of a given system, grouping a set of atoms, which can be a functional chemical group, into a single interaction site called a bead, always trying to maintain the physics of the correct model. This size reduction facilitates larger simulations of larger ones with larger systems, allowing a more efficient exploration of phase space, but losing an atomistic description. Despite the success of the DPD technique, very few groups have used it to obtain detailed information on various aspects of self-assembly of liposomes. Goicochea et al. used DPD for understanding the behavior of biopolymer brushes coating drug-carrying liposomes in an aqueous environment . Shen et al. with the help of DPD simulations, they explored how the tethered PEG polymers will affect the membrane wrapping process of PEGylated liposomes during endocytosis . Yamamoto et al. used DPD simulations to study the spontaneous vesicle formation of amphiphilic molecules in aqueous solution .
Another great advantage of the modeling methods such as DPD is to be able to incorporate their algorithms into powerful supercomputing technologies such as graphic processing units (GPUs). The combination of hardware and software with DPD has allowed us to explore systems at time scales, which were unattainable a few years ago. For instance, Wu et al.  presented a GPU-accelerated DPD with parallel cell-list updating. On the other hand, Nguyen et al.  achieve a speedup of up to 9.5x in the DPD simulations in the well-known LAMMPS code using GPUs. Phillips et al.  also achieve important speedup of DPD using GPUs but the HOOMD-blue code. Other important implementations of DPD in GPUs can be found in Refs [24, 25]. Our group has been a pioneer in the use and application of DPD methodologies using GPUs. We have developed a simulation code called SIMES  (which is an acronym for simulation at mesoscopic time scale), a production DPD software designed and implemented specifically to run on GPUs. SIMES can simulate a broad spectrum of complex systems, including the self-assemblies of liposomes. Readers interested in going deeper into this topic can consult our contributions in the following references [27-30].
Recently, in a previous work of our group, we have studied the stability of structures and the efficiency of the encapsulation of capsaicin, as well as the internal and superficial distribution of capsaicin and chitosan inside of a nanoliposome, constituted by lecithin and coated with a shell of chitosan by DPD simulations . We found that thermodynamic properties, such as the potential mean force, show that the interaction between capsaicin and chitosan polymers is very weak compared to that with lecithin. Besides, an association between capsaicin and chitosan in presence of lecithin is not likely to occur. In this work, we managed to combine the power of the DPD method with the power of the GPUs. The simulations reported have been the longest so far by regarding the encapsulation of capsaicin by liposomal nanoparticles.
Continuing with our research computational about self-assembly of molecular components for the design of bionanomaterials focused on drug delivery applications, the objective of this chapter is to show the advantages offered by the DPD technique to study the self-assembly mechanism of liposomes. As a specific case study, we show the formation of liposomes loaded with capsaicin. Capsaicin, a lipophilic drug, is the pungent vanilloid compound in spicy chili peppers. We chose capsaicin because it is approved as a drug for the treatment of chronic pains . It is also an excellent candidate to be transported by natural polymers such as chitosan . In addition, it is widely recognized in the treatment of urological disorders and the control of satiety and obesity , among other important applications.
2. Dissipative particle dynamic foundations
2.1 Theoretical formulation
The simulation of DPD is one of the most viable methods to study thermodynamic and structural properties of biophysical processes at mesoscopic time and length scales, unlike the simulation with MD, where insights at microscopic time and length scales can achieved. This great advantage of DPD simulation is that the DPD models are built under the coarse-grained scheme, which is obtained from the description of all-atom of a given system. This coarse-graining parameter plays a vital role and has significant impact on the speed of simulation . This feature makes the DPD simulations cheap and very fast computationally speaking.
The theoretical foundations of the DPD method are documented elsewhere. Readers interested can consult the following references [36-40] for more specific details. In this chapter, we will mention in a very summarized way the main mathematical aspects that govern the equations of this methodology. The DPD was introduced by Hoogerbrugge and Koelman in 1992  for simulations of hydrodynamic phenomena. It was subsequently modified by Español and Warren in 1995  who introduced the fluctuation-dissipation theorem in the original algorithm to add the frictional and stochastic forces. Therefore, a coupling between the statistical mechanic laws of the beads and Gibbs canonical ensemble is given, favoring the calculation of thermodynamic properties at longer time and length scales. Similar to molecular dynamics, the time evolution of each DPD particle can be calculated by Newton’s second law:
where , , and are respectively the position, velocity, and momentum vectors of particle i, and is the total interparticle force exerted on particle i by particle j. Specifically, , where a purely repulsive conservative force , a dissipative or frictional force , which represents the effects of viscosity and slows down the particle motion with respect to each other, and a random (stochastic) force , which represents the thermal or vibrational energy of system, is summed to obtain the total force, and these force components can be individually written as
where , , , and . The parameters and are the distance-dependent weight functions for the dissipative force and random force, respectively. The symbols , , and determine the strength of conservative, random, and dissipative forces, respectively. Moreover, in order to generate a correct equilibrium Gibbs-Boltzmann distribution, the dissipative and random forces have to satisfy the following relations:
where is the Boltzmann constant. According to Groot and Warren , we choose a simple form of and as follows:
Technically, DPD simulation differs from MD simulation in two respects. First, the conservative pairwise forces between DPD particles are soft-repulsive, which makes it possible to extend the simulations to longer time scales. Second, the scheme to maintain constant temperature in DPD is implemented in terms of dissipative as well as random pairwise forces such that the momentum is locally conserved, which results in the emergence of hydrodynamic flow effects on the macroscopic scale.
2.2 Mesoscopic model of liposomes
From the molecular point of view, a liposome, which is usually amphiphilic molecules, contains polar heads and hydrophobic hydrocarbon tails. This molecular structure is ideal to build a DPD model, since the polar part can be grouped in a single bead, while the hydrophobic part can be grouped in another bead. The liposome presented in this work is constituted by lecithin, which has two long hydrocarbon chains, is a major component of lipid bilayers of cell membranes and is a natural and biological amphiphile. Lecithin is a natural lipid mixture of phospholipids used frequently for the preparation of various nanosystem delivery vehicles, such as liposomes [9, 10]. According to the molecular structure of lecithin, this can be separated into three different beads that correspond to the head, neck, and tail groups, respectively. A schematic illustration of the fundamental self-assembly of the structure of a liposome adopted in this chapter is shown in Figure 1. In Figure 1A, a detailed molecular description at all-atom level for the lipid is shown; a mapping in the bead representation of the atomistic structure is also shown. Figure 1B is depicting the coarse-grained model of lecithin used in this work with its three fundamental parts. A snapshot of configurational structure of the liposome constructed based on lecithin is shown in Figure 1C. From this representation, the yellow spheres represent the hydrophobic part of the lecithin, while the red spheres represent the hydrophilic part.
On the other hand, we show the coarse-grained model of capsaicin, our target as case of study in drug delivery applications. The molecular structure of capsaicin is also ideal for building a DPD model. Similar to lecithin, the capsaicin can be grouped into three groups or beads, i.e., one bead for the head, neck, and tail groups, respectively. Figure 1A shows the coarse-grained model for capsaicin.
Finally, in Figure 1A, we show the coarse-grained model of water. This model in DPD simulation is already well known; for example, a bead DPD can group up to three water molecules. This model allows to conserve the thermodynamic properties of water.
2.3 Computational details of mesoscopic simulations
Mesoscopic simulations were performed to study the self-assembling of liposomes. The capsaicin was used as an example of drug delivery applications. For this end, the simulations presented in this work were made under canonical ensemble, where N (the total particle number), V (volume), and T (temperature) are kept constant. For this type of simulations at constant temperature, some parameters must be specified. In this work, in all the simulations, we use dimensionless units or DPD units, i.e., the temperature and density were of T = 1.0, ρ = 3.0, respectively, the mass of all particles is . The system studied (liposome in aqueous solution) is constituted by 1000, 2000, 3000, 4000, and 5000 molecules of lecithin; each lecithin is constituted by three beads as shown in Figure 1; the remaining particles form part of solvent or water molecules. A total of 200,000 DPD particles were simulated. The simulation cell has a length of 40 DPD units. In DPD simulations, a key factor to carry them out is the specification of the parameters of interaction between the different components that make up a given system; this is due to the presence of the conservative forces in Eq. (3). The DPD interaction matrix for each type of bead used in this chapter is shown in Table 1.
The nomenclature used in this table goes as follows: the symbols L1, L2, and L3 represent the head, neck, and tail beads of lecithin, while that W symbol represents the bead corresponding the beads of water.
On the other hand, the parameters of intramolecular interactions in lecithin are and for the bonds between two particles , while for angles between three particles, it is taken and . Finally, the cutoff radius is taken as ; the parameters and of the dissipative and random forces are equal to 3.0 and 4.5, respectively, parameters that appear in Eq. (2). The integration time step is , and the simulations are running over the 100 blocks of 10,000 steps every block, having a total of 1 × 106 steps, to reach a total of 4.8 μs. This simulation time is long enough to observe the self-assembly of the liposomes. With a traditional MD, the use of large supercomputing clusters is required to obtain this time scale. The equations of motion are solved using the velocity Verlet algorithm adapted to DPD . The presence of electric charges was not considered in these calculations.
3. Drug delivery applications through numerical experiments
For the study of self-assemblies of liposomes for drug delivery applications, three different processes are presented. First, we analyze lipid self-assembly. Second, self-assembly of capsaicin is presented. Third, we end with the analysis of self-assembly of capsaicin by lecithin. The formation of nanoliposomes was explored by varying the concentration of both lecithin and capsaicin molecules separately. These simulations were analyzed through its density maps. This property is very useful for refining structures or adjusting molecular models; besides, the density maps allow us to spatially inspect the average location of atoms of interest during a simulation .
3.1 Exploring self-assemblies of lipids
To explore the self-assemblies of lipids, the density maps of five different processes were obtained. In these processes, the effect of lipid concentration was explored. Figure 2 shows the snapshots of each of these processes. We observe that in Figure 2A, with a concentration of 1000 lipid molecules, the self-association of lipids is practically negligible. Only a slight formation of small clusters is perceived, similar to the density maps of very diluted solutions. By increasing the concentration of lipids to double, that is, to 2000 molecules as can be seen in Figure 2B, we observe that the association increases with respect to the previous case of 1000 lipid molecules. However, it is not yet a sufficient concentration to detect any indication of the process of self-assembly of lipids. Figure 2C shows the case for a concentration of 3000 lipid molecules. In this system, the formation of very faint regions is observed. The presence of these very light red dots is an indication that the lipids begin to self-assemble. To verify that the lipid molecules are indeed self-assembling, a simulation with 4000 lipid molecules was carried out. The result is shown in Figure 2D. As we expected, the formation of the red dots increases, and they begin more notoriously. To observe the stability in the formation of lipids, a simulation with 5000 molecules was carried out. At this concentration, the lipid association is already clearly observed. A picture of this system is seen in Figure 2E. The red regions show that the self-assembly process has already taken place. These are the optimal conditions in terms of concentration for the process of self-assembly of lipids to be carried out. For higher concentrations of lipids, greater than 5000 molecules, we will find a saturation of the system, which is an amount sufficient for the formation of a nanoliposome. An animation video of each of these systems can be found in the supplementary material section of this work.
3.2 Exploring self-assemblies of capsaicin
We conducted an exploration to study the self-assembly of capsaicin. For this purpose, we carried out four separate numerical experiments. The results of these simulations are found in Figure 3, analyzed through their density maps. These maps show regions of high and low density basically. The regions of low density are the ones that predominate the colors black and purple, while the regions of high density are those that are identified in lighter shades such as red or yellow. Figure 3A shows the density map for a system composed of 250 capsaicin molecules. Due to the low concentration of capsaicin molecules, the black and purple regions predominate, which are located in a density range between 0.02 and 0.04. In the second numerical experiment, we doubled the concentration of capsaicin molecules. Figure 3B shows the density maps for this system composed of 500 capsaicin molecules. The presence of regions in all reddish is notorious, which oscillates between 0.08 and 0.12 in density. However, purple regions are still observed, which indicate that this concentration of capsaicin molecules is not yet enough to have a homogenous system. Figure 3C shows the density map for a system composed of 750 capsaicin molecules. Unlike the density map of Figure 3B, the presence of purple areas is already minimal, while the presence of high-density areas is more consolidated. A higher density of capsaicin molecules is necessary to completely eliminate the purple areas. This can be seen in Figure 3D, which shows the density map for a concentration of 1000 capsaicin molecules. Note in this case how the purple areas have completely disappeared. The entire density map is in a homogenous system with a density that ranges between 0.2 and 0.25. This concentration of capsaicin molecules is optimal for use in the process of formation of nanoliposomes.
3.3 Exploring self-assembly of capsaicin by lecithin nanoliposomes
We have documented that nanoliposomes are effective for the drug delivery process. A strategic requirement for the supply of medicines to a specific site has driven the development of active drug targeting techniques . In this way, to take advantage of the effective potential of the liposome technology and the impact on its formulations, it is necessary to know the optimal conditions in terms of concentration of each of the components that are required for an optimal application of any given drug. To explore the self-assembly of capsaicin by means of lecithin liposomes, we conducted four numerical experiments. The number of capsaicin molecules used in these four experiments was 250, 500, 750, and 1000 corresponding to a concentration of 30, 61, 92, and 123 mM, respectively, and the concentration of lecithin remained constant with 5000 molecules in the four cases. Figure 4 shows the density maps obtained for these four cases. Figure 4A–D shows the density maps for a system loaded with 250, 500, 750, and 1000 capsaicin molecules. For all four we observed a homogeneous formation of the liposome. The high-density regions are representing the region formed by the lipid membrane, which in this case is presented as the region with the highest concentration (in yellow) in the four cases. In addition, in the four cases, they present a redder region in the center of the map; this indicates that the density of the lipids is lower, as expected. Density maps also show that the nanoliposome is stable for 4.8 μs of simulation, since the lecithin molecules do not spread throughout the simulation box nor do they collapse in the aqueous core to form a micelle; this was possible thanks to the previous analysis where we found the optimal concentration of lecithin molecules in the self-assembly process.
In our previous study , we observed the formation of a nanoliposome based on lecithin was carried out with approximately 7000 molecules of lecithin, but in that study, the presence of natural polymers such as chitosan was included, a component that is not present in these simulations.
Derived from the results obtained in the density maps of Figure 4, we proceeded to make an analysis on the encapsulation efficiency (EE). It has been reported that both size and EE are key variables in the preparation of liposomes . However, these variables are difficult to control at the laboratory level even with sophisticated experimental equipment. Computational tools such as DPD are ideal for guiding experiments. The size and EE for the system analyzed in this work are show in Table 2. EE is obtained from the following equation: EE = ((CapsT − CapsF)/CapsT) × 100, where CapsT is the total concentration of capsaicin in the system and CapsF is the free capsaicin in solution. The percentage of EE for the capsaicin obtained in this work was of 96%, which is very close to that reported experimentally of 92% , besides, is very similar to the obtained percentage in our previous study, where the presence of chitosan polymer was included .
|Concentration of capsaicin (Mn)||Size (nm)||± (nm)||EE (%)||± (%)|
In this chapter, simulations at the mesoscopic scale were used to systematically study the self-assembly of liposomes and their potential applications toward the delivery of drugs. The advantages offered by the mesoscopic methods such as DPD over other simulation methods such as MD were exposed; additionally, the acceleration that the DPD simulations have achieved through the new GPUs technologies was also covered. With these technologies, we show simulations in the microsecond scale with a system composed of 150,000 particles in the process of self-assembly of liposomes. We show the optimal conditions for the formation of liposomes in function of the lecithin concentration to achieve a good capsaicin encapsulation. Our percentage of EE obtained is quite acceptable with respect to experimental measurements and simulation studies of similar systems.
This work was supported by CONACYT [Project FON.INST./44/2016] and SIEA-UAEM [Project 4728/2019CIB]. All simulations reported in this work were performed at the Supercomputer OLINKA, located at the Laboratory of Molecular Bioengineering at Multiscale (LMBM), at the Universidad Autónoma del Estado de México.
Conflict of interest
The authors declare that they have no competing interests.