Fabrication and Characterization of β-Cyclodextrin/Mosla Chinensis Essential Oil Inclusion Complexes: Experimental Design and Molecular Modeling

Essential oils (EOs) are primarily isolated from medicinal plants and possess various biological properties. However, their low water solubility and volatility substantially limit their application potential. Therefore, the aim of the current study was to improve the solubility and stability of the Mosla Chinensis (M. Chinensis) EO by forming an inclusion complex (IC) with β-cyclodextrin (β-CD). Furthermore, the IC formation process was investigated using experimental techniques and molecular modeling. The major components of M. Chinensis ‘Jiangxiangru’ EOs were carvacrol, thymol, o-cymene, and terpinene, and its IC with β-CD were prepared using the ultrasonication method. Multivariable optimization was studied using a Plackett-Burman design (step 1, identifying key parameters) followed by a central composite design for optimization of the parameters (step 2, optimizing the key parameters). SEM, FT-IR, TGA, and dissolution experiments were performed to analyze the physicochemical properties of the ICs. In addition, the interaction between EO and β-CD was further investigated using phase solubility, molecular docking, and molecular simulation studies. The results showed that the optimal encapsulation efficiency and loading capacity of EO in the ICs were 86.17% and 8.92%, respectively. Results of physicochemical properties were different after being encapsulated, indicating that the ICs had been successfully fabricated. Additionally, molecular docking and dynamics simulation showed that β-CD could encapsulate the EO component (carvacrol) via noncovalent interactions. In conclusion, a comprehensive methodology was developed for determining key parameters under multivariate conditions by utilizing two-step optimization experiments to obtain ICs of EO with β-CD. Furthermore, molecular modeling was used to study the mechanisms involved in molecular inclusion complexation.


Introduction
Essential oils (EOs) are naturally occurring volatile components isolated from medicinal plants, exhibiting hydrophobic properties and distinct aromas. EOs have been widely used in the pharmaceutical, cosmeceutical, and food industries for their various biological activities, including antioxidant, anti-inflammatory, antibacterial, sedative, antiparasitic, spasmolytic analgesic activities, etc. [1][2][3][4]. It has been shown that EOs contain different kinds of bioactive compounds, particularly phenolic-derived aromatics and terpenes [5]. Furthermore, recent studies revealed that Mosla Chinensis (M. Chinensis) 'Jiangxiangru' has antioxidant and antibacterial properties, which are related to the phenolic monoterpenes, especially carvacrol [6]. M. Chinensis 'Jiangxiangru' is widely grown around the world, such as in Vietnam, India, Japan, and especially in the south of China. Since ancient with molecular modelling and could be an effective tool for understanding processes and analyzing mechanisms involved in IC formation.
Molecules 2022, 27, x FOR PEER REVIEW 3 of 23 gravimetric analysis (TGA), and dissolution experiments. In addition, the binding mode between β-CD and EO was determined by a phase solubility study, and the inclusion mechanism was elucidated by molecular docking and molecular dynamics simulation. In summary, a simple and comprehensive method was developed for rapidly screening key factors under multivariate conditions for preparing β-CD/EO ICs. Experimental design coupled with molecular modelling and could be an effective tool for understanding processes and analyzing mechanisms involved in IC formation.

Plackett-Burman Design (PBD) Study
The preparation of ICs by ultrasonication technique is influenced by a number of factors. Therefore, it was necessary to screen out the critical quality attributes (such as encapsulation efficiency (Y1), loading capacity (Y2), etc.) following the principles of pharmaceutical quality by design [28]. PBD is an efficient tool to select and identify the significant factors from multi-factors by comparing the variation analysis among the experimental design results [29].
In this study, 17 runs of PBD with two levels were introduced to unbiasedly screen the significant variables affecting Y1 and Y2. Table 5 shows the experimental designs and actual response values. Table 1 shows the variables analysis of the regression coefficient and p-value. Generally, a confidence interval greater than 95% or a significance level less The preparation of ICs by ultrasonication technique is influenced by a number of factors. Therefore, it was necessary to screen out the critical quality attributes (such as encapsulation efficiency (Y 1 ), loading capacity (Y 2 ), etc.) following the principles of pharmaceutical quality by design [28]. PBD is an efficient tool to select and identify the significant factors from multi-factors by comparing the variation analysis among the experimental design results [29].
In this study, 17 runs of PBD with two levels were introduced to unbiasedly screen the significant variables affecting Y 1 and Y 2 . Table 1 shows the variables analysis of the regression coefficient and p-value. Generally, a confidence interval greater than 95% or a significance level less than 0.05 (p < 0.05) represents the significant factor [30]. A smaller p-value indicates a greater significance of the correlation coefficient; among the 6 variables, A (p < 0.005), C (p < 0.05), and F (p < 0.05) had a significant influence on Y 1 , and A (p < 0.01) and F (p < 0.05) had a significant effect on Y 2 . As a result, these significant factors (A, C, and F) were investigated for the next optimization stage using CCD [29].

Central Composite Design (CCD) Studies
In the response surface method (RSM), CCD is one of the most commonly used designs, responsible for establishing the model and optimizing the response. The level of CCD optimization factors was designed according to the results of PBD [26,29]. CCD can estimate the interactions between significant factors and optimal level values. The 3D diagrams of the effect of the variables on Y 1 and Y 2 are shown in Figure 2. The optimized model for describing Y 1 and Y 2 were fitted to a second-order polynomial equation as shown in Equations (1) and (2): Equation (1) revealed that factor A positively impacted Y 1 , while the other two factors (C and F) had a negative impact. Equation (2) revealed that all factors (A, C, and F) had a negative impact on Y 2 .
The positive effect of β-CD concentration (A) on Y 1 and Y 2 was attributed to more CD cavities encapsulating the EO with a 1:1 ratio. A study reported that the total polyphenolic content was improved with increasing the β-CD concentration [31]. Moreover, the effect of ultrasonic power (C) on the yield of cholesterol was increased as the power was increased from 0 to 150 W. However, the increasing rate slowed down when the power was increased from 150 to 250 W [32]. Similarly, reports have shown that modest ultrasonic power is beneficial for the formation of ICs [32]. Similarly, in the current study, ultrasonic power had a negative impact and did not show statistically significant changes in Y 1 and Y 2 . It may be due to the fact that the increase in ultrasonic power caused more violent shock waves and high-speed jets, which easily broke the non-covalent bonds between EO and β-CD [33]. In addition, ethanol increased the inclusion capacity of the CD cavity by forming hydrogen bonds with the hydroxyl groups of the CDs [34]. Nevertheless, when the ethanol concentration reaches a certain amount, the competitive inhibition of EO might occur in the β-CD cavity, which was one of the main reasons that the EO/ethanol ratio (F) negatively affected Y 1 and Y 2 .
The ANOVA results of the quadratic model Equations (1) and (2) are presented in Table 2. According to the ANOVA results, the regression analysis of the experimental designs indicated that the linear coefficients (A), interaction coefficient (AF), and quadratic coefficients (A 2 ) were significant (p < 0.05) to Y 1 , and the linear coefficients (A) and quadratic coefficients (A 2 ) were significant (p < 0.05) to Y 2 . The model fitness was assessed by regression coefficient values (R 2 ) and adjusted R 2 [35]. The value of R 2 and adj. R 2 was 0.9906 and 0.9785 in the Y 1 model, respectively, and the value of R 2 and adj. R 2 was 0.9373 and 0.8567 in the Y 2 model, respectively. These values suggested that both calculated models fitted accurately with the experimental data because of the smaller differences between adj. R 2 and R 2 , which were less than 0.2 [35].

Experimental Design Validation
The optimal variables of the encapsulation efficiency (Y1) and loading capacity (Y2) values were obtained by maximization of the intended values. Maximum Y1 and Y2 values were predicted at 85.52% and 9.03%, respectively, under the optimal level of β-CD/EO ratio (8.73/g:g), ultrasonic power (213.56/W), and EO/ethanol ratio (0.50/g:g). The optimized experiments were carried out in triplicate. Table 3 shows the predicted values and specific actual experimental results. The mean values of Y1 and Y2 were 86.17% and 8.92%, and the RD from predicted values were 0.01% and −0.01%, respectively. These results verified that the protocol model had good prediction and reliability. Table 3. Analysis of the significant parameters in the CCD test.

Number
Predicted Value Actual Value Relative Deviation (%) The interaction of variables and optimum levels (A, C, and F) to Y 1 and Y 2 were plotted by constructing the 3D response surface ( Figure 2). Figure 2a shows the response surface plots for Y 1 based on the changes of the β-CD/EO ratio (A) and ultrasonic power (C), while maintaining the EO/ethanol ratio (F) at a constant concentration. The values of Y 1 increased with A and decreased with C. The interaction effect of the β-CD/EO ratio (A), EO/ethanol ratio (F) on Y 1 , and ultrasonic power (C) kept at the middle concentration were shown in Figure 2b. It could be seen that the value of Y 1 was the highest at the highest A and the lowest F. The combined effect of ultrasonic power (C) and EO/ethanol ratio (F) on Y 1 are shown in Figure 2c. The response surface of Y 1 was elliptical, indicating there might be less interaction between C and F [9]. Similarly, Figure 2d-f showed the interaction effects of the β-CD/EO ratio (A), ultrasonic power (C), and EO/ethanol ratio (F) on Y 2 . The results indicated that both A and C had maximum effects on Y 2 .

Experimental Design Validation
The optimal variables of the encapsulation efficiency (Y 1 ) and loading capacity (Y 2 ) values were obtained by maximization of the intended values. Maximum Y 1 and Y 2 values were predicted at 85.52% and 9.03%, respectively, under the optimal level of β-CD/EO ratio (8.73/g:g), ultrasonic power (213.56/W), and EO/ethanol ratio (0.50/g:g). The optimized experiments were carried out in triplicate. Table 3 shows the predicted values and specific actual experimental results. The mean values of Y 1 and Y 2 were 86.17% and 8.92%, and the RD from predicted values were 0.01% and −0.01%, respectively. These results verified that the protocol model had good prediction and reliability. Table 3. Analysis of the significant parameters in the CCD test.

Predicted Value
Actual Value Relative Deviation (%) The surface morphology of the molecule would change in response to interactions between guests (drugs) and CDs (hosts) [9]. The SEM images of β-CD, PM, and β-CD/EO ICs are shown in Figure 3. As shown in Figure 3a, β-CD was composed of smooth blocky particles with irregular shapes. Figure 3b shows the PM image, showing precipitating solids consisting of EO and β-CD with no fixed shape aggregating and adhering to the β-CD surface [36]. The image of β-CD/EO ICs is shown in Figure 3c. The β-CD/EO ICs consisted of irregular shapes and aggregated small particles with smooth surfaces [37]. Compared with β-CD and PM, the size of β-CD/EO ICs showed no visible fractures, cracks, or pores, indicating that the IC was different in morphology formation of a new substance.
ICs are shown in Figure 3. As shown in Figure 3a, β-CD was composed of smooth blocky particles with irregular shapes. Figure 3b shows the PM image, showing precipitating solids consisting of EO and β-CD with no fixed shape aggregating and adhering to the β-CD surface [36]. The image of β-CD/EO ICs is shown in Figure 3c. The β-CD/EO ICs consisted of irregular shapes and aggregated small particles with smooth surfaces [37]. Compared with β-CD and PM, the size of β-CD/EO ICs showed no visible fractures, cracks, or pores, indicating that the IC was different in morphology formation of a new substance.

FT-IR Analysis
FT-IR analysis is frequently used to study the chemical structure and possible interactions within the ICs [37]. The successful formation of ICs was determined by com-

FT-IR Analysis
FT-IR analysis is frequently used to study the chemical structure and possible interactions within the ICs [37]. The successful formation of ICs was determined by comparing the differences among β-CD/EO IC, β-CD, EO, and PM, such as the disappearance, displacement, or reduction of absorption peaks [38]. The FT-IR spectrum of β-CD, EO, PM, and β-CD/EO IC within the range of 400~4000 cm −1 wavelength is shown in Figure 4.
The spectrum of β-CD is presented in Figure 4a. The absorption peaks at 3360 cm −1 and 2925 cm −1 were assigned to O-H and C-H stretching vibrations, respectively [39]. The absorption peak at 1669 cm −1 was assigned to H-O-H stretching. The peaks at 1158 and 1028 cm −1 represented the vibrations of the C-O-C symmetrical and asymmetrical stretching, respectively [40,41]. Due to the presence of multiple components, the EO shows multiple absorption peaks. As shown in Figure 4b, the EO shows absorption bands in the range of 3200~3500 cm −1 , which were related to the presence of -OH stretching. The two peaks at 2875 cm −1 and 2962 cm −1 were the C-H vibrational couplings of -CH 2 . The peaks at 1620, 1584, and 1513 cm −1 were assigned to the C=C stretching of the benzene ring of aromatic substances, such as carvacrol and thymol [42]. Compared with β-CD/EO IC, the absorption peaks of PM at 2959 cm −1 and 1030 cm −1 were redshifted, and the vibration intensity was weakened ( Figure 4c). However, the peak intensity of β-CD/EO IC was slightly changed and was similar to β-CD in the range of 400 cm −1 to 1500 cm −1 . Additionally, there were three peaks with weak stretching at 1516~1621 cm −1 owing to the bands for EO being obscured by strong and broad β-CD bands (Figure 4d) [43]. Furthermore, the peaks of β-CD/EO IC shifted at 3365 cm −1 and weakened at 1165 cm −1 [6,44]. Compared with β-CD/EO IC, the characteristic absorption peaks of EO almost totally disappeared at 500~1600 cm −1 , and there were no new peaks observed. These changes indicated that the β-CD/EO IC had successfully formed via non-covalent bonds, such as hydrogen bonds [45,46]. β-CD/EO IC, the absorption peaks of PM at 2959 cm −1 and 1030 cm −1 were redshifted, and the vibration intensity was weakened (Figure 4c). However, the peak intensity of β-CD/EO IC was slightly changed and was similar to β-CD in the range of 400 cm −1 to 1500 cm −1 . Additionally, there were three peaks with weak stretching at 1516~1621 cm −1 owing to the bands for EO being obscured by strong and broad β-CD bands (Figure 4d) [43]. Furthermore, the peaks of β-CD/EO IC shifted at 3365 cm −1 and weakened at 1165 cm −1 [6,44]. Compared with β-CD/EO IC, the characteristic absorption peaks of EO almost totally disappeared at 500~1600 cm −1 , and there were no new peaks observed. These changes indicated that the β-CD/EO IC had successfully formed via non-covalent bonds, such as hydrogen bonds [45,46].

TGA Analysis
One of the simplest and most commonly used methods for studying the thermal behavior of any formulation or ICs is TG with a continuous increase in temperature [47]. The TGA/DTG curves of β-CD, PM, and β-CD/EO ICs are shown in Figure 5. There were two main mass loss peaks of β-CD, as shown in Figure 5a. The first weight loss was 10.86% from 30 to 100 • C, associated with the evaporation of water molecules in the β-CD cavity [48]. The second weight loss was 67.38% ranging from 270 to 350 • C, and the maximum degradation rate (2.58 mg/min) happened at 325 • C, which was related to the decomposition of β-CD molecules [39,40]. Figure 5b showed three lost peaks of the PM at 96, 135, and 323 • C, which were associated with water evaporation, EO decomposition, and β-CD degradation, respectively [16]. However, the peak of EO volatilization was not observed in the DTG curve of β-CD/EO IC (Figure 5c). At temperatures ranging from 100 to 300 • C, only one decomposition stage of β-CD/EO IC was characterized in the TGA curve with a slight mass loss of 4.42%. In addition, the maximum mass degradation rate was 1.81 mg/min at 322 • C. Furthermore, the DTG curve of β-CD/EO IC showed a peak (1.81 mg/min) at 322 • C lower than β-CD (2.58 mg/min) and PM (3.05 mg/min), indicating that the thermal stability of EO increased and there were no adverse interactions between EO and β-CD [49]. These results suggested that the β-CD/EO IC was effectively fabricated, and the thermal stability was improved, which was consistent with previous studies [50].
to 300 °C, only one decomposition stage of β-CD/EO IC was characterized in the TGA curve with a slight mass loss of 4.42%. In addition, the maximum mass degradation rate was 1.81 mg/min at 322 °C. Furthermore, the DTG curve of β-CD/EO IC showed a peak (1.81 mg/min) at 322 °C lower than β-CD (2.58 mg/min) and PM (3.05 mg/min), indicating that the thermal stability of EO increased and there were no adverse interactions between EO and β-CD [49] . These results suggested that the β-CD/EO IC was effectively fabricated, and the thermal stability was improved, which was consistent with previous studies [50].

In Vitro Dissolution Analysis
The results of dissolution profiles of EO, PM, and β-CD/EO IC are shown in Figure 6. According to the results, the release rate of EO increased gradually over time. The cumulative release amount of EO was 45.42% within 180 min, and the release reached equilibrium at 45 min. The release behavior of PM was similar to EO, and the maximum release was 55.98% within 60 min. On the contrary, the equilibrium time of β-CD/EO IC was 150 min, and the maximum cumulative release was 54.66%. These results indicated that β-CD/EO IC had slower release characteristics than EO and PM, which could prolong the EO's curative effect. Therefore, it can be concluded that the bioactivity of the volatile components can be preserved by forming IC with β-CD.
The dissolution mechanism was associated with the binding ability of EO to β-CD, including nonspecific and specific binding [9]. The nonspecific binding was also named physical adsorption. As shown in Figure 6, a certain amount of EO was released quickly from β-CD/EO IC in the first 30 min, which was attributed to the physical contact formed by nonspecific binding, which is more susceptible to being influenced than the strong binding force in the cavity with continuous agitation [49]. In the latter time, from 30 to 180 min, the release of EO was slow and became more stable, which was related to the EO appropriately encapsulated within the cavity of β-CD through the intramolecular hydrogen bonds and Van der Waals forces. In this study, the dissolution behavior of EO was consistent with the previously reported studies, in which flavonoids had different degrees and rates of dissolution after being encapsulated by β-CD [51].

In Vitro Dissolution Analysis
The results of dissolution profiles of EO, PM, and β-CD/EO IC are shown in Figure 6. According to the results, the release rate of EO increased gradually over time. The cumulative release amount of EO was 45.42% within 180 min, and the release reached equilibrium at 45 min. The release behavior of PM was similar to EO, and the maximum release was 55.98% within 60 min. On the contrary, the equilibrium time of β-CD/EO IC was 150 min, and the maximum cumulative release was 54.66%. These results indicated that β-CD/EO IC had slower release characteristics than EO and PM, which could prolong the EO's curative effect. Therefore, it can be concluded that the bioactivity of the volatile components can be preserved by forming IC with β-CD.
The dissolution mechanism was associated with the binding ability of EO to β-CD, including nonspecific and specific binding [9]. The nonspecific binding was also named physical adsorption. As shown in Figure 6, a certain amount of EO was released quickly from β-CD/EO IC in the first 30 min, which was attributed to the physical contact formed by nonspecific binding, which is more susceptible to being influenced than the strong binding force in the cavity with continuous agitation [49]. In the latter time, from 30 to 180 min, the release of EO was slow and became more stable, which was related to the EO appropriately encapsulated within the cavity of β-CD through the intramolecular hydrogen bonds and Van der Waals forces. In this study, the dissolution behavior of EO was consistent with the previously reported studies, in which flavonoids had different degrees and rates of dissolution after being encapsulated by β-CD [51].

Phase Solubility Study
Phase solubility studies are the most common methods for determining the ability of macromolecules to interact with small molecules [52,53]. Figure 7 shows the results of the phase solubility study, in which the amount of EO increased with the increase of β-CD concentration and the temperature and reached a status of equilibrium. In this study, there were plateaus at 310 and 318 K at phase solubility curves, while the relationship was linear at 328 K, indicating that the temperature would affect the inclusion ratio of β-CD to EO. The possible reason for this phenomenon was that the M. Chinensis EO components were complex, and different components showed a competitive and synergistic relationship with CD, as in our previous study [54]. Generally, a relatively high temperature could enhance the hydrophobic interactions between the molecules [55]. A moderate high temperature was attributed to the increase the complexation of EO in β-CD.

Phase Solubility Study
Phase solubility studies are the most common methods for determining the ability of macromolecules to interact with small molecules [52,53]. Figure 7 shows the results of the phase solubility study, in which the amount of EO increased with the increase of β-CD concentration and the temperature and reached a status of equilibrium. In this study, there were plateaus at 310 and 318 K at phase solubility curves, while the relationship was linear at 328 K, indicating that the temperature would affect the inclusion ratio of β-CD to EO. The possible reason for this phenomenon was that the M. Chinensis EO components were complex, and different components showed a competitive and synergistic relationship with CD, as in our previous study [54]. Generally, a relatively high temperature could enhance the hydrophobic interactions between the molecules [55]. A moderate high temperature was attributed to the increase the complexation of EO in β-CD.

Phase Solubility Study
Phase solubility studies are the most common methods for determining the ability of macromolecules to interact with small molecules [52,53]. Figure 7 shows the results of the phase solubility study, in which the amount of EO increased with the increase of β-CD concentration and the temperature and reached a status of equilibrium. In this study, there were plateaus at 310 and 318 K at phase solubility curves, while the relationship was linear at 328 K, indicating that the temperature would affect the inclusion ratio of β-CD to EO. The possible reason for this phenomenon was that the M. Chinensis EO components were complex, and different components showed a competitive and synergistic relationship with CD, as in our previous study [54]. Generally, a relatively high temperature could enhance the hydrophobic interactions between the molecules [55]. A moderate high temperature was attributed to the increase the complexation of EO in β-CD.

Molecular Docking Analysis
Molecular simulation technology provides an effective method to study the interactions and structures between hosts and guests in a CD system [15,50]. Molecular docking is a powerful way to investigate the inclusion mechanisms at the molecular level and provides the possible inclusion mode [37]. The possible mode is determined by geometric and energy complementarity principles, which ultimately determine the change in binding free energy during the formation of the complex [15]. In the present study, molecular docking was conducted by AutoDock Vina to obtain a visual image of carvacrol with β-CD ( Figure 8).
As the docking process was simulated by AutoDock, the ΔE and RMSD values were used to estimate the docking results. The larger the absolute value of ΔE, the more stable the inclusion structure would be, and the smaller the RMSD, the greater the validity of the confirmation [56]. The results of ΔE and RMSD values are shown in Figure 8a, in which the first pattern configuration was selected for the lowest energy and the largest population [57]. Figure 8b and c show the lowest stable energy conformation combined with carvacrol and β-CD. It could be observed that carvacrol was embedded in the β-CD cavity but was not perpendicularly positioned in the symmetry axis of β-CD to allow the maximum amount of H-bonding (2.4 Å) between the molecules. The docking result was consistent with the changes in the FT-IR spectra of the β-CD/EO IC. Pattern 1 in Figure 8a showed that the best-calculated docking score of carvacrol and β-CD was −4.1 kcal/mol, which indicated a strong binding affinity, including hydrophobicity, polarity, repulsion, entropy, and solvent action [51]. A higher negative docking score indicates a more favorable complexation process [58]. Thus, the hydrophobic force was driven by entropy changes, which was in agreement with phase solubility results. As the docking process was simulated by AutoDock, the ∆E and RMSD values were used to estimate the docking results. The larger the absolute value of ∆E, the more stable the inclusion structure would be, and the smaller the RMSD, the greater the validity of the confirmation [56]. The results of ∆E and RMSD values are shown in Figure 8a, in which the first pattern configuration was selected for the lowest energy and the largest population [57]. Figure 8b and c show the lowest stable energy conformation combined with carvacrol and β-CD. It could be observed that carvacrol was embedded in the β-CD cavity but was not perpendicularly positioned in the symmetry axis of β-CD to allow the maximum amount of H-bonding (2.4 Å) between the molecules. The docking result was consistent with the changes in the FT-IR spectra of the β-CD/EO IC. Pattern 1 in Figure 8a showed that the best-calculated docking score of carvacrol and β-CD was −4.1 kcal/mol, which indicated a strong binding affinity, including hydrophobicity, polarity, repulsion, entropy, and solvent action [51]. A higher negative docking score indicates a more favorable complexation process [58]. Thus, the hydrophobic force was driven by entropy changes, which was in agreement with phase solubility results.

Molecular Dynamics Simulation Analysis
MD simulation studies were carried out to better comprehend the inclusion of carvacrol and β-CD. As shown in Figure 9, the representative conformational changes of the complex structure were intercepted from the simulated initial (0 ns) to stable (19 ns). In the simulation images, carvacrol entered the β-CD cavity via the head side. Carvacrol reached a stable structure through hydrogen bonding at the same location, as shown in the docking analysis. The phenyl ring of carvacrol stacks over the glucopyranose ring of β-CD, demonstrating hydrophobic interactions in which the complex can spontaneously form between carvacrol and β-CD [59].

Molecular Dynamics Simulation Analysis
MD simulation studies were carried out to better comprehend the inclusion of carvacrol and β-CD. As shown in Figure 9, the representative conformational changes of the complex structure were intercepted from the simulated initial (0 ns) to stable (19 ns). In the simulation images, carvacrol entered the β-CD cavity via the head side. Carvacrol reached a stable structure through hydrogen bonding at the same location, as shown in the docking analysis. The phenyl ring of carvacrol stacks over the glucopyranose ring of β-CD, demonstrating hydrophobic interactions in which the complex can spontaneously form between carvacrol and β-CD [59].
Trajectory analysis of the molecular dynamics process was performed to explore the position of carvacrol bounding to β-CD in different timescales ( Figure 10). The RMSD for the carvacrol, β-CD, and their complex were plotted individually, as shown in Figure 10a. According to the RMSD plot analysis, there was a dynamic equilibrium between the guest and host molecules. β-CD was a rigid skeleton with a small RMSD value (0.35 nm). The RMSD (0.56 nm) of the complex was larger than β-CD because the carvacrol molecule was encapsulated in the cavity resulting in structural changes. The binding structure of the complex could be considered stable when the RMSD fluctuates within a small range [48]. Furthermore, Figure 10b shows that the potential energy has no abrupt variation, which indicates a supportive dynamic equilibrium in carvacrol and β-CD [60]. Moreover, one carvacrol was bound to one β-CD, resulting in a steric hindrance effect for the volume of the β-CD cavity and the size of the guest molecules.  Trajectory analysis of the molecular dynamics process was performed to explore the position of carvacrol bounding to β-CD in different timescales ( Figure 10). The RMSD for the carvacrol, β-CD, and their complex were plotted individually, as shown in Figure 10a. According to the RMSD plot analysis, there was a dynamic equilibrium between the guest and host molecules. β-CD was a rigid skeleton with a small RMSD value (0.35 nm). The RMSD (0.56 nm) of the complex was larger than β-CD because the carvacrol molecule was encapsulated in the cavity resulting in structural changes. The binding structure of the complex could be considered stable when the RMSD fluctuates within a small range [48]. Furthermore, Figure 10b shows that the potential energy has no abrupt variation, which indicates a supportive dynamic equilibrium in carvacrol and β-CD [60]. Moreover, one carvacrol was bound to one β-CD, resulting in a steric hindrance effect for the volume of the β-CD cavity and the size of the guest molecules.

Preparation of ICs
β-CD/EO ICs were prepared under different preparation conditions using the ultrasonication method. Based on previous research [39], β-CD/EO ratio (A/g:g), ultrasonic time (B/min), ultrasonic power (C/W), temperature (D/ • C), reaction water (E/mL:g), and EO/ethanol ratio (F/g:g) were selected as the influencing factors. Initially, 30 mL of deionized water was measured with a graduated cylinder, and then a certain amount of β-CD was weighed over a ME303E balance (Mettler-Toledo Instruments Co., Ltd., Shanghai, China). Next, different proportions of EO and ethanol were added slowly to the above solution, followed by ultra-sonication at different temperatures (40, 50, or 60 • C), powers (180, 270, or 360 W), and times (15, 30 or 45 min). Finally, the white β-CD/EO ICs were obtained after allowing the solution to stand at a low temperature for 24 h and then drying it in an oven at 40 • C for 72 h. The β-CD/EO ICs were stored in a desiccator with granular silica gel at 25 • C for subsequent use.
The physical mixture (PM) was prepared at the optimal conditions based on the optimization results of the central composite design, in which the mass of EO and β-CD were weighted out at a specified proportion (8.73:1) and mixed homogeneously for subsequent use.
Encapsulation efficiency (Y 1 ) and loading capacity (Y 2 ) of EO in ICs were measured using the alcohol extraction method. About 10 mg of β-CD/EO IC was extracted with 10 mL of ethanol for 30 min in an JY92-IIDN ultrasonic bath (Ningbo Xinyi Ultrasonic Equipment Co., Ltd., Ningbo, China), and then filtered by 0.22 µm filter membrane. The content of Y 1 and Y 2 was calculated using a UV-spectrophotometer (UV-2600, Japan) at a wavelength of 277 nm, and the concentration of EO was calculated using the established standard curve equation (Y = 0.0147X (µg/mL) − 0.0292, R 2 = 0.9992). Y 1 and Y 2 were calculated using Equation (3) and (4) [61], respectively: where m 1 refers to the mass of encapsulated EO, m 2 is the mass of total EO added initially, and m 3 is the total mass of β-CD/EO IC.

Step 1: PBD for Screening Key Variables
The PBD experiments were designed using JMP software (version 10.0.0, JMP Statistical Discovery), which generated 15 runs based on 6 variables (A, B, C, D, E, and F). The experiment value was set for each factor at two levels: high level (+1) and low level (−1). Moreover, the center point was set at zero level (0) [20]. Table 4 shows the matrix and corresponding variables of the PBD. Table 5 shows the specific experimental sequence and coding as well as the experimental designs and actual response values. The test results and regression coefficients were analyzed by analysis of variables (ANOVA) to select the significant factors (p < 0.05) for the next step of the experimental design. Table 4. Levels of the factors tested in PBD.

Step 2: CCD Optimization
The experimental values with significant factors for Y 1 and Y 2 were analyzed and optimized by CCD. Y 1 is the main response value that significantly influences the inclusion outcome [58]. In CCD, significant factors were set at five levels (−α, −1, 0, +1, and +α, where α = 1.287), and non-significant factors were maintained at the zero level [27,62]. The details of the variable labels of CCD are shown in Table 6. The specific codes and the response results of experiments are shown in Table 7 as well as the experimental designs and actual response values. The second-order polynomial equation was fitted on the results of the CCD optimization experiment. The relationship between independent variables and the response value was calculated using Equation (5) [20,27].
where Y, β 0 , β i , β ii , and β ij represented the experimental response, intercept term, linear coefficient, square coefficient, interaction and coefficient, respectively. X i and X j were independent factors.

Experimental Design Validation
According to the experimental results of CCD, the binomial regression model was used to optimize the response values (encapsulation efficiency (Y 1 ) and loading capacity (Y 2 )) based on the optimization module. The maximum willingness mode was utilized to predict the range of optimization parameters. The ANOVA tests between predicted and actual results were conducted to verify whether there were significant differences. Each actual test was conducted in triplicate. Then, the relative deviation (RD) [63] was used to verify the reliability of the optimization parameter combination, and the RD was calculated using Equation (6): Relative Deviation (%) = Predicted value − Actual value Predicted value * 100 (6)

SEM Determination
The morphology of β-CD, PM and β-CD/EO IC was investigated by FEI Quanta 250 (Oxford Instruments, Czech) scanning electron microscope. Before measurement, each sample was sputtered with gold plating to increase the electrical conductivity and imaged under 15 kV acceleration voltage. The SEM images were collected under a certain magnification [42].

FT-IR Determination
Infrared spectra of the samples were obtained at room temperature using Fouriertransform infrared (FT-IR) spectroscopy (PerkinElmer, UK). β-CD, PM, and β-CD/EO IC with KBr at a ratio of 1:100 were ground and mixed thoroughly together, and two drops of EO were attached to the surface of the KBr tablet, in which the samples were tested at the wavenumber range of 400~4000 cm −1 and the scanning speed was 16 cm −1 /s. The FT-IR spectral data were recorded after the baseline had been smoothed and corrected by the built-in software [64].

TGA Determination
The thermal properties of β-CD, PM, and β-CD/EO IC were investigated by Exstar TG/DTA 6300 TG analyzer (SII Nano, Japan). In each case, about 5~10 mg of sample was heated from 30 to 700 • C under a 200 mL/min nitrogen atmosphere with a heating rate of 10 • C/min. The thermal behavior of samples was analyzed by thermogravimetry (TG) and derivative thermogravimetry (DTG) curves [64].

In Vitro Dissolution Study
The in vitro dissolution tests were carried out to determine the dissolution behavior of β-CD/EO IC using a ZRS-8G intelligent dissolution apparatus (Tiandatianfa Technology Co. Ltd., Tianjin, China). A certain amount of EO (about 37.30 mg), PM, and β-CD/EO IC (containing an equivalent amount of EO) were separately added to 500 mL phosphate buffer solution (pH = 6.8) at 37 ± 0.5 • C, and then, the rotation speed was set at 50 rpm [65,66]. At predetermined time intervals, 5 mL of the dissolution solution was taken out and centrifuged at 4000 rpm for 5 min. Subsequently, the fresh buffer with the same temperature and volume was replenished to keep the volume constant. The EO concentration in the solution was assayed at 277 nm using a UV spectrophotometer.

Phase Solubility Study
Phase solubility study is widely used to evaluate the interactions between the guest molecule and host molecule at different concentrations [67]. Briefly, the excess EO (300 µL) was added to 5 mL aqueous solution with a series of β-CD solutions ranging from 0 to 10 mmol/L. Moreover, the β-CD/EO suspensions were placed into a heating kettle with a stirring device for 48 h to achieve equilibrium at three temperatures (310, 318, and 328 K). Then, it was filtered through 0.22 µm filters, and the EO absorbance was determined at a wavelength of 277 nm via UV-spectrophotometer. The amount of EO in the solution was plotted against β-CD quantity [68].

Molecular Docking Study
To determine the most possible conformation between EO and β-CD, the interactions were explored using carvacrol, which is the main component found in EO. The complex between carvacrol and β-CD was generated, and the docking fraction was calculated by AutoDock Vina (version 1.5.7) [69]. In the first step, the crystal structure of β-CD was obtained from Cambridge Crystallographic Data Centre (CCDC) and the CCDC ID was 1,235,577 [48]. The 3D structure of carvacrol was downloaded from PubChem (CID 10,364). Initially, the β-CD was set as rigid acceptor molecules, and the carvacrol was used as the ligand molecule that was allowed for flexible twisting [70]. The three-dimensional grid box with 40*40*40 Å size was created and calculated by AutoGrid. Then, 10 simulations with other default parameter settings were performed by Autodock Vina based on the Lamarckian genetic algorithm (LGA) [71].
The docking results were estimated according to binding energy (∆E) parameters and root mean square deviation (RMSD) values. ∆E represented the energy released, which was mainly used to judge the tightness of the inclusion structure. The RMSD was a deviation statistic that evaluated the stability of the simulation system, which represented the structural changes and atomic fluctuations at the initial position [72].

Molecular Dynamics Simulation Study
Molecular dynamics (MD) simulation was performed to understand the mechanism of the carvacrol encapsulation in the β-CD. MD calculations were performed with the Gromacs molecular dynamic simulation package (version 2021.1) [59]. The visual molecular dynamics (VMD) software was used for all molecular visualizations [73]. The RMSD of all atoms was calculated based on the initial conformations. The Antechamber program was used to generate carvacrol and β-CD force field parameters, and the tLeap (Amber 18.0 Module) was used to generate the molecular topology.
The CHARMM27 force field parameters were used for both molecules, and the system was solvated using the TIP3P model [13]. The system containing a total of 4000 water molecules, 1 molecule of carvacrol and 1 molecule of β-CD, was used for the simulation. The steepest descent method was used to minimize the system's energy, and the NPT ensemble was performed with 100 ps. The simulation was run for 30 ns using a 2 fs time step for the production. The V-rescale method was used for temperature coupling, and the Parrinello-Rahman method was used for pressure coupling. The particle mesh ewald (PME) method was used for long-range electrostatic interactions. A 14 Å cut-off was employed for Van der Waals force (VdW) and short-range coulombic interactions.

Conclusions
In summary, two-step experimental designs and molecular dynamics simulation were utilized to investigate the formation and stabilization mechanisms of ICs. Three critical variables, which have an impact on the formation of β-CD/EO IC were obtained from six parameters by PBD and optimized the levels using CCD, including the β-CD/EO ratio, ultrasonic power, and EO/ethanol ratio. The optimum response values of encapsulation efficiency (86.17%) and loading capacities (8.92%) were obtained at 8.73:1 (g:g) of β-CD/EO ratio, 213.56 (W) of ultrasonic power, and 1:2 (g:g) of EO/ethanol ratio. The results of SEM and FT-IR verified that the β-CD/EO IC was prepared successfully. TGA showed that the thermal stability of EO was increased obviously, and the dissolution experiment illustrated that the EO had a sustained-release function after encapsulation. Molecular docking and dynamic simulation revealed a dynamic equilibrium with a molar ratio of 1:1 between carvacrol and β-CD. Therefore, this study suggests that a two-step experimental design manifests a promising potential to screen out key parameters and the potential for developing new technologies combining experiments with molecular modeling.