An in Silico Approach to Reveal the Nanodisc Formulation of Doxorubicin

Molecular dynamic behaviors of nanodisc (ND) formulations of free doxorubicin (DOX) and DOX conjugated lipid prodrug molecules were investigated by molecular dynamics (MD) simulations. We have unveiled how formulation design affects the drug release profile and conformational stability of ND assemblies. Our simulation results indicate that free DOX molecules loaded in the ND system experienced rapid dissociation due to the unfavorable orientation of DOX attached to the lipid surface. It is found that DOX tends to form aggregates with higher drug quantities. In contrast, lipidated DOX-prodrugs incorporated in ND formulations exhibited sufficient ND conformational stability. The drug loading capacity is dependent on the type of lipid molecules grafted on the DOX-prodrug, and the drug loading quantities in a fixed area of NDs follow the order: DOX-BMPH-MP > DOX-BMPH-TC > DOX-BMPH-PTE. To gain further insight into the dynamic characteristics of ND formulations governed by different kinds of lipidation, we investigated the conformational variation of ND components, intermolecular interactions, the solvent accessible surface area, and individual MSP1 residue flexibility. We found that the global conformational stability of DOX-prodrug-loaded ND assemblies is influenced by the molecular flexibility and lipidated forms of DOX-prodrug. We also found that the spontaneous self-aggregation of DOX-prodrugs with increasing quantities on ND could reduce the membrane fluidity and enhance the conformational stability of ND formulations.

Introduction

Doxorubicin (DOX), a potent anthracycline cytotoxic drug, has been routinely used as a frontline chemotherapeutic agent in the treatment of various cancers. The anticancer mechanisms of DOX are known to involve multiple mechanisms, including intercalation of DOX into DNA double helix (Coldwell et al., 2008), inhibitions of several molecular targets (topoisomerase enzymes I and II, DNA and RNA polymerase) (Momparler et al., 1976; Foglesong et al., 1992; Marinello et al., 2018), generation of reactive oxygen species (ROS) (Davies and Doroshow, 1986), upregulation of C6 ceramide level (Ji et al., 2010), induction of autophagic cell death (Velez et al., 2011), and immunostimulatory effects (Mattarollo et al., 2011).

However, its clinical utility is restricted by dose-dependent toxicity on noncancerous cells of major organs, including the heart, liver, kidney, and brain (Chatterjee et al., 2010; Octavia et al., 2012; Tacar et al., 2013). Due to these unspecific side effects of conventional DOX therapy, various drug delivery systems (DDSs) such as liposomes, micelles, dendrimers, and inorganic nanoparticles have been developed to enhance therapeutic efficacy minimize the adverse effects of DOX. However, the optimized biocompatibility and pharmacokinetics of these DDSs are critically dependent on PEGylation technology. This approach potentially results in several negative aspects, including restriction of cellular uptake, anti-PEG immune response, non-degradability, and heterogeneity in production. These unfavorable impacts have plagued the clinical translation of most DDS, and hence this has motivated the discovery of alternative PEG-free DDSs.

Endogenous and synthetic lipoproteins, such as high-density lipoprotein (HDL)-like nanodiscs, have rapidly evolved as a promising DDS for small molecules, peptides, and nucleic acids because of their intrinsic long in vivo half-lives and passive targeting properties (Kuai et al., 2017; Kuai et al., 2018; Kadiyala et al., 2019). Therefore, nanodisc drug formulation has attracted significant attention in the field of nanomedicine (Lacko et al., 2015; Wang et al., 2018; Tang et al., 2021). Several studies regarding nanodisc-based DOX delivery that employed different encapsulation strategies have shown increased efficacy and reduced side effects in anticancer therapy (Murakami et al., 2010; Yuan et al., 2013; Wang et al., 2014; Kuai et al., 2018). However, nanodisc (ND) formulation of DOX approached by different encapsulation designs can lead to variable drug-loading capacity (DLC), distinct drug release profiles, altered structures, and surface charges, as well as stability of the ND system. These characteristics are critical for DOX-DDS development. Although conventional trial and error approaches may provide sufficient guidance on formulation development, this process is usually sophisticated, time-consuming, and costly (Boruah et al., 2019). Besides, the rational design of DDS formulation depends on a detailed understanding of molecular principles; the mechanistic insights at the molecular level are hard to obtain through experimental research (Ramezanpour et al., 2016). Theoretical studies based on computational modeling methods enable the initial screening of formulation variables to predict potential conditions to accelerate the experimental process (Mehta et al., 2019). In this regard, molecular dynamics (MD) simulations offer a powerful tool to inspect structural, molecular interactions, and dynamic features of drugs and carriers (Bunker and Rog, 2020). Taking into account these considerations, it appears desirable to initially predict the optimum ND formulation of DOX through computational methods.

Here we present an in silico approach that adopted the ND modeling and MD simulation techniques to investigate the effect of different formulation designs on DOX loaded ND DDS (Figure 1). Free DOX and lipidated-DOX prodrugs are designed to access the formulation variables on the structural and dynamical behavior of ND assemblies. The present study focuses on the influence of different lipidated forms of DOX and the drug-to-lipid ratio (D/L ratio) on drug release profile, the structural and dynamic properties of ND-DOX formulations. Moreover, the mechanism of the interaction between DOX and the ND lipid bilayer is poorly understood. This drove us to investigate the behavior of free DOX and lipidated-DOX prodrug molecules in the DPPC bilayer ND through MD simulations.

FIGURE 1

FIGURE 1. A flow chart summarizes the individual steps for the MD simulation of ND systems.

Materials and Methods

Preparation of the Nanodisc System

Membrane scaffold protein 1 (MSP1) and 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) were assembled into NDs using CHARMM-GUI (Qi et al., 2019). The diameter of all nanodiscs was 9.8 nm. For preparing four different ND formulation systems: 1. Free DOX, 2. DOX-BMPH-PTE (DOX-N-β-maleimidopropionic acid hydrazide-1,2-Dipalmitoyl-sn-Glycero-3-Phosphothioethanol), 3. DOX-BMPH-MP (DOX-N-β-maleimidopropionic acid hydrazide-2-Mercaptoethyl palmitate), and 4. DOX-BMPH-TC (DOX-N-β-maleimidopropionic acid hydrazide- Thiocholesterol), four corresponding structural analogs (1. Sitosteryl Glucoside, 2. DPPC, 3. Hexadecanoate, and 4. Campesteryl Glucoside) were selected from CHARMM-GUI for subsequent positional labeling. The selected structural analogs were randomly inserted into the ND system with a predefined set of proportions (10, 25, and 50%). Then DOX and Lipidated-DOX were superimposed to the labeled analog using flexible docking, and the labeled analog was subsequently removed. To eliminate the impact of the wrong docking poses on simulation accuracy, only the pose where the DOX sugar moiety is at the hydrophilic interface, and the molecule as a whole is wrapped in the pocket of the membrane layer is preserved. Each step in the process is presented in a flowchart (Figure 2). The control system (empty ND) consists of two MSP1 and 190 DPPCs, and the concrete composition of other ND formulation systems is described in Table 1. Each MSP1 contained 200 residues and a total of 400 residues in each ND system. The system was protonated at 298 K, the salt concentration was set to 0.145 M, and the pH was set to 7.35. All systems were solvated in cubic boxes using TIP3P water and neutralized by adding 0.145 M NaCl. Energy minimization for up to 10,000 steps was performed using the gradient descent algorithm to optimize the system.

TABLE 1

TABLE 1. The composition of each ND formulation system.

FIGURE 2

FIGURE 2. MD simulation results of free-DOX incorporated DPPC ND formulation. The RMSD trajectory of DPPC, MSP1, and ND complex as a function of simulation time in nanoseconds (ns) for (A) the empty ND, (B) 10% DOX-ND, (C) 25% DOX-ND, and (D) 50% DOX-ND. (E) The RMSD trajectory of DOX molecules with increasing proportions (10, 25, and 50%) as a function of simulation time in ns. The initial and final snapshot of the simulation for each DOX-ND formulation: (F) 10% DOX-ND, (H) 25% DOX-ND, and (J) 50% DOX-ND. The interactions of the DOX molecules with the ND differed depending on the proportions of DOX: (G) 10% DOX, (I) 25% DOX, and (K) 50% DOX.

Parameter Settings of Molecular Dynamics Simulations

0.25 ns of NVT equilibration and 1.625 ns of NPT equilibration were performed. System temperature and pressure were maintained at 298 K, and at 1 bar using the Berendsen algorithm, respectively. The Nose-Hoover and Parrinello-Rahman algorithms were used to maintain temperature and pressure in the production simulation (Martyna et al., 1992). The constraints of heavy atoms in the system were gradually released until free in the equilibration stage. The simulation with the time step of 2 fs lasting 100 ns was performed for each system. The CHARMM36 force field was applied to all simulations (Huang and MacKerell Jr, 2013). The LINCS algorithm was used to constrain hydrogen bonds, and Coulomb interactions were calculated using particle-mesh Ewald (Hess et al., 1997). The root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), solve accessible surface areas (SASA), and hydrogen bonds of the system were calculated in the trajectory analysis. All structures in the trajectory were clustered based on the RMSD of MSP1, and the structure in the largest cluster with the lowest RMSD relative to all other structures was defined as the central structure. All simulations were performed using GROMACS 2020.5 accelerated by CUDA modules, and all programs ran on AMD Ryzen™ Threadripper™ PRO 3995WX and The GeForce RTX^™ 3080 Ti graphics card (Van Der Spoel et al., 2005).

Results and Discussion

Dynamic Behaviour of DOX Molecules in DPPC ND Formulation

We initially investigated the MSP1/DPPC ND formulations of free DOX molecules through 100 ns all-atom MD simulations. Three different proportions of DOX ranging from 10 to 25% and 50% of total lipids were randomly placed into the ND bilayer. Considering the structural stability of ND and molecular motion of DOXs in the DPPC membrane, the root-mean-square deviation (RMSD) of DOX, DPPC, MSP (MSP1), and ND complex is determined in Figures 3A–E as a function of time. In the case of an empty ND system (Figure 3A), DPPC and MSP1 reach equilibrium in the first 12 ns of the 100 ns trajectories. Both DPPC and MSP show similar trends with an average RMSD (after 50 ns of simulation) of 1.00 and 1.20 nm, respectively. The global motions of the ND complex are plotted by the RMSD values comprising of DPPC and MSP in the trajectory, with RMSD quickly reaching a plateau of ∼1 nm after about 12 ns of simulation. In contrast to the empty ND system (control), the 10% DOX-incorporated ND assembly was unstable during the 100 ns of simulation (Figure 3B), with RMSD values undergoing a progressive increase of approximately 2.5-fold higher than control at the end of the simulation. However, The average RMSD value (after 50 ns of simulation) of DPPC and MSP1 is 1.18 and 1.22 nm, respectively. This indicates the significant RMSD increment of the ND complex is attributed to the random motion of DOXs.

FIGURE 3

FIGURE 3. RMSD trajectories and snapshots of center clusters for various types of lipidated-DOX prodrug ND formulations with increasing drug proportions (10, 25, and 50%). (A) ND formulation of DOX-BMPH-PTE. (B) ND formulation of DOX-BMPH-MP. (C) ND formulation of DOX-BMPH-TC.

To further study the role of DOX content in ND stability, ND with elevated DOX content of 25 and 50% were simulated under the same conditions. The structure of the ND complex (25% DOX) was increasingly unstable throughout the simulations, with the highest RMSD above 3.5 nm (Figure 3C). In contrast, the ND complex with higher DOX content (50%) generally shows lower RMSD with a maximum value of ∼2.6 nm (Figure 3D). However, the RMSD values of DPPC display the same increasing trend in 25 and 50% DOX of ND complex by maintaining an average of 1.31 and 1.35 nm during the second half of 50 ns, respectively. This implies the loss of constraints or reduced density between DPPC molecules. Moreover, all simulations show that an increase in the ratio of DOX results in little effect on the RMSD of MSP1. In this regard, the global instability mainly comes from the dissociation of DOX from the ND complexes. Interestingly, the RMSD profiles of 10 and 25% of DOX molecules in ND formulations show a similar trend with increasing RMSD values along the trajectories, while a higher proportion of DOX (50%) in the ND formulation leads to improved global stability, as reflected in lower RMSD pattern of the DOX molecules (Figure 3E). The simulation results suggest the self-aggregation propensity of DOX molecules.

In order to gain further insight into the dynamic behavior of DOX molecules, the first and last snapshots of MD trajectories are illustrated, and DOX-DPPC/MSP1 interactions are analyzed (Figures 3F–K). The drug leakage evaluation shows that 55% of DOX quantities were dissociated from the ND (10% DOX) throughout 100 ns MD simulations, suggesting variable interactions between DOX and DPPC ND bilayer (Figure 3F). However, the degree of DOX dissociation behavior is dropped to 25 and 15%, with increases in DOX-loading content from 25 to 50%, respectively (Figures 3H,J), indicating multiple network interactions were created via homomolecular and heteromolecular interactions. In the case of the low content of DOX (10%), DOX can interact with the DPPC molecules because of its amphiphilic structure. The positively charged amino sugar moiety can readily interact with the negatively charged phosphate group of DPPC through electrostatic interactions. The anthraquinone (AQ) ring of DOX can insert into the DPPC bilayer via hydrophobic interaction. In the case of the mid content of DOX (25%) (Figure 3I), additional interactions were formed between the DOX and polar or charged residues in MSP1. Specifically, Arg, Lys, and Glu residues frequently interact with the amino sugar and AQ moieties of DOX molecules. This results in the relocation of DOX molecules toward the ND rim, which is also reflected in the decrease of average RMSD values of MSP1 (Figures 3A–D). In addition to DOX-MSP1 interactions, the high content of DOX (50%) can form a giant network of DOX molecules on the interface of ND via both hydrogen bonding and hydrophobic π-π interactions (Figure 3K). This self-aggregation led to an arrested molecular motion and improved association of DOX with the ND. Therefore, we consider that the design of ND formulation involved free DOX molecules should be primarily concerned with the increased self-aggregation tendency at higher DOX concentrations (Fonseca et al., 1997; Fülöp et al., 2013). Such spontaneous homomolecular interactions may distort the orientational-dependent interactions between DOX and ND assembly.

The Impact of Drug Conjugation on ND Formulation

Three existing types of DOX lipidation are evaluated by a similar MD simulation approach. MD simulation results of ND formulation containing increasing quantities (10, 25, 50%) of BMPH-PTE, BMPH-MP, and BMPH-TC conjugated pH-sensitive DOX prodrugs and their corresponding central structure with the lowest RMSD score in the cluster are illustrated in Figure 4. In the group of ND formulations of DOX-BMPH-PTE (Figure 4A), no DOX leakage was observed during the MD simulations. Visual inspection on the ND complexes showed that DOX-BMPH-PTE molecules were tightly assembled with DPPC NDs. The RMSD profile of each component trajectory remains generally stable after 50 ns of simulations. Lower average RMSD values of DOX-BMPH-PTE and MSP1 are observed with an increase in DOX-BMPH-PTE proportion (Figure 4A). This implies the loss in the degree of ND conformational freedom and plausibly relative to the concentration-dependent DOX self-aggregation and DOX-MSP1 interactions.

FIGURE 4

FIGURE 4. MD simulation analysis of hydrogen bonds formed between DOX prodrug molecules and MSP1 protein in ND complexes. Dynamic measurements of hydrogen bond numbers in ND formulations of (A) DOX-BMPH-PTE (10%), DOX-BMPH-MP (10%), and DOX-BMPH-TC (10%). (B) DOX-BMPH-PTE (25%), DOX-BMPH-MP (25%), and DOX-BMPH-TC (25%). (C) DOX-BMPH-PTE (50%), DOX-BMPH-MP (50%), and DOX-BMPH-TC (50%). Dynamic measurements of hydrogen bonds per drug molecule in ND formulations of (D) DOX-BMPH-PTE. (E) DOX-BMPH-MP. (F) DOX-BMPH-TC. Electrostatic interactions between DOX prodrug molecules and MSP1 protein in ND formulations of (G) DOX-BMPH-PTE. (H) DOX-BMPH-MP. (I) DOX-BMPH-TC.

Similar dynamic behavior was observed in ND formulations of DOX-BMPH-MP and DOX-BMPH-TC. It is worth noting that each ND assembly possessed the same initial surface area, and the packing capacity is dependent on the surface area per lipid molecule. In this regard, the DOX-BMPH-PTE molecule occupies the largest area (63.0 Ų) on the ND among the simulated prodrugs, followed by DOX-BMPH-TC (40 Å²) and DOX-BMPH-MP (30 Å²). Therefore, ND formulated DOX-BMPH-MP appears to have the highest packing capacity. However, a parallel comparison between DOX-BMPH-PTE and DOX-BMPH-MP ND formulations shows that more DOX-BMPH-MP loading quantities do not seem to govern the stability of the MPS1 and ND complex significantly (Figures 4A,B). But the single acyl-chain characteristics of MP may subsequently influence membrane curvature. As the large head groups (DOX) of lipid clusters group together, DOX-BMPH-MPs impose a positive curvature on both sites of the ND membrane. Notably, the MD trajectory of DOX-BMPH-TC ND formulation exhibits the lowest RMSD values and fluctuations (except the case with a drug loading of 25%) than the formerly discussed ND formulations, which can be typically observed in the case with a 50% drug loading ratio (Figure 4C). We speculated that this significant RMSD improvement is primarily attributed to the cholesterol domain of DOX-BMPH-TC, which can stabilize the ND by increasing the order of the lipid packing, enhancing the rigidity, and reducing the fluidity of the membrane. Moreover, DOX-BMPH-TC with fewer rotatable bonds than DOX-BMPH-PTE and DOX-BMPH-MP can also result in reduced structural flexibility and lower RMSD.

The Electrostatic Interactions between DOX-Prodrug molecules and MSP1 Proteins

In MD simulations of DOX-prodrug ND formulations, we observed that the number of hydrogen bonds between DOX-prodrug molecules and MSP1 proteins is proportional to the loading quantities in all types of ND formulations (Figures 5A–C). The DOX-BMPH-PTE ND formulation consistently has the highest average counts of hydrogen bonds (after 50ns of simulation) across different loading quantities (10%, H_n = 16.1; 25%, H_n = 29.5; 50%, H_n = 45.5), followed by a consistent pattern for DOX-BMPH-MP ND formulation (10%, H_n = 13.1; 25%, H_n = 18.2; 50%, H_n = 35.8), and received the lowest average counts for DOX-BMPH-TC ND formulation (10%, H_n = 10.3; 25%, H_n = 11.5; 50%, H_n = 17.8). This regularity difference between the different types of ND formulations could be explained by the fact that the number of rotatable bonds and the structural flexibility of DOX-prodrug molecules are in descending order: DOX-BMPH-PTE > DOX-BMPH-MP > DOX-BMPH-TC (Figure 1). The loss of conformational freedom in DOX-prodrug molecules results in limited DOX-MSP1 interactions, and hence the less number of hydrogen bonds. This implies that increasing in conformational restriction of drug molecules may have a reduced influence on MSP1 dynamics.

FIGURE 5

FIGURE 5. A schematic illustration of the ND formulation designs for free DOX and lipidated DOX-prodrugs.

The fluctuations of hydrogen bonds per-DOX interaction with MSP1 over the simulation time show that hydrogen bonds per-DOX are significantly reduced with increasing DOX quantities (Figures 5D–F), and the minimum limits are achieved when the proportions of DOX rise to above 25%. In accordance with the center clusters of DOX-BMPH-PTE ND, we found abundant electrostatic interactions between DOX prodrugs and MSP1 proteins, including salt bridges and hydrogen bonds (Figures 5G–I). Similar interactions are observed in DOX-BMPH-MP ND central structures, whereas in DOX-BMPH-TC ND central structures, only a few hydrogen bonds were formed with MSP1 proteins. Collectively, Arg, Lys, Asp, and Glu residues of MSP1 proteins have commonly participated in hydrogen bonding and salt-bridge interactions with DOX-prodrug molecules.

Solvent Accessible Surface Area and Stability of DOX-Prodrug Loaded ND Assemblies

Solvent accessible surface area (SASA) of DOX-prodrug loaded ND assemblies could be an important factor in stability studies. The SASA per prodrug molecule of different ND formulations is shown in Figures 6A–C. We noticed that a smaller head moiety of lipidated DOX-prodrugs can lower SASA per prodrug molecule. Besides, comparable trajectories are observed for 10 and 25% of both DOX-BMPH-MP ND and DOX-BMPH-TC ND (Figures 6B,C), which suggests rigidity of the DOX-BMPH-MP and DOX-BMPH-TC head moieties. It is notable that the SASA per prodrug molecule appears to be dropped at high proportions of DOX-prodrug quantity (50%), indicating that the tendency of decreasing in SASA could be governed by the quantity-dependent self-aggregation of DOX (Fülöp et al., 2013). Various interactions such as hydrogen bonding, π-π, π-alkyl stacking, and hydrophobic interaction are frequently involved in the aggregation networks of DOX. The visual inspection of central structures found that the presence of six-membered daunosamine sugar of DOX molecules plays a dominant role in forming hydrogen-bonding networks (Figures 6D–F). Other substructures, including aglycone moiety, phosphate groups, glycerol groups, and maleimide groups, also facilitated the intermolecular and intramolecular electrostatic interactions.

FIGURE 6

FIGURE 6. MD simulation analysis of the solvent accessible surface area (SASA) and self-aggregation behavior of DOX-prodrug ND formulations. Time evolution of the SASA per prodrug molecule during 100 ns of MD simulation for ND formulations of (A) DOX-BMPH-PTE. (B) DOX-BMPH-MP. (C) DOX-BMPH-TC. The aggregation networks of DOX molecules in ND formulations of (D) DOX-BMPH-PTE. (E) DOX-BMPH-MP. (F) DOX-BMPH-TC. Overall Stability and Flexibility of MSP1 Protein Assembled with Different DOX-ND Formulations.

To observe and compare the motion of MSP1 residues, we performed the root mean square fluctuation (RMSF) analyses for both the empty ND (reference structure) and different DOX-ND Formulations (Figure 7). In all cases, the most flexible regions are the C-terminal residues of the MSP1 (Figures 7A–D). As the DOX loading increased, each type of formulation followed a similar tendency with varying degrees of reduction in global RMSF values (Figure 7E). It is observed that residues of the ND formulation of DOX-BMPH-TC are more rigid than other ND formulations. This is in agreement with RMSD results and indicates that there were relatively strong interactions between the lipids/Drugs and MSP1 proteins, as well as lower fluidity of the membrane, which stabilized the DOX-BMPH-TC ND in a more compact conformation than others. Nevertheless, the rigidity of residues could be improved by loading a higher proportion of lipidated-prodrugs except for free DOX ND formulation.

FIGURE 7

FIGURE 7. Root Mean Square Fluctuation (RMSF) analyses MSP1 protein assembled with different DOX-ND formulations. Each MSP1 contained 200 residues, and the abscissa in the line charts represents the average value of RMSF in the same residue of two MSP1 proteins. RMSF of four systems compared with MSP1 residues in the empty ND: (A) MSP1 residues in the free DOX-ND formulations. (B) MSP1 residues in DOX-BMPH-PTE ND formulations. (C) MSP1 residues in DOX-BMPH-MP ND formulations. (D) MSP1 residues in DOX-BMPH-TC ND formulations. (E) Comparison of the fluctuations of each residue between four ND formulations and the empty ND with variants colored and sized according to its RMSF contribution. The coloring-thickness scale varies from deep red and thin (low fluctuations) to dark blue and thick (high fluctuations).

Conclusion

The molecular dynamic behavior of both free DOX and DOX conjugated lipid prodrug molecules loaded in ND formulations was investigated in depth by a series of all-atom MD simulations. We found that the interactions between free DOX molecules and the ND system are labile due to the strong electrostatic interaction between the amino sugar moiety of DOX and the phosphate head group of DPPC. These hydrogen bonds and intrinsic steric hindrance altered the ideal binding mode of DOX and resulted in an unfavorable orientation of DOX attached on the lipid surface rather than the vertical insertion of the AQ ring. Such amphiphilic characteristics of DOX may enable rapid membrane penetration and diffusion (Toroz and Gould, 2019). The stability of DOX binding modes to the ND membrane could be dependent on lipid composition (Siani et al., 2022). Through the comparison between three types of ND formulations of lipidated-DOX prodrugs, we observed that lipidation approaches could provide sufficient conformational stability for prodrug loading and different lipid used for lipidation led to variable loading capacity. ND formulation of DOX-BMPH-MP may be selected for more drug loading quantities as compared to other formulations. Our results also demonstrated that the global conformational stability of DOX-prodrug-loaded ND Assemblies is influenced by the molecular flexibility and lipid composition of DOX-prodrug molecules. Thus, ND formulation of DOX-BMPH-TC may be considered for better dynamic stability. We have also demonstrated that the increase in proportions of DOX-prodrug within a specific range (e.g., 10–50%) could induce DOX self-aggregation and hydrogen bonding with MSP1 residues, incidentally lowering the membrane fluidity and hence improving the conformational stability of ND formulations.

The main limitation of our study lies in the assumption that ND formulations of DOX satisfy the discoidal double-belt model and can be properly prepared under experimental conditions. Although this study mainly focused on comparing ND formulations with one another, simulated time could be another limitation that affects the accuracy. Besides, the effect of drug loading quantity on ND solubility was not within the scope of the study. Nevertheless, our findings offered a molecular insight into how lipidation of DOX influences the dynamic performance of ND formulations. The obtained information may guide the practical design and optimization of the ND formulation of DOX. The in silico approach adopted here exhibits future potential as a virtual screening or evaluation method to develop ND formulations.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

JW conceived the idea, designed the study, directed the project. JW, XC, and ZHC developed the ND formulation dynamic simulation method and wrote the manuscript. XC, ZHC, and DX performed all the experiments and analyzed the data. MW, ZHW, XC and ZHC evaluated the computational modeling for formulation design and supervised the experimental procedure. DX, YLV, YLI, SL, WX, YM, ZIC, and TC provided technical support and revised the manuscript. All authors contributed text to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (82001887 to JW and 81900105 to ZW), in part by a horizontal project of Sun Yat-sen University (K21-75110-007), Science and Technology Innovation Committee of Shenzhen Municipality (JCYJ20210324115003009), and the Guangdong Basic and Applied Basic Research Foundation (2019A1515110326).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor ML declared a shared affiliation, though no other collaboration, with all authors at the time of the review.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Boruah, N., Sarma, H., and Sharma, H. (2019). Computational Formulation Development and Drug Delivery. Tamil Nadu: i-manager Publications, 191–195.

Google Scholar

Bunker, A., and Róg, T. (2020). Mechanistic Understanding from Molecular Dynamics Simulation in Pharmaceutical Research 1: Drug Delivery. Front. Mol. Biosci. 7, 604770. doi:10.3389/fmolb.2020.604770

PubMed Abstract | CrossRef Full Text | Google Scholar

Chatterjee, K., Zhang, J., Honbo, N., and Karliner, J. S. (2010). Doxorubicin Cardiomyopathy. Cardiology 115 (2), 155–162. doi:10.1159/000265166

PubMed Abstract | CrossRef Full Text | Google Scholar

Coldwell, K. E., Cutts, S. M., Ognibene, T. J., Henderson, P. T., and Phillips, D. R. (2008). Detection of Adriamycin-DNA Adducts by Accelerator Mass Spectrometry at Clinically Relevant Adriamycin Concentrations. Nucleic Acids Res. 36 (16), e100. doi:10.1093/nar/gkn439

PubMed Abstract | CrossRef Full Text | Google Scholar

David Foglesong, P., Reckord, C., and Swink, S. (1992). Doxorubicin Inhibits Human DNA Topoisomerase I. Cancer Chemother. Pharmacol. 30 (2), 123–125. doi:10.1007/BF00686403

PubMed Abstract | CrossRef Full Text | Google Scholar

Davies, K. J., and Doroshow, J. H. (1986). Redox Cycling of Anthracyclines by Cardiac Mitochondria. I. Anthracycline Radical Formation by NADH Dehydrogenase. J. Biol. Chem. 261 (7), 3060–3067. doi:10.1016/s0021-9258(17)35746-0

CrossRef Full Text | Google Scholar

Fonseca, M. J., van Winden, E. C. A., and Crommelin, D. J. A. (1997). Doxorubicin Induces Aggregation of Small Negatively Charged Liposomes. Eur. J. Pharmaceutics Biopharmaceutics 43 (1), 9–17. doi:10.1016/S0939-6411(96)00018-5

CrossRef Full Text | Google Scholar

Fülöp, Z., Gref, R., and Loftsson, T. (2013). A Permeation Method for Detection of Self-Aggregation of Doxorubicin in Aqueous Environment. Int. J. Pharmaceutics 454 (1), 559–561. doi:10.1016/j.ijpharm.2013.06.058

CrossRef Full Text | Google Scholar

Hess, B., Bekker, H., Berendsen, H. J. C., and Fraaije, J. G. E. M. (1997). LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 18 (12), 1463–1472. doi:10.1002/(sici)1096-987x(199709)18:12<1463::aid-jcc4>3.0.co;2-h

CrossRef Full Text | Google Scholar

Huang, J., and MacKerell, A. D. (2013). CHARMM36 All-Atom Additive Protein Force Field: Validation Based on Comparison to NMR Data. J. Comput. Chem. 34 (25), 2135–2145. doi:10.1002/jcc.23354

CrossRef Full Text | Google Scholar

Ji, C., Yang, B., Yang, Y.-L., He, S.-H., Miao, D.-S., He, L., et al. (2010). Exogenous Cell-Permeable C6 Ceramide Sensitizes Multiple Cancer Cell Lines to Doxorubicin-Induced Apoptosis by Promoting AMPK Activation and mTORC1 Inhibition. Oncogene 29 (50), 6557–6568. doi:10.1038/onc.2010.379

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadiyala, P., Li, D., Nuñez, F. M., Altshuler, D., Doherty, R., Kuai, R., et al. (2019). High-Density Lipoprotein-Mimicking Nanodiscs for Chemo-Immunotherapy against Glioblastoma Multiforme. ACS Nano 13 (2), 1365–1384. doi:10.1021/acsnano.8b06842

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuai, R., Ochyl, L. J., Bahjat, K. S., Schwendeman, A., and Moon, J. J. (2017). Designer Vaccine Nanodiscs for Personalized Cancer Immunotherapy. Nat. Mater 16 (4), 489–496. doi:10.1038/nmat4822

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuai, R., Yuan, W., Son, S., Nam, J., Xu, Y., Fan, Y., et al. (2018). Elimination of Established Tumors with Nanodisc-Based Combination Chemoimmunotherapy. Sci. Adv. 4 (4), eaao1736. doi:10.1126/sciadv.aao1736

PubMed Abstract | CrossRef Full Text | Google Scholar

Lacko, A. G., Sabnis, N. A., Nagarajan, B., and McConathy, W. J. (2015). HDL as a Drug and Nucleic Acid Delivery Vehicle. Front. Pharmacol. 6. doi:10.3389/fphar.2015.00247

PubMed Abstract | CrossRef Full Text | Google Scholar

Marinello, J., Delcuratolo, M., and Capranico, G. (2018). Anthracyclines as Topoisomerase II Poisons: From Early Studies to New Perspectives. Ijms 19 (11), 3480. doi:10.3390/ijms19113480

PubMed Abstract | CrossRef Full Text | Google Scholar

Martyna, G. J., Klein, M. L., and Tuckerman, M. (1992). Nosé-Hoover Chains: The Canonical Ensemble via Continuous Dynamics. J. Chem. Phys. 97 (4), 2635–2643. doi:10.1063/1.463940

CrossRef Full Text | Google Scholar

Mattarollo, S. R., Loi, S., Duret, H., Ma, Y., Zitvogel, L., and Smyth, M. J. (2011). Pivotal Role of Innate and Adaptive Immunity in Anthracycline Chemotherapy of Established Tumors. Cancer Res. 71 (14), 4809–4820. doi:10.1158/0008-5472.CAN-11-0753

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehta, C. H., Narayan, R., and Nayak, U. Y. (2019). Computational Modeling for Formulation Design. Drug Discov. Today 24 (3), 781–788. doi:10.1016/j.drudis.2018.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Momparler, R. L., Karon, M., Siegel, S. E., and Avila, F. (1976). Effect of Adriamycin on DNA, RNA, and Protein Synthesis in Cell-free Systems and Intact Cells. Cancer Res. 36 (8), 2891–2895.

PubMed Abstract | Google Scholar

Murakami, T., Wijagkanalan, W., Hashida, M., and Tsuchida, K. (2010). Intracellular Drug Delivery by Genetically Engineered High-Density Lipoprotein Nanoparticles. Nanomedicine 5 (6), 867–879. doi:10.2217/nnm.10.66

PubMed Abstract | CrossRef Full Text | Google Scholar

Octavia, Y., Tocchetti, C. G., Gabrielson, K. L., Janssens, S., Crijns, H. J., and Moens, A. L. (2012). Doxorubicin-induced Cardiomyopathy: From Molecular Mechanisms to Therapeutic Strategies. J. Mol. Cell Cardiol. 52 (6), 1213–1225. doi:10.1016/j.yjmcc.2012.03.006

CrossRef Full Text | Google Scholar

Qi, Y., Lee, J., Klauda, J. B., and Im, W. (2019). CHARMM‐GUI Nanodisc Builder for Modeling and Simulation of Various Nanodisc Systems. J. Comput. Chem. 40 (7), 893–899. doi:10.1002/jcc.25773

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramezanpour, M., Leung, S. S. W., Delgado-Magnero, K. H., Bashe, B. Y. M., Thewalt, J., and Tieleman, D. P. (2016). Computational and Experimental Approaches for Investigating Nanoparticle-Based Drug Delivery Systems. Biochim. Biophys. Acta (Bba) - Biomembranes 1858 (7, Part B), 1688–1709. doi:10.1016/j.bbamem.2016.02.028

CrossRef Full Text | Google Scholar

Siani, P., Donadoni, E., Ferraro, L., Re, F., and Di Valentin, C. (2022). Molecular Dynamics Simulations of Doxorubicin in Sphingomyelin-Based Lipid Membranes. Biochim. Biophys. Acta (Bba) - Biomembranes 1864 (1), 183763. doi:10.1016/j.bbamem.2021.183763

CrossRef Full Text | Google Scholar

Tacar, O., Sriamornsak, P., and Dass, C. R. (2012). Doxorubicin: an Update on Anticancer Molecular Action, Toxicity and Novel Drug Delivery Systems. J. Pharm. Pharmacol. 65 (2), 157–170. doi:10.1111/j.2042-7158.2012.01567.x

CrossRef Full Text | Google Scholar

Tang, Z., Xiao, Y., Kong, N., Liu, C., Chen, W., Huang, X., et al. (2021). Nano-bio Interfaces Effect of Two-Dimensional Nanomaterials and Their Applications in Cancer Immunotherapy. Acta Pharmaceutica Sinica B 11 (11), 3447–3464. doi:10.1016/j.apsb.2021.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Toroz, D., and Gould, I. R. (2019). A Computational Study of Anthracyclines Interacting with Lipid Bilayers: Correlation of Membrane Insertion Rates, Orientation Effects and Localisation with Cytotoxicity. Sci. Rep. 9 (1), 2155. doi:10.1038/s41598-019-39411-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Der Spoel, D., Lindahl, E., Hess, B., Groenhof, G., Mark, A. E., and Berendsen, H. J. C. (2005). GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 26 (16), 1701–1718. doi:10.1002/jcc.20291

CrossRef Full Text | Google Scholar

Velez, J. M., Miriyala, S., Nithipongvanitch, R., Noel, T., Plabplueng, C. D., Oberley, T., et al. (2011). p53 Regulates Oxidative Stress-Mediated Retrograde Signaling: a Novel Mechanism for Chemotherapy-Induced Cardiac Injury. PLoS One 6 (3), e18005. doi:10.1371/journal.pone.0018005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, B., Yuan, Y., Han, L., Ye, L., Shi, X., and Feng, M. (2014). Recombinant Lipoproteins Reinforce Cytotoxicity of Doxorubicin to Hepatocellular Carcinoma. J. Drug Target. 22 (1), 76–85. doi:10.3109/1061186X.2013.839687

CrossRef Full Text | Google Scholar

Wang, J., Wang, A. Z., Lv, P., Tao, W., and Liu, G. (2018). Advancing the Pharmaceutical Potential of Bioinorganic Hybrid Lipid-Based Assemblies. Adv. Sci. 5 (9), 1800564. doi:10.1002/advs.201800564

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, Y., Wang, W., Wang, B., Zhu, H., Zhang, B., and Feng, M. (2013). Delivery of Hydrophilic Drug Doxorubicin Hydrochloride-Targeted Liver Using apoAI as Carrier. J. Drug Target. 21 (4), 367–374. doi:10.3109/1061186X.2012.757769

CrossRef Full Text | Google Scholar