Interactions between β-cyclodextrin as a carrier for anti-cancer drug delivery: a molecular dynamics simulation study

Abstract A series of molecular dynamics simulations were performed on 5-fluorouracil (5-Fu), Alendronate (Ald), and Temozolomide (TMZ) anticancer drugs in the presence and absence of β-cyclodextrin (βCD) as a carrier. Thermodynamic investigations showed that the van der Waals interaction energy was dominant in loading all drugs inside the βCD cavity. The sum of the interaction energies illustrated that the highest affinity was related to Ald (−136.5 kJ/mol), which in turn was due to the presence of bulky and charged atoms of phosphorus and oxygen, although TMZ (−115.92 kJ/mol) showed a very high affinity as well. At the same time, the hydrogen bond analysis also represented that Ald had the most hydrogen bond (1.97) with the highest half-life (3.13 ps) with βCD. Investigation of the root mean fluctuation (RMSF) indicated that all the drugs had a relatively rigid structure and maintain this rigidity during loading in the βCD cavity, and in the meantime, Ald was slightly more flexible than 5-Fu and TMZ. The area of the primary hydroxyl rim decreased in all drug-containing systems, which in turn was caused by the attractive interaction of drugs with oxygens in the primary hydroxyl rim. Especially for those drugs that were able to penetrate to the end of the primary hydroxyl rim of the βCD, that means TMZ and 5-Fu. Meanwhile, due to the lack of Ald penetration to the end of the primary hydroxyl rim, the area change in the Ald-containing system was less than in the two others. Communicated by Ramaswamy H. Sarma


Introduction
Nanotechnology provides a platform for researchers to overcome physiological barriers to treating different types of cancer.Hence, a wide range of nanoparticles containing organic and inorganic materials have evolved as nanocarriers (Cho et al., 2008;Ehdaie, 2007;Ferrari, 2005;Lavan et al., 2003;Luk & Zhang, 2015;Mart� ın del Valle et al., 2009;Park, 2013).Novel drug delivery systems encompass new approaches to the delivery of drugs that is now taking benefit of nanoparticles technology and it is generally designed to improve the pharmacological and therapeutic profile of a drug molecule and nanoparticle-based drug delivery systems have considerable potential in the treatment of numerous illnesses (Ganji et al., 2013;Niknam et al., 2022;Rasoolidanesh et al., 2021;Rezvani et al., 2013).Among the carriers used in drug delivery, especially in cancer drug delivery systems, we can mention all kinds of cyclodextrins.Due to this fact that the administration of drugs through different routes especially of those, which are poorly soluble and belong to class II or IV of the biopharmaceutical classification system, represents a major challenge.Remarkably, that most of the cytotoxic anticancer drugs belong to the BCS class IV which comprises substances with both low solubility in aqueous fluids and low apparent permeability (Gidwani & Vyas, 2015;Menezes et al., 2019;Sancho et al., 2016).Cyclodextrins are sufficient to overcome certain forms of disadvantages associated with anticancer drugs, including limited aqueous solubility (hydrophobicity), degradation in digestive fluids, insufficient stability in vitro (shelf life), low bioavailability, short in vivo stability (half-life).The lack of efficient treatment has created the need to develop and implement new technology based on the combined strategy of cyclodextrin complex and nanotechnology with the aim of making the treatment more useful and acceptable (Shelley & Babu, 2018).To improve/ increase solubility and stability, increase bioavailability and dissolution, reduce toxicity, and modify physico-chemical properties, many anticancer drugs along with cyclodextrins and their derivatives have been considered as an inclusion complex (Canbolat et al., 2014;Dos Passos Menezes et al., 2017;Lima et al., 2016;Poulson et al., 2021).The main types of cyclodextrins include 6, 7, or 8 glucose units named a, b, and c, which are crystalline, homogeneous, and non-hygroscopic materials (Vyas et al., 2008).The general structure of these molecules is in the form of a hollow cone, which consists of several glucose units that are covalently connected by an oxygen atom, and due to the lack of free rotation on the surface of the bonds between glucopyranose units, they cannot have a cylindrical shape, so they can be called showed a short conical structure with a hydrophobic cavity in the center and a hydrophilic outer surface (Boonyarattanakalin et al., 2012).In the meantime, b-cyclodextrin (bCD) seems to be the best natural cyclodextrin due to its cheapness, availability, and capacity to form complexes with a wide range of substances, pore size, and more effective pharmaceutical complexes (Crini et al., 2018).In fact, like other cyclodextrins, the characteristic feature of the bCD is its ability to form compounds with various substances through guest-host interactions.The core of this structure can play a role in trapping or encapsulating other substances, including anticancer drugs.The type of solvent, chemical structure and physicochemical properties of CD-drug, and thermodynamic interactions between the guest (CD) and host (drug) could be considered as the factors influencing the formation of suitable host-guest complexes.For example, it has well known the number of CD-CD hydrogen bonds in the bCD was more than the aCD and cCD.As a result, selfaccumulation occurs in the bCD more than in other types (Budhwar, 2018;Mixcoha et al., 2014).
Many studies confirm the use of nanocarrier complex with various anticancer drugs including temozolomide (TMZ) (Bazzazzadeh et al., 2020), alendronate (Ald) (Zhu et al., 2014), and 5-fluorouracil (5-Fu) (Wang et al., 2020), as a potential nanosystem for drug delivery in the treatment of cancers.The study of the favorable interactions of TMZ, Ald, and 5-Fu with organometallic carriers (MOFs) was investigated experimentally (Gu & Meng, 2021;He et al., 2021;L� azaro & Forgan, 2019) and theoretically (Boroushaki et al., 2022) and found that these drugs tend to be loaded on the metal sites of MOFs, that are the adsorption centers of the guest molecules.Recently, it has been well shown that 5-FU tends to absorb stronger with gold-fullerene than C60 (Sabet et al., 2022).In recent years, experimental and theoretical studies of the dynamics and thermodynamics of drug-CD have grown and especially could be attributed to the investigation conducted by Bilensoy et al., in the formulation of a type of vaginal gel from the complex of b-CD/5-fluorouracil in treatment for diseases related to HPV such as cervical cancer or genital warts with a lower dose (Bilensoy et al., 2007).Furthermore, an experimental/theoretical study was done by Biernacka et al., to enhance the drug delivery approach by bCD/Ald inclusion complex in the treatment of bone cancer can be mentioned (Biernacka et al., 2014).Overall, there are a wide set of available techniques in principle allow not only to detect of the host-guest complex formation but also to study of its structure and dynamics of formation/decomposition.Among these, molecular dynamics (MD) simulations are powerful tools that could provide valuable complements to experiment, with information about the details of investigating the drug concentration on the diffusion rate of drug molecules into lipid bilayer membrane (Sodeifian & Razmimanesh, 2019), interactions between drugs-bCD, protein-ligands, peptides-membranes, and any other bio-molecules (Aghazadeh et al., 2020;Ganjali Koli & Azizi, 2016;Koli & Azizi, 2019;Pereva et al., 2016;Roosta et al., 2021).In this research, to get an insight based on dynamics and thermodynamics as well as provide a detailed understanding of molecular mechanisms of the interactions between the bCD (as a carrier) with TMZ, Ald, and 5-Fu drugs, a series of MD simulations have been done.This study, could facilitate the selection of a suitable carrier for these drugs and gain an initial assessment compared to other carriers used with them (Boroushaki et al., 2022).

Simulation methods
The molecular structures of TMZ, Ald, 5-Fu, and the bCD are shown in Figure 1.Simulations were carried out for the seven different systems.The first system which comprises only one bCD and water molecules as a reference system.The next three systems consisted of only one of these drugs and water molecules.The other three systems included one of each of these drugs, b-CD, and water molecules, as explained in Table S1 in detail.The initial positions of drugs and bCD were adjusted so that the drugs were placed at a distance of 8 Å above the secondary hydroxyl rim of bCD (The rim of the wider side that contains O2 and O3 atoms), and at the center of the box, as shown in Figure 2, where box lengths in the XYZ dimensions for all systems were 5 � 5 � 5 nm.The topology and coordinates files based on GROMOS 54a7 force field (Malde et al., 2011;Stroet et al., 2018) were used to describe the neutral form of all molecules, and validation tests for the force field parameters have been performed in another study (Mojdehi et al., 2021).All of these systems were hydrated by 4000 water molecules and with periodic boundary conditions in all three dimensions.For water molecules, the extended simple point charge model (SPC/E) was used (Berendsen et al., 1987).In all systems, steepest descent energy minimization was used to remove undesirable atomic contacts and equilibrated the water molecules within systems ( � Zilinskas, 2006).At first, positions of drugs and bCD were restrained and equilibration was conducted in the NVT ensemble for 1 ns and then followed by equilibration in the NPT ensemble for 9 ns and box lengths in the XYZ dimensions for all bCD/Drug/Water and Drug/Water containing systems were about �4.95 � 4.95 � 4.95 and �4.94 � 4.94 � 4.9 nm, respectively.The range of box volume changes was between 120 (for Drug/Water containing systems) and 122 nm 3 (for bCD/Drug/Water containing systems).After the equilibration steps, all simulations were run for 100 ns from their starting conditions and coordinates of the atoms using the GROMACS 5.1.1 package (Abraham et al., 2015).All the bond lengths were constrained with the LINCS algorithm (Hess et al., 1997).The temperature was set at 300 K through the V-rescale thermostat with a coupling time of 0.1 ps (Evans & Holian, 1985).The pressure was held constant at 1 bar by compressibility 4.5 � 10 À 5 bar À 1 and by coupling the simulation cell to a Parrinello-Rahman barostat, with a coupling time constant of 2 ps (Parrinello & Rahman, 1980).To integrate Newton's equations of motion, the leap-frog algorithm with time step 2 fs was applied (Hockney et al., 1974).The long-range Coulomb interaction was taken into account using the particle mesh Ewald (PME) (Darden et al., 1993) technique and the cut-off distance for Coulomb interactions and van der Waals (vdW) interactions were both set to be 1.2 nm.

Results and discussion
The results of this research are divided into three parts as follows: 1. Thermodynamic properties of the simulated systems; examining the interaction energies of different components, as well as the hydrogen bonds formed between the components.2. Dynamic properties and structural changes of water inside the bCD; investigation of the mobility of components in the simulated systems, studying the drug loading mechanism using the results obtained from the distance of the drugs from the bCD along with comparing the probability of the presence of drugs from the bCD cavity considering the radial distribution functions (RDFs) of water in different simulated systems.
Examining the mobility of drugs using the mean square displacement (MSD) of drugs and bCD in simulated systems.3. Investigating the structural drift, flexibility, and structural changes by calculating the root mean square deviation (RMSD) of drugs and bCD, the root mean square fluctuation (RMSF) of drugs and bCD, and finally changes in different areas of the bCD rims.

Interaction energies
The study of the interaction energies between the various components of simulated systems shows the thermodynamic affinity of drugs to water and bCD.The interaction energy values between the components are summarized in Table 1.Examining the 5-FU-containing systems shows that the dominant interactions between water/drug are electrostatic (Coul) (À 100.16 and À 132.26 kJ/mol in the bCD-containing system and the system without bCD, respectively), although the van der Waals (vdW) interaction in the bCD/5-Fu system (À 20.67 kJ/mol) was also somewhat effective.As such, from a thermodynamic point of view, the solubility of 5-FU was controlled by Coul interaction energies.The reason is due to the presence of highly electronegative F and O atoms, Figure 1.
On the other hand, by considering the final configurations of the simulated systems (Figure 2), it was found that the drug fully loaded into the bCD and concurrently the vdW and Coul interactions energies of the water/drug were greatly reduced.
But finally, the loading of the drug into the bCD is supported and strengthened by the vdW interactions inside the bCD, so that the vdW and Coul interaction energies of the bCD/drug were equal to À 56.46 and À 6.04 kJ/mol,  respectively.Examining the reference system was shown that both vdW and Coul interaction energies (À 274.15 and À 262.40 kJ/mol, respectively) between bCD/water have a very significant contribution to the solubility of bCD.In the bCD þ 5-FU þ Water system, although the amount of both interaction energies has decreased (Coul and vdW energies between bCD/water were À 234.90 and À 238.10 kJ/mol, respectively) the amount of interaction energies remains significant, which shows that bCD keeps the tendency of solvation in water even after the drug has loaded in it, and its interactions with 5-FU and water were thermodynamically favorable.
In the Ald þ Water system, because of the large number of O and P atoms in the Ald structure as well as the structural symmetry of Ald, which increases the interactions, the thermodynamic tendency to solvate was considerable (Coul and vdW energies between Ald/Water were À 328.77 and À 12.63 kJ/mol, respectively).
Table 1 shows that the Coul interaction energy between Ald/Water in the bCD þ Ald þ Water system is still negative and significant (À 169.44 kJ/mol).
Examining the interaction energies between bCD/Ald demonstrated that vdW interaction energy was dominant (À 100.05 kJ/mol) and determined the thermodynamics of loading the drug into the bCD cavity, although the contribution of Coul interaction was also significant (À 36.49kJ/mol).It seems that the bulky O and P atoms in the Ald structure strengthen the vdW interactions and increase the thermodynamic tendency of this drug to enter the bCD cavity, and at the same time, the negative charge of O atoms strengthens the Coul interactions with bCD/Ald.The vdW interaction energy between bCD/water was shown to a further decrease compared to the same interactions in the bCD þ 5-Fu þ Water system, but the thermodynamic tendency of solvation in water was still preserved.
In the TMZ þ Water system, in comparison to two similar drug-containing systems, the Coul interaction energy between TMZ/Water (À 113.76 kJ/mol) was somewhat weaker due to less number of O atoms in the TMZ structure than two other drugs, as well as the absence of an electronegative atom such as F (as in 5-FU).On the other hand, due to the presence of numerous N atoms inside and outside the ring (Figure 2), the amount of vdW interaction was significant and far more (À 68.25 kJ/mol) than the similar case for Ald and 5-Fu.The TMZ/Water Coul and vdW interactions were decreased in the bCD þ TMZ þ water system (À 74.22 and À 20.03 kJ/mol, respectively), which was more severe for the vdW interaction energy.Such as two other drugs, bCD/ TMZ interaction was thermodynamically favorable, and the main contribution was related to the vdW interaction (À 103.83 kJ/mol).
Finally, it could be concluded that the thermodynamics of the system including 5-FU, Ald, and TMZ drugs in water was controlled by Coul interactions, which in the case of TMZ, the vdW contribution was also significant.But in the bCDcontaining systems, the loading of drugs into the bCD cavity was controlled and supported by the vdW interactions, which can be attributed to the presence of bulky atoms such as O and P in the structure of drugs and bCD itself.This strongly strengthens the possibility of the vdW interactions between bCD/drug.Also, in all systems, although the interaction energy between bCD/water decreases to some extent, the solvation of bCD in water was still thermodynamically favorable, and both Coul and vdW interactions have a significant contribution to the solubility of bCD in water.These results are in good agreement with the experimental findings so the voltammetric and spectroscopy studies showed the improvement of the TMZ solubility and formation of stable complexes.Phase solubility diagrams showed an increase in TMZ solubility and 1:1 stoichiometry of the complexes formed (Krzak & Bilewicz, 2020).Also, the weaker interaction energy between 5-Fu/bCD could be considered as thermodynamic evidence, following the experimental results, regarding the faster release of 5-Fu (Bilensoy et al., 2007).

Hydrogen bond
Considering the ability to establish hydrogen bonds (H-bond) between different components in the simulated systems, the competition in establishing H-bond can reflect the thermodynamic tendency of drugs to bCD, to load or solvation in water.Looking at the H-bond between water and the drug, it was found that all drugs show a lower H-bond with water in the presence of bCD, Table 2.This finding was in line with the reduction of water/drug interaction energy in the bCD þ drug þ water containing systems compared to the just water/drug containing systems, which was because of the tendency of drugs to interact with bCD through the penetration into the bCD cavity that keeps them away from water.Meanwhile, the quality of the H-bond (lifetime of the Hbond) was greatly increased because they do not have space to escape from the bCD cavity during the loading of drugs.By examining the number of water molecules in the range of 0-0.5 nm from the center of bCD in all the b-CD-containing systems, it can be seen that the number of water molecules were pushed out of the bCD cavity, and this amount was more noticeable in the case of 5-FU, which could be attributed to its small molecular structure and the creation of less spatial hindrance in the center of the cavity, as well as its greater hydrophilic property (caused by the presence of F and O atoms), and therefore it maintains more water molecules.While this has not occurred in the case of Ald and TMZ because they are both bulky and Ald is also expected to expel more water molecules, Ald despite its steric hindrance as well as larger structure, compared to TMZ, due to the presence of more hydrophilic atoms (P and O) expels fewer water molecules and keeps relatively more water molecules than TMZ in the center of the cavity, Table S2 and Figure 1.Then, by examining the number of H-bond between b-CD/drug (Table 2), it was revealed that the number of H-bond in the Ald-containing system was more than two other drugs (1.973) along with higher quality (lifetime 3.13 ps).
This result is supported by the experimental findings so that the specific host-guest complexations between Ald and bCD were easily detected by electrospray ionization mass spectrometry (ESI-MS) in the gas phase, which suggests that hydrogen bonds were even maintained upon vaporization and ionization (Biernacka et al., 2014).
The number of H-bond between bCD/bCD in both Ald and 5-FU-containing systems was not changed largely (6.19 and 6.11, respectively) compared with the reference system (6.04), but it was increased in the presence of TMZ (6.56).Since in the TMZ-containing system, the number of H-bond between bCD/bCD was increased and due to this fact that the interaction energy of bCD/water in the presence of TMZ was somewhat more positive than the other two drugs (Table 1), it can be concluded that the solubility of bCD in the bCD þ TMZ þ Water system decreases more than the two other systems.This increase in the number of H-bond between bCD/bCD is since among the atoms of TMZ structure there have been fewer suitable candidates for making H-bond with bCD, and therefore, it strengthens the possibility of hydrogen bonding between bCD/bCD.Finally, by studying the number of H-bond between bCD/water, it can be seen that in all the drug-containing systems, the number of H-bond was decreased in comparison with the reference system (19.558).However, in the 5-FU-containing system, the number of H-bond between bCD/water was slightly more (17.163)than in the two other drug-containing systems.The reason for the higher H-bond between bCD/water in the 5-FU-containing system could be related to the greater dipole moment of bCD in this system compared to the two other drug-containing systems (the dipole moment of bCD in the 5-FU, Ald, and TMZ-containing systems were 5.093, 4.868 and 4.918 respectively).Because of the rigid structure of bCD, the H-bond quality (lifetime) of bCD/water was very close to that of bCD/bCD in the bCD þ water system.

Distances and drug loading mechanism
One of the important dynamic characteristics of drug molecules in bCD is the loading mechanism in the bCD cavity (Mojdehi et al., 2021), to this end, initially, the distances between the center of mass (COM) drugs and bCD were obtained in different systems.As shown in Figure 3, 5-Fu loaded later than the two other drugs (after around 5 ns from the beginning of the simulation).
5-Fu has the lower molar mass and moment of inertia (6.49 a.m.u.nm 2 ) than the two other drugs (the value of the moment of inertia for TMZ and Ald was 16.72 and 21.05 a.m.u.nm 2 , respectively).Hence, 5-Fu was available to water molecules for a longer period of time and therefore, the possibility of displacement of 5-Fu was more than two other drugs at the first few nanoseconds, so it takes more time for the simulation system to reach equilibrium and the result of the interactions causes loading of the drug, so the average value of the distance of 5-Fu from the COM of bCD was about 0. 19 nm. Figure 2 shows that 5-Fu loading into the bCD was occurred from the F atom side and from the secondary rim toward the primary rim.In the following, it was seen that the Ald loaded about 3 ns later than TMZ, and average distances from the center of bCD were 0.17 and 0.12 nm, respectively.It could be attributed to the presence of bulky O and P atoms in the structure of Ald (Figure 1), which causes proper interaction of these atoms (interaction of hydrophilic atoms) with water, and also was observed that Ald loaded from the O and P atoms side into the middle of the bCD cavity, and its hydrocarbon chain (which includes the amino group) was located out of the secondary hydroxyl rim of bCD towards the water bulk.Also, the vdW repulsions between the O and P atoms in the Ald structure with the O atoms in the middle of the bCD structure (O1) make the drug unable to get approach the primary hydroxyl rim of bCD and be fully loaded, this in turn results in Ald being further away from the COM of bCD compared to the two other drugs.Finally, due to the lack of hydrophilic atoms in the TMZ structure as well as the larger hydrophobicity than two other drugs, TMZ was loaded into the bCD cavity (which has hydrophobic properties) much closer and at a shorter  distance of 0.12 nm.The orientation of N, O, and C atoms towards of the primary hydroxyl rim and a shorter distance from the center of bCD causes complete and faster loading of TMZ than the two other drugs.The probability distribution of the distances, Figure 4, showed that the maximum possible values were very close to the average distances of the drugs from the COM of bCD which indicates the fact that the distribution around the average values was significant.As a result, the average distance values had the most repetitions in the simulation and were statistically acceptable and reliable.Finally, the position of drugs in the equilibrium distance from the COM of bCD during the time intervals when the system was in equilibrium is well shown in Figure S1.
To better understand the effect of drug loading, the analysis of several radial distribution functions (RDF), as well as the number of water molecules in the spheres with different radii from the center of the bCD cavity were considered and found that loading of all drugs was accompanied with the outflow of water molecules, as seen in Figure 5 and Table S2.
By considering the RDF, it is seen that the presence of water molecules in the center of bCD was quite evident for the reference system (the presence of a significant peak at a distance of �0.2 nm) and was in good agreement with the results of another simulation study (Gebhardt et al., 2018).As can be seen, when the drugs were entered into the bCD cavity, the first peak disappeared and the presence of water molecules were significantly reduced so that no peak was observed in the distance of 0-0.9 nm for the drug-containing systems.This finding was consistent with the number of water molecules in the spheres made in the center of the bCD cavity.The number of water molecules in the reference system between 0-0.5 nm from center of the bCD cavity was �9 molecules, which with the entry of drug molecules, a significant number of water molecules pushed out of the cavity.Thus, the greatest outflow of water molecules was due to the loading of TMZ (�7.5), and the least was related to the loading of 5-Fu (�5.5).The reason that 5-Fu expels a smaller number of water molecules could be related to its small structure, which provides a suitable space for water molecules inside the bCD cavity.Because of the larger molecular structure, Ald was placed in the middle of the bCD cavity and could not penetrate throughout the bCD cavity.Thus, Ald provided more free volume of the bCD cavity for water molecules.For TMZ, as shown in Figure 2, the bCD cavity was fully occupied by TMZ and reduced the presence of water molecules.By studying the total amount of Coul and vdW interaction energies between bCD/water in the reference system and bCD/drug/water-containing systems, it was found that the thermodynamics of the system were proportional to the number of water molecules in the innermost of the bCD cavity.In all the drug-containing systems, the sum of the Coul and vdW interaction energies of bCD/water became somewhat more positive (À 473, À 461, and À 452 kJ/ mol for 5-Fu, Ald, and TMZ, respectively) than the reference system (À 536.55 kJ/mol), but they were still enough significant negative.

Mean square displacement (MSD)
To study the movement of drugs and bCD in various simulated systems, the mean square displacement (MSD) was analyzed, as shown in Figures 6 and 7. Comparing the movement of drugs shows that the movement in all bCD-containing systems was much lower than in the water/drug-containing systems, which is caused by the loading of the drugs inside bCD, which severely limits the movement of the drugs.In the meanwhile, in the water/drug-containing systems, the movement of drugs follows their moment of inertia (as mentioned in the previous section).So that the largest movement was related to 5-Fu and the least movement was related to Ald.
When examining the movement of drugs in the bCD-containing systems, the movement was observed to be greatly reduced.Therefore, when the drugs were loaded into the bCD cavity, the movement was no longer a function of the moment of inertia or even influenced by the thermodynamic properties of their interaction energies with water, but it was caused by the fact that the drugs were loaded into the bCD cavity and drugs were out of reach of the water molecules.Hence, the movements of drugs were greatly reduced.
Figure 7 shows the movement of bCD in different systems.As one can see, the movement of bCD in the presence of 5-FU decreased more than the two other drugs.Thermodynamic results (Table 1) confirm the movement behavior of bCD in the presence of drugs.So that the sum of the Coul and vdW interaction energies between bCD/water in the presence of 5-FU was more negative than the two other drugs, it had less movement due to a more favorable thermodynamic tendency to solve in water.The interaction energies of b-CD/water in the presence of Ald are higher than TMZ, and it seems that the closeness of their interaction energies makes the movement behavior of bCD in the presence of these two drugs almost similar and close to each other.

Flexibility
The root mean square fluctuations (RMSF) are defined as the fluctuation of every single atom around its average position.
When a simulation is equilibrated, that is, when the structure of interest fluctuates around a stable average conformation, it makes sense to compute the fluctuations of each subset of the structure (each atom, for example) relative to the average structure of the simulation, the RMSF (Mart� ınez, 2015).

Root mean square fluctuations (RMSF) of drugs.
Figure 8 illustrates the flexibility of each drug atom in different simulated systems.
The flexibility of 5-Fu in both 5-Fu-containing systems was similar and did not change significantly due to being only one small rigid ring with F and O atoms.
Although TMZ has two 6-and 5-sided rings connected and rigid, due to the existence of a linear part including an amino group, it is expected to have some flexibility.But as one can see, whether in water or in the presence of b-CD, there were no significant changes in the flexibility of atoms.There is a large linear part in the molecular structure of Ald (atoms 1 to 7).Ald's atoms were more flexible than the two other drugs, but it still maintains their rigidity to a large extent.Examining the Ald's atoms shows that atoms number 1, 2, and 3, which correspond to the N and H atoms of the amino chain, have been exposed to water bulk more than others and had more flexibility than others in both Ald-containing systems.Changing the flexibility of atom number 8, which is connected to the quaternary carbon as well as the symmetrical phosphate groups on both sides, was significant to some extent.In the symmetrical phosphate groups, atoms number 11, 12, 14, and 17 which correspond to the oxygens connected to the P show the most changes in flexibility.According to the loading mechanism of drugs and by taking into account the distance from the center of the bCD cavity as well as the placement of the Ald inside the bCD cavity that occurred from the side of the secondary rim, it can be seen that the bulky atoms of the phosphate groups (due to being hydrophilic) during loading inside the secondary rim affected by contact with water molecules.Since Ald has  been loaded into the bCD cavity from the side of phosphorus and oxygens atoms, it seems that these atoms need a rotation and rearrangement, and most changes were related to the flexibility of these atoms.

Root mean square fluctuations (RMSF) of b-CD.
The flexibility of bCD was studied based on the type of bCD monomer atoms (Figure c1).Two separate groups, including the oxygen atom group and the carbon atom group, were evaluated separately, Figure 9.In oxygen atoms, the greatest flexibility was associated with O6 and in carbon atoms was C6.However, the O6 and C6 display slight changes in the different systems.The C6 atoms have been connected to the O6 atoms from the top side and to the 6-membered ring on the bottom side (Figure c1).Hence, the flexibility of C6 atoms in different systems was almost unchanged.The C6 atoms have been located on the primary hydroxyl rim of the bCD cavity (the region of the bCD where O6 atoms were located) and have been exposed the most to the water molecules and loaded drugs, therefore, had shown the greatest flexibility.However, the lack of significant changes in the flexibility of these atoms in the different systems indicates the fact that these atoms are not affected by drugs, and their flexibility was more related to the nature of bCD itself.The same behavior was seen in the case of O6.Even though there was space for movement, the O6 atoms did not show a noticeable change.Examining the secondary hydroxyl rim (the region of bCD where O2 and O3 atoms were located), the most changes in the flexibility were seen upon the O2 and O3 atoms.This could be attributed to the Coul and vdW repulsions between the oxygen atoms in the secondary hydroxyl rim (O2 and O3) as well as the tendency of drugs to enter into the bCD cavity during loading, which could be associated with the flexibility changes and displacement of the oxygen atoms in the secondary hydroxyl rim.Ald contains two bulky P and O atoms that have a large vdW radius and when it enters through the secondary rim causing a  strong repulsion with O2 and O3 atoms.Therefore, by increasing the mobility of these atoms, flexibility is increased.The changes caused by the presence of Ald were more than the two other drugs.After Ald, it seems that the O2 and O3 oxygen atoms had the most flexibility changes in presence of TMZ.TMZ has two rings and two bulky N and O atoms.The amino chain was fully loaded into the bCD cavity, and the flexibility changes of O2 and O3 atoms were more than the reference system (bCD in water).Finally, 5-FU was unable to make significant changes in O2 and O3 flexibility due to its small structure that has only one ring during loading.Investigations in the middle of the bCD cavity show that O1 atoms had the greatest flexibility in the presence of TMZ.Because TMZ was completely loaded in the bCD cavity and affected the flexibility of O1.But Ald was loaded over the bCD cavity and did not fill the bCD cavity, so it had less effective than TMZ to change the flexibility of O1.
Eventually, 5-FU has the least effect due to its small structure.In the analysis of the ring that contains the C1-C5 carbon atoms, C5 is very rigid because it has located on the side of the primary rim.On the other hand, it has been fixed by the O6 and C4 atoms of the ring and is not shown special flexibility changes.Its flexibility could be related only to the nature of bCD itself, as the previous study confirms this finding (Mixcoha et al., 2014).
C2 and C3 atoms, which connected to the secondary rim of bCD through O2 and O3 atoms, during loading when the drugs enter the cavity exposed to strong vdW repulsion of drug atoms on the secondary hydroxyl rim of b-CD.The effects of these interactions could be seen especially at the moment of loading and entering Ald on C2, C3, O2, and O3 atoms, which causes more flexibility and displacement of these atoms.C1, C4, and O1 atoms, which are located in the middle of the bCD structure, were investigated and found that the mentioned atoms were affected by the strong repulsion caused by entering the drug during loading and made them more flexible.The presence of TMZ and Ald was noticeably changing the flexibility, but 5-FU did not cause a specific change.

Stability and structural changes
To obtain information about the stability of the structures during the simulation time, the amount of structural drift of bCD as well as the drugs using the Root mean square deviation (RMSD) parameter was considered.The RMSDs (bCD and drugs) from the start of the simulation (t ¼ 0), were shown in Figures S2 and S3.As can be seen, the average RMSD value for bCD in the reference system was �0.11 nm, which was in good agreement with the previous MD simulation study (Mixcoha et al., 2014).Examining all the simulated systems shows that after a short time from the beginning of the simulation, the RMSD value approaches an average value, and around this point, it shows slight fluctuations during the simulation time.This indicates that the simulation systems were well equilibrated.Also, in the case of the bCD, the average RMSD values were close to each other in different simulated systems so it confirms that the bCD structure is rigid in all systems.
As mentioned earlier, bCD is a cyclic oligosaccharide, which is observed in the form of a short cone based on the orientation of the hydroxyl groups in its primary and secondary hydroxyl rims (Khuntawee et al., 2015).bCD has particular importance due to its rim size and special use in biological fields, especially as mechanism (host-guest).So, a more detailed examination of the surface of the bCD rims during the entry of guest molecules (drugs) could provide a deeper insight into the understanding of the mechanism as well as the structural changes of bCD (Mojdehi et al., 2021).Also, the investigation of these findings may be related to the mechanism of drug release, in which the open and closed states of the secondary rim of the bCD molecule may play an important role at different stages (Nutho et al., 2014).For this purpose, the area of rims associated with the primary and secondary hydroxyl rim was examined and analyzed.It was observed that the area associated with the primary hydroxyl rim in the reference system had a numerical average of 0.9 nm 2 and with the entry of drugs shows a sharp decrease, Figure 10.The greatest reduction was related to the 5-FU because the penetration of the 5-FU was from the side of the F atom and causes complete penetration of this drug to the end of the primary hydroxyl rim, (as previously in the RMSF and loading mechanism section stated).It seems that the atoms of 5-FU have an effective attractive interaction with the O6 hydroxyl group and caused the area of this rim to become smaller (0.75 nm 2 ).Also, TMZ has affected this area through the favorable interaction of N and O atoms with the O6 hydroxyl group (0.77 nm 2 ).Finally, the bulky P and O atoms in the molecular structure of Ald, as mentioned before, prevented the complete penetration of this drug into the bCD cavity.So, although it has reduced the area of this rim, it is due to fewer effective interactions of the atoms with hydroxyl groups of O6 causing a decrease and less change in this area compared to the presence of two other drugs, so the average value of 0.79 nm 2 was obtained.
In further investigations, the surface of the secondary hydroxyl rim was not shown a significant changes when drugs enters into the bCD cavity.So, the average value of this area for the reference system was 1.32 nm 2 , which decreases to an average of 1.28 nm 2 when drugs enters into the bCD cavity.Although, there was a decrease, this reduction is not noticeable and is almost the same for all drugs.Finally, by examining the middle area that includes O1 atoms in the bCD, it can be seen that unlike the primary and secondary hydroxyl rim areas (which were more flexible due to the greater freedom of their hydroxyl groups).Due to the strong bonds of O1 hydroxyl groups, the middle area of bCD in different simulated systems was not changed, and the average value (0.83 nm 2 ) of this area for all the simulated systems was almost the same.

Conclusion
In this study, the loading of anti-cancer drugs in the bCD was evaluated from the dynamic and thermodynamic point of view.The results of the thermodynamic analysis showed that all three drugs had a thermodynamic affinity to the bCD.Both vdW and electrostatic interactions were effective in the bCD solubility in water, in the presence and absence of drugs, while the electrostatic interaction between drug/ water was mostly dominant.In all drug-containing systems, hydrogen bonds between water/bCD were somewhat reduced but still significant.In general, due to the negative interaction energies and the establishment of significant hydrogen bonds, the bCD could be considered a suitable carrier for these three drugs (especially Ald).Examining the loading of drugs illustrated that TMZ was loaded earlier and 5-Fu later than the two other drugs.Loading of Ald took place in the middle of the bCD cavity and its amino group was placed outside and toward the water bulk, while 5-Fu and TMZ penetrated to the end of the primary hydroxyl rim.The analysis of the RDFs showed that the loading of drugs into the bCD cavity was accompanied by the exit of water molecules from the center of the bCD cavity.The movement of drugs in the system without the bCD is followed by their moment of inertia, while after loading inside the bCD cavity, it was completely independent of it.Most of the changes in the flexibility of bCD atoms were caused by the presence of Ald, and among these, atoms C1-C4 and O2-O3 showed the most flexibility.

Figure 2 .
Figure 2. Snapshots of the initial structures (a), the final structure from the side view (b), and the final structure from the top view (c) of the simulated systems.

Table 1 .
Analysis of interaction energies in different simulated systems.

Figure 3 .
Figure 3. Average distance of TMZ, Ald and 5-FU drugs from b-CD in different simulated systems.

Figure 4 .
Figure 4. Probability distribution of distance of TMZ, Ald and 5-FU drugs from b-CD in different simulated systems.

Figure 5 .
Figure 5. Radial distribution function (RDF) of Oxygen of water in different simulated and reference systems.

Figure 6 .
Figure 6.Mean square displacement of drugs in different simulated systems.Figure 7. Mean square displacement of b-CD in different simulated systems and reference system.

Figure 7 .
Figure 6.Mean square displacement of drugs in different simulated systems.Figure 7. Mean square displacement of b-CD in different simulated systems and reference system.

Figure 8 .
Figure 8. Root mean square fluctuations of drugs in different simulated and reference systems.

Figure 9 .
Figure 9. Root mean square fluctuations of b-CD in simulated systems.

Figure 10 .
Figure 10.Small aperture surface area in different simulated and reference systems.

Table 2 .
Lifetime and number of hydrogen bonds in different simulated systems.