The four actomyosin complexes in different nucleotide states (nucleotide-free, ATP, ADP·Pi and ADP) were minimized, solvated and equilibrated at a temperature resembling the experimental one. After about 600 ps of equilibration dynamics, the root mean square deviation (RMSD) of the Cα atom position with respect to the starting structures reached the stable value of 0.57±0.03 nm (mean±SD), 0.26±0.02 nm, 0.83±0.08 and 0.24±0.02 nm, respectively in the case of nucleotide-free, ATP, ADP·Pi and ADP complexes.

Table 1 provides the conformational details of the complexes as calculated by using the Procheck software [29]; the individual stereochemical quality of the structures is meant as a measure of the reliability of the equilibrated models. The fully, additionally and generously allowed regions indicate the permitted backbone conformations for psi/phi angle pairs for each residue of the protein. The high percentage of residues belonging to the allowed regions supports the good reliability of the obtained results.

In the literature, experimental investigations of the interaction forces between actin and myosin at the single molecule level are usually done using two different types of set-ups. The former consists in the actin filament attached to a probe held by a force transducer and presented to myosin molecules bound to a substrate [30–32]. The latter assumes that actin is bound to a surface, while the motor protein, attached to plastic or glass beads held by force transducers is presented to the counterpart filament [33]. The most of the experimental apparatuses frequently used as force transducers are atomic force microscopes (AFM), and optical tweezers; the proteins’ interactions are estimated using different types of photodiode detectors.

In the present work we used a MD approach to investigate the characteristics of the forced unbinding of myosin motor domain from an actin monomer. This study was motivated by experimental measurements of the unbinding process [33]; both in single-molecule experiments and in computer simulations, the proteins’ unbinding is obtained by applying external forces that are large enough to overcome the potential energy barriers. A correct estimation of the interaction properties of the actomyosin complex depends on the choice of the appropriate parameters for the SMD simulations. Therefore, a preliminary sensitivity analysis was performed to define the optimal pulling parameters to be used for the simulation of the actomyosin complex unbinding.

The equilibrated structure of the actomyosin complex in ADP·Pi state was used to evaluate the influence of the pulling rates (ν) on the molecular behaviour of the system during the forced unbinding simulation. Three pulling rates – 0.1 nm/ps, 0.01 nm/ps and 0.001 nm/ps – were tested using a spring constant of 3000 kJ mol⁻¹ nm⁻². The obtained interaction energy vs. intermolecular distance profiles (Figure 2a) are qualitatively similar in the three cases; the minimum interaction energy decreases when the value of ν is lowered. Moreover, at higher pulling rates (0.1 and 0.01 nm/ps), a plateau of minimum interaction energy was reached when the two proteins were approached; at a pulling rate of 0.001 nm/ps, a minimum value of the interaction energy was obtained at an intermolecular distance of about 6 nm.

In agreement with the results of the experimental unbinding forces of different biological complexes [34, 35], our results (Figure 2b) suggest that the lower the pulling velocities are, the smaller the forces needed to move the myosin apart from the actin become; decreasing the velocity from 0.1 to 0.001 nm/ps reduces the value of the maximum pulling force from 7492 pN to 1023 pN. The force required to unbind biological complexes strongly depends on the pulling velocity. Merkel et al. [34] studied both avidin-biotin and streptavidin-biotin unbinding behaviour by means of dynamic force microscopy. Different loading rates were tested over six orders of magnitude (from 0.05 to 60000 pN/s) and the unbinding force increased from 5 to 170 pN. In addition, atomic force microscopy experiments on P-selectin revealed that the force required to unbind the ligand-receptor complex has a logarithmic dependence on the pulling velocity [35].

Even if the lowest maximum pulling force (1023 pN) obtained by performing SMD simulation with a pulling rate of 0.001 nm/ps is still one order of magnitude higher than the experimental one [34], it is weak enough and changes very slowly in time to induce minor distortions in the molecular structures during pulling experiment simulation.

Generally, SMD simulations meant to induce conformational changes or proteins’ unfolding are very sensitive to the pulling rate [27]; when high pulling rates are being used, unphysical trajectories are yielded. Unlike unfolding SMD simulations, the set-up adopted for the characterization of the actomyosin complex allows a proper estimation of proteins’ behaviour during unbinding as long as proteins distortions are limited [27].

The pulling velocity used in the SMD simulations is inversely related to the unbinding time; at low velocities the molecular system has more time to thermally equilibrate and consequently oscillations were recorded due to the continuous rearrangements of the proteins (Figure 2b); this phenomenon resembles the biological behaviour of the proteins.

Moreover, since we deal with globular-shaped proteins, the large size of the solvated models (about 127000 atoms) imposes a pulling rate of 0.001 nm/ps, comparable with the one adopted in previous computational studies on different complexes, which successfully matched the experimental data [36, 37] but still higher than the pulling rates used regularly in AFM experiments [38], which are in the nanometres per milliseconds range. Indeed, considering that MD simulations were performed in parallelization conditions (Intel Core2 CPU 6300 1.86 GHz), the CPU time was 0.2 ns per day. Decreasing the pulling rate of one order of magnitude (from 0.001 nm/ps to 0.0001 nm/ps) would increase the total CPU time required to evaluate the interaction properties of each state of the actomyosin complex, from 25 days up to more than 250 days due to the high oscillations of the system.

Taking into account the abovementioned observations, the lowest tested pulling velocity (0.001 nm/ps) was used to perform the SMD simulations for the analysis of the different nucleotide states of the actomyosin complex.

The effect of the spring constant on the mechanical behaviour of myosin when it is moved apart from the actin monomer was evaluated through SMD simulations at a constant rate of 0.001 nm/ps with three different k values, 300, 3000 and 10000 kJ mol⁻¹ nm⁻².

When a stiff spring is used to pull apart the myosin from the actin monomer, large amplitude fluctuations of the pulling force profile can be observed, and some filtering might be considered in order to reveal the local events. Softer springs exert a less restricted action on the molecular system; in this case local displacements are detected (Figure 3a). Our results in terms of force oscillations agree with the observations of Lorenzo and Bisch [39]. When thermodynamic systems are subjected to external oscillatory forces, the amplitude of the displacement (δ_x) due to the thermal energy fluctuation can be estimated by means of the expression: δ_x= (k_BT/k)^1/2, where k_B is the constant of Boltzmann and T is the absolute temperature (T = 300K). Since the amplitude of the force oscillations is δF = -k δ_x, then the following relationship is obtained: δF = -(k_BTk)^1/2 [39]. The results of the SMD simulations performed with different k (Figure 3a) comply with these theoretical models. For instance, in case of k=3000 kJ mol⁻¹ nm⁻² (4.98 N/m), the mean value of the amplitude of the force oscillations measured during SMD simulation is about 180 pN (Figure 3a), while the estimated value obtained by theoretical models is equal to 145 pN.

The simulation time required to pull the myosin motor domain 2 nm apart from actin was 4.8 ns when using k =300 kJ mol⁻¹ nm⁻², while in case of higher spring constants (k = 3000 kJ mol⁻¹ nm⁻² and k = 10000 kJ mol⁻¹ nm⁻²) the simulation time decreased of about 27% (Figure 3b).

The best compromise in terms of limiting the amplitude of oscillations introduced in the force values and reducing the simulation time, was obtained for a spring constant of 3000 kJ mol⁻¹ nm⁻².

Using the optimal pulling parameters defined previously (ν = 0.001 nm/ps, k = 3000 kJ mol⁻¹ nm⁻²), the interaction properties of the actomyosin complex in different nucleotide-dependent conformations were estimated via SMD simulations, starting from the equilibrated model of each state. We monitored the intermolecular distance and the interaction energy between the two proteins of the actomyosin complex and we evaluated the interaction forces and binding stiffness (Figure 4 and Table 2) by means of Eq. 3 and Eq. 4 (see Experimental Section).

Our results predict that the conformational changes occurring in the motor domain of the myosin during a cycle of ATP hydrolysis markedly influence the interaction forces and the energy profiles of the actomyosin complex.

The interaction force vs. intermolecular distance curves (Figure 4) are approximations, valid only around the minimum value of the forces, but are not representative for systems behaviour in general. Nevertheless the maximum interaction forces estimated for each nucleotide-dependent actomyosin complex (F_max, Table 2) are a quantitative estimation, which reflect the coupling behaviour shown by the interaction energy profiles in Figure 4 and by the non-covalent interactions (van der Waals and Coulomb energy terms, Figure 5) acting between the actin monomer and the myosin in the complex. The rigor conformation stands out for a depth falling of the interaction energy curve when myosin motor domain approaches the actin monomer, while the energy curves calculated for ATP, ADP·Pi and ADP states have shallow minima when the proteins are moved toward. This behaviour is representative of a weak coupling between two proteins and enables their relative motility.

In all the studied cases the non-bonded energy term lowers with the decrease of the intermolecular distance until the minimum interaction energy is achieved. As shown in Figure 5, at the equilibrium distance (intermolecular distance corresponding to the energy minimum) the contribution of the electrostatic attraction to the net interaction energy is dominant when the actomyosin complex is in the nucleotide-free and ADP states (62% and 65%, respectively), in agreement with Liu et al. [26]; on the other hand, in the ATP and ADP·Pi states the van der Waals term dominates (75% and 55%, respectively). Consequently, the Pi release from the complex is coupled with a decrease of the van der Waals contribution (from 271 to 147 kJ/mol) and a growth of electrostatic contribution (from 219 to 276 kJ/mol) to the total interaction energy term.

In addition, the number of the inter-chain hydrogen bonds (Hb) vs. intermolecular distance (Figure 5) is in agreement with the affinities between the proteins during a cycle of ATP hydrolysis (Figure 4); in fact complexes characterized by higher interaction forces are connected by a higher number of Hb. A detailed analysis of Hb connecting the interfaces of the equilibrated actomyosin complexes in different nucleotide-dependent states gives also valuable qualitative indications about the affinities at the atomic level within the complex during an ATP hydrolysis cycle. The hydrogen bridging of myosin (M) and actin (A) occurring when the nucleotide-binding cleft of the motor protein is open, is mainly due to the following Hb: Lys448 (M) - Glu93 (A), Lys550 (M) - Met47 (A), Tyr552 (M) - Glu93 (A), Asp553 (M) - Lys50 (A). When the nucleotide-binding cleft is closed, the proteins’ interfaces are mainly bridged by 5 Hb: Lys530 (M) - Met45 (A), Lys530 (M) - Glu58 (A), Arg546 (M) - Glu100 (A), Phe547 (M) - Arg96 (A), Lys549 (M) - Tyr92 (A). In both cases the hydrogen bridging of the actomyosin complex involves the residues of the motor protein’s primary binding site. The low number of Hb connecting the complex indicates that the electrostatic and van der Waals terms are primarily responsible for the actomyosin interaction.

Our findings are in agreement with the experimental observations by Rayment et al. [20]. The in-vitro motility assay, combined with X-ray crystallography and molecular biology techniques, shows that myosin in the rigor state (nucleotide-free) is usually locked tightly to the actin filament (interaction force F_max=0.83 nN, Table 2; electrostatic attraction is dominant, Figure 5), while the affinity between the two proteins decreases when the ATP nucleotide is bound to the large cleft of the myosin head (F_max=0.27 nN, Table 2; van der Waals contribution is dominant, Figure 5); after ATP hydrolysis (F_max=0.46 nN, Table 2; van der Waals contribution is dominant, Figure 5), the release of the Pi from the complex triggers the power stroke and generates the myosin(ADP)-actin state (F_max=0.60 nN, Table 2; electrostatic attraction is dominant, Figure 5), which is the intermediate complex that leads to the rigor state after the ADP release [20].

The interaction strength modulation is supported by the conformational changes of the myosin motor domain. In strongly bound states, the actin-binding cleft is closed and 50K domain of the myosin interact with actin monomer, while in weakly bound states the actin-binding cleft is open and only the lower 50K domain establish weak stereospecific interactions with the actin [40]. These conformational changes of the myosin determine large variations of the interaction surface of the actomyosin complex. The interaction surface of both myosin motor domain and actin monomer varies with the type of the adenine nucleotide attached to the nucleotide binding pocket of the motor protein. As suggested by Holmes et al. [40], our results (Table 2) indicate that the interaction surface (S) of the actomyosin complex halves when passing from rigor state (833.32 Ų) to weakly bound state (402.08 Ų). Concerning the nucleotide-free complex, the interaction surface is comparable to the one estimated by Liu et al. [26] when performing MD simulations on the complex with two actin monomers. This similarity gives indications about the reduced contacts between myosin and the adjacent actin monomer.

With regard to experimental data from literature to be used for comparison, most of the investigations on the single actomyosin complex are focused on the displacement and the force generated during a cycle of ATP hydrolysis [32,41]. Other experiments performed at single molecule level give indication about the value of the unbinding force [30–33]. Nakajima et al. [33] monitored the interaction events occurring on single myosin head and actin by means of scanning force microscopy and measured an unbinding force of 14.8 pN [33], in good agreement with the value obtained by Nishizaka et al. [30, 31] via optical tweezers. These values are lower than the ones obtained in the present study by means of SMD. This discrepancy is possibly due to the differences in the set-up adopted in the MD approach and in the experimental procedure. In SMD simulations the myosin motor domain was pulled through an unbinding direction perpendicular to the longitudinal axes of the actin filament while in experimental set-ups the direction of the imposed loads varied within a wide range. The peeling procedure used in Nakashima’s and Nishizakas’ experiments [30, 31, 33] likely cause a distortion of the actomyosin interface altering their conformational match. Small deformations of proteins’ structures could reduce significantly the binding interface and consequently the binding strength.

Moreover, as shown previously, the pulling force and consequently the interaction force highly depend on the pulling velocity. Since the unbinding force is smaller when the external load is applied slowly, we may expect a slight further decrease of the interaction force when performing SMD simulations with pulling rates lower than 0.001 nm/ps.

A significant difference between MD simulations [42] and measurements by AFM experiments [43] was observed for the avidin-biotin molecular system. As pointed out by Izrailev et al. [42], the different values of the interaction forces arise from the different time scales of the simulated and experimental tests. In the AFM experiment performed by Florin et al. [43] the unbinding time was ~1 ms and the measured force was about 160 pN, while in the case of MD [42] the simulation time was about 1 ns and the obtained unbinding forces were significantly larger (450 pN). In the case of AFM and similar manipulation techniques, the unbinding time is usually of the order of a millisecond, long enough to allow the proteins to thermally equilibrate and to overcome thermal barriers. In the case of molecular systems of the same size of the one considered in this study, MD simulations can cover only a timescale of the order of nanoseconds; hence, the unbinding event occurs in a dissipative regime by applying large forces, which are able to overcome relevant energy barriers and induce a sufficiently rapid dissociation [42].

The X-ray crystallographic models of both open and closed forms of the myosin motor domain and of the actin monomer were obtained from the Brookhaven Protein Data Bank. Four different nucleotide-dependent binding states (nucleotide-free, ATP, ADP·Pi and ADP) of the actomyosin complex were built. With regard to myosin, we used Dictyostelium discoideum data; the nucleotide-free state (when the motor protein is tightly locked to actin filament) was modelled starting from 1FMV.pdb [21], while the weak binding states (ATP, ADP) were obtained starting from 1FMW.pdb [21] and 1MND.pdb [2], respectively. The myosin(ADP·Pi) complex was built starting from the equilibrated model of 1FMW.pdb. In order to model actin, we used 1ATN.pdb [18], which is an atomic model of the complex between rabbit skeletal muscle actin and bovine pancreatic DNase I. The myosin models were docked on actin monomer using 1O1B.pdb [19] as template configuration. In particular, the Least Squares Fitting (LSQ-fitting) algorithm was used to set the relative positions of the proteins and to properly fit the interfaces.

Each complex was prepared for the MD simulations by removing all the water molecules embedded in the structures and all the chemical compounds (chloride, magnesium and tetrafluoroaluminate ions) added during the crystallization procedure for the stabilization of the proteins and not physiologically present in biological systems.

Energy minimizations and MD calculations were performed by means of the Gromacs 3.3 software, using the GROMOS96 43a1 force field. All simulations were carried out in parallelization on an Intel(R) Core(TM)2 CPU 6300 1.86 GHz processor-based Linux running machine.

To ensure that the investigated actomyosin complexes were free from strongly repulsive non-bonded contacts or geometric distortions inconsistent with the potential energy function, a preliminary in vacuum energy minimization of each system was performed. The energy minimization was carried out using a sequence of two algorithms (steepest descent and limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-bfgs).

The actomyosin models were then solvated in rectangular boxes (17 nm x 6 nm x 13 nm) with about 38460 Single Point Charge (SPC) water molecules, bringing the total size of the systems to about 126764 atoms each. In order to reduce the computational costs associated with water molecules calculations, the box dimensions were chosen as small as possible to get a water shell surrounding the protein complex with a thickness of at least 1 nm (except for x direction, Figure 1). The box dimension in the x direction was instead set large enough to perform SMD simulations. The total charge of the solvated protein complexes was neutralised by adding Na⁺ counter ions. Each system was then prepared for the MD steps by performing a further energy minimization.

In order to gradually increase the system temperature from 0 to 300K, proteins and solvent were separately coupled to a Berendsen thermostat for 100 ps, and a time constant of 0.1 ps was used to produce a weak coupling between temperature and atom velocities [44]. To finally achieve the starting models for the SMD simulations, each system was further equilibrated for 1 ns.

MD simulations were performed using the leap-frog integration with a time-step of 2 fs and the atomic coordinates of the whole system were recorded every ps. The van der Waals interactions were cut off beyond 1 nm. The Particle-Mesh-Ewald (PME) method [45] was used to allow for the calculation of the long-range electrostatic interactions using a Fourier spacing of 0.12 nm and a fourth order interpolation.

The investigation of the interaction forces that characterise the actomyosin complex in different nucleotide-dependent states was carried out by performing SMD simulations that started from the equilibrated model of each complex. During the simulation the actin monomer was fixed by blocking the spatial coordinates of the Cα atoms, while the myosin motor domain was moved closer and further by means of a virtual spring. One end of the spring (R) was attached to the center of mass of the myosin molecule (Figure 1), while the other end (S) was moved along the x direction at a constant rate v defined as: (1)ν=dxS(t)dtwhere x_S (t) represents the coordinate in x direction of the free end of the spring as a function of time, t.

The pulling force acting on the myosin motor domain was calculated as: (2)FP(t)=k(xS(t)-xR(t))where x_R (t) is the coordinate of the end of the spring attached to the myosin center of mass (R), and k is the value of the spring constant.

In order to properly set the parameters of the SMD, the final equilibrated structure of the actomyosin complex in the ADP·Pi state was used to perform an analysis of sensitivity to the pulling rate (v) and to the spring elastic constant (k). Concerning the pulling rate tests, the k value was set at 3000 kJ mol⁻¹ nm⁻², while the myosin was moved with velocities of 0.1, 0.01 and 0.001 nm/ps. Three other SMD simulations with a v of 0.001 nm/ps and a k of 300, 3000 and 10000 kJ mol⁻¹ nm⁻² were performed to test the influence of the spring constant.

For each nucleotide-dependent actomyosin complex, the interaction energy (V_int) was calculated as the sum of the Coulomb and the van der Waals energy terms and sampled regularly together with the intermolecular distance (r) during the simulations. The intermolecular distance was calculated as the distance between the centers of mass of the myosin and actin monomer.

The interaction energy vs. intermolecular distance (r) data points were fitted using a third order polynomial function. Consequently, the interaction force (F_int) was estimated as the first derivative of the obtained curve: (3)Fint=-dVintdr

The binding stiffness (k_int) was evaluated as the second order derivative of the interaction energy with respect to the intermolecular distance, calculated at the energy minimum distance (r_Vmin): (4)kint=d2Vintdr2|r=rV min