## 2 Materials and Methods
### In Silico Solubility Study: COSMO-RS  
The in silico solubility of main antioxidant compounds (carnosol, carnosic, and rosmarinic acid) and VACs (eucalypt, camphor, \(\upalpha\)-pinene, borneol, \(o\)-cymene, limonene, and terpinen-4-ol) from rosmenary in the mixture of optimal oil solvent and its derivatives was studied by the more sophisticated COSMO-RS software package from COSMO_logic_, Leverkusen, Germany. The chemical structure of all solvents and solutes studied was mutually transformed by JChemPaint version 3.3 (GitHub Pages, San Francisco, CA, USA) into the simplified molecular input line entry syntax (SMILES) notations, which could be further simulated to their three-dimensional \(\upsigma\)-surface modelling (Figure S1, Supplementary Materials) by embedded Turbomole (TmoleX, version 7.1), where green to yellow codes the weakly polar surfaces, blue represents electron-deficient regions (\(\delta^{+}\)), and red codes electron-rich regions (\(\delta^{-}\)). This 3D information on the repartition of charge density on the molecular surface can be reduced to \(\upsigma\)-profiles and \(\upsigma\)-potentials, which could be used for subsequent prediction of the interactions between oily solvents and bioactive compounds based on a quantum-chemical approach. COSMO-RS is known as a modern and powerful method with the combination of quantum chemical considerations (COSMO) and statistical thermodynamics (RS) for molecular description with predicted thermodynamic properties and further solvent prescreening without any experimental data [18]. As depicted in Figure 2, COSMO-RS generally involves two major steps. All molecules are embedded into virtual conductors simulated in the first step by the COSMO model, where the molecule induces a polarization charge density (\(\upsigma\)) on its surface that is a good local descriptor of the molecular surface polarity. Therefore, molecules during the quantum calculation are converged to their energetically optimal states in the conductor with respect to electron density and geometry. Density functional theory (DFT) with triple zeta valence polarized basis set (TZVP) were used as the standard quantum chemical method throughout this study. The statistical thermodynamics calculation was used in the second step. This polarization charge density was used for the quantification of the interaction energy of pair-wise interacting surface segments concerning electrostatics and hydrogen bonding. The 3D distribution of the polarization charges on the surface of each molecule was converted into a surface composition function (\(\sigma\)-profile) that provided the information about molecular polarity distribution. The thermodynamics of molecular interactions based on the \(\sigma\)-profile obtained was then calculated to chemical potential of the surface segment (\(\sigma\)-potential) using the COSMO_thermX_ program (version C30 release 13.01). The \(\sigma\)-potential describes the likeliness of the solute compound interacting with solvents according to their polarities and hydrogen bonds, where the part of the negative charge of molecules (hydrogen bond acceptor) was located on the right side with positive \(\sigma\)-profile values, whereas the part of the positive charge (hydrogen bond donor) was located on the left side with negative \(\sigma\)-profile values. The solvent's affinity for polarity surface could be interpreted according to the \(\sigma\)-profile and \(\sigma\)-potential in order to better understand the solute-solvent interaction in a mixture state. In addition, COSMO_thermX_ also calculated the relative solubility of COSMO-RS simulated volatile and non-volatile bioactive compounds in various oil solvents and mixtures with their derivatives based on the logarithm of the solubility in mole fractions.