**4. Visualizing NCI using a charge density topology**

The topological features of the electronic charge density have played an important part in several schemes aiming to provide a route to understanding molecular structure from first principles, that is, without resorting to any form of empirical or approximate models (other than the physically justified approximations underlying the electronic structure methods employed). Methods such as the theory of Atoms in Molecules (AIM) [92] and the Electron Localisation Function (ELF) [93–95] can provide great insight into features of the density. Both AIM and the ELF are useful for studying strong interactions such as covalent or ionic bonding but are not as useful in the analysis of very weak interactions of the type that are important for subjects as diverse as protein structure, drug design, catalysis, materials selfassembly, and many more. To address this, Johnson et al. [91, 96] introduced the NCI analysis method. This approach is based on the electronic charge density, *ρ*, and its derivatives. The first derivative *ρ* enters into the expression for the reduced density gradient, *s*,

$$S = \frac{1}{(2(3\ \pi^2))^{4/3}} \frac{|\nabla \rho|}{\rho^{4/3}} \tag{4}$$

readily obtained with first principles electronic structure methods, such calculations remain too computationally expensive for the large systems of interest in biological or materials science applications. In such cases it is possible to approximate the charge density with a sum of atomic densities providing a pro-molecular density, *ρpro*. This approximate density does not include relaxation of the atomic densities as would happen in self-consistent calculations for bonded systems. Fortunately, the largest deviation of *ρ*pro from the fully relaxed density occurs in covalent bonding regions and has a minimal effect on the (low-density) NCI regions. Substituting *ρpro* into the equations above has been found to have have very little impact on the resulting NCI analysis. For purposes of interfacing with the results of molecular dynamics simulations, pro-molecular densities have the additional advantages that they are readily computed for both finite and extended/periodic systems and the fact that this approach requires orders of magnitude less computational time than the method based on first prin-

Modeling Soft Supramolecular Nanostructures by Molecular Simulations

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

17

A recent development of the NCI method that makes use of the pro-molecular route to building the charge density of the system under study is the Independent Gradient Method (IGM) of Lefebvre et al. [90] As with the NCI method, IGM employs in addition to *ρ* quantities related to the first and second derivatives of the density. However, in IGM the reduced density gradient, *s*, is replaced with the descriptor *δginter*, defined as the difference between the

*ginter* = |∇*ρIGM*, *inter*|−|∇*ρ*| (6)

*δginter > 0* indicates the presence of NCI and the magnitude of the descriptor at a point in space gives an indication of the strength of the interaction. *ρIGM, inter* is obtained from sums over all the atomic densities in the different fragments (here for two fragments A and B in the

> <sup>=</sup> <sup>|</sup><sup>∑</sup> *i*=1 *NA* \_\_\_*<sup>i</sup> x*|<sup>+</sup> |∑ *i*=1 *NB* \_\_\_*<sup>i</sup>*

The IGM approach has the attractive feature that *δginter* allows for a more facile comparison of the strength of the weak interactions than the quantity *s* in the original NCI [90]. This is due to the fact that the presence of significant non-covalent interactions is indicated when *s* → *0* making interpretation difficult whereas the magnitude of the IGM quantity *δginter* increases in

An illustrative example of the use of the NCI and IGM methods is given in **Figures 3** and **4**. The host-guest complex shown is between curcurbit[7]uril and the dodecyl serine-based monomeric cationic surfactant 12SerTFAc [97] (to make the example more accessible to the reader only a single configuration is shown). Calculation of the NCI and IGM descriptors was

As can be seen from the molecular structure images, both the *s* and *δginter* surfaces provide qualitatively similar information and display the general distribution of the most significant

*x*<sup>|</sup> (7)

first derivatives of the charge densities for the total system and the fragments

\_\_\_ *x*) *IGM*,*inter*

performed using IGMPlot version 1.0 (http://kisthelp.univ-reims.fr/IGMPLOT/).

ciples electronic structure calculations.

x-direction)

(

proportion to the strength of the interaction.

In regions of low *ρ* such as the density tails far from a molecule the gradient remains large and so *s* displays high values. However, in regions of low *ρ* between atoms where weak NCI occur the gradient (and therefore *s*) drop to zero.

In order to differentiate between stabilizing attractive interactions and those that are unfavorable and destabilize the system it is necessary to analyze the second derivative (Laplacian) of the charge density *<sup>2</sup> ρ*. Decomposition of *<sup>2</sup> ρ* into the three eigenvalues representing the axes of maximal variation

$$
\nabla^2 \rho = \lambda\_1 + \lambda\_2 + \lambda\_{\mathcal{Y}} \text{ ( $\lambda\_1 \le \lambda\_2 \le \lambda\_3$ )}\tag{5}
$$

provides information on the nature of the interactions at a given point in space. *λ<sup>2</sup>* displays negative values in stabilizing/bonding regions (charge is flowing into this region indicating a local build-up of *ρ*) while in destabilizing/repulsive regions it is positive (charge flowing out, indicating a local depletion of *ρ*). Plotting *s* against *sign(λ<sup>2</sup> )ρ* will therefore permit the identification of NCI in regions where s and *sign(λ<sup>2</sup> )ρ* → *0*. Although these quantities can be readily obtained with first principles electronic structure methods, such calculations remain too computationally expensive for the large systems of interest in biological or materials science applications. In such cases it is possible to approximate the charge density with a sum of atomic densities providing a pro-molecular density, *ρpro*. This approximate density does not include relaxation of the atomic densities as would happen in self-consistent calculations for bonded systems. Fortunately, the largest deviation of *ρ*pro from the fully relaxed density occurs in covalent bonding regions and has a minimal effect on the (low-density) NCI regions.

of regions of low charge density corresponding to stabilizing/destabilizing NCI, based on the analysis of the electronic charge density of the interacting molecules and the respective gradients. Although the original NCI method of Johnson and co-workers [91] provide a very similar qualitative analysis of NCI, the reduced density gradient, s, used to identify the interaction types is a dimensionless quantity and therefore difficults the evaluation of the respective strengths. The new IGM allows for a quantitative comparison of the strength of NCI

The topological features of the electronic charge density have played an important part in several schemes aiming to provide a route to understanding molecular structure from first principles, that is, without resorting to any form of empirical or approximate models (other than the physically justified approximations underlying the electronic structure methods employed). Methods such as the theory of Atoms in Molecules (AIM) [92] and the Electron Localisation Function (ELF) [93–95] can provide great insight into features of the density. Both AIM and the ELF are useful for studying strong interactions such as covalent or ionic bonding but are not as useful in the analysis of very weak interactions of the type that are important for subjects as diverse as protein structure, drug design, catalysis, materials selfassembly, and many more. To address this, Johnson et al. [91, 96] introduced the NCI analysis method. This approach is based on the electronic charge density, *ρ*, and its derivatives. The

first derivative *ρ* enters into the expression for the reduced density gradient, *s*,

*ρ*. Decomposition of *<sup>2</sup>*

out, indicating a local depletion of *ρ*). Plotting *s* against *sign(λ<sup>2</sup>*

provides information on the nature of the interactions at a given point in space. *λ<sup>2</sup>*

(2(3 *π*<sup>2</sup>))1/3

In regions of low *ρ* such as the density tails far from a molecule the gradient remains large and so *s* displays high values. However, in regions of low *ρ* between atoms where weak NCI occur

In order to differentiate between stabilizing attractive interactions and those that are unfavorable and destabilize the system it is necessary to analyze the second derivative (Laplacian)

negative values in stabilizing/bonding regions (charge is flowing into this region indicating a local build-up of *ρ*) while in destabilizing/repulsive regions it is positive (charge flowing

\_\_\_\_ |∇*ρ*|

, which corresponds directly to the charge

*<sup>ρ</sup>*4/3 (4)

*ρ* into the three eigenvalues representing the

, (*λ*<sup>1</sup> ≤ *λ*<sup>2</sup> ≤ *λ*3) (5)

*)ρ* → *0*. Although these quantities can be

displays

*)ρ* will therefore permit the

through the calculation of the IGM descriptor, δ<sup>g</sup>

*s* = \_\_\_\_\_\_\_\_ <sup>1</sup>

the gradient (and therefore *s*) drop to zero.

∇<sup>2</sup> *ρ* = *λ*<sup>1</sup> + *λ*<sup>2</sup> + *λ*<sup>3</sup>

identification of NCI in regions where s and *sign(λ<sup>2</sup>*

of the charge density *<sup>2</sup>*

axes of maximal variation

**4. Visualizing NCI using a charge density topology**

density gradient(s) in real-space.

16 Molecular Dynamics

Substituting *ρpro* into the equations above has been found to have have very little impact on the resulting NCI analysis. For purposes of interfacing with the results of molecular dynamics simulations, pro-molecular densities have the additional advantages that they are readily computed for both finite and extended/periodic systems and the fact that this approach requires orders of magnitude less computational time than the method based on first principles electronic structure calculations.

A recent development of the NCI method that makes use of the pro-molecular route to building the charge density of the system under study is the Independent Gradient Method (IGM) of Lefebvre et al. [90] As with the NCI method, IGM employs in addition to *ρ* quantities related to the first and second derivatives of the density. However, in IGM the reduced density gradient, *s*, is replaced with the descriptor *δginter*, defined as the difference between the first derivatives of the charge densities for the total system and the fragments

$$\delta \mathbf{g}^{\text{inter}} = \left| \nabla \rho^{\text{IGM}, \text{ inter}} \right| - \left| \nabla \rho \right| \tag{6}$$

*δginter > 0* indicates the presence of NCI and the magnitude of the descriptor at a point in space gives an indication of the strength of the interaction. *ρIGM, inter* is obtained from sums over all the atomic densities in the different fragments (here for two fragments A and B in the x-direction)

$$\left(\frac{\delta\rho}{\delta x}\right)^{\text{IGM},\text{totr}} = \left| \sum\_{i=1}^{N\_e} \frac{\delta\rho\_i}{\delta x} \right| + \left| \sum\_{i=1}^{N\_e} \frac{\delta\rho\_i}{\delta x} \right| \tag{7}$$

The IGM approach has the attractive feature that *δginter* allows for a more facile comparison of the strength of the weak interactions than the quantity *s* in the original NCI [90]. This is due to the fact that the presence of significant non-covalent interactions is indicated when *s* → *0* making interpretation difficult whereas the magnitude of the IGM quantity *δginter* increases in proportion to the strength of the interaction.

An illustrative example of the use of the NCI and IGM methods is given in **Figures 3** and **4**. The host-guest complex shown is between curcurbit[7]uril and the dodecyl serine-based monomeric cationic surfactant 12SerTFAc [97] (to make the example more accessible to the reader only a single configuration is shown). Calculation of the NCI and IGM descriptors was performed using IGMPlot version 1.0 (http://kisthelp.univ-reims.fr/IGMPLOT/).

As can be seen from the molecular structure images, both the *s* and *δginter* surfaces provide qualitatively similar information and display the general distribution of the most significant

of) systems such as those arising from MD simulations it is likely that the future will see the

Modeling Soft Supramolecular Nanostructures by Molecular Simulations

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

19

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

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,

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

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

naphthalene promoted an increase in the binding constant (up to 100-fold), irrespective of the

kJ.mol−1), pointing out the important effects of the substituents. Substitution of

the main driving forces responsible for the formation of this type of complexes [99].

2− (**Figure 5**), and ferrocenyl alkanethiols with both free and gold-surface

newer IGM method being used much more extensively than the original NCI one.

**5. Applications**

and proteins.

such as ClO<sup>4</sup>

128 to 2.1×10<sup>4</sup>

− and SO<sup>4</sup>

**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> )ρ* (−0.03 ≤ *sign(λ<sup>2</sup> )ρ* ≤ 0.03 a.U.). Blue = stabilizing, red = destabilizing and green = weak interactions. (right) scatter plot of *s* and *sign(λ<sup>2</sup> )ρ* values.

**Figure 4.** Host-guest complexation between curcurbit[7]uril and 12serTFAc sufactant. (left) IGM *δginter* = 0.01 a.U.Isosurface colored by value of *sign(λ<sup>2</sup> )ρ* (−0.03 ≤ *sign(λ<sup>2</sup> )ρ* ≤ 0.03 a.U.). Blue = stabilizing, red = destabilizing and green = weak interactions. (right) scatter plot of *δginter* and *sign(λ<sup>2</sup> )ρ* values.

host-guest interactions. Similarly, the color-coding provided by the value of *sign(λ<sup>2</sup> )ρ* can be 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 spikes in the non-covalent interaction region −0.025 ≤ *sign(λ<sup>2</sup> )ρ* ≤ 0.025 in a very similar manner 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 of) systems such as those arising from MD simulations it is likely that the future will see the newer IGM method being used much more extensively than the original NCI one.
