**3. Materials and Methods**

The interaction of PX with the β-CD NS model was studied using molecular mechanics (MM) and molecular dynamics (MD) methods at the constant temperature of 300 K. The MM and MD simulations are performed using the consistent valence force field (CVFF) [33], and the Materials Studio and Insight II/Discover packages [34] adopting the same simulation protocol as proposed in yjr previous work [20–23]: after the initial energy minimization, MD runs were carried out until the equilibrium state was achieved followed by final geometry optimizations of the final configuration and of some conformations assumed during the MD run. The PX structure (see Scheme 1) was generated using the Insight II/Discover Module Builder and the final most stable geometry (Figure 3) was obtained after an MD run lasting for 1 ns, with the final optimization of different conformations assumed by the system during the MD run. At first, the possible β-CD/PX inclusion compounds in a 1:1 and in a 2:1 stoichiometry were studied. Then, the interaction between the β-CD NS model and PX was considered. The structure of the β-CD NS model reported by Raffaini [30] was generated by linking 8 β-CD (Model 2 in [30]) with pyromellitic dianhydride (PMA): each β-CD carries two PMA linking agents bound to a primary hydroxyl at diametrically opposite sides of the macrocycle, linking also two topologically distant β-CDs. The drug molecules were initially placed in a random distribution cell far from to the NS surface in order to study the interactions without assuming any a priori inclusion complexes. The interaction was studied considering 4, 8, 16, and 40 PX molecules randomly distributed around the NS model in a 220 Å cubic cell with a varying CDs/PX ratio. All simulations were performed in implicit water using a distance-dependent dielectric constant with periodic boundary conditions. All energy minimizations were carried out up to an energy gradient < 4 × 10−<sup>3</sup> kJ mol−<sup>1</sup> Å−<sup>1</sup> . The MD simulations were performed in an NVT ensemble (canonical ensemble, see Abbreviations) at a constant temperature (300 K) controlled using a Berendsen thermostat. The integration of dynamical equations was carried out with the Verlet algorithm using a time step of 1 fs and the instantaneous coordinates were periodically saved for further analysis. The MD runs in implicit water lasted for 10 ns in Section 2.1 and 100 ns for the NS systems in Section 2.2. During the MD runs, the time evolution of the potential energy and of the van der Waals components was calculated in order to monitor the significant conformational changes.
