Utilizing Molecular Simulation to Examine Nanosuspension Stability

1. Introduction

Previous estimates suggest that approximately 40% of active pharmaceutical ingredients are poorly soluble in water. [1] To improve the bioavailability of such compounds, many formulations utilize drug nanosuspensions. Nanosuspensions possess a high ratio of surface area to volume, which greatly aids drug dissolution and increases saturation solubility. [2,3,4,5,6,7,8,9,10,11] Further, they can be formed with minimal amounts of potentially toxic solvents and can be delivered through oral, [12] dermal, [13] pulmonary, [14] or intravenous administration, [15] which demonstrates both safety and convenience. [16]

Such nanosuspensions can be formed through two classes of methods, broadly considered to be top-down methods and bottom-up methods. Top-down methods utilize wet milling to produce nanoparticles from larger drug crystals, [17,18,19,20] while bottom-up methods precipitate nanocrystals out of a drug solution. [21,22,23] In both cases, it is often necessary to stabilize the nanosuspension through the use of excipients. These excipients adsorb to the drug crystal and prevent aggregation or Ostwald ripening, which can degrade nanosuspension stability. [24] Drug-excipient interactions are highly specific, and different excipients work better for different drugs; however, discovering the optimal excipient is often done through experimental screening. [25,26,27,28,29,30,31] Gaining insight into the interactions that drive excipient adsorption to drug nanoparticles could help rationalize formulation design.

Computational modeling offers a promising approach for examining the interactions that dictate nanosuspension stability. [32,33,34] Electronic structure methods such as density functional theory can be used to most accurately evaluate drug-nanoparticle interactions. [35,36] At the cost of some accuracy relative to electronic structure methods, molecular dynamics (MD) simulations can study larger systems and explore longer timescales, making them advantageous in a variety of drug nanosuspension studies. Indeed, MD simulations have been used to calculate rate of dissolution of different drug crystal polymorphs, [37] determine the drug-excipient binding free energy, [38] and to examine the structural features of excipient adsorption to a drug. [39,40,41,42,43,44] The results of these studies suggest that computational modeling may provide crucial insight into how excipient-drug interactions dictate nanosuspension stability.

The aim of this work was to utilize MD simulations to elucidate the structural features that dictate nanosuspension stability. Previous work from Ferrar et al. synthesized drug nanoparticles via a top down milling technique, characterized their stability through dynamic light scattering experiments, and noted that excipient amphiphilicity seemed important for nanosuspension stabilization. [39] To rationalize these previous expeirmental results, we created a computational framework to systematically investigate the same drug-excipient combinations. We examined the molecular features that drive excipient adsorption to drug nanoparticles. Next, we explored the correlation between experimental nanosuspension stability and computational properties for a variety of drug-excipient combinations. Our results suggested that the polar fraction of solvent-accessible surface area dictates nanosuspension stability. To test this model, we then applied it to a system not previously studied, GDC-0810, which is a former breast cancer candidate compound, and can be crystallized in both a neutral and salt form. Our model predicted that the salt form would have a higher fraction of polar surface area, and, thus, more easily form stable nanosuspensions, which we confirmed by creating nanosuspensions of GDC-0810 salt through flash nano precipitation (FNP). We further explained these computational results by demonstrating that our predictions can be made without simulation by using only the hydrogen bonding and charge properties of each the drug and the excipient. Overall, our work connects a microscopic quantity, polar fraction of solvent-accessible surface area, to a macroscopic quantity, nanosuspension stability, and thus creates a pathway for rational formulation design.

2. Methods

2.1. Molecular dynamics simulation

To examine nanosuspension stability, all-atom molecular dynamics simulations were performed on a training set of three different drug molecules (naproxen, indomethacin, and itraconazole) with six different excipients (sodium dodecyl sulfate (SDS), sodium octyl sulfate (SOS), sodium deoxycholate (SDC), polyethylene glycol (PEG), polypropylene glycol (PPG), and poloxamer), as characterized by previous experiments, reported in Ferrar et al. and in Table S1–S3. [39] After characterizing our training set, we then performed simulations with two different crystal forms of GDC-0810, [45] one of which is neutral and one of which is ionic. The GDC-0810 test simulations were performed with all six excipients that were utilized in the training set and Tween 80, for a total of seven excipients. Each simulations began by identifying the crystal structure of the drug. Naproxen (ID = COYRUD), indomethacin (ID = INDMET), and itraconazole (ID = TEHZIP) were extracted from the Camridge Structural Database, [46] while the GDC-0810 crystal structure was determined experimentally in a neutral (ID = gened) and ionic (ID = genea) form. [45] Each crystal was then replicated to a cube of approximately 35 Å per side using the periodic boundary conditions builder function in MOE 2020.9. [47] Next, GROMACS 2020.4 editconf was used to center the crystal in a square simulation box with side lengths of 110 Å. Excipients were added randomly in solution to the simulation box using GROMACS insert-molecules. [48] The number of excipients added was chosen to match experiments:drug (wt%) of 10%, and the exact number of drugs and excipients used in our simulations are provided for the training set (Table S4) and test set (Table S5).

Systems built by GROMACS were then loaded into Schrödinger 2020.4 Maestro. [49] Simulations were prepared using Desmond Simulation Builder, which solvated the simulation box with side length of 110 Å with simple point charge (SPC) water, added sodium ions when they were necessary to balance the charge of excpients, and prepared the simulation to run in the OPLS3e force field. [50] Each setup was relaxed using Desmond Minimization protocol, which performed 100 ps of Brownian dynamics at 10K with a 1fs timestep. We then ran production Desmond Molecular Dynamics simulations [51] under constant number, particle, and temperature (NPT ensemble) for 200 ns at 300 K and 1.01325 bar. Restraints were placed on the drug crystal, with a force constant of 5 kcal ${\mathrm{mol}}^{-1}$ ${\AA }^{-2}$. Data was saved every 100 ps, and the first 100 ns were excluded from all analysis for equilibration.

2.2. Simulation analysis

Analysis of each simulation began by calculating a molecular contact matrix using MDAnalysis 1.1.1. [52,53] Contacts between two molecules, either drug or excipient, were included if any heavy atoms were within 4.5Å. This contact matrix then served three purposes: (1) we computed the atomistic contact maps, or frequencies that excipient atoms and drug atoms were in contact, (2) we computed the probability distribution of whole molecule contacts, and (3) we output PDB files of the drug crystal and all excipients or counterions that are in contact with the drug crystal, which thus excludes water and any excipient or counterion not in contact with the crystal. We then input these PDB files into MOE 2020.09 to calculate the surface properties of our excipient-drug complex, including the fraction of polar surface area (${\mathrm{FASA}}_{\mathrm{p}}$), accessible surface area (ASA), radius of gyration (${\mathrm{R}}_{\mathrm{g}}$) and dipole. Fitting to linear regression and least absolute shrinkage and selection operator (LASSO) were performed using MATLAB 2019b. [54] After calculating each initial linear regression model using all data available, we predicted its accuracy using k-fold cross validation. We repeated the k-fold division of our data 100 times, and report the average and standard deviation of the classification accuracy for the data excluded from model fitting. Molecular properties of drugs and excipients were calculated using sdfCalcProps.csh, a script available through the chemalot open source package. [55]

2.3. Flash nano precipitation of GDC-0810

Nanoparticles with a diameter of 225 nm were prepared through FNP with a multi-inlet vortex mixer (MIVM) as previously described. [56] GDC-0810 was dissolved at 15mg/mL in 0.4 eq NaOH and 1mL was added to a 1mL syringe. Equal volume of water was added to 3 additional 1mL syringes. Syringes were placed on the FNP and the plate was compressed at an equal rate across the syringes. After passing through the MIVM, the nanoparticles entered a dilution chamber of 14 mL of water.

Nanoparticles were characterized by dynamic light scattering (Wyatt DynaPro Plate Reader III). Samples were diluted 10 fold then 30 µL was placed in a 384 black wall, clear bottom polystyrene microplate (Corning). Particle diameter by intensity was reported from 15 acquisitions.

3. Results

3.1. Hydrophobic interactions drive excipient adsorption to drug crystals

While nanosuspensions present a promising approach to improve the solubility and dissolution of active pharmaceutical ingredients, choosing the correct formulation to stabilize a chosen drug remains a difficult task. Traditional approaches have focused on more experimental characterization methods, but recent efforts have shown computational modeling may aid in the design of stable nanosuspensions. [32,33,34] To continue in this direction, we applied molecular dynamics simulation to quantify drug-excipient interactions and find underlying physical characteristics that can predict nanosuspension stability.

We started by characterizing a set of systems studied previously by Ferrar et al. [39] These previous experiments had characterized nanoparticle size for three drug molecules (naproxen, indomethacin, and itraconazole) with six different excipients (sodium dodecyl sulfate (SDS), sodium octyl sulfate (SOS), sodium deoxycholate (SDC), polyethylene glycol (PEG), polypropylene glycol (PPG), and poloxamer). For each of these 18 possible drug-excipient combinations, we performed molecular dynamics simulations. For each simulation, a drug crystal was placed at the center of an aqueous simulation box, with excipient molecules dispersed randomly throughout the aqueous phase. Over time, the excipients adsorb to the drug crystal (Figure 1), with differences being observed for each drug-excipient pair (Figure S1).

To quantify the specific interactions between drugs and excipients, we computed drug-excipient contact maps for naproxen (Figure 2), indomethacin (Figure S2), and itraconazole (Figure S3). These atomistic contact maps give insight into the parts of the drug and excipient that are most likely and least likely to form contacts between the drug and excipient. In all cases, we find that excipient-drug interactions are driven by hydrophobicity. As a result, drugs are contacted through exposed hydrophobic areas, such as the aromatic rings on naproxen and indomethacin or the terminal alkane group of itraconazole. Further, evidence of hydrophobicity driven interactions can also be seen on the excipients. Alkane groups of SDS and SOS are frequently in contact with the drug particle, but the sulfate group is exposed to water. Similarly, SDC binds primarily through cyclohexane groups, but exposes alcohol and carboxylic acid groups. Finally, polymers such as PEG, PPG, and poloxamer interact preferentially through carbon atoms, while oxygen atoms are less likely to form contacts with the drug crystal. Our simulations also show that the middle of the polymer has less contacts with the crystal than the ends of the chain, likely as a result of self-association or interactions with water. Viewed in total, these contact maps show that excipients interact with drugs through hydrophobic regions, and leave more polar portions of the excipient or drug exposed.

While atomistic pictures of these contact probabilities provides an interesting example of drug-excipient interactions, we can also quantify interactions between the excipient and drug nanocrystal in terms of molecule to molecule contacts, or coordination number. These molecular contacts describe how many neighboring molecules of a given type a reference molecule may have. By examining the frequency of the number of drug molecules per excipient (Figure 3), the number of drug molecules per drug (Figure S4), the number of excipient molecules per drug (Figure S5), and the number of excipient molecules per excipient (Figure S6), we can better understand the ways in which drugs and excipients interact. We divide the data into surfactants (SDS, SOS, and SDC) and polymers (PEG, PPG, and poloxamer) because the polymers have a larger molecular weight and, therefore, have a higher number of drug molecules per excipient.

The number of drug molecules per surfactant reveals trends in the binding modes of surfactants with drug molecules (Figure 3A,C,E). In particular, SDC generally has a lower average coordination number than SOS or SDS. Previous simulations observed that the alkyl chains of SOS and SDS can penetrate the surface of the drug crystal, while SDC tends to lie on the surface of the crystal. [39] This mechanism of adsorption likely results in less drugs bound to SDC than SOS or SDS. Further, SOS and SDC frequently have 0 drug molecules per surfactant. This observation indicates that SOS and SDC have a lower affinity for the drug molecule than the more hydrophobic SDS, which results in SOS and SDC molecules that do not adsorb to the drug crystal.

The number of drug molecules per polymer also demonstrates the role of excipient hydrophobicity in adsorption. In comparison to the other polymers, PPG tends to have a lower number of drug molecules per polymer (Figure 3B,D,F) and higher number of polymer molecules per polymer (Figure S6). This indicates that being partially soluble may be important for a polymer to prevent excipient self-association and allow the polymer to interact with a higher fraction of the drug crystal. Taken together with the surfactants, these trends demonstrate that the interactions between drug and excipient are intricately controlled by balancing the excipient’s affinity for drug, the excipient itself, and the aqueous phase.

3.2. Atomistic modeling of nanosuspension stability

After examining how excipients interact with drug molecules, we wondered if we could relate our findings back to experimental nanosuspension stability by analyzing the excipient-drug complex. To accomplish this task, we isolated the drug crystal and all excipients and ions in contact with the crystal, and calculated the accessible surface area (ASA), radius of gyration (${\mathrm{R}}_{\mathrm{g}}$), fraction of the total surface area that is polar (${\mathrm{FASA}}_{\mathrm{p}}$) and dipole moment for each drug-excipient pair (Figure 4A-C, Figure S7A). We normalized ASA and ${\mathrm{R}}_{\mathrm{g}}$ by that of the original crystal to remove effects arising from the differences in size of each drug crystal.

With various simulated descriptors in hand, we then aimed to compare our results to experimental measures of nanosuspension stability. As an ideal formulation will minimize the ratio of excipient:drug, we used the minimum weight percentage of excipient needed to form a stable nanosuspension (min wt%) as our target metric, and performed additional experiments to complement the data set from Ferrar et al. (Table S1). [39] Of the 16 test drug-excipient combinations, 8 formed stable nanosuspensions, for which we recorded the min wt% (Figure 4D, Table S2). Next, we calculated the Pearson correlation coefficient between min wt% and ${\mathrm{FASA}}_{\mathrm{p}}$ (Figure 4E), ${\mathrm{ASA/ASA}}^0$ (Figure S7B), ${\mathrm{R}}_{\mathrm{g}}$/${\mathrm{R}}_{\mathrm{g}}^0$ (Figure S7C), and dipole moment (Figure S7D). We found only ${\mathrm{FASA}}_{\mathrm{p}}$ is highly correlated to the min wt%. Further, ${\mathrm{FASA}}_{\mathrm{p}}$ is inversely correlated with min wt%, indicating a higher polar fraction of solvent-accessible surface area from our simulation will experimentally form a stable nanosuspension at lower excipient:drug ratios. This result is easily rationalized, for more polar particles are more soluble in water, and thus less likely to undergo destabilization from aggregation or Ostwald ripening. [24]

While our modeling supports the hypothesis that a higher polar fraction of solvent-accessible surface area results in a more stable nanosuspension at lower concentrations of excipients, we still lacked a quantitative way to predict nanosuspension stability. To make our predictions more quantitative, we simplified our data by dividing it into two groups. In the first group we label as “successes", milling at certain excipient:drug ratios results in a stable nanosuspension. In the other group which we label as “failures", milling does not result in a stable nanosuspension, regardless of the excipient:drug ratio. Each group consists of eight drug-excipient combinations (Table S3). We then used a logistic regression to calculate the probability of getting a “success" based on the ${\mathrm{FASA}}_{\mathrm{p}}$, and found

$$log\left(\frac{P_f}{P_s}\right)=6.14-26.1\ast {\mathrm{FASA}}_{\mathrm{p}},$$

(1)

where $P_f$ is the probability of not successfully forming a nanosuspension, and $P_s$ is the probability of successfully forming a nanosuspension (Figure 4F, n=16). To examine the accuracy of predictions from this logistic regression model, we recomputed our model using 100 independent rounds of 4-fold cross validation and found a classification accuracy of $0.74\pm 0.21$, indicating that our model should be applicable to data beyond the training set. Our model has a negative coefficient in front of ${\mathrm{FASA}}_{\mathrm{p}}$, again suggesting that a higher polar fraction of solvent-accessible surface area is more likely to successfully form a stable nanosuspension.

4. Discussion

4.1. Simulated polar fraction of solvent-accessible surface area successfully predicts experimentally determined nanosuspension stability

Comparisons between our MD simulations and experiments on naproxen, indomethacin, and itraconazole suggest that a higher polar fraction of solvent-accessible surface area results in a more stable nanosuspension. To test this hypothesis, we studied compound GDC-0810, a former breast cancer candidate compound. GDC-0810 was selected as it can be crystalized in a neutral and ionic form (Figure 5A,B). [45] As we hypothesized an ionic crystal would be more polar than a neutral one, GDC-0810 offered an opportunity to examine the effect of polarizability on effect nanoparticle stability within different states of the same compound. For each the neutral and the ionic crystal, we performed simulations with the same six excipients used in our training data, plus Tween 80 (Figure 5C), resulting in 14 new simulations.

We immediately observed that the binding mode of the excipient varies greatly by the crystallization state of the drug. The atomistic contact maps find that SDS adsorbs to different regions of the neutral crystal than those of the ionic crystal (Figure 5D-F). This trend is seen with all excipient simulations (Figure S9,S10). Similarly, molecular contacts vary greatly with the charge state of the drug. Most significantly, nearly all excipients molecules remain bound to the neutral crystal throughout the entire simulation, resulting in a very low probability of having 0 drug molecules per excipient. However, in the ionic crystal, a significant number of excipients are in solution, as indicated by a much higher frequency at 0 drug molecules per excipient (Figure 5G-J). These differences suggest that most excipients have a lower affinity for the ionic crystal than the neutral crystal. Similar trends for the drug molecules per drug (Figure S11), the drug molecules per excipient (Figure S12), and the excipient molecules per excipient (Figure S13) demonstrate that the ionic crystal and neutral crystal have different structural features.

Next, we aimed to connect knowledge of property differences within our training set to the example of GDC-0810. We first calculated the surface properties for both the ionic and neutral crystal (Figure 6, Figure S14). ${\mathrm{FASA}}_{\mathrm{p}}$ is much higher for the ionic crystal than the neutral crystal, indicating that the ionic crystal is more likely to form stable nanosuspensions. To quantitatively predict nanoparticle stability, we inserted the ${\mathrm{FASA}}_{\mathrm{p}}$ into Equation (1). Our model predicts that the ionic crystal was likely to form stable nanosuspensions under any conditions, while the neutral crystal may form a stable nanosuspension with more polar excipients like SOS or SDS (Table 1).

Based on our MD simulations with different excipients, we hypothesized that ionic GDC-0810 may form nanosuspensions without excipients. Indeed, predictions of nanoparticle stability using the ${\mathrm{FASA}}_{\mathrm{p}}$ of the crystal without excipients (Figure 6) suggested that the GDC-0810 salt (${\mathrm{P}}_{\mathrm{s}}$=0.871), but not the protonated form (${\mathrm{P}}_{\mathrm{s}}$=0.161) is likely to form a stable nanosuspension. To evaluate this prediction, we utilized flash nano precipitation [56] to create nanosuspensions of GDC-0810 salt. Even without excipients, we found that the GDC-0810 salt formed nanosuspensions, which we characterized with dynamic light scattering. Our analysis revealed that the nanoparticles were stable, with an average size of 244.6 nm and % PD (percent polydispersity) of 26.8% upon formation. Further, the nanoparticles remained stable and had an average size of 225 nm and % PD of 49.7% after incubation for over 3 years (Figure 7). These results suggest that the salt nanosuspension would be a good choice of formulation for this poorly soluble compound.

4.2. Leveraging two-dimensional properties to predict nanosuspension stability

The results with GDC-0810 support our previous theory that a higher polar fraction of solvent-accessible surface area results in a more stable nanosuspension. In principle, one could run molecular dynamics simulation and test various combinations to discover the optimal excipient for a given drug. However, such simulations are computationally expensive, which motivated us to find underlying two-dimensional properties of drugs and excipients that can predict the simulated values of ${\mathrm{FASA}}_{\mathrm{p}}$.

We then hypothesized that hydrogen bond acceptors, hydrogen bond donors, and charged molecules contribute most to ${\mathrm{FASA}}_{\mathrm{p}}$. Using these properties for both the drug (Table S6) and excipient (Table S7), we used the least absolute shrinkage and selection operator (LASSO) to fit a linear regression model for ${\mathrm{FASA}}_{\mathrm{p}}$ with all systems we studied (Figure 4C and Figure 6). LASSO is a technique that performs both variable selection and regression, improving the accuracy and interpretability of the resulting regression model (Figure S15). The model that minimizes the mean squared error (n=32) takes the expression

$${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}=13.9\frac{H{A}_{\mathrm{drug}}}{M{W}_{\mathrm{drug}}}+27.0\frac{H{D}_{\mathrm{drug}}}{M{W}_{\mathrm{drug}}}+2.3\frac{H{A}_{\mathrm{exc}}}{M{W}_{\mathrm{exc}}}+4.3\frac{H{D}_{\mathrm{exc}}}{M{W}_{\mathrm{exc}}}+118\frac{|q_{\mathrm{drug}}|}{M{W}_{\mathrm{drug}}}+10.6\frac{|q_{\mathrm{exc}}|}{M{W}_{\mathrm{exc}}}-0.108,$$

(2)

where ${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}$ is the predicted value of ${\mathrm{FASA}}_{\mathrm{p}}$, $M{W}_{\mathrm{drug}}$ is the molecular weight of the drug molecule, $H{A}_{\mathrm{drug}}$ is the number of hydrogen bond acceptors in the drug molecule, $H{D}_{\mathrm{drug}}$ is the number of hydrogen bond donors in the drug molecule, $q_{\mathrm{drug}}$ is the charge of each drug molecule, $M{W}_{\mathrm{exc}}$ is the molecular weight of the excipient molecule, $H{A}_{\mathrm{exc}}$ is the number of hydrogen bond acceptors in the excipient molecule, $H{D}_{\mathrm{exc}}$ is the number of hydrogen bond donors in the excipient molecule, and $q_{\mathrm{exc}}$ is the charge of each excipient molecule. This model fits the calculated properties from simulations well, resulting in a high correlation between the 2D property-predicted and simulated ${\mathrm{FASA}}_{\mathrm{p}}$ (Figure 8). Further, hydrogen bond acceptors, hydrogen bond donors, and charge all contribute to a higher value of ${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}$. Interestingly, we also observe that hydrogen bond donors contribute more than hydrogen bond acceptors, a common trend in drug discovery settings. [57] This model suggests that excipients with more hydrogen bond donors, hydrogen bond acceptors, and net charge will result in more stable nanosuspensions. Finally, the charge state of drug molecules contributes significantly to ${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}$, which suggests that salt forms of drugs are more likely to form stable nanosuspensions.

To test whether our model of two-dimensional properties is capable of predicting nanoparticle stability, we applied the same tests on our ${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}$ that we applied on ${\mathrm{FASA}}_{\mathrm{p}}$. First, we examined the correlation between ${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}$ and the min wt% (Figure 9A). Unlike the simulated values of ${\mathrm{FASA}}_{\mathrm{p}}$, the ${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}$ shows little correlation with the min wt%. This deviation could originate from a variety of factors, such the binding mode of a given excipient or crystal size and shape.

Despite our modeled ${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}$ not correlating with min wt%, it does show some success at predicting nanosuspension stability. Following our procedure for ${\mathrm{FASA}}_{\mathrm{p}}$, we labeled a drug-excipient combination as a “success" if it formed a stable nanosuspension upon milling at certain excipient:drug ratios and labeled it as a “failure" if it did not form stable nanosuspensions upon milling at any excipient:drug ratios. For this binary data, we performed a logistic regression to calculate the probability of getting a “success" based on the ${\mathrm{FASA}}_{\mathrm{p}}$, and found

$$log\left(\frac{P_f}{P_s}\right)=13.5-58.0\ast {\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}},$$

(3)

where $P_f$ is the probability of not successfully forming a nanosuspension, and $P_s$ is the probability of successfully forming a nanosuspension (Figure 9B, n=16). Like our previous regression model in Equation (1), we recomputed this updated model using 100 independent rounds of 4-fold cross validation and found a nearly-identical classification accuracy of $0.74\pm 0.20$, indicating that our model should be applicable to data beyond the training set. The expression itself is also similar to that found in Equation (1), and suggests that a higher values of ${\mathrm{FASA}}_{\mathrm{p}}^{\mathrm{pred}}$ result in more stable nanosuspensions. Finally, we apply Equation (3) to predict the probability of “success" for GDC-0810 nanoparticles (Table 2). Our model predicts that any excipient would be able to form a stable nanosuspension with the salt form of GDC-0810, which aligns with our experimental observation that the salt of GDC-0810 can form a nanosuspension without excipient.

Overall, our model suggests that drug-excipient combinations with a higher polar fraction of solvent-accessible surface area will form the best nanosuspensions. This property can be maximized by increasing the charge and hydrogen bonding properties of excipients or by using salt forms of drug molecules. However, this simplified model should be used with some amount of caution. First, to test the ability of excipients to adsorb to stable nanoparticles, we placed restraints on the position of drug atoms. These extra potentials prevent crystal rearrangements, which would be overestimated in unrestrained simulations, because computational feasibility required using smaller nanoparticles than those predicted experimentally. As such, effects arising from crystal deformation may not be fully captured by our model. Further, it is possible that a sufficiently polar drug-excipient combination will be soluble in solution, and not form a stable nanosuspension. We begin to observe such trends with the ionic form of GDC-0810, where less excipients adsorb to the drug crystal in our simulations. These effects may also manifest themselves in drug solubility, where potentially soluble salt isoforms may prevent the formation of stable nanosuspensions. Such considerations are not taken into account in our model, which may mean predictions from Equation (1) and Equation (3) might overestimate the stability of some nanosuspensions. As our work here primarily addressed explored hydrophobic drug-excipient combinations in which aggregation and Ostwald ripening are the primary concerns, future studies may want to employ similar methods to test increasingly polar excipients and elucidate the balance between stable nanosuspensions and solubility.

5. Conclusions

In this study, we utilized molecular dynamics simulations to evaluate the characteristics that determine nanosuspension stability. We began by examining how hydrophobic interactions drive how excipients adsorb to naproxen, indomethacin, and itraconazole. Then, we showed that the polar fraction of solvent-accessible surface area from our simulated drug-excipient nanoparticle was correlated with the minimum ratio of excipient:drug necessary to form stable nanosuspensions. Inspired by this relationship, we created a logistic regression model that predicts the probability of forming a stable nanosuspension from the polar fraction of solvent-accessible surface area from simulation. To test this model, we studied GDC-0810, which can be crystallized in different isoforms. Our simulations revealed that the crystal isoform had a large impact on excipient adsorption. When we studied the fraction of polar surface area, we found that the ionic crystal form is far more polar, and, therefore, is more likely to form stable a stable nanosuspension than the neutral crystal. These results aligned with experimental observations where we observed stable nanosuspensions of GDC-0810 without any excipient. We were then able to extend our model to two-dimensional properties and discovered our trends can be explained by hydrogen bonding and charge properties. While this simplified model was unable to resolve the more subtle differences given by the minimum weight percentage, it could predict the probability of forming a successful nanosuspension. Overall, we utilized a combination of modeling and targeted experiments to discover that salt forms of drug crystals and more polar excipients can improve nanosuspension stability.