1. Introduction

Since the 1970s, there has been a rapid increase of interest in the industrial application of cyclodextrins (CDs) ^[1]. This growth was associated with an unambiguous confirmation of the nontoxicity and considerable price decrease of CDs. CDs were earlier considered as “toxic”, as mistakenly ascribed by French in 1957 ^[2]. Fortunately, Szejtli postulated the lack of toxicity, a view that was thoroughly examined and, finally, widely accepted [3,4]^[3][4]. Since then, the price of CDs has rapidly dropped; they currently cost as low as USD 5 per kilogram ^[5].

CDs are primarily used in pharmaceutical formulations due to their unique properties ^[6]. Through the formation of host–guest complexes, they increase the solubility of poorly soluble drugs and protect substances against external factors, such as light, humidity, and heat. Moreover, CDs could mask unpleasant smells or flavors of drugs, which is especially important in formulations dedicated to children. Currently, more than 100 original drugs are manufactured with CDs as excipients [7,8,9]^[7][8][9]. Interestingly, CDs play a role in fighting the COVID-19 pandemic. For example, CDs can be found in formulations of remdesivir ^[10] (an antiviral drug administered to treat COVID-19) and Johnson & Johnson’s Jannsen (J&J/Janssen) single-shot COVID-19 vaccine ^[11]. CDs are also utilized in the production of face masks, which are widely used to help slow the spread of the coronavirus ^[12].

The desirable properties of CDs in the pharmaceutical field can be explained at the molecular level ^[13]. Cyclodextrin (CD) molecules resemble a “doughnut” ring, in which small, non-polar substances can be entrapped. The external fragments of CD molecules are polar due to the presence of hydroxyl groups. When a non-polar substance (e.g., a poorly soluble drug) enters the molecular hole of cyclodextrin, the formed host–guest complex is polar (at outside) and, therefore, is more soluble than the separated guest molecule.

Solid complexes of active pharmaceutical ingredients (APIs) with CDs are amorphous, in major cases ^[14]. Therefore, one can conclude that they are characterized by high disorder and dynamics. Hence, molecular dynamics (MD) simulations seem to be an ideal tool to investigate their structure and properties.

2. Molecular Modeling of CD Host–Guest Complexes—Theoretical and Practical Aspects

While designing and preparing the CD complexes, the formation of the covalent bond between host and guest molecules is usually avoided in order to enable the guest to freely escape from the CD cavity. Therefore, CD host–guest complexes are stabilized by non-covalent interactions, such as H-bonds and van der Waals forces (vdW). In most cases of such complexes, a few energetically-similar minima, found as the creation of one bond, usually requires the breaking of another of similar strength. This is also in agreement with the experimental results—the solid complexes of CDs are usually amorphous, as the guest molecule can exist in various conformations and in various poses. Therefore, modeling static structures of CD complexes through geometry optimization, without considering their dynamics, can be a source of inaccuracy, and significant differences between the experimentally and computationally obtained results. Hence, as will be proven in the next section, MD simulation is a perfect method to study CD complexes.

For MD simulations of large molecular complexes, such as ligand–protein, molecular mechanics (MM) methods are commonly used. On the contrary, when MD simulations are performed on relatively small molecules, it is usually at the quantum mechanics (QM) level of theory, which significantly increase the accuracy of calculations, but also their computational costs. In terms of the sizes of the modeled objects, CD complexes are somewhere in between. While geometry optimization calculations on the static structures of CD complexes are, nowadays, performed mostly at the QM level, usually by the means of DFT [25]^[15], the MD simulations are still being performed at the MM level. However, to increase the accuracy of the calculations, while still maintaining their reasonable computational costs, multiple solutions, such as dedicated force fields and special sampling methods, have been developed, and will be described below.

In the reviewed studies, a pre-requisite for each MD simulation involved the creation of a complex between the CD and a guest molecule. To the authors’ best knowledge, so far, there have been no published attempts to simulate this process other than through the docking procedure. In other words, we found no work in which the process of the formation of the complex from substrates, CD, and guest molecules, would be studied using MD simulations. This was probably due to the long simulation needed to observe this reaction. The literature review shows that the most often used docking software are open-source AutoDock Vina [49]^[16] and commercial CDOCKER from the Accelrys Discovery Studio Package [50]^[17]. Surprisingly, the authors of some of the reviewed works underestimate the docking step and do not describe it with sufficient details. The proper initial orientation of the guest inside the CD cavity is particularly important to the linear molecules that may form two substantially different complexes: parallel and antiparallel. During classical MD simulations, it is almost impossible to observe the transition between those two complexes, as it would require complex dissociation, guest reorientation, and subsequent association of the guest and CD.

Coarse-grained (CG) MD are quite different from the above-mentioned all-atom MD simulations. In the CG MD approach, molecules are not represented by individual atoms, but by coarse-grained sites, approximating groups of atoms, such as whole amino acid residue. By decreasing the degrees of freedom, much longer simulation times can be studied at the expense of molecular details [72]^[18].

3. Application of the Molecular Dynamics Simulations for Systems Including Cyclodextrins—The Most Important and Interesting Cases

Two main purposes for CDs application in the research have emerged: (i) CDs as drug carriers and (ii) CDs as extracting agents. As this article is meant to provide guidance for those wanting to perform MD simulations of CDs, we constructed this review from the applicability perspective of CDs.

The main driving force of the guest inclusion in CDs are vdW interactions.

Such reoccurring results, as well as the number of articles on the MD–CD topic, emphasize the usefulness of this calculation technique, in terms of CD analysis.

The third topic is composed of other applications of CDs, for instance CDs as capping agents for gold nanoparticles (AuN) [177]^[19] or recognizing agents for N -methylfulleropyrrolidine regioisomers [179]^[20]. In those cases, multiple CDs are applied, regarding their roles. In the mentioned fullerene-derivatives study, one or two γ-CD molecules were used to separate the isomers. Whereas for AuN, γ-CD served as stabilizing agents; as an example: 1007 AuN and 60 CDs were simulated at once, in a way that they completely surrounded the guest.

4. Conclusions

As shown above, the number of works presenting the results of MD simulations on CD host–guest complexes have rapidly increased since the early 2010s. Moreover, the applied methods are becoming more sophisticated at increasing the accuracy of such calculation, to extend their application.

For example, while in the oldest works on this topic, the standard FFs were applied, currently, CD-dedicated FFs, such as Glycam06, are commonly used. Further, in the literature, examples of hybrid QM/MM methods, widely applied in protein–ligand interaction models, could be found. Similarly, solvent treatment methods are being improved. Initially, the implicit methods were used; however, they were found inaccurate and, thus, were replaced by explicit methods, such as TIP3P or even TIP4P.

MD of CD complexes are no longer simulated solely to analyze the RMSD or RMSF, but also for post-MD calculations, to better assess the binding and, thus, thermodynamic stability. Hence, GBSA or PBSA methods are commonly used. Moreover, very recent applications of umbrella sampling and steered (biased) MDs prove that state-of-the-art methods of MD could be useful in studying and designing CD complexes.

The versatility of MD simulations allows for studying all kinds of complexes. Both native and substituted CDs (as hosts) and a whole range of APIs (as guests) are being frequently modeled this way, in 1:1 and 2:1 molar ratios. Moreover, the application of CDs as extracting agents could also be evaluated by the means of MD simulations.

As shown in this review, MD simulations could be used to predict the structure, solubility, and stability of CD complexes and, thus, used to explain, at the molecular level, the experimental results. In addition, such simulations are now being used at the stage of the design of CD complexes, preceding their experimental preparations. With the anticipated progress in MD simulations, this second application is expected to become even more popular.