Solvent-Dependent Conformational States of a [2]Rotaxane-Based

Nov 4, 2013 - We contribute to this subject by looking into the solvent influence on the conformational states of a synthetic molecular rotor, namely,...
1 downloads 11 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCC

Solvent-Dependent Conformational States of a [2]Rotaxane-Based Molecular Machine: A Molecular Dynamics Perspective Naresh K. Jena‡ and N. Arul Murugan*,† †

Division of Theoretical Chemistry and Biology School of Biotechnology, Royal Institute of Technology, SE-10691 Stockholm, Sweden ‡ Department of Physics, Stockholm University, Roslagstullsbacken 21, Albanova, SE-10691 Stockholm, Sweden ABSTRACT: Motion is an essential and fundamental feature of any living organism. The evolved organisms have developed sophisticated and perfect machineries and highly delicate mechanisms to carry out directional and coordinated movements which eventually lead to motion at the macroscopic length scale. By mimicking these natural machineries, attempts to design and synthesize similar molecular motors are made in relevance to their applications in drug delivery, data storage, and molecular sensing. It is highly desirable to establish the rules for controlling the conformational states of molecular motors by tuning some of the external variables which can be used for the design strategies. We contribute to this subject by looking into the solvent influence on the conformational states of a synthetic molecular rotor, namely, diketopyrrolopyrrole (DPP) based [2]rotaxane, using the force-field molecular dynamics approach. We study this system in three different solvents, and we report a strong solvent dependence in the population of three different translational isomers. In chloroform solvent we report the dominant population of the 2-P isomer which is in excellent agreement with experimental results based on H NMR spectra (Org. Lett. 2013, 15, 1274). However, there is a striking difference seen in the population of translational isomers in DMSO solvent, and we attribute these features to negligence of solvent hydrogen bonding induced upfield and downfield effects in the interpretation of experimental proton NMR spectra. In addition, we also report a solvent-polarity-induced fully unstretched to folded conformational transition in the [2]rotaxane system. On the basis of the molecular mechanics Poisson−Boltzmann (and generalized Born) surface area approach, we identify the driving force for the formation of the supramolecular guest−host [2]rotaxane system. Finally, we calculate the relative binding free energies for the macrocycle at different binding sites of the DPP skeleton using the molecular dynamics simulations performed for the macrocycle−rotaxane system in water solvent which suggests the increased stability of the 2-O isomer in polar solvent. explore their intricate way of functioning.10 Such efforts not only provide molecular level understandings of relevant biological processes but also result in charting new paradigms in terms of their potential applications. In [2]catenanes, circumrotation of the two interlocked rings controllable by oxidation has been exploited to generate the next generation of electronic displays based on the electrochromic response of these molecules.11 The underlying principle, in this particular example, has been the fact that the interchange between different types of host−guest interactions produce variable charge transfer and generation of RGB colors from the same molecule depending on the site-specific interactions.11 In a major breakthrough reported by Green et al.,12 they have demonstrated a monolayer of bistable [2]rotaxane molecules to be the data storage device in the design of a futuristic memory device patterned at an unmatched density of 1011 bits per square centimeter. Ideally, in all these molecular electronics

1. INTRODUCTION The incessant functioning of life processes in any living organism hinges upon, at the molecular level, seamless working of a particular class of molecules, among others, which carry out organized movements of molecular components with unprecedented control and precision. Examples of these molecular motors, as they are commonly known, include kinesin, myosin, titin, ATP synthase, etc.1−4 The complexity inherent in these molecules comes from the fact that they function in response to several external stimuli (pH variation, changes in solvent nature, temperature, radiation, etc.),5 and the relative motion of the motor components is highly regulated to accomplish some kind of meaningful cellular activity. These classes of molecules have inspired chemists to design artificial molecular machines,6 commonly referred to as mechanically interlocked molecules (MIMs).6−9 The most common examples of MIMs are [2]catenanes and [2]rotaxanes.8,9 The former consists of two mechanically interlocked rings, whereas in the latter there is a macrocyclic ring held onto a dumbbell-shaped molecule where the axis is terminated by two bulky groups at both ends serving as stoppers for the ring. In recent years, there have been consistent efforts in designing these classes of molecules and to © 2013 American Chemical Society

Received: July 3, 2013 Revised: November 1, 2013 Published: November 4, 2013 25059

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

skeleton (refer to Figure 1(a)) and their solvent dependence; (ii) to study the intermolecular hydrogen bonding that is

applications, shuttling between two distinct on/off states of the MIMs is of central importance.13 An unique feature of MIMs is bistability which comes from the interconversion between two distinct coconformations.13 Over the past few years, the template-directed synthesis of these molecules has been well explored in the solution phase.13 Typically, the two components of MIMs are held together by noncovalent interactions. In the ground state these two coconformations can be at equilibrium, and their distribution is decided by the free energy difference.13 A range of factors like structural flexibility, nature of the solvent, pH of the media, electronic environment, radiation, etc. can affect the distributions of coconformers. Hence, it is highly imperative to understand their dynamical behavior in solution, not only from the point of view of garnering deeper fundamental insights but also to tune them for a desired application. Keeping in mind the fact that the functioning of these molecular machineries sensitively depends on several factors, it is always challenging for the experimentalists to elucidate the structure−function relation with respect to the external environment like the effect of solvent. In a significant development by Rijs et al., they were able to provide detailed structural and binding information of isolated jet-cooled [2]rotaxanes based on gas-phase IR spectroscopic studies.14 Furthermore, a very recent attempt by Rijs et al. has added more understanding into the structure and binding interactions of two-station rotaxanes from their investigations based on high-resolution IR spectroscopy.15 The difficulty to extract such a wealth of information, however, is manifolded in the solution phase. Hence, it is immensely imperative to supplement experimental rationalizations with theoretical/computational approaches. In this context, molecular dynamics (MD) simulation is a particularly efficient tool to investigate their dynamical behavior. The group of Goddard and co-workers has done pioneering work on the atomistic simulations of [2]rotaxanes at several interfaces such as air/ water or on Au(111), and they were able to validate their findings with experiments.12−14 Such atomistic level details are crucial for computer-guided modeling of new-generation nanoscale devices. In another related study, Mendoza et al.16 have attempted to understand the adsorption of a [2]rotaxane on Au(111) and Ag(111) surfaces based on both MD simulations and experiments. In a more recent investigation by Zazza et al.,17 MD simulations and density functional theory (DFT) based studies provide valuable information regarding the hydrogen bonding (H-bond) interactions and conformational flexibility of a [2]rotaxane molecule in acetonitrile solution. Very recently, we came across the design of a novel diketopyrrolopyrrole (DPP) based [2]rotaxane by Raju et al.18 In this particular molecule there is the presence of an orthogonal hydrogen bond which sensitively depends on the solvent polarity.18 Additionally, this molecule showed excellent optical sensing properties of anions, particularly the F−.18 A remarkable naked-eye color change from pink to pale-green is reported when the macrocycle binds to the F− anion which is only specific to this anion due to the restricted dimension of the binding cavity of the macrocycle which is not suitable for other anions. On the basis of the 13C and 1H NMR spectra, the relative position of the macrocycle on the DPP rotaxane has been established.18 These interesting features prompted us to model and investigate the solvent-dependent dynamics of this [2]rotaxane system. Following this line of thought, the objectives of this manuscript are (i) to study the relative positions of the macrocycle on the molecular axle of the DPP

Figure 1. Molecular structure of the DPP skeleton (a) and macrocycle (b) of the [2]rotaxane (c) system. Part (a) also shows the labeling for different sites of the DPP skeleton based on which translational isomerism in the [2]rotaxane system is defined. Part (b) shows the hydrogen atoms of macrocycle that are involved in hydrogen bonding with O1, N1, O2, O3, and N2 sites of the DPP skeleton.

responsible for the solvent-dependent macrocycle contact positions (refer to Figure 1(b)); (iii) to study the conformational flexibility of this [2]rotaxane structure and the solvent effect; and (iv) to study the driving force for the formation of interlocked structures and to know whether these are thermodynamically stable. To address the solvent dependence of macrocycle contact positions and shuttle dynamics, we carry out extensive long time scale molecular dynamics simulations for the macrocycle (ring)−DPP skeleton (rod) complex (shown in Figure 1(c)) in different solvents. To understand the stability of the complexes, we carry out binding free energy calculations using molecular mechanics/Poisson−Boltzmann (and generalized Born) surface area approaches.19 The details of the two computational approaches employed are given in some detail in the next section. From our results, apart from being able to successfully reproduce the key experimental observations of Raju et al.,18 we have attempted to provide a generalized picture on the solvent-dependent shuttle dynamics and population of different translational isomers. These atomistic level findings are likely to stimulate further interest.

2. COMPUTATIONAL DETAILS 2.1. Molecular Dynamics Simulation. We have carried out molecular dynamics simulations for DPP-stoppered [2]rotaxane in three different solvents, namely, chloroform, dimethyl sulfoxide (DMSO), and water. The experimental 25060

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

for this analysis. Eventually, using this approach one computes the enthalpy change obtained in the association process. The free energy is computed using the end products where the binding free energy is obtained as the difference between the free energy of complex and individual molecules. The free energy of the complex is defined using the following expression

NMR spectra have been reported for this rotaxane complex only in the former two solvents.18 Our aim was to understand the behavior of this rotaxane system over the whole range of solvent polarity from chloroform to water. Independently, the DPP skeleton and macrocycle were optimized in three different solvents using the polarizable continuum to incorporate the specific solvent effect. Density functional theory at B3LYP/631G* level as implemented in Gaussian0920 has been used for these optimizations. Further, the DPP skeleton (refer to Figure 1(a)) associated with the macrocycle (refer to Figure 1(b)) has been optimized using the MO6/6-31G* basis set in chloroform solvent described as a polarizable continuum. The optimized geometry for the whole DPP rotaxane−macrocycle complex is shown in Figure 1(c) which has been used as the input geometry for molecular dynamics simulations in all three solvents. The charges fitted to reproduce molecular electrostatic potential using the CHELPG21 procedure have been used to get the atomic charges for individual units of the DPP skeleton and macrocycle. In particular, the set of charges used for these two units was based on the calculation in that specific solvent. The general amber force field22 has been used for both rotaxane units and chloroform and DMSO solvents to describe their nonbonded interaction. However, for water solvent the TIP3P force field23 has been employed. Overall three different simulations were carried out, each one corresponding to DPP rotaxane in a specific solvent. However, the starting structure for the macrocycle−DPP complex corresponds to the optimized geometry of [2]rotaxane in chloroform solvent (as shown in Figure 1(c)). It is interesting to note that, it has been reported experimentally the macrocycle binds to the pyridine site of the DPP rotaxane and is involved in two intermolecular hydrogen bonds. The simulations were carried out in orthorhombic boxes. In particular, for DPP rotaxane in chloroform solvent the box dimensions were 109, 86, and 84 Å and included 4043 solvent molecules. For the DMSO solvent the box lengths were 112, 81, and 88 Å which included 3204 DMSO molecules. Finally, in the case of water solvent, the simulation box dimensions were 108, 84, and 82 Å, and it included 19 093 water molecules. All three simulations were carried out in an isothermal−isobaric ensemble using Amber11 software.24 Primarily, equilibration runs were performed until the densities of all three [2]rotaxane−solvent systems reach convergence. The time step for integration of equation of motion was 1 fs, and the total time scale for simulation was kept around 50 ns. A cutoff of 10 Å was maintained for the calculation of nonbonded interactions. Instantaneous configurations after the equilibration run were stored for the purpose of analysis. 2.2. Binding Free Energy Calculation. Binding free energy calculations were carried out to understand the stability of the DPP rotaxane and macrocycle complexes. In the present case, the rotaxane has stoppers at both ends, and the interlocked structure has been generated during the synthesis. The internal ring diameter of the macrocycle is smaller than the size of the bulky groups of the DPP rotaxane, making the complex kinetically stable. However, in the situation where there are no stoppers at the end of rotaxane, the complex stability will play a major role in the functionality of the rotaxane. For this reason we have carried out the binding free energy calculations using molecular mechanics/Poisson− Boltzmann (and generalized Born) surface area approaches.19 Usually the binding free energy calculations are carried out in water solvent, and so we have specifically chosen water solvent

Gcomplex = UMM + GGB/PB + Gnp − TS

(1)

where Gcomplex is the free energy of the complex; UMM is the interaction energy of the complex; GGB/PB is the polar contribution of the solvation free energy; Gnp is the nonpolar contribution to the solvation free energy; and S is the entropy. Similar to the free energy of the complex, the free energies of macrocycle and rotaxane are defined, and binding free energies are obtained as a difference in these values. These two approaches, namely, the Poisson−Boltzmann and generalized Born, differ with respect to the polar and nonpolar terms of solvation free energy. All three terms, namely, the interaction energy and nonpolar and polar terms in the solvation free energy, are computed using the MMPBSA software.25 However, the entropy of the complex and its components are computed separately using normal-mode analysis as available in the NMODE module of MMPBSA software.25 In fact, the absolute free energies of complex and the components and binding free energy are computed for various snapshot configurations as obtained from the MD (the solvent molecules are stripped out in this analysis). In particular, the results on free energies presented here are based on few hundreds of such snapshots from the MD tajectory. 2.2.1. Relative Binding Free Energies of Macrocycle at Different Sites. The aforementioned calculation can only describe the energetics for the formation of the macrocycle− rotaxane complex system. However, to understand the population of different translational isomers (described by the relative position of macrocycle hydrogen bonding sites with respect to different hydrogen bonding acceptor sites of rotaxane), this is not sufficient. Either free energy calculation for all possible translational isomers separately should be carried out or the free energy profile across a reaction coordinate connecting different translational isomers has to be estimated. The latter one also provides information about the barrier heights involved in the interconversion between different translational isomers. However, in the usual approach a single degree of freedom is identified as a reaction coordinate, and sets of simulations are carried out for different predefined values along the reaction coordinate.26 Finally, the potential energy surface is obtained by combining the results from these simulations.26 This approach usually relies on computing the potential of mean force (pmf) within the framework of the umbrella sampling method.27 In one particular example, Zhang et al.26 have applied this approach to compute free energy profiles for a molecular machine based on cyclodextrin. However, in the present case, this approach cannot be employed since the transition from translational isomer 2-O to 2-R not only involves translation but also involves rotation of the macrocycle by 180° (refer to Figure 1(a)). So, we have carried out the binding free energy calculations for a macrocycle when it is bound to different hydrogen bonding acceptor sites of rotaxane (i.e., O1, N1, O2 sites as shown in Figure 1(a)). Since the conformational flexibility of the rotaxane will also strongly influence the energetics associated with different isomers, in this set of calculations the rotaxane has been treated as a rigid body, while the macrocycle and solvents were treated 25061

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

tional flexibility of the rotaxane itself and see whether it has any significant solvent dependence. Finally, we will address the stability of such rotaxanes based on the binding free energy calculations and also discuss the relative binding free energies of the macrocycle at different hydrogen bond acceptor sites of the DPP skeleton. Primarily, we have analyzed the most probable contact position of the macrocycle with the DPP skeleton and based on this established the population of different translational isomers in different solvents. The different binding sites and labeling of the specific atoms involved in the intermolecular hydrogen bonding between the macrocycle and DPP skeleton are presented in Figure 1(a) and 1(b). In this particular molecule, the favorable interactions that hold the macrocycle to the rotation axis are the hydrogen bonding between the acidic H atoms of the macrocycle ring with several N or O atoms present on the DPP skeleton of the rotaxane. Keeping this in mind, we have computed the distance distribution between the pyridine nitrogen (N) atom of the macrocycle and the nitrogen and oxygen atoms of the DPP skeleton (O1, N1, O2, O3, N2), and the results are shown in Figure 3. For simplicity, we will

using flexible molecular models. Overall, three sets of simulations were carried out for each of the translational isomers, namely, 2-O, 2-P, and 2-R in water solvent. The simulation protocol is the same as described in Section 2.1 except that in this set of simulations the rotaxane is treated as a rigid body. The representative snapshot configurations for these three sets of simulations are as shown in Figure 2(a), 2(b), and 2(c), respectively, depicting the translational isomers, namely, 2-O, 2-P, and 2-R.

Figure 3. Distance distributions of macrocycle pyrimidine nitrogen and selected oxygen and nitrogen sites (sites 1−5 as shown in Figure 1(a)) of the DPP skeleton in chloroform (a), DMSO (b), and water (c) solvents.

Figure 2. Four different translational isomers explored in this article. It appears like the molecular flexibility of the DPP skeleton controls the population of different isomers. The former three isomers referred to in the figure as (a), (b), and (c) (representing 2-O, 2-P, and 2-R isomers, respectively) are seen in the flexible DPP skeleton simulations, while the fourth isomer 2-S, (d), is observed in the rigid DPP skeleton simulations.

refer to different contact sites of DPP as sites 1−5 as shown in Figure 1(a): two of these are nitrogen sites, and the remaining three are oxygen sites. Interestingly, in the case of chloroform (refer to Figure 3(a)) and DMSO (refer to Figure 3(b)) solvents, we find that the first three sites are relevant since these are involved in the short length contacts with the macrocycle. In the case of water solvent, almost all the contact distances between the macrocycle and DPP sites appear within a distance r < 10 Å (see Figure 3(c)). By comparing Figures 3(a) and 3(b), it is clearly evident that in chloroform solvent the macrocycle pyridine nitrogen atom is mostly in closer contact with site 2 (N1) which is in complete agreement with experimental reports which are based on proton NMR spectra.

3. RESULTS AND DISCUSSIONS In this section, we will discuss the solvent polarity effect on the macrocycle contact position with that of the DPP skeleton, and we aim to understand the solvent-dependent distribution of different translational isomers by analyzing the intermolecular hydrogen bonding patterns. We will also address the conforma25062

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

In the original experimental paper,18 this translational isomer was referred as 2-P (refer to Figure 2(b)). The sites 1 and 3 appear as first neighboring sites. However, in DMSO and water solvent, the macrocycle pyridine nitrogen atom is in closer contact with both sites 1 and 2 which is seen from the two distribution curves having peak position at r = 3.5 Å (see Figure 3(b)). We also find that it is located closer to site 3 with a nonzero population. Overall, our results only agree partly with the experimental reports which predicted the population of a translational isomer referred to as 2-R, and this refers to a situation where the DPP macrocycle is in closest contact with site 3 (O2) of DPP. However, the force-field-based simulation studies rather predict that in DMSO solvent there exist three translational isomers, namely, 2-P, 2-R, and a newer one (we refer to this as 2-O), with their population decreasing in the order site 1 > site 2 > site 3 (refer to Table 1). The values in

Most importantly, our simulations are able to capture these intricate features, although the experimental observation of only one translational isomer can be thought to be representing only one facet of the overall dynamical phenomena. There is no experimental report in the case of water solvent for the possible comparison of our results in this solvent. It is also worth noticing that the population of any of the translation isomers does not exceed 30%, suggesting that the rotaxane is not stabilized much due to intermolecular hydrogen bonding between the macrocycle and DPP skeleton. This is not the case in chloroform solvent where the population of translational isomer 2-P is almost 100%, suggesting that the rotaxane system is stabilized by intermolecular hydrogen bonding between the macrocycle and DPP. Considering the nonpolar nature of chloroform and the inability to make hydrogen bonding, these results look obvious. However, in DMSO solvent there appears to be competition between the macrocycle and DMSO solvent to make intermolecular hydrogen bonding with DPP which eventually contributes to the reduction of different translational isomers. Interestingly, in water solvent, the macrocycle pyridine nitrogen is in closer contact with both sites 1 and 2 and contributes to the population of translational isomers 2-O and 2-P, respectively, to 48% and 91%. Interestingly, the translation isomer 2-R is also populated up to 15%. Another striking observation is that in DMSO solvent the distribution curves (refer to Figure 3(b)) have a larger width when compared to chloroform (refer to Figure 3(a)) and water (refer to Figure 3(d)), suggesting that the macrocycle covers a larger translational distance along the axis. It is surprising to see that in water solvent the distribution curves are very narrow. Overall, the translational mobility of the macrocycle along the axis decreases in the order DMSO > chloroform > water which shows a nonmonotonic behavior with increasing solvent polarity. It is important to understand the solvent dependence of the macrocycle contact with different polar sites along the translational axis of the DPP skeleton. By looking into the hydrogen bond making ability of the solvents studied, it can be easily proposed that the interplay between intermolecular hydrogen bonding between the macrocycle and DPP skeleton with solvents can be a key parameter for the observation of different translational conformers observed. In the case of chloroform, since there is no possibility for hydrogen bonding between the macrocycle−solvent and DPP−solvent, the rotaxane−macrocycle complex is stabilized by hydrogen bonding between the macrocycle hydrogen atoms (NH1 and NH2) with site 2. However, an interesting question arises why in this case it does not make any intermolecular hydrogen bonding through site 1 as in the case of DMSO solvent. In the case of DMSO and water solvents, the intermediate states for the macrocycle−rotaxane complex (between the two translational conformers, namely, 2-R and 2-P) can be stabilized through intermolecular hydrogen bonding with solvents. This eventually reduces the kinetic barrier for the conformational transition between conformers 2-O and 2-P in the case of DMSO and water solvents and so making the conformational jumps a frequent event. However, in the case of chloroform solvent, there exists a significant kinetic barrier since the intermediate states between the conformers 2-P and 2-R cannot be stabilized through hydrogen bonding with solvent. It is interesting to refer to a study where a similar solvent polarity driven decrease in the activation energy for conformational jumps has been discussed in the case of adenosine.28

Table 1. Average Population of Different Translational Isomers solvent

2-O

2-P

2-R

chloroform DMSO water

9 26 48

98 12 91

9 1 15

Table 1 were computed by integrating the distance distribution curves in Figure 3 up to a distance value, r < 4 Å. Each of these distribution curves is normalized to unity, and the integration of this curve over a range of distances gives the probability of finding the macrocycle in this distance range. However, the population of different translational isomers is presented as a percentage which is obtained by multiplying the integrated value by a factor 100. The macrocycle and DPP system are said to be in the associated state when the distance is less than 4 Å . Moreover, here the translational isomers 2-O, 2-P and 2-P, 2-R are not mutually exclusive, while the isomers 2-O and 2-R are mutually exclusive. The reason is that the macrocycle can occupy a site where it is mutually hydrogen bonded to both sites 1 and 2. However, the DPP skeleton has to undergo conformational change such that O1 and N2 sites rearrange to one side. In the starting configuration as seen in Figure 1(a) these two sites are arranged in opposite directions. As can be seen from Table 1, we predict a larger population of translational isomer 2-O, while the experimental reports suggest a 2-P isomer. The reason for this can be straightforward. The experimental report is based on the proton NMR spectra where the population of different translational isomers is based on the upfield or downfield effect of DPP hydrogens. A downfield effect seen in the hydrogens near sites 1−5 has been attributed to the increase in the population of specific isomers. However, in the case of solvents which can involve in hydrogen bonding with the DPP, hydrogens can also contribute to such effects, and so the population of different isomers has to be calculated carefully. The ambiguity seen between the experimental and our results has to be attributed to the dynamical nature of the system. It is often difficult to understand the structure of a dynamical system using NMR data which only provide an averaged picture. Moreover, if we carefully analyze the sites 1 and 3, we will notice that these two oxygen atoms are almost in similar chemical environment as both are connected to a methylene group (−CH2−) and phenyl groups in either side. Hence, in a dynamical scenario involving solvents, the macrocycle has the possibility to be present at both the sites. 25063

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

To understand the role of intermolecular hydrogen bonding in stabilizing the translational isomers 2-O, 2-P, and 2-R, we have analyzed the possible hydrogen bonding contacts between the macrocycle and DPP, and the results are shown in Figure 4.

Figure 4. Distance distribution of NH1···O1, NH2···O1, NH1···N1, NH2···N1, NH1···O2, and NH2···O2 contacts in chloroform (a), DMSO (b), and water (c) solvents.

Figure 5. Contour plots of hydrogen bonding contacts in chloroform, DMSO, and water solvents.

For each solvent, there are six different distance distribution curves: The distances used for the analysis are computed between the two macrocycle NHs with the three sites, namely, 1, 2, and 3. As can be seen in chloroform solvent (refer to Figure 4(a)), the short hydrogen bonding distances are due to macrocycle NH (NH1 and NH2) with site 2 of DPP. It can also be seen that there is nonzero population of translational isomers stabilized due to contacts of NH1 and NH2 hydrogens with oxygen sites (i.e., sites 1 (O1) and 3 (O2)). The distance distribution curves for both NH···N contacts look almost similar, suggesting that the energetics of both intermolecular hydrogen bondings are comparable which looks obvious by considering the symmetry of the macrocycle. In the case of DMSO (Figure 4(b)) solvent, the short distance hydrogen bondings are due to the four possibilities such as (i) macrocycle NH1···O1, (ii) macrocycle NH2···O1 of site 1, (iii) macrocycle NH1···N1, and (iv) macrocycle NH2···N1. It is interesting to notice that the contacts of NH1 and NH2 with that of site 3 (O2) are not very significant in the case of DMSO. However, these contacts (i.e., macrocycle NH1···O2 and macrocycle NH2···O2) also contribute significantly in the case of water solvent (refer to Figure 4(c)). Figure 5 shows the contour plots of g(R1,R2) where R1 refers to the distance of macrocycle hydrogens (NH1 and NH2) with site- (O1), while R2 refers to the distance with that of site 2 (N1). As can be seen from Figure 5, in chloroform solvent, the high density spot in this contour plot is associated with intermolecular hydrogen bonding between the macrocycle and site 2 and corresponds to the dominant population of the 2P translational isomer. However, in the case of water solvent the high density spot is associated with the hydrogen bondings between site 1 and site 2 and so reveals the dominant population of 2-O and 2-R isomers. On the contrary, in DMSO

solvent, the conformer corresponding to the larger density does not have any intermolecular hydrogen bonding which again confirms the results on the populations of different translational isomers as presented in Table 1. Finally, we aimed to address the conformational flexibility of the DPP in different solvents, and for this we have computed the end-to-end distance (i.e., the distance between the atoms referred to as E1 and E2 (refer to Figure 1(a))). The results are shown in Figure 6. Interestingly, the average end-to-end distance (Rend‑to‑end) decreases in different solvents in the order

Figure 6. End-to-end distance distribution of [2]rotaxane in different solvents. 25064

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

Figure 7. Representative (a) fully stretched conformer, (b) fibril-like folded conformer, and (c) trigonal-like folded conformer. The stretched conformers were often seen in chloroform solvent, while the folded conformers were often observed in water solvent.

peak beyond this distance. The representative conformers for these two cases are shown in Figure 7(b) and 7(c) along with a fully stretched conformer (refer to Figure 7(a)) in chloroform solvent. To identify the driving force for the formation of such rotaxanes, the binding free energy calculations were carried out in water solvent, and the results are shown in Tables 2 and 3. The macrocycle (guest), DPP (host), and supramolecular free energies and the difference in binding free energies are shown. On the basis of the values for Ggas and Gsolv terms in Table 2(a) and 2(b), it appears that the supramolecule formation is mostly driven by the intermolecular interaction between the macrocycle and DPP, and solvation free energy contributions act against the complex formation. In particular, to the binding free energy, the van der Waals interaction contributes dominantly (≈ −40 kcal/mol) when compared to the electrostatic interaction (≈ −8 kcal/mol). This clearly suggests that the macrocycle−rotaxane complex system is stabilized by a hydrophobic-like interaction rather than a hydrophilic-type interaction. It also appears that the solvation free energy acts against such complex formation, and such contribution amounts to 11−14 kcal/mol. Overall, the binding free energy predicted from both MM/GB(PB)SA approaches amounts to

chloroform > DMSO > water, and the values correspond to 20.1, 15.2, and 9.6 Å, respectively. Moreover, there appear to be two and three distinct conformers in the case of water and DMSO, respectively, while the chloroform structure is dominated by only a single conformer. It is really striking to see that the Rend‑to‑end can be reduced to 7.5 Å in water solvent which is three times smaller than in the case of chloroform solvent. This is only possible in that the DPP skeleton folds in water solvent from a fully stretched form as in chloroform solvent. This also explains why the distance distribution curves (as seen in Figure 3(c)) are narrower in water solvent. It is also worth noticing that such elongated and folded conformers have been reported for (bam)succ-[C12 alkyl spacer]-ni-oxidized molecular rotaxane in methyl cyanide solution from a previous molecular dynamics study.17 However, in the present case such stretched-to-folded conformational transition in the [2]rotaxane is driven by the solvent polarity. We have analyzed the conformers responsible for the two peaks appearing in the end-to-end distance distribution in water solvent. The conformers with Rend‑to‑end less than and greater than 10.2 Å have been analyzed. A misfolded fibril-like conformer29 appears to be responsible for the peak appearing for a distance less than 10.2 Å, and a trigonal-like conformer seems to be responsible for the 25065

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

Finally, we have investigated the relative binding free energies of macrocycles in the hydrogen bond acceptor sites of the DPP skeleton (i.e., sites 1, 2, and 3). The results are presented in Table 4(a) and 4(b) and, respectively, correspond to the generalized Born and Poisson−Boltzmann approaches. Only the binding free energies of macrocycle and DPP complex are shown, and individual subsystem free energies and complex free energy are not shown. The results shown here in the table differ significantly when compared to results in Table 2(a) and Table 2(b) which clearly explains the role of conformational flexibility of the DPP skeleton. The results shown in Table 4(a) and 4(b) use trajectories of three sets of simulations with a rigid DPP skeleton. To verify that the binding free energies calculated are for the specific translational isomer (such as 2O, 2-P, and 2-R), we have computed the NH1(or NH2)···X distance distributions (refer to Figure 8) where X refers to O1,

Table 2. Binding Free Energy (in kcal/mol) Computed from the (a) MM/GBSA and (b) MM/PBSA Approach energy component (a) bond angle dihed EVDWAALS EEL EGB ESURF Ggas Gsolv total (b) bond angle dihed EVDWAALS EEL EPB ECAVITY Ggas Gsolv total

complex

guest

host

difference

54.93 165.24 65.82 −69.32 −203.73 −44.11 9.39 12.93 −34.71 −21.77

39.29 134.44 44.45 −22.78 −29.24 −35.49 8.30 166.15 −27.19 138.95

15.64 30.80 21.36 −6.50 −166.50 −23.79 4.93 −105.20 −18.86 −124.06

−0.00 −0.00 0.00 −40.03 −7.98 15.18 −3.84 −48.01 11.34 −36.67

54.93 165.24 65.82 −69.32 −203.73 −45.94 8.00 12.93 −37.94 −25.00

39.29 134.44 44.45 −22.78 −29.24 −37.23 7.43 166.15 −29.80 136.35

15.64 30.80 21.36 −6.50 −166.50 −25.87 3.64 −105.20 −22.22 −127.43

−0.00 −0.00 0.00 −40.03 −7.98 17.15 −3.07 −48.01 14.08 −33.93

Table 3. Entropy (in kcal/mol) Calculations for the [2]Rotaxane entropy

complex

guest

host

translational rotational vibrational total

14.45 13.18 154.81 182.46

14.17 12.84 110.44 137.46

13.29 11.45 41.04 65.79

difference

Figure 8. Distance distribution of NH1 (or NH2)···O1, NH1 (or NH2)···N1, and NH1 (or NH2)···O3 contacts. Respectively, the first translational isomer is referred to as 2-O, and the second and third are referred to as 2-P, and 2-R. Note the bimodal distribution in the case of 2-R where in fact the second peak corresponds to a new translational isomer referred to as 2-S.

−20.79

−34 to 37 kcal/mol. The complex formation results in the decrease of microstates available and otherwise for the individual host−guest subsystems, and so the entropy acts against the complex formation as suggested from the computed value of −21 kcal/mol (refer to Table 3).

N1, and O2 sites of the DPP skeleton, respectively. As can be seen for the first two cases the unimodal distance distribution curve, it suggests that there are no conformational jumps, and during the entire simulation the complex remains in a specific

Table 4. Binding Free Energy (in kcal/mol) Computed from the (a) MM/GBSA and (b) MM/PBSA Approacha energy component (a) EVDWAALS EEL EGB ESURF Ggas Gsolv GTotal (b) EVDWAALS EEL EPB ECAVITY Ggas Gsolv GTotal a

2-O

2-P

2-R

2-S

−27.63(0.13) −1.46(0.08) 6.28(0.03) −2.96(0.01) −29.09(0.15) 3.32(0.03) −25.77(0.14)

−22.15(0.17) −2.49(0.15) 6.91(0.06) −2.38(0.02) −24.64(0.18) 4.53(0.06) −20.11(0.15)

−20.67(0.20) −1.77(0.15) 6.23(0.07) −2.20(0.01) −22.43(0.26) 4.03(0.05) −18.41(0.23)

−21.06(0.22) −1.80(0.10) 6.03(0.07) −2.16(0.02) −22.86(0.22) 3.87(0.06) −18.99(0.22)

−27.63(0.13) −1.46(0.08) 6.75(0.04) −2.11(0.01) −29.09(0.15) 4.64(0.04) −24.45(0.14)

−22.15(0.17) −2.49(0.15) 7.29(0.10) −1.58(0.03) −24.64(0.18) 5.71(0.04) −18.93(0.15)

−20.67(0.20) −1.77(0.15) 6.50(0.08) −1.25(0.01) −22.43(0.26) 5.25(0.05) −17.18(0.23)

−21.06(0.22) −1.80(0.10) 5.88(0.07) −1.15(0.02) −22.86(0.22) 4.73(0.06) −18.13(0.22)

The numbers in the parentheses refer to the standard error. 25066

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

with a flexible molecular model the scenario can be different. Overall, this set of calculations clearly demonstrates the importance of solvent polarity in the binding free energy of different translational isomers and explains the reason for the solvent-polarity-induced conformational population.

translational isomeric form (i.e., as 2-O or 2-P). However, in the third case, the distribution curve is rather bimodal suggesting the possibility of conformational transition to another isomer which we refer to as 2-S. A careful analysis of the structure suggests that in this isomer there is hydrogen bonding between NH1 (or NH2) atoms of the macrocycle and the O3 site of the DPP skeleton (refer to Figure 1(a)). A snapshot configuration for this translational isomer is shown in Figure 2(d). The reason for the conformational transition between the 2-R and 2-S can be understood easily by looking into the skeleton structure which shows that these two sites are arranged in the same side of the skeleton. It is also interesting to note that the macrocycle travels around 4.8 Å during this transition. During a flexible body simulation, due to the dihedral rotation this distance can be reduced significantly, and so one should see much faster conformational jumps. However, there is another scenario where the folding of the DPP skeleton may make this transition infrequent. In fact, this folding should be the reason why we are not observing a 2-S isomer in a set of flexible body simulations. So, the rigid and flexible body simulations provide some useful information for an effective design of molecular machines. The conformational flexibility of the DPP skeleton is going to play a major role in dictating the possible conformational states allowed for the macrocycle. The rigidity of the skeleton and presence of favorable binding sites in it (in the same side) make the shuttling between several translational isomers smooth which is central for any desired applications. This can be accomplished by structural modifications of the skeleton. So, we suggest the DPP skeleton to be modified with alternative single double bonds where the hydrogens are replaced with functional groups (that can serve as a hydrogen bond acceptor) to stabilize different conformational states. Alternatively, as we have learned from our ongoing discussions, the solvent polarity can also be a deciding factor. Solvent which favors more flexibility in the skeleton structure (particularly leading to the folded state) may not be a suitable choice for designing a molecular machine. Below, we will discuss the binding free energies of different translational isomers (including 2-S). Independent of the binding free energy calculation method employed (refer to Table 4(a) and 4(b)), the relative binding free energy decreases in the order site 1 > site 2 > site 4 ≈ site 3. The binding energies of different translational isomers obtained from force-field-based modeling and so their relative population should be carefully interpreted since in some cases the difference in relative binding strengths may be at the order of errors involved in the derivation of such force fields. However, in the current study the differences in binding free energies of isomers 2-O, 2-P, and 2-R are in the order of 1.7− 7.4 kcal/mol, while the standard errors are in the order of 0.14−0.23 kcal/mol. It is important to remember that these calculations correspond to water solvent. Even though the Ggas term will not be affected in different solvents, the Gsolv can be altered significantly. We have also tested this within the MM/ PBSA approach by varying the dielectric constant of external media between 5 and 80 where the lowest value corresponds to chloroform solvent and the highest value corresponds to water solvent. In this case, the solvation free energy varied from 4.2 to 6.2 kcal/mol. Even though such change in solvation free energy cannot alter the position of 2-O, it can change the energetics of translational isomers 2-R, 2-P, and 2-S and so the population of these isomers in different solvents. However, this discussion is based on the rigid molecular model of the DPP skeleton, and

4. CONCLUSIONS Force-field molecular dynamics simulations were carried out to understand the solvent-dependent distribution of translational isomers of the DPP [2]rotaxane system. Among the three possible isomers, 2-P dominates in the case of chloroform solvent which is in complete agreement with the experimental reports. However, in DMSO solvent the simulations predict the 2-O isomer to dominate the population, while the experimental report predicts the 2-R isomer to be dominant. We attribute this ambiguity to the interpretation of proton NMR spectra where the upfield and downfield shifts in protons cannot be just discussed based on the guest−host complex alone since a polar solvent like DMSO which has a tendency to form hydrogen bonding also contributes to this effect. In water solvent we predict comparable populations of 2-O and 2-P translational isomers. Moreover, we also find that the [2]rotaxane system exhibits a dramatic solvent-induced conformational transition from a fully stretched to folded state. On the basis of the binding free energy calculations, the driving force for the macrocycle−DPP, a guest−host supramolecular complex formation, have been discussed. The binding free energy calculations for macrocycle on different hydrogen bonding acceptor sites of the DPP skeleton reveal that the conformational flexibility of the DPP skeleton and polarity of the media can alter the binding free energies of the corresponding translational isomers and so can change their relative population. This explains the solvent polarity induced difference in translational isomer distribution in this [2]rotaxane system. On the basis of the simulations with a rigid and flexible molecular model for the DPP skeleton we suggest a rigid DPP skeleton to serve as a better molecular switch. In summary, our report can provide fundamental insights into the dynamical behavior of molecular motors in response to solvent polarity and thereby can also supplement toward a decent understanding of their biological counterparts. Additionally, it is likely to give more credence to the growing importance of MIMs in the fabrication of nanoscale devices in general.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS N.K.J. gratefully acknowledges the Carl Tryggers Foundations for the funding and also Michael Odelius, Fysikum, Stockholm University, for the support. This work was supported by a grant from the Swedish Infrastructure Committee (SNIC) for the project “Multiphysics Modeling of Molecular Materials”, SNIC 023/07-18.

■

REFERENCES

(1) van den Heuvel, M.; Dekker, C. Motor Proteins at Work for Nanotechnology. Science 2007, 317, 333−336.

25067

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068

The Journal of Physical Chemistry C

Article

(23) Jorgensen, W. L.; Madura, J. D. Quantum and statistical mechanical studies of liquids. 25. Solvation and conformation of methanol in water. J. Am. Chem. Soc. 1983, 105, 1407−1413. (24) Case, D. et al. Amber 11; University of California: San Francisco. 2010. (25) Miller, B. R.; McGee, T. D.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E. MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 3314− 3321. (26) Zhang, Q.; Tu, Y.; Tian, H.; Zhao, Y.-L.; Stoddart, J. F.; Ågren, H. Working Mechanism for a Redox Switchable Molecular Machine Based on Cyclodextrin: A Free Energy Profile Approach. J. Phys. Chem. B 2010, 114, 6561−6566. (27) Torrie, G. M.; Valleau, J. P. Monte Carlo free energy estimates using non-Boltzmann sampling: Application to the sub-critical Lennard-Jones fluid. Chem. Phys. Lett. 1974, 28, 578−581. (28) Murugan, N. A.; Hugosson, H. W. Solvent Dependence of Conformational Distribution, Molecular Geometry, and Electronic Structure in Adenosine. J. Phys. Chem. B 2009, 113, 1012−1021. (29) Murugan, N. A.; Olsen, J. M. H.; Kongsted, J.; Rinkevicius, Z.; Aidas, K.; Ågren, H. Amyloid Fibril-Induced Structural and Spectral Modifications in the Thioflavin-T Optical Probe. J. Phys. Chem. Lett. 2013, 4, 70−77.

(2) Sweeney, H.; Houdusse, A. Structural and Functional Insights into the Myosin Motor Mechanism. Annu. Rev. Biophys. 2010, 39, 539−557. (3) Lv, S.; Dudek, D.; Cao, Y.; Balamurali, M.; Gosline, J.; Li, H. Designed biomaterials to mimic the mechanical properties of muscles. Nature 2010, 465, 69−73. (4) Itoh, H.; Takahashi, A.; Adachi, K.; Noji, H.; Yasuda, R.; Yoshida, M.; Kinosita, K. Mechanically driven ATP synthesis by F1-ATPase. Nature 2004, 427, 465−468. (5) Wang, Q.-C.; Qu, D.-H.; Ren, J.; Chen, K.; Tian, H. A Lockable Light-Driven Molecular Shuttle with a Fluorescent Signal. Angew. Chem., Int. Ed. 2004, 43, 2661−2665. (6) Balzani, V.; Credi, A.; Raymo, F.; Stoddart, J. Artificial Molecular Machines. Angew. Chem., Int. Ed. 2000, 39, 3348−3391. (7) Lehn, J. M. Supramolecular Chemistry: Concepts and Perspectives; Wiley-VCH: New York, 1995. (8) Amabilino, D. B.; Stoddart, J. F. Interlocked and Intertwined Structures and Superstructures. Chem. Rev. 1995, 95, 2725−2828. (9) Raymo, F. M.; Stoddart, J. F. Interlocked Macromolecules. Chem. Rev. 1999, 99, 1643−1664. (10) Berná, J.; Leigh, D. A.; Lubomska, M.; Mendoza, S. M.; Pérez, E. M.; Rudolf, P.; Teobaldi, G.; Zerbetto, F. Macroscopic transport by synthetic molecular machines. Nat. Mater. 2005, 4, 704−710. (11) Ikeda, T.; Stoddart, J. Electrochromic materials using mechanically interlocked molecules. Sci. Technol. Adv. Mater. 2008, 9, 014104. (12) Green, J. E.; Wook Choi, J.; Boukai, A.; Bunimovich, Y.; Johnston-Halperin, E.; DeIonno, E.; Luo, Y.; Sheriff, B. A.; Xu, K.; Shik Shin, Y.; et al. A 160-kilobit molecular electronic memory patterned at 1011 bits per square centimetre. Nature 2007, 445, 414− 417. (13) Fahrenbach, A. C.; Bruns, C. J.; Cao, D.; Stoddart, J. F. GroundState Thermodynamics of Bistable Redox-Active Donor–Acceptor Mechanically Interlocked Molecules. Acc. Chem. Res. 2012, 45, 1581− 1592. (14) Rijs, A. M.; Compagnon, I.; Oomens, J.; Hannam, J. S.; Leigh, D. A.; Buma, W. J. Stiff, and Sticky in the Right Places: Binding Interactions in Isolated Mechanically Interlocked Molecules Probed by Mid-Infrared Spectroscopy. J. Am. Chem. Soc. 2009, 131, 2428−2429. (15) Rijs, A. M.; Kay, E. R.; Leigh, D. A.; Buma, W. J. IR Spectroscopy on Jet-Cooled Isolated Two-Station Rotaxanes. J. Phys. Chem. A 2011, 115, 9669−9675. (16) Mendoza, S. M.; Whelan, C. M.; Jalkanen, J.-P.; Zerbetto, F.; Gatti, F. G.; Kay, E. R.; Leigh, D. A.; Lubomska, M.; Rudolf, P. Experimental and theoretical study of the adsorption of fumaramide [2]rotaxane on Au(111) and Ag(111) surfaces. J. Chem. Phys. 2005, 123, 244708. (17) Zazza, C.; Mancini, G.; Brancato, G.; Sanna, N.; Barone, V. Neutral molecular shuttle in acetonitrile dilute solution investigated by molecular dynamics and density functional theory. Comput. Theor. Chem. 2012, 985, 53−61. (18) Raju, M. V. R.; Lin, H.-C. A Novel Diketopyrrolopyrrole (DPP)-Based [2]Rotaxane for Highly Selective Optical Sensing of Fluoride. Org. Lett. 2013, 15, 1274−1277. (19) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; et al. Calculating Structures and Free Energies of Complex Molecules: Combining Molecular Mechanics and Continuum Models. Acc. Chem. Res. 2000, 33, 889−897. (20) Frisch, M. J. et al. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (21) Breneman, C.; Wiberg, K. Determining Atom-Centered Monopoles from Molecular Electrostatic Potentials. The Need for High Sampling Density in Formamide Conformational Analysis. J. Comput. Chem. 1990, 11, 361−373. (22) Wang, J.; Wolf, R.; Caldwell, J.; Kollman, P.; Case, D. Development and testing of a general AMBER force field. J. Comput. Chem. 2004, 34, 1157−1174. 25068

dx.doi.org/10.1021/jp406576m | J. Phys. Chem. C 2013, 117, 25059−25068