Diffusivity of α-, β-, γ-cyclodextrin and the inclusion complex of β-cyclodextrin: Ibuprofen in aqueous solutions; A molecular dynamics simulation study

Cyclodextrins (CDs) are widely used in drug delivery, catalysis, food and separation processes. In this work, a comprehensive simulation study on the diffusion of the native α-, β- and γ-CDs in aqueous solutions is carried out using Molecular Dynamics simulations. The effect of the system size on the computed self-diffusivity is investigated and it is found that the required correction can be as much as 75% of the final value. The effect of the water force field is examined and it is shown that the q4md-CD/TIP4P/2005 force field combination predicts the experimentally measured self-diffusion coefficients of CDs very accurately. The self-diffusion coefficients of the three native CDs were also computed in aqueous-NaCl solutions using the Joung and Cheatham (JC) and the Madrid-2019 force fields. It is found that Na+ ions have higher affinity towards the CDs when the JC force field is used and for this reason the predicted diffusivity of CDs is lower compared to simulations using the Madrid-2019 force field. As a model system for drug delivery and waste-water treatment applications, the diffusion of the β-CD:Ibuprofen inclusion complex in water is studied. In agreement with experiments for similar components, it is shown that the inclusion complex and the free β-CD have almost equal self-diffusion coefficients. Our analysis revealed that this is most likely caused by the almost full inclusion of the ibuprofen in the cavity of the β-CD. Our findings show that Molecular Dynamics simulation can be used to provide reasonable diffusivity predictions, and to obtain molecular-level understanding useful for industrial applications of CDs.

knowledge, there are only two studies reporting diffusion coefficients of CDs calculated using MD simulations to date. Naidoo et al. [48] performed pulse-field-gradient spin-echo nuclear magnetic resonance (PGSE NMR) measurements and MD simulations to determine the diffusion coefficients of α-, β-, and γ -CDs. In that study, the CSFF and SPC/E force fields for the CDs and water were used, respectively. The diffusion coefficients of the CDs were measured at 1 mM concentration in D 2 O. Naidoo et al. [48] concluded that considering the different solvents used in the simulations and experiments, the simulation results are in a good agreement with the measurements. Tang and Chang [37] presented an MD simulation approach for calculating the binding free energy of host-guest complexes based on the diffusion coefficients of the host and the guest. Using this approach, Tang and Chang computed binding free energies of seven different guest molecules with β-CD using the q4md-CD/GAFF/TIP3P combination of force fields.
In this work, a comprehensive study of the diffusivity of CDs in aqueous solutions is reported. Such a study is largely lacking and is important for not only supplementing and further guiding experiments, but also for providing the necessary molecular understanding for the design and/or optimization of applications in which mass transport is crucial, e.g., controlled drug delivery. As it has been shown that the effect of the system size on the computation of self-diffusivity can be significant in MD simulations [49,50] , the effect of the system size on the self-diffusion of CDs is investigated. To model CDs, we use the q4md-CD force field which has been shown to be very accurate in predicting the crystal properties of CDs [51] . Since it has been shown in literature that the choice of force field to represent the solvent has a significant effect on the computed transport properties of the solutes [52][53][54] , the diffusion coefficient of CDs are calculated using four different water force fields. After identifying the best performing solute (i.e. water) force field, the effect of NaCl concentration on the diffusion of CDs is also investigated using two different NaCl force fields. Finally, as a model system for a drug-delivery application, the diffusion coefficient of the ibuprofen: β-CD inclusion complex is computed.
The remainder of this article is structured as follows. The following section contains details about the force fields and the MD simulation scheme used. In Section 3 (i.e., Results and Discussion), the effect of the system size is investigated, the computed selfdiffusion coefficients of CDs in different aqueous solutions are presented, and the diffusion of the ibuprofen: β-CD inclusion complex is discussed. The conclusions of this work are provided in Section 4.

Force fields
All bonded and non-bonded force field parameters of CDs were taken from the q4md-CD force field [51] . This force field is shown to be very accurate in reproducing the experimentally measured structure of CDs, as well as the binding properties of CDs with a wide variety of guest molecules [33,34,51,55] . Water was modelled using TIP3P [56] , SPC/E [57] and TIP4P/2005 [58] force fields, as well as Bind3P [34] . The latter has been recently developed by Yin et al. [34] by reparameterizing the Lennard-Jones size and energy parameters of the TIP3P water model. The Bind3P force field is fitted to experimentally measured binding free energy and enthalpy of CD with various guest molecules [34] . For all waterion simulations, TIP4P/2005 [58] was used combined with the Joung and Cheatham [59,60] (JC) and the Madrid-2019 force fields [61] for Na + and Cl − . This choice was based on the recent study by Döpke et al. [62] who investigated the performance of different ion force fields combined with the TIP4P/2005 force field in reproducing experimental properties of salt solutions such as diffusion coefficients, ion hydration free energies and hydration radius. The JC force field has full ionic charges (i.e., +1, −1) while the Madrid-2019 force field has scaled atomic charges (i.e., +0.85, −0.85). Moreover, the Madrid-2019 [61] force field uses specific terms for the Lennard-Jones interactions between the ions and water molecules instead of the general combining rules. Ibuprofen is represented as a fully flexible molecule using the General AM-BER force field [63] . The partial charges of Ibuprofen are computed with the Restrained Electrostatic Potential (RESP) method at the 6-31G * level of theory, using the R.E.D III.52 script [64] with the Gaussian09 RevB.01 software package [65] . To check the reliability of the derived charges, the electrostatic potentials (ESPs) generated by the fitted charges and DFT at the 6-31G * level are compared. The relative root mean square between the ESPs in atomic units are found to be smaller than 0.02. The derived charges correspond to the simulation of a molecule at 0 K in vacuum. This approach is used since it is the recommended charge derivation method for the GAFF force field [63] . Recently, Schauperl et al. [66] reported the RESP2 method in which the partial charges are calculated from a combination of gas-and aqueous-phase charges. Schauperl et al. [66] showed that RESP2 charges are more accurate than the traditional RESP in reproducing experimentally measured properties (e.g. hydration free energy, heats of vaporization) of organic liquids [66] . The Lennard-Jones and electrostatic interactions are truncated at 9 Å . Analytic tail corrections are included in the pressure and energy calculations. The Particle-Mesh Ewald (PME) method [67] is used for the long range electrostatic interactions. In all simulations, fixed atomic charges are used and polarization effects are neglected according to the original parametrization of the used force fields. All force field parameters are provided as Supplementary material in the form of GROMACS parameter files.

Simulation details
MD simulations using the GROMACS 2018.2 [68] software package were performed for the diffusivity computations, while LAMMPS (version 16 Feb. 2016) [69] was used for computing viscosities. To integrate the equation of motion, the leap-frog algo-rithm with 2 fs timestep is used. The initial configurations are obtained by solvating a CD molecule in water. For the aqueous electrolyte systems, additional water molecules are used for the solvation of the CD molecule, and afterwards these are randomly replaced by the appropriate number of Na + and Cl − ions. All initial configurations are created using the built-in GROMACS tools [68] . The equilibration scheme used is as follows: (1) energy minimization of the system using the steepest descent method, (2) 50 ps MD simulation at 50 K at constant volume, (3) 1 ns MD simulation in which the temperature is ramped up from 50 K to the production temperature in the NVT ensemble, and (4) 2 ns MD simulation in the NPT ensemble at the production temperature and pressure. This thorough equilibration scheme is based on the inclusion complexation studies by Gilson and co-workers [35,36] . For consistency, this approach is used in all our simulations. The temperature is regulated using the modified velocity rescaling thermostat of GROMACS [70] , which is able to properly reproduce the canonical ensemble. The coupling constant used in the thermostat is 0.1 ps. The pressure was regulated using the Parrinello-Rahman barostat [71] with a coupling constant of 2 ps. In all simulations, unless stated otherwise, a CD molecule is solvated in 30 0 0 water molecules.
Self-diffusion coefficients, D i ,self , are calculated using the Einstein relation [72] : where t is the time, N i is the number of molecules of type i , and r j,i is the position vector of the j th molecule of species i . The brackets ... represent an ensemble average. D i ,self is computed by the slope and intercept of the long-time mean squared displacement (MSD) curve at time-scales where the slope of MSD as a function of time is 1 in a log-log plot [72,73] . Because of the use of one solute molecule and the slow dynamics of CDs due to their size, long MD simulations in the range of 600 ns -1 μs were needed for obtaining a linear regime with slope equal to 1 in the log-log MSD vs time curve. The self-diffusion coefficients reported in this work are calculated as the average of four independent simulations starting from different initial configurations. Each initial configuration was created by carrying out the equilibration procedure using different random seeds. Each MD simulation was carried out on 12 -24 CPUs, achieving a performance of 100-160 ns/day. Self-diffusion coefficients computed from MD simulations depend on the size of the simulation box [40,[74][75][76][77] . Yeh and Hummer [74] derived an analytic correction term which should be added to the self-diffusivities computed in MD simulations to obtain the self-diffusion coefficients in the thermodynamic limit [74] ( D ∞ self ). The YH correction term is: where ξ is a constant equal to 2.837297, k B is the Boltzmann constant, T is the temperature, L is the size of the simulation box, and η is the shear viscosity of the solvent computed from MD simulations. It is important to note here that viscosities exhibit no finite-size effects [76] . The viscosity of water and of aqueous electrolyte solutions required for correcting our diffusivity results are not available for all the force fields used and all conditions considered in the current study. To this purpose, the shear viscosities of all systems are calculated using the OCTP plugin [73] in LAMMPS, an easy-to-use tool which computes transport properties on-thefly during the simulation. The reported viscosities are calculated as the average of four independent 20 ns long simulations. All details about the OCTP plugin can be found in the original paper [73] . All self-diffusivities reported in this work are corrected for system size effects using Eq. (2) , unless it is stated otherwise. The computed viscosities for all systems are listed in Tables S1 and S2 of the Supplementary material. In this work, the statistical error in the self-diffusion coefficients and viscosities are calculated as the standard deviation of the mean from the four independent simulations.

System size effects
The effect of the system size on the diffusivity of pure fluids or mixtures computed with MD simulations has been investigated in detail by several authors [40,49,74,75,78] . In Fig. 2 , the computed self-diffusion coefficients of β-CD with and without the YH correction at 298.15 K are shown as a function of the simulation box length. The TIP4P/2005 water force field is used. As expected, the finite size self-diffusivities scale linearly with the size of the simulation box, and thus, the application of the YH correction yields the self-diffusivity value in the thermodynamic limit. However, it is important to note that the magnitude of the required correction is substantial. As shown in Fig. 2 , by using 1500 water molecules the computed self-diffusion coefficient (i.e. 1.9 ×10 −10 m 2 s −1 ) is almost equal to the required correction (i.e. 1.8 ×10 −10 m 2 s −1 ). For smaller systems the correction is even higher. For a system size of 500 water molecules, the YH correction is approximately 76% of the final self-diffusion coefficient value. In the study by Moultos et al. [49] , the effect of the system size on the self-diffusion coefficient of pure CO 2 and CH 4 was studied. The authors show that the magnitude of the YH correction for these systems ranges from ca. 6 to 15% of the final self-diffusion coefficient, depending on the system size. For the systems investigated in this study, in which the size of the solute is considerably larger than the size of the solvent (i.e. CDs in water), correcting the self-diffusion coefficient of the solute for system size effects is absolutely crucial because the correction term can be in many cases larger than the computed diffusivity value in MD. Similarly to the results obtained for the β-CD, the required system size correction for the αand γ -CDs is also substantial as shown in Figures S1 and S2 of the Supplementary material. Hereafter, all of the reported self-diffusion coefficients in tables and figures are corrected for system size effects using using Eq. (2) .  [56] , Bind3P [34] , SPCE [57] , and TIP4P2005 [58] force fields for water, respectively. The self-diffusivities are corrected for finite size effects using Eq. (2) . The experimentally measured self-diffusion coefficients are shown with open black symbols [79][80][81] . The maximum statistical error in the computed self-diffusion coefficients is 10%. The error bars are smaller than the symbol size.

Effect of the force field on the diffusivity of cyclodextrins in water
In Fig. 3 , the computed and experimentally measured selfdiffusion coefficients of the α-, β-, and γ -CDs in water are shown as a function of temperature for the range 298.15-312.15 K at a pressure equal to 1 bar. The experimental results are taken from the studies of Ribeiro et al. [79][80][81] . As can be seen, by using the SPC/E and TIP4P/2005 water models the computed selfdiffusion coefficients are in a reasonable agreement with the experimental results. The self-diffusion coefficients predicted using the TIP4P/2005 water force field are the most accurate, showing a mean average deviation from the experimental measurements equal to ca. 8 %. From Fig. 3 , it can be also observed that when the TIP3P and Bind3P water force fields are used, the predicted diffusivities are higher than the experimental measurements by a factor of 3 to 4, for the whole temperature range. In the case of TIP3P water model, this overestimation is expected since the predicted diffusivity and viscosity of pure water with this model are substantially over-and underestimated, respectively [82] . The Bind3P force field was recently developed by Yin et al. [34] , who refitted the Lennard-Jones size and energy parameters of the TIP3P water model to better predict the experimentally measured binding free energies and enthalpies of host-guest systems. It was shown that compared to TIP3P, the Bind3P force field systematically improves the prediction of the experimental binding free energies and enthalpies of CDs with different guest molecules [34] . To the best of our knowledge, transport properties for the Bind3P water model are not reported in literature, and thus it is not known how accurately this new force field is in predicting transport properties of pure water. In Table S1 of the Supplementary material, it can be seen that the viscosity of pure water calculated with TIP3P is Fig. 4. The average number of hydrogen bonds (HBs) formed between the (a) α-CD, (b) β-CD and (c) γ -CD with water as a function of temperature at P = 1 bar. The symbols in blue, orange, green, and red colors indicate the use of TIP3P [56] , Bind3P [34] , SPC/E [57] , and TIP4P/2005 [58] water force fields, respectively. The maximum statistical error in the computed number of HBs is 10%. The error bars are smaller than the symbol size. on average 12% higher than with Bind3P. This relatively small difference in viscosities justifies the similar poor performance of the two force fields in predicting the diffusivity of the CDs in water. In general, it is evident that the predicted self-diffusion coefficients of the CDs strongly depend on the accuracy of the force field used for modelling the solvent (i.e. water). This finding is in-line with the literature reporting MD simulations of self-diffusivity in aqeuous solutions [49,52] .
To obtain a better understanding of the effect of the water force fields on the physical properties of the CDs which can influence the diffusion, the average number of water molecules residing in the cavity of the CD during the diffusion process, the number of intra-and intermolecular hydrogen bonds between the CD and water, and the size of the small and large rims of the CDs are studied. The rims refer to the two openings in the CD molecule as can be seen in Fig. 1 b. In Fig. 4 , the average number of hydrogen bonds formed between the water molecules and the α-, β-, and γ -CD are shown as a function of temperature for the four water force fields used. The hydrogen bond analysis is carried out by using the built-in hbond GROMACS tool. In all simulations, the criterion for the formation of a hydrogen bond is a cut-off distance of 3.5 Å between the donor and acceptor atoms, and a cut-off angle of 30 • between the donor-hydrogen-acceptor atoms [51,83,84] . The difference in the average number of intermolecular hydrogen bonds formed using the different force fields is lower or equal to one in all cases investigated. In Figure S3 of the Supplementary material, the average number of intramolecular hydrogen bonds between the hydroxyl groups of the α-, β-, and γ -CD is shown when using the four different water force fields. Similarly to the intermolecular hydrogen bonds, there is no significant difference observed. Thus, the significant difference between the self-diffusion coefficients computed by using the different water models is not reflected in the number of inter-or intramolecular hydrogen bonds in the system.
In Figure S4 of the Supplementary material, the average number of water molecules residing inside the cavity of the α-, β-, and γ -CDs are shown as a function of temperature. Since the number of water molecules inside the cavity of the CDs continuously changes as the CD diffuses (i.e., water molecules continuously entering and leaving the cavity), the number of water molecules is calculated as a time average. In the studies by Cézard et al. [51] , Rodriguez et al. [85] , and Shikata et al. [86] the residence time of water molecules in the CD cavity was investigated by MD simulations and experimental techniques. Both the computed and the measured residence times are in the range of 20 -75 ps. This shows that the water molecules are continuously flowing in and out of the CD cavity with no stable water clusters being formed. In our study, a water molecule is considered to be inside the cavity of the α-, β-, and γ -CD if the distance between the centers of mass of the CD and the water is smaller than 0.34, 0.4, and 0.46 nm, respectively. This definition is adapted from the work of Zhang et al. [83] . Similarly to our findings for hydrogen-bonding, the average number of water molecules inside the cavity of the CDs predicted by using the different water models is almost identical. In Figures S5 and S6 of the Supplementary material, the diameter of the upper and lower rim of the α-, βand γ -CDs are shown as a function of temperature.
The rim sizes of the CDs do not show any significant difference for the different tested force fields. From the analysis presented above, it is evident that the shape and size of the CDs is almost unaffected by the choice of the force field of the solvent. Moreover, the large differences in the self-diffusion coefficients predicted by using different water models (see Fig. 3 ) are not reflected in the number of hydrogen bonds  [59] and Madrid-2019 [61] NaCl force fields, respectively. For water the TIP4P2005 [58] water force field is used. The selfdiffusivities are corrected for finite size effects using Eq. (2) . The maximum statistical error in the computed self-diffusion coefficients is 10%. The error bars are smaller than the symbol size.
(see Fig. 4 ) or in the number of water molecules residing in the CD cavity (see Figure S4 of the Supplementary material). These findings strongly suggest that the observed differences in the selfdiffusion coefficients of the CDs are mainly caused by the ability of the various water force fields to accurately predict the density and transport coefficients of the pure solvent. Thus, it is not surprising that the force field combination q4md-CD/TIP4P/2005 yields the most accurate results for the diffusivity of CDs in water, since TIP4P/2005 predicts the relevant properties of pure water very accurately [58,82] .

Effect of NaCl concentration
Since in most applications CDs are dissolved in aqueous solutions containing salts, we investigate the effect of NaCl concentration in water on the diffusion of CDs using two different ion force fields (i.e. JC and Madrid-2019) combined with the TIP4P/2005 water force field. The TIP4P/2005 force field was chosen since it was found to be the most accurate one combined with the q4md-CD force field (see Fig. 3 ). The two ion force fields are selected based on the study by Döpke et al. [62] who showed that the JC and Madrid-2019 ion force fields combined with TIP4P/2005 can reproduce properties including ion hydration free energy, hydration radius reasonably well. In Fig. 5 , the computed self-diffusion coefficients of the α-, β-, and γ -CDs in the aqueous electrolyte solutions are shown as a function of the NaCl molarity from 0 up to 1 mol L −1 at 298.15 K. The diffusivity of all CDs decreases with the increase in NaCl concentration, regardless of the ion force field used. This is most likely caused by the preferential interaction of the ions with the hydroxyl groups of the CDs altering their hydration shell. Moreover, the viscosity of the solvent increases with the NaCl concentration (see Table S2 of the Supplementary Material), which also contributes to the slower diffusion of the CDs. In Fig. 5 , it can also be observed that for NaCl molarities up to 0.4 mol L −1 , both ion force fields yield the same self-diffusivities (within the statistical uncertainty) for all three CDs. For higher ion concentrations (i.e., M NaCl > 0.6 mol L −1 ), the computed self-diffusion coefficients of αand β-CD using the two different ion force fields start to deviate from each other, showing that simulations using the Madrid-2019 force field always predict higher values. The difference between the two force field combinations becomes smaller as the size of the CD molecule increases (i.e from α to γ ). In the case of γ -CD, both force fields practically predicts the same diffusion coefficients (see Fig. 5 c).
To obtain a better understanding of the cause for these differences, the radial distribution functions (RDFs) of the Na + and Cl − ions in respect to the center of mass of the CDs are calculated. Since the highest difference in the diffusion coefficients occurs at the highest NaCl concentration, we calculate the RDFs for the simulations at M NaCl = 1.0 mol L −1 . In Fig. 6 a, b, and c, the radial distribution function of the α-, β-, and γ -CDs with the Cl − ions are shown for the JC and Madrid-2019 force fields for M NaCl = 1.0 mol L −1 . JC Cl − ions have a slightly higher affinity to interact with the CDs, expressed by a higher first peak on the RDF, than the ions modelled using the Madrid-2019 force field. Fig. 6 b and c show that for the βand γ -CDs a small peak appears at very short distances (ca. at d = 1.5 Å ) if the Madrid-2019 force field is used. These peaks suggest that, although it is rare, the Cl − ions can enter into the cavity of the CDs during the simulation. In Fig. 6 d, e, and f, the radial distribution function of α-, βand γ -CDs with Na + ions are shown for the two force field combinations, at M NaCl = 1.0 mol L −1 . Na + ions modelled with the JC force field have a higher affinity to interact with the CDs compared with the Madrid-2019 Na + ions. The difference between the height of the peaks in the RDFs is decreasing as the CDs get larger. This means that the affinity of Na + ions to interact with the CDs becomes more similar. These   [58] force field was used to model water. The green, orange, and blue colors indicate the Cl − , Na + and water, respectively. findings explain the differences in the computed self-diffusivities caused by the size of the CD and the use of the different ion force fields shown in Fig. 5 .
In Fig. 7 , density contour plots of water, sodium, and chloride around the α-, β-, and γ -CDs are shown. The contour surfaces show the locations where the density of the compound is equal to the set contour level. The contour levels of water, Na + , and Cl − are 1.05, 1.5 and 1.3 times the bulk densities of the corresponding compounds, respectively. In-line with the RDFs shown in Fig. 6 , Na + ions modelled with the JC force field (see Fig. 7 a, c and e) oc- Table 1 Computed self-diffusion coefficients of the βCD:Ibuprofen inclusion complex (the two orientations can be seen in Fig. 1 ) and free βCD in water at T = 298 . 15 K and P = 1 bar. The q4md-CD [51] force field was used for modelling the β-CDs, the TIP4P2005 [58] for water, and the General AMBER force field [63] was used for ibuprofen. The self-diffusivities are corrected for finite size effects using Eq. (2) . σ is the error with a 95% confidence interval. cupy a larger area around the CDs (corresponding to higher peaks in the RDF) than the Na + using the Madrid-2019 force field (see Fig. 7 b, d and f). From the contour plots it can also be seen that the preferred location of the Na + ions is near the hydroxyl groups of the larger rim of the CDs when the JC force field is used. When the Madrid-2019 model (see Fig. 7 b, d and f), Na + ions are located further from the center of mass of the CDs but still lie in the vicinity of the hydroxyl groups of the large rims. Cl − ions, regardless of the ion force field used, are always located outside of the CD, close to the center of the large rim.

Diffusion of the β-CD: Ibuprofen inclusion complex
Ibuprofen is an anti-inflammatory drug which is widely used worldwide for treating fever, pain and inflammation. [87] Due to its common use, ibuprofen is also a frequently detected contaminant in surface and ground waters [88][89][90][91] . Since β-CD can be potentially used both as a drug delivery agent [2,3,6] and in waste water treatment applications [92][93][94] , the mass transport of β-CD:Ibuprofen inclusion complex is an important aspect in both applications. In our MD simulations, the two possible orientations of ibuprofen inside the β-CD cavity are considered. In Orientation 1, the carboxyl group of ibuprofen is located near the small rim of the β-CD. In Orientation 2, the ibuprofen molecule is rotated by 180 degrees to place the carboxyl group near the large rim of the β-CD. A schematic representation of the orientations is shown in Fig. 1 c and d. The computed self-diffusion coefficients of the β-CD:Ibuprofen inclusion complex for the two orientations, and the free β-CD at 298.15 K are shown in Table 1 . The values of the selfdiffusivities of the two inclusion complexes are very close to each other (i.e. within 5%), and ca. 5-10% lower than the one of free β-CD.
Although no experimental measurements or molecular simulation computations of the self-diffusivity of the β-CD:Ibuprofen inlcusion complex in aqeous solution are reported, Santos et al. [30] presented Taylor dispersion experiments for the free β-CD and the β-CD:Caffeine inclusion complex. These experiments revealed that the diffusivity of the free β-CD is only slightly higher than the one of β-CD:caffeine inclusion complex (3.17 ×10 −10 m 2 s −1 and 3.05 ×10 −10 m 2 s −1 , respectively). Santos et al. [30] concluded that the inclusion complex has practically the same diffusivity with the free β-CD. The authors justified this finding by the fact that caffeine molecule is fully included in the cavity of the β-CD. Although these experiments were carried out using a different guest molecule compared to our study, a comparison between the two can be justified because the size of caffeine is comparable to the size of ibuprofen. As shown in Table 1 , the molecular simulation results are in full agreement with the experimental findings, showing no significant difference in the diffusivity of the free β-CD with the inclusion complex in water. In Figure S7 in the Supplementary material, the RDFs of water in respect to the free β-CD and the inclusion complex in both orientations are shown. In the case of free β-CD, the RDF peak at a distance of approx. 2 Å indicates that water molecules enter the cavity, however, in the case of the inclusion complex, no peak in the RDF is observed. This means that ibuprofen fully occupies the cavity of the β-CD preventing water molecules from entering.

Conclusions
In this work, the self-diffusion coefficients of α-, β-, γ -CDs and the β-CD:Ibuprofen inclusion complex in aqueous solutions are computed by means of MD simulation and compared with the respective experimental measurements. The effect of the system size on the self-diffusion of CDs is investigated. It is shown that the YH correction can be 76% of the final self-diffusion coefficient value. This suggests that the use of system size corrections for mixtures with compounds which are considerably different in size is absolutely necessary even if relatively large system sizes are used in the MD simulations. The self-diffusion coefficients of the native CDs in water are computed using four different water force fields at temperatures ranging from 298.15 to 312.15 K, at a pressure of 1 bar. The q4md-CD force field by Cézard et al. [51] was used to model CDs. We found that by using the TIP4P/2005 force field, the experimentally measured self-diffusion coefficients can be computed very accurately. The physical properties of CDs which can influence their diffusivity, such as the number of hydrogen bonds formed within the CDs and between the CDs and water, the average number of water molecules in the cavity of the CDs, and the size of the lower and upper rims of the CDs are also investigated to obtain a better understanding of effect of the water force field on these properties. Our results show that the effect of the water force field on these properties is almost negligible. This suggests that the large differences in the predicted self-diffusion coefficients are not caused by structural differences (such as the size and shape) of the CDs when different water force fields are used but mainly depend on the performance of these force fields in predicting the density and transport properties of the pure solvent. To reveal the effect of the salt concentration of the solvent on the diffusion of CDs, the self-diffusion coefficient of the three native CDs with two different ion force fields are computed as a function of the NaCl concentration at 298.15 K and 1 bar. It was found that the diffusivity of all CDs decreases with the increasing NaCl concentration, regardless of the ion force field. It is also shown that the Na + ion modelled by the JC force field has a higher affinity toward the CDs than the Madrid-2019 force field. This difference is likely the reason of the slower diffusivity of the αand β-CDs when the JC force field is used. Finally, as a model system for drug delivery and also for waster treatment application of β-CD, the diffusion of the β-CD:Ibuprofen inclusion complex in water is studied. It is shown that self-diffusion coefficient of the inclusion complex is approximately 5-10% lower than the free β-CD's. These almost equal diffusion coefficients are most likely caused by the almost full inclusion of ibuprofen in the cavity of the β-CD. Our findings indicate that MD simulations can provide reasonable predictions for supporting experiments, and the necessary molecular level understanding which is useful for industrially relevant applications of CDs.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.