Force probe simulations using an adaptive resolution scheme

Molecular simulations of the forced unfolding and refolding of biomolecules or molecular complexes allow to gain important kinetic, structural and thermodynamic information about the folding process and the underlying energy landscape. In force probe molecular dynamics (FPMD) simulations, one pulls one end of the molecule with a constant velocity in order to induce the relevant conformational transitions. Since the extended configuration of the system has to fit into the simulation box together with the solvent such simulations are very time consuming. Here, we apply a hybrid scheme in which the solute is treated with atomistic resolution and the solvent molecules far away from the solute are described in a coarse-grained manner. We use the adaptive resolution scheme (AdResS) that has very successfully been applied to various examples of equilibrium simulations. We perform FPMD simulations using AdResS on a well studied system, a dimer formed from mechanically interlocked calixarene capsules. The results of the multiscale simulations are compared to all-atom simulations of the identical system and we observe that the size of the region in which atomistic resolution is required depends on the pulling velocity, i.e. the particular non-equilibrium situation. For large pulling velocities a larger all atom region is required. Our results show that multiscale simulations can be applied also in the strong non-equilibrium situations that the system experiences in FPMD simulations.


Introduction
Force spectroscopy is a standard experimental technique to investigate unfolding pathways, details of the energy landscape and the mechanical properties of single biomolecules and molecular complexes [1][2][3]. Usually, one end of the molecule is fixed in space and an external force is applied to the other end. This force is either constant (force clamp) or changes linearly in time (force ramp, FR). In the first case one observes the extension of the system as a function of time while in the FR protocol the force measured at the pulling device is recorded as a function of the extension [4]. The information extracted from these types of experiments not only allows * Author to whom any correspondence should be addressed.
to determine (un)folding rates, but also important properties of the folding landscape like the position of transition states, transition path times or the existence of stable intermediates [5,6].
Molecular dynamics (MD) simulations are routinely used to investigate conformational transitions in soft matter systems and in particular the folding and unfolding of biomolecules like peptides, proteins or RNA [7,8]. The mechanical properties of biological systems can be studied with atomistic resolution using the techniques of FPMD simulations (also called steered MD simulations) [9,10]. In most cases, however, there is a gap of up to five orders of magnitude in the time scales of such simulations and experimental realizations of force spectroscopy [11]. Only recently it has become possible to match the time scales of FPMD simulations and experiments performed by using a high-speed atomic force microscope to study the unbinding of a streptavidin-biotin complex [12]. In general, however, it is challenging to reach the experimentally relevant long time scales using atomistic FPMD simulations. One reason lies in the need of a rather large simulation box and the resulting large number of solvent molecules. Additionally, a large number of (un)folding trajectories are required to allow for a meaningful statistical analysis of the results. In order to speed up FPMD simulations, some of the well established techniques of coarse graining have successfully been applied to study the mechanical folding pathways of proteins and of RNA [13,14]. However, typical coarse grained models employ simplified interaction potentials and reduced numbers of particles. Therefore, the atomistic details of the formation and rupture of noncovalent bonds responsible for conformational transitions cannot be resolved. Markov state models allow to study the kinetics of conformational transitions on long time scales using dynamical information from short atomistic simulation runs [15,16] and they have successfully been applied to extent the dynamical range of FPMD simulations [17,18]. Also methods that are developed directly to increase the efficiency of FPMD simulations are available [19,20].
In most cases the primary interest of FPMD simulations lies in the study of the mechanical (un)folding kinetics of the solute and the dynamics of the solvent molecules only plays a minor role. Therefore, mixed resolution schemes that treat the solute in an all-atom (AA) manner and the solvent in a coarse-grained (CG) way should be applicable. There are different methods to set up mixed resolution schemes. One hybrid method that is particularly well suited to treat systems in which no exchange between particles treated with different resolution takes place uses the definition of virtual interaction sites [21]. These virtual sites are positioned at the center of mass of a group of atoms in the AA part of the system. The CG forces acting on the virtual sites are then distributed uniformly among the neighboring atoms to achieve the coupling between the different parts of the system. We have applied this methodology to the special non-equilibrium situation encountered in FPMD simulations and we have found that the scheme is applicable in principle but the accuracy is not comparable to the one in equilibrium situations [22]. Other approaches to enable simulations with mixed resolution are (among others) the method introduced by Izvekov and Voth [23] and the adaptive resolution scheme (AdResS) developed by Praprotnik et al [24] and by Krekeler et al [25]. This method has recently been applied to the computation of solvation free energies and it proved to yield reliable results in a non-equilibrium situation [26]. In the present paper, we use the AdResS methodology for a CG description of the solvent allowing the solute to be described in an AA manner in FPMD simulations. The AdResS is based on a partitioning of the simulation box into a region with AA resolution, one with CG resolution and a crossover or hybrid region which allows for particle exchange between the regions of different resolution. Thus, in our study, the AA region consists of the solute and some surrounding solvent molecules while the solvent that is not in the immediate neighborhood of the solute is treated in a CG manner. As in FPMD simulations the system is strongly driven out of thermal equilibrium, it is not clear a priori that methods originally developed for equilibrium simulations are also applicable in these situations.
As in earlier investigations, we will perform FPMD simulations using a calix [4]arene catenane dimer in mesitylene solvent as a model system. Apart from the mentioned hybrid simulations, we have investigated this system in AA simulations and have found that its reversible unfolding kinetics can well be understood in terms of a simple two-state model [27][28][29]. In equilibrium, the two calix [4]arene 'cups' form a complex stabilized by a ring of 16 hydrogen bonds (H-bonds). Complete dissociation of the cups is prevented by a set of four intertwined aliphatic loops consisting of 14 methylene groups each. We will use mesitylene as a solvent because due to its aprotic nature it does not form H-bonds that interfere with the intramolecular H-bonds between the two calix [4]arene monomers.
The paper is organized as follows. In the next section, the computational details are presented including a brief recapitulation of the AdResS methodology and the presentation of results of equilibrium AdResS simulations of the calix[4]arene dimer system. We then compare the results of FPMD simulations performed employing AdResS to those of AA simulations and close with some concluding remarks.

All-atom simulations
All AA simulations were performed using the GROMACS 2018.4 program package employing the OPLS-AA force field [30][31][32]. We used a stochastic dynamics integrator [33] at a temperature of 298 K with a friction constant of 0.1 ps. All bonds were constrained using the LINCS algorithm [34] allowing for a time step of 2 fs. Short range electrostatic and van der Waals interactions were computed using a cut-off of 1.2 nm. The long range Coulomb interactions was treated using the reaction field method with a relative dielectric constant of 2.4. For the van der Waals interactions, we applied a dispersion correction [35]. For the AA simulations the neighbor list was updated after 25 simulation steps and for the CG simulations a pair list with a cutoff of 1.37 nm was used. We used Cartesian periodic boundary conditions in all simulations. We performed an energy minimization starting with a (7.5 nm) 3 cubic simulation box containing one calix [4]arene catenane dimer and 1780 mesitylene molecules, cf figure 1. After this, the system was equilibrated in the canonical ensemble for 200 ps. Then it was coupled to a Berendsen barostat [36] with a time constant of 0.5 ps and an isothermal compressibility of 8.26 × 10 −5 bar −1 . The box size was determined to be (7.49 nm) 3 for a pressure of 1 bar and this was used in all AA simulations. All production runs were performed in the canonical ensemble using this box size. We mention that this box size is larger than the box size we used in earlier investigations ((5.8 nm) 3 ) [37,38]. However, the larger box size used here is also used for the AdResS simulations and therefore a direct comparison is possible.  [4]arene dimer along with the definition of the reference group and the pulled group used in the setup of the FPMD simulations. These groups are defined as the center of mass of the methoxy-carbon atoms at the narrow rim of one calix [4]arene monomer. The reference group is fixed in space and the pulled group is moved along the vector connecting the two groups using a harmonic potential. The end-to-end distance r ee is defined as the distance between the two groups. The aliphatic loops are omitted for clarity.

Coarse-grained potentials
The CG potentials for the solvent molecules were computed using the iterative Boltzmann inversion (IBI) method [39]. IBI relates the free-energy of a pair of particles to the logarithm of the radial distribution function (RDF) to obtain a potential of mean force (PMF) as a function of the distance r between the particles [39,40]. Starting from the expression for the PMF, U(r) = −k B T ln(g(r)), with g(r) denoting the RDF and k B the Boltzmann constant, one obtains the effective pair potential iteratively. The initial PMF is estimated from the reference RDF, U(r) CG 0 = −k B T ln(g(r) ref ) and the iteration cycle is defined by The mesitylene molecules were treated as spheres and the RDF of the center of mass was used in the determination of the PMF. In the iterative scheme we used a cut-off for the potential of 1.2 nm and each simulation run had a duration of 200 ps (the first 20 ps were omitted for equilibration). The effective pair potential was smoothed using cubic splines and we iterated the procedure 325 times without pressure correction. The simulations were performed using the program package VOTCA [41,42] (VOTCA 1.3) and the reference RDF was obtained from an AA simulation (25 ns) of a pure mesitylene system in a box with size (7.49 nm) 3 . The results for the RDF are presented in figure 2. Without showing the results here, we note that the resulting PMF does not exhibit an observable minimum, cf the supporting information (https://stacks.iop.org/JPCM/33/194005/mmedia). It is well known that coarse graining speeds up the dynamics relative to AA simulations [43,44]. We quantified this effect by a measurement of the diffusion coefficients which are given by D CG = 11.8 × 10 −10 m 2 s −1 and D AA = 1.7 × 10 −10 m 2 s −1 . Note that these values differ from the corresponding ones given in reference [22] (D CG = 8.8 × 10 −10 m 2 s −1 and D AA = 7.2 × 10 −10 m 2 s −1 ). We attribute this to the different box sizes, different treatment of the long range electrostatic interactions and different integrators.

Adaptive resolution scheme
For the AdResS methodology it is important to balance the chemical potential in the AA region and the CG region. This is accomplished by the computation of the thermodynamic force that is calculated iteratively according to and works on the CG part in the hybrid region allowing the exchange of molecules between the different regions [25,45]. Here, M is the mass of the molecules, κ T is a constant conceptionally related to the isothermal compressibility and ρ i (x) the density in the ith iteration step. The iteration starts with the initial density and a vanishing thermodynamic force. We have used a spherical setting with different radii of the AA region and a constant slab thickness of the hybrid region of s Hy = 1.2 nm. The range of the thermodynamic force was enlarged by 0.2 nm in order to incorporate molecules that are only partly placed in the hybrid region. We varied the radius of the AA region, r AA , in order to study the dependence of the results on this choice. The values used are r AA = 1.6, 0.8 and 0.4 nm. Since the average end-to-end distance in equilibrium is about 1.4 nm, this means that for the smallest value of r AA parts of the calix [4]arene dimer are located in the hybrid region. The center of the calix [4]arene dimer was always kept in the center of the simulation box ensuring that the distance between the pulled group and the reference group to the border of the CG region stayed the same. We used a force capping methodology in the hybrid region [46] with a maximum force of F max = 2.5 × 10 5 kJ (mol −1 K −1 ). This ensures that the forces between two particles at the border between the CG region and the hybrid region do not increase too much. Force capping is applied because the CG potential depends only on the distance of the center of mass of the mesitylene molecules and is independent of the relative orientation of the molecules. Therefore, the distance between parts of the molecules can become very small. All simulations had a duration of 1 ns with the first 400 ps omitted due to equilibration. We performed 88 iterations for r AA = 1.6 nm, 149 for r AA = 0.8 nm and 244 for r AA = 0.4 nm in order to obtain flat density profiles throughout the simulation box. (For the smaller AA regions a flat density profile in the hybrid region is more important.) In figure 3, we plot the density as a function of the distance from the center of the simulation box which coincides with the center of mass of the calix [4]arene dimer. Since the calix [4]arene dimer resides in the center of the box, the solvent density vanishes for small values of r center . For larger distances, the density follows the one of the AA simulation to a very good approximation. Due to the structure of the calix [4]arene dimer the distribution of the solvent density in the immediate neighborhood of the solute is not isotropic. This is the reason for the appearence of the hump like structure for r center ∼ 1 nm. For some further information regarding the thermodynamic force we refer to the supporting information.

Equilibrium simulations
The AdResS is well known to reproduce the AA results in equilibrium simulations for a number of different situations [25]. In order to assure that it also works for our particular system, we performed equilibrium simulations and monitored the most important structural features. In figure 4(a) we show the end-to-end distance r ee as a function of the simulation time and its probability distribution for a 50 ns AA simulation. We additionally present the distributions for AdResS simulations using various values for the radius of the AA region, r AA , in figure 4(b). It can be seen that the average value of r ee is almost independent of the size of the AA region. Only for the smallest value of r AA = 0.4 nm, it differs from the AA value by about 1.5%. We mention that for this small AA region parts of the calixarene capsules are located in the hybrid region. However, since the center of the dimer is placed in the AA region, the stabilizing UU-bonds are described correctly. From the similarity of the widths of the distributions we conclude that also the fluctuations are very well sampled by the AdResS simulations using the given parameters. We have also monitored the number of UU-bonds stabilizing the closed structure. In all simulations one observes that most of the time (more than 50%) the maximum number of 16 UU-bonds are formed and there are quite frequent fluctuations in which one or two bonds open. However, there are hardly any significant differences between the AA simulation and the AdResS simulations. We thus conclude that the AdResS simulations give a good representation of the AA results.

Simulation setup
All FPMD simulations presented in this work were performed using the FR protocol. We used two modes, a pull mode and a relax mode where after pulling the dimer into the open conformation the pulling direction is inverted and all other parameters remain the same. The reference group was fixed in space and a harmonic potential was applied to the pulled group, where the groups are defined as in figure 1. The force measured at the spring is given by Here, K is the spring constant, V the pulling velocity, and z(t) denotes the deviation of the position of the pulled group from its initial value. Note that the extension is defined as x = V · t and the so-called loading rate is given by the product of the force constant and the pulling velocity, μ = K · V.
In FPMD simulations the system is driven out of equilibrium and it is pulled through the solvent. In our earlier simulations employing virtual sites [22], the coarse graining of the solvent did not have a dramatic impact on the equilibrium properties of the system, but the rupture forces were significantly reduced relative to those obtained from AA simulations. Therefore, in the present study, we investigate the dependence of the results of AdResS simulations on the size of the AA region, the pulling velocity and the stiffness of the pulling device.
As mentioned above, we used a box of size of (7.49 nm) 3 for all AA simulations and the AdResS simulations although we have found in earlier studies of the calix[4]arene catenane dimer system that a box length of 5.8 nm is sufficient for all FPMD simulations performed so far [37,38]. A box of this size appears to be a good compromise between the two extreme scenarios that have to be considered in the case of FPMD simulations. For quasistatic pulling, r ee approximately follows the pulling protocol r ee x max while for fast pulling one expects r ee x max where x max denotes the extension reached at the end of the pulling simulation. For a typical value of x max ∼ 4 nm a box of size (5.8 nm) 3 is large enough. Such a box contains one dimer and 700 mesitylene molecules with 21 atoms per molecule. In order to estimate the relevant number of atoms in an AdResS simulation, we compute the average number of mesitylene molecules in the respective regions by assuming a homogeneous density for simplicity. For instance, using r AA = 1.6 nm, we have about 90 mesitylene molecules in the AA region, 380 in the hybrid region and 1315 in the CG region. As a result of this naive estimate, we have to treat explicitly only approximately half of the number of solvent particles as compared to the AA simulations. Using noncubic boxes and nonspherical AA regions in the AdResS simulations this ratio can even further be reduced. Furthermore, adaptive schemes to determine the size of the AA region are also expected to be very effective. However, a direct estimate of the computational efficiency of the AdResS simulations as compared to the AA FPMD simulations is not possible with the present preliminary implementation of AdResS.
For a brief presentation of the results of AA simulations, in figure 5 we show examples of the most important observables as a function of the extension, x = V · t, for a pulling simulation with K = 1 N m −1 and V = 1 m s −1 (i.e. μ = 1 N s −1 ). Shown are the end-to-end distance r ee (red), the force measured at the spring attached to the pulled group, F, (black) and additionally the number of UU-bonds (blue) and UE-bonds (green). The transition from the closed state to the open state that takes place at an extension of roughly 2.6 nm is observable in all quantities. The end-to-end distance increases almost linearly from the equilibrium value of 1.4 nm to about 1.55 nm and jumps at the transition to 2.1 nm. The force also increases linearly until there is a rip in the force versus extension curve (FEC) after which it increases again. If the molecular energy landscape is assumed to be harmonic, one can extract the corresponding stiffness from the slope of the FECs [29,37,47]. It is also evident that the number of UU-bonds (cf figure 1) slowly decays for small extensions and then abruptly drops to zero at the transition point. At this point, the UE-bonds stabilizing the open state are formed. Due to the non-equilibrium nature of the pulling procedure the maximum number of 8 UE-bonds is not reached. We mention that the behavior in the relax mode simulations is quite similar for the chosen parameters albeit with a finite hysteresis. For very fast pulling the transition becomes irreversible in the sense that in the relax mode the closed state is no longer reached [28].

Force versus extension curves
In figure 6 we show examples of FECs as obtained from AA simulations and from AdResS simulations with r AA = 1.6 nm for K = 1 N m −1 and V = 1 m s −1 . The simulations were always performed in the same way. First a pulling simulation was performed until the extension V · t reached a value of 4 nm and then a relax mode simulation followed. Therefore, for large extensions the two curves are on top of each other. The hysteresis demonstrating the non-equilibrium nature of the FPMD simulations is manifested by the different extensions at which the respective transitions take place. Only for quasistatic pulling, i.e. very small pulling velocities, the results of the pull mode and the relax mode simulations coincide. It is evident that the FECs for the AA simulations and the AdResS simulations are very similar for the parameters used and the chosen realization. In order to test the applicability of the AdResS methodology when applied to FPMD simulations in more detail, we used two values for the loading rate, μ = 1 N s −1 and μ = 10 N s −1 . These are quite high values for μ and therefore the system is driven strongly out of equilibrium. Furthermore, we used different values of the pulling parameters as given in table 1. As the conformational transition from the closed state to the open state and vice versa are stochastic processes, the rupture forces and the rejoin forces will vary in different realizations of the FPMD simulation using the same parameters. Therefore, in order to provide a meaningful statistical analysis of the results, we performed 300 simulations for each set of parameters.  We first consider the FECs as they have been shown for one example in figure 6. In figure 7 we show averaged FECs [47,48] obtained from pull mode simulations. The shaded areas are meant to represent the width of the distributions (the second moment). It is evident that there is a very good agreement between the results of the AA simulations and of the AdResS simulations for r AA = 1.6 nm. For r AA = 0.8 nm, the system appears to be softer in the sense that the rupture force is smaller. The slopes of the averaged FECs, however, are the same meaning that the molecular stiffness is unaltered. This changes for the smallest AA region used, r AA = 0.4 nm. Here, the rupture event appears at still smaller force and has more the form of a broad shoulder than of a rip. However, the most prominent difference to all other simulations is given by the fact that the slope in the open state is smaller than in the other cases and furthermore for large extensions, x 3.5 nm" the force decreases again indicating that the system becomes unstable. This can be understood from the fact that the end-toend distance exceeds 2.2 nm for these extensions, cf figure 5. Thus, a substantial part of the calix [4]arene dimer enters the hybrid region and the aliphatic loops are destabilized because the relevant interactions are not considered in the AdResS protocol. Due to this failure of AdResS for such a small value of r AA when applied to FPMD simulations we will no longer consider simulations performed using r AA = 0.4 nm. In figure 8 we show averaged FECs in the pull mode and the relax mode for various choices of the parameters used in the protocol. The general features are well reproduced qualitatively by all AdResS simulations independent of the pulling velocity and the loading rate. The hysteresis shows the known dependence on V and also an increase of the rupture force with the loading rate is observed [37,47]. In particular, the results of the AA simulations and of the AdResS simulations for r AA = 1.6 nm coincide very well for all sets of parameters in both modes. The finite size of the AA region becomes important for r AA = 0.8 nm. In this case, the mean rupture forces are smaller and the mean rejoin forces are somewhat larger than in case of the AA simulations. In case of the pull mode simulations the differences of the mean rupture forces exceed the widths of the distributions. Additionally, the differences apparently are more pronounced for the larger loading rate.

Rupture force distributions
In order to discuss these findings in more detail, we consider the distributions of the characteristic forces obtained from the individual FECs as shown in figure 6. In the present paper we will consider the mean force between the maxima and the minima in the transition region for both cases, pull and relax mode FECs, cf figure 6. These forces behave somewhat different from the maximum rupture force and the minimum rejoin force but capture all important features [37,38]. Furthermore, and more important in the present context, for reversibly binding systems like the calix [4]arene dimer investigated here, F rupt and F rejoin will converge to the equilibrium value, F eq , for slow pulling [47]. Figure 9 shows these distributions for the four values of the force constants and the concomitant pulling velocities as listed in table 1. The mean rupture forces (vertical bars in figure 9) slightly decrease as a function of the force constant and the mean rejoin forces slightly increase. This indicates that they will merge for vanishing pulling velocity.
As pointed out above, the AdResS simulations using r AA = 1.6 nm yield results almost identical to the AA simulations. In case of the AdReS simulations with r AA = 0.8 nm the distributions of rupture forces are shifted to smaller forces while the F rejoin distributions tend to be shifted to larger forces. Thus, in these AdResS simulations the hysteresis (cf figure 8) is less pronounced and this apparently has its origin in the impact of the enhanced solvent mobility in the CG region. We mention that the number of simulations performed (300 for each set of parameters) is not large enough for a detailed discussion of the shape and the widths of the distributions. Furthermore, the observed differences appear to depend only weakly on the force constant and more strongly on the loading rate. Only for the slowest pulling velocity of V = 0.125 m s −1 the difference in the rupture forces is somewhat smaller than for faster pulling. At the same time the F rejoin distribution becomes rather broad. We interpret this finding as resulting from the fact that for very slow pulling, in the relax mode simulation the system resides in the open state for quite a long time and r ee decreases only slowly. Therefore, the higher solvent mobility in the CG region (and also in the hybrid region) apparently gives rise to some further softening of the dimer.
The most prominent features that can be observed in the force distributions are also reflected in the mean forces. In figure 10, we present the mean rupture and the mean rejoin forces as a function of the pulling velocity. An increase of the rupture force and a slight decrease of the rejoin force is observed for both loading rates. The differences of the results for r AA = 0.8 nm to the other simulation results discussed above are evident and they are somewhat more pronounced for μ = 10 N s −1 than for μ = 1 N s −1 . As mentioned above, the results regarding the differences appear to depend mainly on the loading rate. The average difference in the values of F rupt is on the order of 400 pN for μ = 1 N s −1 and it is about 650 pN for μ = 10 N s −1 . The differences in the mean rejoin forces is almost the same for both loading rates (approximately 150 pN).
The decrease of the rupture force for the AdResS simulations with r AA = 0.8 nm indicates that the increased solvent mobility (D CG /D AA 7) is important. If we assume that the dynamics of the calix [4]arene dimer is dominated by activated barrier crossing, the celebrated Bell model [49] yields a logarithmic dependence on the diffusion coefficient, F rupt ∼ log (μ/D) [50]. However, the ratio between the rupture forces should not vary with loading rate. The larger differences for μ = 10 N s −1 possibly are due to the fact that the regime of drift motion might be entered for fast pulling [37]. In this regime the rupture force behaves as F rupt (V → ∞) ∼ (μ/D) [50,51] and therefore depends somewhat stronger on μ than in case of activated dynamics.
The situation apparently is somewhat different if the relax mode simulations are considered. The differences in the mean rejoin forces are almost independent of V and of μ. The increased solvent mobility experienced by the dimer in the open state seems to have a smaller impact on the rejoin forces.

Analysis of the H-bond network
We now discuss some aspects of the dynamics of the H-bond network and some properties of the energy landscape. As mentioned above, the network of UU-bonds stabilizes the closed state and these bonds open at the transition into the open state where the network of UE-bonds is formed, cf figure 5. Since the transition is a stochastic process and takes place at different extensions depending on the pulling velocity, the number of H-bonds also strongly depends on the extension. On the other hand, the different conformational states are characterized by the corresponding ranges of the end-to-end distance. Therefore, in the following we consider the mean number of H-bonds as a function of the average end-to-end distance of the dimer, r ee , and present the results in figure 11. This representation has been shown to depend only weakly on the pulling protocol [29,37]. In the pull mode, the average number of UU-bonds, #UU , decreases from almost 16, the maximum possible number, to zero. The average number of UE-bonds, #UE , reaches a maximum at r ee ∼ 2 nm and for further stretching they decrease again. A similar behavior is observed in the relax mode simulations with the only difference that the maximum number of UE-bonds formed is somewhat larger than in the pull mode. In all cases, the AdResS simulations using r AA = 1.6 nm yield results that are in excellent agreement with those obtained from AA simulations. For simulations employing the smaller AA region, r AA = 0.8 nm, #UU behaves very similar to the other results while #UE exceeds the values obtained from the other simulations in the region of the maximum indicating a slightly higher stability of the open state in this case as it would be expected for a somewhat smaller loading rate. However, these differences still are very small and within the width of the distributions. It is interesting to note that #UU is the same as for the AA simulations also in the relax mode simulations and the enhanced stability of the open state plays almost no role for the structure of the closed state.
As in our earlier work on the unfolding kinetics of the calix [4]arene dimer, we use the average number of H-bonds for a characterization of the energy landscape of the system. In particular, we define the minima of the closed state (q C ), the open state (q O ) and the transition state (q T ) in the following way [37]: and present the results in figure 12. For both loading rates, the results are independent of the force constant and coincide within the statistical error. We conclude from these findings that the energy landscape of the dimer is hardly altered in the AdResS simulations.

Conclusions
We have presented a detailed investigation of the applicability of the AdResS methodology to FPMD simulations. As an example we chose a well studied system, a calix [4]arene dimer, that is known to undergo conformational transitions that can be described by a two-state model. Additionally, the system is small enough to allow for the direct comparison with AA simulations. We studied the system for different sizes of the AA region and a fixed slab thickness of the hybrid region. In equilibrium, the AdResS is known to work with excellent results when compared to AA simulations and we confirmed this finding with the simulations of our system consisting of one dimer in mesitylene solvent.
Since AdResS has been developed for equilibrium simulations, it is not clear to which extent it can be applied to nonequilibrium situations. We therefore studied the performance of AdResS as applied to FPMD simulations for high pulling velocities and large loading rates. We found that all AdResS simulations using an AA region with a radius of r AA = 1.6 nm yield results that are basically identical to the results obtained with AA simulations. This holds for all quantities that characterize the kinetics of the conformational transition in our model system and also the structural features like the number of the H-bonds formed in each of the states. It is important to point out that a box size comparable to the AA region in the AdResS simulations is much too small to give reliable results in AA simulations, in particular in the open state of the dimer.
The results of AdResS simulations employing a smaller AA region (r AA = 0.8 nm) worked very well in equilibrium but gave rise to deviations from the AA results in FPMD simulations. This holds in particular when the rupture forces and rejoin forces are considered. Here, the system appears softer and the distance to equilibrium for a given loading rate seems to be reduced. We attribute this to the higher solvent mobility in the hybrid region and in the CG region. The observed differences apparently increase with increasing loading rate meaning that for larger μ = K · V a larger AA region is required and the results hardly depend on the stiffness K or the pulling velocity V separately. When structural features like the number of H-bonds stabilizing the conformational states of the dimer are considered, the AdResS simulations gave a very good representation of the results of the AA simulations. The same holds for the characteristic features of the energy landscape as deduced from the properties of the H-bond network.
The model system considered in the present work is very simple in the sense that the (un)folding pathway is determined by two states only. Furthermore, the solvent is of an aprotic nature and therefore does not interfere with the H-bonds stabilizing the conformations of the dimer. However, the solvent mobility has proven to play a crucial role even in this well defined situation. Therefore, the size of the AA region in applications of the AdResS methodology to FPMD simulations has to be chosen with care. This holds in particular if polar solvents are considered since in that case electrostatic interactions become of major importance. A study of this effect is under way.
In conclusion we have demonstrated that the AdResS methodology can be applied successfully to perform FPMD simulations. The number of particles treated in an AA manner can be reduced considerably but the current preliminary implementation of AdResS for FPMD simulations does not allow for a quantitative estimate of the computational efficiency of the methodology.