**5. Applications**

host-guest interactions. Similarly, the color-coding provided by the value of *sign(λ<sup>2</sup>*

*)ρ* values.

spikes in the non-covalent interaction region −0.025 ≤ *sign(λ<sup>2</sup>*

*)ρ* (−0.03 ≤ *sign(λ<sup>2</sup>*

interactions. (right) scatter plot of *δginter* and *sign(λ<sup>2</sup>*

colored by value of *sign(λ<sup>2</sup>*

18 Molecular Dynamics

seen to predict that the interactions are largely of the weak van der Waals type with the exception of a single blue region denoting the hydrogen bond between the serine hydroxyl group and the host. A significant difference between the *s* and *δginter* surfaces is that the latter also contains information on the strength of the interactions through the volumes of the isosurface regions at a given point. The IGM approach therefore provides an important extra source of information on the nature of the interactions. The intrinsic difference between the two methods is even more clearly seen in the adjacent scatter plots. In the IGM plot (**Figure 4**, right) peaks of different magnitude corresponding to different strengths of interaction can be seen which give rise to the differing isosurface volumes. The NCI plot (**Figure 3**, right), however, displays all

**Figure 4.** Host-guest complexation between curcurbit[7]uril and 12serTFAc sufactant. (left) IGM *δginter* = 0.01 a.U.Isosurface

**Figure 3.** Host-guest complexation between curcurbit[7]uril and the dodecyl serine-based monomeric cationic surfactant

12SerTFAc. (left) reduced density gradient, *s* = 0.4, isosurface colored by value of *sign(λ<sup>2</sup>*

Blue = stabilizing, red = destabilizing and green = weak interactions. (right) scatter plot of *s* and *sign(λ<sup>2</sup>*

*)ρ* ≤ 0.03 a.U.). Blue = stabilizing, red = destabilizing and green = weak

and makes interpretation much less easy. For this reason, despite very similar computational requirements for both methods making them equally suitable for analysis of large (ensembles

*)ρ* can be

*)ρ* ≤ 0.03 a.U.).

*)ρ* values.

*)ρ* ≤ 0.025 in a very similar manner

*)ρ* (−0.03 ≤ *sign(λ<sup>2</sup>*

This section presents some examples in which MD or MC simulations are used for establishing the contributions of NCI (e.g. electrostatic, hydrophobic and hydrogen-bond interactions) to the free energy, in systems in which these are the main source of recognition, assembly, and stability. These include, for instance, host-guest complexes, polyelectrolytes and proteins.

MD and free-energy calculations have provided insight on relevant aspects including the effect of substituents, the role of solvation, and rationales for the conformation and thermodynamic characterization of the systems under investigation. Among several studies on the topic, available in the literature, only a few are briefly reviewed. For instance, MD simulations in water, with PMF estimation based on both TI and constraint force methods, have been used for inspecting the thermodynamic properties of binding of some inorganic ions, such as ClO<sup>4</sup> − and SO<sup>4</sup> 2− (**Figure 5**), and ferrocenyl alkanethiols with both free and gold-surface confined β-cyclodextrins [14, 98]. In general, the association process between the host and guest molecules is analyzed along the reaction coordinate defined by the distance between the centers of mass of both host and guest, along a reference coordinate (e.g. z-axis). In these particular systems, electrostatic interactions and hydrogen bonding have played a major role on the binding process. For instance, simulations have elucidated the association mode of sulfate anion with free and grafted β-cyclodextrin (see **Figure 5**). For the latter system, a small minimum in the PMF profile (**Figure 5**, top), positioned at a larger distance relative to the cavity of surface-grafted cyclodextrin, suggested the formation of a noninclusion complex. The ion binds to the host by forming hydrogen bonds with the free-cyclodextrin portal. Also relevant is the individual energy contributions to the cyclodextrin-ion interaction, which is governed by electrostatic interactions. However, this favorable electrostaction contribution is not sufficient to compensate desolvation of the anion, considering the respective hydration energy [14].

The release and transport of drugs mediated by cyclodextrin-based carriers, have also been investigated [99] systematically using MD and free-energy calculations, showing the preferred inclusion modes of such drugs for cyclodextrins. One example (see **Figure 6**), refers to the binding of amphotericin B (AmB), which possesses two sites, within the respective prolonged macrolide, with higher binding affinity to γ-cyclodextrin. The decomposition of the PMFs into free-energy contributions have suggested that van der Waals and electrostatic interactions are the main driving forces responsible for the formation of this type of complexes [99].

Recently, some of the authors [5] have demonstrated the relevance of non-included moieties in the stability constants of several cyclodextrin-based complexes. The binding constants for naphtalene derivatives forming complexes with β-cyclodextrin were calculated (ranging from 128 to 2.1×10<sup>4</sup> kJ.mol−1), pointing out the important effects of the substituents. Substitution of naphthalene promoted an increase in the binding constant (up to 100-fold), irrespective of the

**Figure 5.** (Top) schematic representation of the association pathway (left) between SO<sup>4</sup> 2− and a β-cyclodextrin grafted to a gold-surface, used for the construction of the PMF profile (right). The free-energy profile characterizing the association process, in water, of the ion into the free β-cyclodextrin is also shown. (Bottom) illustrative spanshot showing the presence of hydrogen bonding (right) as the ion approximates the cavity of the host, and (left) graphical distribution of distances and angles correspondonding to the hydrogen bonds. Reproduced from Ref. [14] with permission from the Royal Society of Chemistry.

of both guest and host cavity. One of the most interesting observation was that chloride was hermetically sealed inside the host cavity, as a result of a concerted action of both conforma-

**Figure 6.** (A) Schematic representation of the two inclusion pathways along which AmB approximates to β- or γciclodextrins. (B) PMF profiles for the inclusion of AmB with β- and γ-cyclodextrins considering the two orientations. (C) Equilibrium configuration of the complex between AmB and β-cyclodextrin (orientation I) followed by a representation of the hydrogen bonds formed between host and guest molecules and the respective averaged number (O⋯H⋯O angle >135° and O⋯O distance <3.5 Å), as a function of the simulation time. (D) Decomposition of the PMFs into van der Waals and electrostatic host-guest and host-solvent contributions, for the complexes between AmB and β- and γ-cyclodextrin in orientation I, and AmB and γ-cyclodextrin in orientation II. Reproduced from Ref. [99] with permission from the

Modeling Soft Supramolecular Nanostructures by Molecular Simulations

http://dx.doi.org/10.5772/intechopen.74939

21

Solvent contributions to the thermodynamics of binding have been scrutinized [13], emphasizing the importance of using explicit models for the accurate calculation of binding thermodynamics. A detailed analysis on the energetics of the host-guest complexation between β-cyclodextrin and several model drugs (e.g. puerarin, daidzin, daidzein, and nabumetone) have demonstrated that the flexibility of the binding partners and solvation-related enthalpy and entropy changes must be included explicitly for the precise estimation of thermodynamic parameters involved in molecular association. In this study, it was also demonstrated that implicit models are not suitable to provide detailed information on how free-energy is decom-

The dimerization of cyclodextrins has also been recognized as a relevant step in the construction of nanostructured materials. Only a few studies (see e.g. refs [15, 48].) involving umbrella

tion change and desolvation (see **Figure 7**) [26].

American Chemical Society.

posed into enthalpy and entropy [13].

nature of the substituent, the latter comprising small hydrophobic and hydrophilic, including charged, groups. Enthalpic and entropic contributions were separated in order to estimate their weight in the free-energy. Also, the preferred orientation of the guest molecules within the cyclodextrin was established.

In a completely different system [26], the same strategy was used for exploring the ion caging ability of bambusurils in aqueous media (**Figure 7**). Specifically, new insights on the conformation, hydration and energy changes involved in the complexation process between the water-soluble benzoyl-substituted bambus[6]uril and chloride ion were provided. The structural features of three single bambusuril derivatives, and the relative energy contributions to the formation of the benzoyl-substituted bambus[6]uril-chloride complex were computed. The estimated standard free-energy of binding and the respective binding constant were found to be −11.7 kJ.mol− 1 and 112, respectively. Binding occurred with complete desolvation

**Figure 6.** (A) Schematic representation of the two inclusion pathways along which AmB approximates to β- or γciclodextrins. (B) PMF profiles for the inclusion of AmB with β- and γ-cyclodextrins considering the two orientations. (C) Equilibrium configuration of the complex between AmB and β-cyclodextrin (orientation I) followed by a representation of the hydrogen bonds formed between host and guest molecules and the respective averaged number (O⋯H⋯O angle >135° and O⋯O distance <3.5 Å), as a function of the simulation time. (D) Decomposition of the PMFs into van der Waals and electrostatic host-guest and host-solvent contributions, for the complexes between AmB and β- and γ-cyclodextrin in orientation I, and AmB and γ-cyclodextrin in orientation II. Reproduced from Ref. [99] with permission from the American Chemical Society.

of both guest and host cavity. One of the most interesting observation was that chloride was hermetically sealed inside the host cavity, as a result of a concerted action of both conformation change and desolvation (see **Figure 7**) [26].

nature of the substituent, the latter comprising small hydrophobic and hydrophilic, including charged, groups. Enthalpic and entropic contributions were separated in order to estimate their weight in the free-energy. Also, the preferred orientation of the guest molecules within

a gold-surface, used for the construction of the PMF profile (right). The free-energy profile characterizing the association process, in water, of the ion into the free β-cyclodextrin is also shown. (Bottom) illustrative spanshot showing the presence of hydrogen bonding (right) as the ion approximates the cavity of the host, and (left) graphical distribution of distances and angles correspondonding to the hydrogen bonds. Reproduced from Ref. [14] with permission from the

2− and a β-cyclodextrin grafted to

**Figure 5.** (Top) schematic representation of the association pathway (left) between SO<sup>4</sup>

In a completely different system [26], the same strategy was used for exploring the ion caging ability of bambusurils in aqueous media (**Figure 7**). Specifically, new insights on the conformation, hydration and energy changes involved in the complexation process between the water-soluble benzoyl-substituted bambus[6]uril and chloride ion were provided. The structural features of three single bambusuril derivatives, and the relative energy contributions to the formation of the benzoyl-substituted bambus[6]uril-chloride complex were computed. The estimated standard free-energy of binding and the respective binding constant were found to be −11.7 kJ.mol− 1 and 112, respectively. Binding occurred with complete desolvation

the cyclodextrin was established.

Royal Society of Chemistry.

20 Molecular Dynamics

Solvent contributions to the thermodynamics of binding have been scrutinized [13], emphasizing the importance of using explicit models for the accurate calculation of binding thermodynamics. A detailed analysis on the energetics of the host-guest complexation between β-cyclodextrin and several model drugs (e.g. puerarin, daidzin, daidzein, and nabumetone) have demonstrated that the flexibility of the binding partners and solvation-related enthalpy and entropy changes must be included explicitly for the precise estimation of thermodynamic parameters involved in molecular association. In this study, it was also demonstrated that implicit models are not suitable to provide detailed information on how free-energy is decomposed into enthalpy and entropy [13].

The dimerization of cyclodextrins has also been recognized as a relevant step in the construction of nanostructured materials. Only a few studies (see e.g. refs [15, 48].) involving umbrella

**Figure 7.** Typical conformation of the inclusion complex between one chloride ion (purple sphere) and the dodeca (4-carboxybenzyl)bambus[6]uril, in water, sampled during the MD run, at 298 K. (Right) Scheme of the dissociation pathway, along the z-component (ξ<sup>z</sup> ), between host and guest. The steering force is applied on the chloride ion for generating the initial configurations for the umbrella sampling procedure. Reproduced from Ref [26] with permission from Elsevier.

sampling simulations and PMF calculations have been carried out for investigating the binding affinity of dimers in the presence/absence of a guest molecule, and the preferred orientation of interglucopyranose hydrogen bonds, at the cyclodextrin portals. These include β-cyclodextrin monomers and dimers in aqueous and nonaqueous media [48]. Polar solvents with stronger hydroden bond accepting abilities can easily disrupt intermolecular hydrogen bonds, decreasing the stability of the dimers. In the same environment, higher binding affinities are achieved if guest molecules are included in the channel-type cavity of such dimers. These results allowed concluding that the formation of dimers is solvent-dependent and guest-modulated [48]. In another related study [15], it was shown that the cooperative binding of cyclodextrin cavities to guest molecules can facilitates the dimerization process, favoring the overall stability and assembly of nanostructures. Different dimerization modes yielded different binding strengths. This has been demonstrated in systems consisting of isoflavone drug analogs used as drug templates, and cyclodextrins (**Figure 8**). A detailed quantification of host, guest and solvent contributions to the thermodynamics of binding has revealed that head-to-head dimerization promotes the most stable complexes, which can be used as building blocks for template-stabilized nanostructures. Desolvation of cyclodextrin dimers and entropy changes upon complexation also affect significantly the cooperative binding [15].

polyelectrolytes, based on MC and MD, have been precluded by the respective chain lengths and the need of large concentrations of compacting agents. The conformational and energetic changes in DNA delivery systems have been studied computationally by some of the authors [101, 102]. For instance, a coarse-grained model was proposed [101] for exploring structural and thermodynamics aspects of confining polyions in spherical cavities, with different linear charge densities and with counterions of different valences. The free-energy of the confined polyion and counterions were estimated as a function of the sphere radius, from a dilute solution corresponding to the largest sphere. A positive free-energy difference was found for all systems and was associated to the compression. This means that the system containing the polyion with the largest linear charge density and monovalent counterions displays the largest resistance to being compressed. The penalty is thus largest for the polyion with the largest linear charge density and monovalent counterions. The replacement of the monovalent counterions with trivalent ones induced a compactation of the polyion, as a consequence of the stronger electrostatic polyion-counterion attraction. This leads to a nearly full compensation

permission from the American Chemical Society.

**Figure 8.** (Top) schematic representation of the association coordinate for the 2:1 complex between a β- cyclodextrin dimer and daidzein. (Bottom) PMF profiles and representative configurations of the 2:1 complexes between the specific BHHP arrangement of β- cyclodextrin dimer and three different model drugs. Reproduced from Ref. [15] with

Modeling Soft Supramolecular Nanostructures by Molecular Simulations

http://dx.doi.org/10.5772/intechopen.74939

23

of the ideal and electrostatic contributions to the free-energy of their confinement.

Other interesting example of application relates to protein-ligand interactions [103–107], which are typically pursued by significant conformational and energetic changes. These

Other strategies have been adopted resorting to both MD and MC simulations. As an example of application [100], free-energy calculations and lattice chain MC simulations have been used for understanding the formation of polyrotaxanes (**Figure 9**), including the identification of the most favorable conformations of cyclodextrin molecules in a polyrotaxane and the quantification of dimerization free energies, related to different spatial arrangements of two consecutive cyclodextrin molecules. It has been suggested, that the binary association is controlled by the formation of hydrogen bonds between two adjacent molecules, promoting an overall structural stabilization [100].

The thermodynamic cost of confining polyelectrolytes spherical cells of different radii can also be quantified by free-energy calculations, resorting to TI [101]. Simulation studies of

Modeling Soft Supramolecular Nanostructures by Molecular Simulations http://dx.doi.org/10.5772/intechopen.74939 23

sampling simulations and PMF calculations have been carried out for investigating the binding affinity of dimers in the presence/absence of a guest molecule, and the preferred orientation of interglucopyranose hydrogen bonds, at the cyclodextrin portals. These include β-cyclodextrin monomers and dimers in aqueous and nonaqueous media [48]. Polar solvents with stronger hydroden bond accepting abilities can easily disrupt intermolecular hydrogen bonds, decreasing the stability of the dimers. In the same environment, higher binding affinities are achieved if guest molecules are included in the channel-type cavity of such dimers. These results allowed concluding that the formation of dimers is solvent-dependent and guest-modulated [48]. In another related study [15], it was shown that the cooperative binding of cyclodextrin cavities to guest molecules can facilitates the dimerization process, favoring the overall stability and assembly of nanostructures. Different dimerization modes yielded different binding strengths. This has been demonstrated in systems consisting of isoflavone drug analogs used as drug templates, and cyclodextrins (**Figure 8**). A detailed quantification of host, guest and solvent contributions to the thermodynamics of binding has revealed that head-to-head dimerization promotes the most stable complexes, which can be used as building blocks for template-stabilized nanostructures. Desolvation of cyclodextrin dimers and entropy changes

**Figure 7.** Typical conformation of the inclusion complex between one chloride ion (purple sphere) and the dodeca (4-carboxybenzyl)bambus[6]uril, in water, sampled during the MD run, at 298 K. (Right) Scheme of the dissociation

generating the initial configurations for the umbrella sampling procedure. Reproduced from Ref [26] with permission

), between host and guest. The steering force is applied on the chloride ion for

Other strategies have been adopted resorting to both MD and MC simulations. As an example of application [100], free-energy calculations and lattice chain MC simulations have been used for understanding the formation of polyrotaxanes (**Figure 9**), including the identification of the most favorable conformations of cyclodextrin molecules in a polyrotaxane and the quantification of dimerization free energies, related to different spatial arrangements of two consecutive cyclodextrin molecules. It has been suggested, that the binary association is controlled by the formation of hydrogen bonds between two adjacent molecules, promoting

The thermodynamic cost of confining polyelectrolytes spherical cells of different radii can also be quantified by free-energy calculations, resorting to TI [101]. Simulation studies of

upon complexation also affect significantly the cooperative binding [15].

an overall structural stabilization [100].

pathway, along the z-component (ξ<sup>z</sup>

from Elsevier.

22 Molecular Dynamics

**Figure 8.** (Top) schematic representation of the association coordinate for the 2:1 complex between a β- cyclodextrin dimer and daidzein. (Bottom) PMF profiles and representative configurations of the 2:1 complexes between the specific BHHP arrangement of β- cyclodextrin dimer and three different model drugs. Reproduced from Ref. [15] with permission from the American Chemical Society.

polyelectrolytes, based on MC and MD, have been precluded by the respective chain lengths and the need of large concentrations of compacting agents. The conformational and energetic changes in DNA delivery systems have been studied computationally by some of the authors [101, 102]. For instance, a coarse-grained model was proposed [101] for exploring structural and thermodynamics aspects of confining polyions in spherical cavities, with different linear charge densities and with counterions of different valences. The free-energy of the confined polyion and counterions were estimated as a function of the sphere radius, from a dilute solution corresponding to the largest sphere. A positive free-energy difference was found for all systems and was associated to the compression. This means that the system containing the polyion with the largest linear charge density and monovalent counterions displays the largest resistance to being compressed. The penalty is thus largest for the polyion with the largest linear charge density and monovalent counterions. The replacement of the monovalent counterions with trivalent ones induced a compactation of the polyion, as a consequence of the stronger electrostatic polyion-counterion attraction. This leads to a nearly full compensation of the ideal and electrostatic contributions to the free-energy of their confinement.

Other interesting example of application relates to protein-ligand interactions [103–107], which are typically pursued by significant conformational and energetic changes. These

changes often occur on time scales that make direct atomic-level simulations useless. One possible alternative relies on the assumption that the change in the binding free-energy resulting, for e.g. from a mutation of a ligand bounded to a protein, complies a linear response scheme [104], with parameters estimated from training sets of protein-ligand complexes, and used for predicting binding affinities of new ligands. Another strategy is the computational design of the ligand in free and bound states. In this type of systems, the emergence of some reliable implicit solvation models and classical statistical-mechanical approaches have been part of the solution. These include, for instance, the use of the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) model [46]. A theoretical study, recently published [108], encompasses MD and free-energy calculations for evaluating the structural and thermodynamic signatures involved in the interaction of siRNA molecules, bearing chemical modifications at 3'-overhang, with the PAZ domain of human Argonaute (**Figure 10**). In these systems, a reduction of the complex binding affinity has been observed upon siRNA modification, being explained by the hampering of hydrogen bond formation in the active site of PAZ. Analyses of free-energy, achieved from simulations based on the Generalized Born and surface area continuum solvation (MM-GBSA) method, and of hydrogen bonds, have provided a complete identification of the most relevant residues for PAZ-siRNA interaction (see **Figure 10**). This data will contribute to the improved design of synthetic nucleotide analogues, circumventing

Modeling Soft Supramolecular Nanostructures by Molecular Simulations

http://dx.doi.org/10.5772/intechopen.74939

25

Molecular simulations, including free-energy calculations is still a fertile ground for research, from both theoretical and applied perspectives. Among the various challenges in the field of free-energy calculations, the integration of thermodynamics and kinetics through the concomitant determination of potentials of mean force and diffusivities in biased simulations, will benefit from substantial efforts. Although there is a general consensus that free-energy calculations are an important tool of computational chemistry and structural biology, thriving beyond academic walls, namely in industrial environments, some concerns still remain, dealing with the question to whether or not those methods provide convincing and reliable

The authors acknowledge the Fundação para a Ciência a Tecnologia (FCT), Portuguese Agency for Scientific Research, through the Projects n. 016648 POCI-01-0145-FEDER-016648 and COMPETE POCI-01-0145-FEDER-007440. The Coimbra Chemistry Centre is also supported by FCT through the Projects PEst-OE/QUI/UI0313/2014 and POCI-01-0145-FEDER-007630. Tânia F.G.G. Cova, Sandra C. C. Nunes and Andreia. F. Jorge acknowledge, respectively, the PhD and post-doctoral research Grants SFRH/BD/95459/2013, SFRH/BPD/71683/2010 and

some of the intrinsic siRNA drawbacks.

**6. Concluding remarks**

**Acknowledgements**

SFRH/BPD/104544/2014, assigned by FCT.

answers.

**Figure 9.** (A) Structural arrangement of a polyrotaxane formed formed by a poly(ethylene glycol chain included in multiple α-cyclodextrins, (B) definition of the reaction coordinate, ξ, in wich both the geometrically restrained (cyan) and free (red) cyclodextrins are presented. (C) Three possible conformations (HH, HT and TT) showing different spatial arrangements of two consecutive α-cyclodextrin molecules. (D-F) Free-energy profiles corresponding to the dimerization of α-cyclodextrins on the poly(ethylene glycol chain, for the HH, HT and TT conformations. (G-I) Individual cyclodextrin-cyclodextrin, cyclodextrin-water, and cyclodextrin-thread energy contributions for the total free-energy, in the three conformations. Reproduced from Ref. [100] with permission from the American Chemical Society.

**Figure 10.** Representative conformation of the complex between PAZ domain of the human argonaute 2 (hAgo2) and the chemically modified siRNAs at 3′ overhang, obtained from the last nanoseconds of the MD production runs. Panel A illustrates the relative positioning of domains and interdomain linker of the hAgo2 (using the crystal coordinates of PDB ID: 4F3T) and the positioning of the siRNA antisense strand. Panels **B** and **C** show the most representative structures of the complexes formed between PAZ domain and siRNAs, modified with phosphorothioate thymidine (PS) and L-threoninol-thymine (THR), respectively. PS and THR correspond to the second-last residues of the modified siRNAs while PS3' and THR3' correspond to the last one. Hydrogen-bonding interactions between the siRNA nucleotides and the PAZ amino acid residues are represented in black dash lines. The phosphorus, oxygen, nitrogen, hydrogen and sulfur atoms are colored in yellow, red, blue, white, and green respectively. Panels **D** and **E** correspond to the occupancies of the most prominent hydrogen bonds formed between the PAZ binding pocket and 2-nt modified nucleotides with PS and THR, respectively. The analysis of the instantaneous hydrogen bond formation was obtained using an in-house algorithm. Adapted from Ref. [108].

changes often occur on time scales that make direct atomic-level simulations useless. One possible alternative relies on the assumption that the change in the binding free-energy resulting, for e.g. from a mutation of a ligand bounded to a protein, complies a linear response scheme [104], with parameters estimated from training sets of protein-ligand complexes, and used for predicting binding affinities of new ligands. Another strategy is the computational design of the ligand in free and bound states. In this type of systems, the emergence of some reliable implicit solvation models and classical statistical-mechanical approaches have been part of the solution. These include, for instance, the use of the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) model [46]. A theoretical study, recently published [108], encompasses MD and free-energy calculations for evaluating the structural and thermodynamic signatures involved in the interaction of siRNA molecules, bearing chemical modifications at 3'-overhang, with the PAZ domain of human Argonaute (**Figure 10**). In these systems, a reduction of the complex binding affinity has been observed upon siRNA modification, being explained by the hampering of hydrogen bond formation in the active site of PAZ. Analyses of free-energy, achieved from simulations based on the Generalized Born and surface area continuum solvation (MM-GBSA) method, and of hydrogen bonds, have provided a complete identification of the most relevant residues for PAZ-siRNA interaction (see **Figure 10**). This data will contribute to the improved design of synthetic nucleotide analogues, circumventing some of the intrinsic siRNA drawbacks.
