For the first part of the study, we used micelles self-assembled from BS, PL, or a combination of the two in a 4:1 ratio. The motivation was to understand the individual contributions of these two major components of the intestinal fluid to the solubilization of each of the four molecules as well as any synergetic effects. Initially, all the systems contained more molecules than those specified in Table 1. By monitoring the equilibration of micelle aggregation kinetics, it was found that the aggregation number of NaTC molecules in a BS-only micelle was approximately eight, and the number of DLiPC molecules stably residing in a PL micelle was approximately 14.

When both bile salts and phospholipids were present at the beginning of the equilibration process, the resulting micelles contained approximately 24 NaTC and six DLiPC molecules (i.e., the 4:1 ratio was maintained upon micelle formation in the simulations). In addition, the exchange of individual NaTC and DLiPC molecules continued even after 100 ns; however, the ratio in the micelles did not change significantly. The micelles were prolate ellipsoids with eccentricity values between 0.12 and 0.4 (see Table 1). Pure PL aggregates were generally more spherical, whereas the BS aggregates were the least spherical.

The number of molecules and sizes of the self-assembled PL and mixed BS-PL micelles were significantly larger than those of BS-only micelles. The average gyration radii of the mixed BS-PL self-assembled micelles were found to be 1.5, 1.4 nm, and 1.02 nm. Owing to the hydrophobic effect during the self-assembly process, the tails of the DLiPC molecules tended to aggregate at the center of the micelles to minimize contact with water. Sodium taurocholate molecules are more complex in terms of lipophilicity, with multiple sites across the molecules containing hydrophobic and hydrophilic groups. This makes them prefer the surface of the aggregate, creating a shell around the DLiPC core, with the DLiPC headgroups oriented towards the outside of the micelle. Coalescence events between individual micelles were also observed during the simulations; however, these aggregates would split into smaller colloids within nanoseconds.

As discussed above, the inhomogeneous distribution of the BS and PL components makes it difficult to compare the free energy profiles of different molecules because the local environment might differ during pulling (and US simulations) for each compound. Moreover, the shape of the micelles may be affected by the small drug or surfactant molecule. For example, sodium caprate (a fatty acid) has a hydrophilic head in contact with water in the randomly assembled micelles. Its interactions with the colloid are, to a great extent, dependent on the coverage of the hydrophobic tail by the micelle molecules. In the case of the BS-only aggregate, this tail was only partially inserted into the micelle. The free energy of a sodium caprate at the BS micelle surface was 14.5 kJ/mol lower than that of the bulk (Fig. 2). In the PL micelle, the free energy difference between the surface and bulk is ~ 21 kJ/mol, which is approximately equal to that in the mixed BS-PL micelle case, indicating no obvious synergetic effect for a mixed BS and PL micelle for sodium caprate. In contrast, danazol was solubilized twice as well in PL compared to BS micelles, but even better in the mixed micelle. The total effect of the mixed micelle on the solubilization capacity of danazol is stronger than the sum of the individual contributions of NaTC and DLiPC (~ 53 kJ/mol vs ~ 16.5 kJ/mol and 33 kJ/mol respectively). Thus, danazol tightly interacts with both BS and PL molecules and inserts almost completely into the mixed micelles. Prednisolone is more solubilized by the phospholipid molecules. The change in the free energy when prednisolone was pulled from the mixed micelle was slightly lower than that for the smaller PL aggregate (~ 26 versus ~ 24 kJ/mol). The free energy of solubilization of prednisolone in the BS micelle is of the same order as that of sodium caprate (~ 15 kJ/mol). Finally, probucol, the most hydrophobic compound, is almost equally well solubilized within BS and PL micelles (− 54 kJ/mol and − 55 kJ/mol), and in the mixed micelle, the free energy is almost equal to the sum of the two (− 96 kJ/mol). The drug molecule is in the most favorable energy state once it is completely screened from contact with water.

Combined free energy profiles and radial distribution of select beads. Free energy differences (\(\Delta G\)) are expressed in kJ/mol. The positions of various groups of the small molecules are depicted in dark and light green and are mapped in the lowest panel

From the results above, it seems that both danazol and probucol are better solubilized in the micelles than prednisolone and sodium caprate.

Prednisolone solubilization in PL aggregates seems to be more advantageous compared to BS micelles. Danazol resides on the surface of the colloid upon approaching the membrane, and it is possible that the molecule would transfer to the membrane. Probucol was well-solubilized by both BS and PL micelles without significant differences in affinity. However, as probucol is entirely covered with NaTC and DLiPC molecules, any transfer from the micelle to the cell membrane towards systemic circulation would require either breakage of the micelle or fusion between it and the bilayer. In addition, because the results suggest that a probucol molecule would prefer to occupy the center of the micelle, aggregation, and as a consequence, crystallization might potentially occur within the colloid.

Free energy profiles obtained in this way for colloids randomly assembled around a small molecule thus provide valuable information, but they might be insufficient if energy barriers are present on the way from the bulk to the FE minimum state. This may be a more realistic scenario for drug molecules released from a dosage form in the small intestine. Thus, a profile spanning from the center of the colloid to the bulk may be needed to verify the properties of the free energy well within the colloid.

Simulations with preorganized micellular structures were performed for pure phospholipids and mixed BS-PL micelles. We were unable to construct a preassembled BS micelle, as those colloids tended to fall apart during the initial simulation phases. A substantial increase in the size of the colloids (Table 1) was necessary to ensure symmetry with respect to the molecular orientation and composition of the colloid.

The free energy profiles were compared for both types of micelles. The free energy profiles of the self-assembled and pre-organized pure PL micelles were qualitatively consistent (Fig. 3a). In all cases, except for sodium caprate, the energy wells were deeper in the larger, pre-organized micelles. For poorly soluble drugs, it makes sense that a larger colloid would screen drug molecules better from aqueous contact. Danazol and prednisolone occupied the PL headgroup region, whereas probucol moved to the center of the colloid after overcoming a small energy barrier at 2–2.5 nm from the interface with water.

Comparison of free energy profiles of self-assembled and pre-organized PL and BS-PL colloids. The top panel shows a schematic representation of the steered MD simulations, in which the small molecules were pulled from the initial position, either at the center of the pre-organized micelle or at the surface of the self-assembled aggregates

For the BS-PL micelles, manual organization into separated PL and BS layers leads to differences in free energy compared to randomly assembled micelles (Fig. 3b). Although the rank order was the same, the energy minima for probucol and danazol were lower (178 kJ/mol and 129 kJ/mol compared to 99 kJ/mol and 51 kJ/mol, respectively), and the positions of the minima were also shifted. Probucol is solubilized to the highest degree, in the same way as for the self-assembled micelles, while caprate, on the other hand, does not show a significant free energy minimum at the outer layer of the mixed pre-organized micelle. There is also a potential barrier at a depth of 2 nm from the water-micelle interface, which indicates the loss of contact with water for the caprate head group.

An alternative way to study micelle-drug interactions is to perform free energy calculations only at specific points of interest along the reaction coordinate. To investigate this approach, we used pre-organized micelles and determined the free energy at three different points: in the micelle center, at the unbiased simulation equilibrium position, and in the bulk (two points were placed on the graph for clarity, at zero and one nm away from the outermost contact with the micelle; only one measurement in bulk was made). As before, the zero point was set as a reference to correspond to the bulk water.

For all four small molecules, there was good (and somewhat expected) quantitative agreement between the free energy change values measured at the chosen points and the profiles observed from the US simulations (Fig. 4). The largest discrepancy was observed for probucol, where the values of \(\Delta G\) were higher than those predicted by the US. Nevertheless, this approach seems to work well for the comparison of small molecule profiles. The limitation of the method, as compared to US simulations, is the lack of information on whether any intermediate energy barriers are present. For example, for danazol and probucol, no information can be obtained from the FE calculations about the region with small local free energy maxima between the FE points. This can be compensated, to some extent, by introducing more measurement points along the reaction coordinate.

Mapping of the selected FE point (FEP) calculations in the bulk, at the interface, and in the center of the pre-organized PL colloids onto the corresponding umbrella sampling (US) profiles. The profiles were normalized to the free energy change in bulk. Four points are shown, with the bulk being depicted with two points: right at the interface and at a distance of one nm from it in the bulk. Standard deviations of the FE calculations do not exceed the size of the corresponding markers

Another solution to speed up the simulations is to use coarse-grained resolution. Therefore, we tested how much faster the simulations are with CG compared to AA, while at the same time examining the differences in the results as a consequence of the reduced level of detail associated with coarse-graining.

The self-assembled BS-PL mixed micelles were chosen, to include interactions with both types of molecules, and also let the small drug/surfactant molecule freely find its optimal position within the micelle. The simulations in this section focused on sodium caprate, danazol, and probucol, and the resulting free energy profiles observed from the CG US simulations can be seen in Fig. 5.

Comparison of the coarse-grained (CG) and all-atom (AA) model profiles for sodium caprate (a), danazol (b), and probucol (c). The combined CG and AA visual representations of the micelles and the equilibrated position of the API within (marked with dashed red ellipsoids), d, e, and f, respectively. In the CG representation, the danazol molecule is completely covered, whereas in AA it has contact with water molecules

Sodium caprate and probucol showed a good match between AA and CG profiles. However, in the danazol case, the results appeared to be significantly different between AA and CG, both about the position and depth of the energy minimum. The energy minimum for danazol in the AA simulations is ~ 51 kJ/mol, with the position being approximately 1.65 nm from the surface of the water. The corresponding free energy change for the CG simulations is more than double that, 112 kJ/mol, and the position is deeper at a distance of 2.8 nm from the micelle surface.

To further understand the AA and CG behavior, we analyzed the SASA of the APIs while in their equilibrium position in the micelle (specifically, only the surface exposed to water was calculated) using two probe sizes, 0.26 nm and 0.4 nm (Table 2). The ratio of this surface exposed SASA to the total surface area of the API was then calculated to find out the portion of the drug surface in contact with water. The data showed an expected trend between the APIs. The contact area with water was minimal for probucol, the most hydrophobic drug, covering 0–4.4% of the surface. Sodium caprate, as a surfactant, had the highest percentage, between 15.6–31.6% of its surface area in contact with water. Danazol had 3.8–4.2% of its surface at the interface with the aqueous phase. Interestingly, the SASA differences between the all-atom and coarse-grained simulations were the least pronounced for danazol, the API that had the most different free energy profiles between the AA and CG resolutions (Fig. 5). This illustrates the complexity of the interplay between the model molecules and the validation required to perform reasonable simulations.

CG simulations have significantly higher performance, but extensive validation, beyond octanol–water partitioning might be required. One of the possibilities might be the combination of the AA and CG techniques, where key points from the AA simulations can serve as a check for the CG free energy profiles. Alternatively, the CG free energy profile can be used as an approximate map or reference for the targeted point-wise free energy calculations.

The computational time required to run the simulations varies significantly and depends on hardware, software, and simulation parameters. However, two characteristics can provide an estimate of the computational costs for simulations with different protocols described in this study. One is the number of windows required to determine the free energy profile, and another is the average time for running such a window. We summarised these data in Table 1 in Online Appendix A. As can be seen, CG and FEP approaches are more efficient. Moreover, multiple screenings might be required with analysis time, as the spring constant for umbrella sampling needs to be gradually adjusted over several iterations, covering the entire space along the reaction coordinate. Nevertheless, CG mapping and validation might be a time-consuming procedure too. Even though several automated tools for topology building are available, it is highly recommended to ensure the validity of the molecules and fine-tune them in many cases. Therefore, fast screening might be organized with scripts and little validation, but a more thorough analysis would be required to introduce the accurate CG topologies. Table S1 is only shown to give an estimate of the time required to run simulations with different protocols used in this study. As can be seen from the last row, one the hardware we used, one can run a complete series of AA simulations for self-assembled colloids of mixed BS and PL in 48 h, whereas the corresponding CG series can be finished within 8.8 h. FE point calculations can be finished within 5.3 h for three points (API in water, outmost energy well, and in the center of the colloid). However, for the later approach, CG simulations should be first run to define the approximate position of the energy well.