Computational Study of the Binding Mechanism of Actin-Depolymerizing Factor 1 with Actin in Arabidopsis thaliana

Actin is a highly conserved protein. It plays important roles in cellular function and exists either in the monomeric (G-actin) or polymeric form (F-actin). Members of the actin-depolymerizing factor (ADF)/cofilin protein family bind to both G-actin and F-actin and play vital roles in actin dynamics by manipulating the rates of filament polymerization and depolymerization. It has been reported that the S6D and R98A/K100A mutants of actin-depolymerizing factor 1 (ADF1) in Arabidopsis thaliana decreased the binding affinity of ADF for the actin monomer. To investigate the binding mechanism and dynamic behavior of the ADF1–actin complex, we constructed a homology model of the AtADF1–actin complex based on the crystal structure of AtADF1 and the twinfilin C-terminal ADF-H domain in a complex with a mouse actin monomer. The model was then refined for subsequent molecular dynamics simulations. Increased binding energy of the mutated system was observed using the Molecular Mechanics Generalized Born Surface Area and Poisson–Boltzmann Surface Area (MM-GB/PBSA) methods. To determine the residues that make decisive contributions to the ADF1 actin-binding affinity, per-residue decomposition and computational alanine scanning analyses were performed, which provided more detailed information on the binding mechanism. Root-mean-square fluctuation and principal component analyses confirmed that the S6D and R98A/K100A mutants induced an increased conformational flexibility. The comprehensive molecular insight gained from this study is of great importance for understanding the binding mechanism of ADF1 and G-actin.


Introduction
Actin is a highly conserved protein. As one of the most abundant proteins in most eukaryotic cells, actin plays important roles in cellular functions such as endocytosis, organelle movement, cell division, cell mobility, and maintenance of cell shape [1][2][3]. Those functions are influenced by rapid transitions between monomeric (G-actin) and filamentous (F-actin) states regulated by a large number of actin-binding proteins (ABPs) in the cell, such as capping proteins and severing proteins. Actin-capping proteins bind to actin and enhance filament depolymerization, and actin-severing proteins enhance fragmentation. Actin-depolymerizing factor (ADF)/ cofilin proteins are actin-severing proteins, and are highly conserved among eukaryotes [4][5][6][7][8]. They bind to both monomeric and filamentous actin, and play vital and complicated roles in actin dynamics by controlling the rate of filament polymerization and depolymerization [9]. The ADF/cofilin proteins are involved in primary filament depolymerization, and facilitate actin turnover by severing actin filaments and increasing the rate of dissociation of actin monomers from the pointed ends of actin filaments [8,[10][11][12]. The activity of ADF/cofilin proteins is tightly controlled in response to various cellular activities. In plants, the activity of ADF is regulated by several factors such as N-terminal phosphorylation and pH [13][14][15][16].
In the genus Arabidopsis, there are 11 expressed members of the ADF family that are grouped into four ancient subclasses [17]. AtADF1 belongs to subclass I, and is strongly expressed in an extensive range of tissues including flowers, seedlings, roots, and mature leaves [10,17,18]. In 2000, the crystal structure of ADF1 from Arabidopsis thaliana was determined by Bowman et al.; it was the first ADF/cofilin structure from the plant kingdom to be determined [19]. However, the structure of an ADF in complex with actin was not determined until 2008 when Paavilainen et al. reported the crystal structure of the twinfilin C-terminal ADF-H domain in a complex with a mouse actin monomer (PDB ID: 3DAW), which indicated that the ADF-H domain binds to G-actin with the long α-helix inserted into the hydrophobic cleft between subdomains 1 and 3 of actin [20]. By then, numerous crystal structures of ADF/coffin from different organisms had been determined [21][22][23][24][25][26][27][28].  reported that in A. thaliana ADF1 is predominantly phosphorylated by AtCDPK6 at serine 6, which prevents ADF1 from binding to actin. The subsequent mutation experiment demonstrated that the S6D and R98A/K100A mutants of ADF1 in A. thaliana decreased the binding affinity of the ADF for both actin monomers and filaments [29,30]. Others have explored the mechanism of interaction between ADF/cofilin and actin using computational methods. Wriggers et al. built a structure model of an ADF/cofilin-G-actin complex based on the crystal structure of the actin-gelsolin segment-1 complex by docking and molecular dynamics (MD) simulations [31][32][33][34][35][36]. Sept et al. then studied the association rate of actin monomers bound with ADF based on this model using the Brownian dynamics method [32]. The molecular interaction mechanism between cofilin and actin filaments has also been investigated using all-atom MD simulations, coarse-grained MD simulations, and normal mode analysis, which provided insight into the overall mechanism how ADF/cofilin binding influences the structure and mechanical properties of actin filaments [33][34][35][36]. However, detailed understanding of the direct molecular interactions between ADF and G-actin, and the dynamic behavior after ADF1 mutation in A. thaliana, were still lacking.
Here, we modeled an AtADF1 actin monomer complex structure based on the crystal structure of AtADF1 and the twinfilin C-terminal ADF-H domain in a complex with an actin monomer. We then refined the features of the model in order to perform MD simulations. Molecular Mechanics Generalized Born Surface Area and Poisson-Boltzmann Surface Area (MM-GB/PBSA) methods were used to calculate the binding free energy between ADF1 and actin. Per-residue decomposition and computational alanine scanning analyses were performed to determine the residues that make decisive contributions to ADF1 and actin binding affinity. Root-mean-square fluctuations (RMSFs) and principal component analysis confirmed that S6D and R98A/K100A mutations increased conformational flexibility. Our study provides a comprehensive molecular insight into the binding mechanism of ADF1 and G-actin, and gives an explanation for the reduced binding affinity of mutated ADF1 to G-actin at atomic level. Moreover, such information provides novel clues for further mutation experiments on the ADF/cofilin family.

Construction of simulation systems
The ADF-actin complex was built using A. thaliana ADF1 and actin1, based on the twinfilin C-terminal ADF-H domain in a complex with a mouse actin monomer (PDB ID: 3DAW). The crystal structure of AtADF1 (PDB ID: 1F7S) was obtained from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank [37] and aligned to the twinfilin C-terminal ADF-H domain. The missing residues (2-ANA-4, 129-MDLDVFRSRAN-139) were also built based on this structure by Modeller [38]. The sequence of actin1 (accession number: AEC09427.1) was downloaded from the National Center for Biotechnology Information (NCBI) [39]. A homology model of actin1 was built using the actin monomer of 3DAW as a template (protein sequence identity = 88%) by Modeller [38].

Molecular dynamics simulation
The model was then refined to eliminate bad contacts by using Amber 14. The force field parameters for Mg 2+ ion and ATP were downloaded from the Amber parameter database [40,41]. A standard AMBER ff03.r1 force field was assigned to the protein [42,43]. In the experimental work, the purification for ADF1 and actin were performed under different pH conditions (pH = 7 for ADF1 and pH = 8 for actin) [29]. The ionizable residues were predicted with the same protonation states under both conditions based on the pKa values calculated by the H+ + server [44] and PROPKA3.1 [45,46], which were consistent with the default set of AMBER. The protonation state of ionizable residues was set at the default value for pH 7. Systems were neutralized by adding Na + ions. The built complex was solvated in a rectangular box filled with TIP3P water molecules [47], maintaining a 10-Å distance between any solute atom and the boundary. In the first stage, the water molecules and ions were relaxed by restraining the whole protein and ligands (5000 cycles of steepest descent and 2000 cycles of conjugate gradient minimizations); second, the side chains of the protein were relaxed by restraining the backbone of the proteins and ligands (5000 cycles of steepest descent and 2000 cycles of conjugate gradient minimizations); third, the whole system was relaxed without any restrains (5000 cycles of steepest descent and 5000 cycles of conjugate gradient minimizations). After minimization, the model was gradually heated from 0 to 300 K within 50 ps with the backbone of proteins restrained (500 kcal/mol/Å 2 ) in the NVT ensemble. Then, the model was relaxed within 2.55 ns from 500 to 0 kcal/mol/Å 2 in the NPT ensemble. The final equilibration phase lasted 1 ns without restraints.
All molecular dynamics simulations were performed using the GPU version of the PMEMD engine provided with Amber 14. First, the initial model was subjected to a 40 ns MD simulation to reach equilibrium. A snapshot was then extracted from the equilibrium stage as the wild type (WT) and the mutation was carried out. Four systems, ADF1-WT, ADF1-S6D, phosphorylated form ADF1-S6 phos and ADF1-R98A/K100A in a complex with actin, were subjected to simulations. The detailed protocol was mentioned above. MD simulations were conducted in the WT, ADF1-S6D, ADF1-S6 phos and ADF1-R98A/K10A systems at 100 ns to produce trajectories, respectively.
The covalent bonds to hydrogen atoms were constrained using the SHAKE algorithm, and the Particle Mesh Ewald (PME) method [48] was employed to calculate long-range electrostatic interactions. The cut-off for van der Waals interactions was set to 10 Å. The time step used for the simulations was set to 2 fs. The atom coordinates were saved every 10 ps for subsequent analysis.

Binding free energy calculation
The MM-GB/PBSA methods were applied to calculate the binding free energies between ADF1 and actin [49][50][51]. For each system, 500 snapshots were collected from the 10 ns of the trajectory after equilibration with 20 ps intervals. The binding free energy between ADF1 (receptor) and actin (ligand) was calculated as follows: G complex , G receptor , G ligand are the free energies of complex, receptor, and ligand, respectively. According to MM-PBSA and MM-GBSA theory, the binding free energy (ΔG binding ) is composed of two parts: the gas phase molecular mechanical (MM) energy (ΔG gas ) and the solvation free energy (ΔG solv ). It can also be decomposed into three terms: the molecular mechanical energy term (ΔE MM , a sum of the changes of ΔE int , ΔE ele , and ΔE vdw ), the solvation energy term (ΔG solv ), and the vibrational entropy term (TΔS). ΔE int , ΔE ele , and ΔE vdw are given as internal energy contribution, electrostatic, and van der Waals interaction terms, respectively. ΔE int is canceled between ligand, receptor, and complex by using a single trajectory strategy, and it can significantly reduce the noise in most cases. The change of solvation energy (ΔG solv ) comprises the polar component (ΔG GB/PB ) and the nonpolar component of the desolvation energy (ΔG SA ). The polar solvation contribution (ΔG GB/PB ) can be calculated using the Poisson-Boltzmann (PB) and Generalized Born (GB) equation. The dielectric constant for solvent was set to 80 and for solute was set to 1, 2 or 4, respectively. The nonpolar component of the desolvation energy (ΔG SA ) can be estimated using Eq 6, where ΔA represents the change of the solventaccessible surface area (SASA) of the system calculated using the LCPO algorithm [52], and the fitting coefficients γ and β were set to 0.0072 kcal/molÁÅ 2 and 0 in GB [53], and 0.00542 kcal/ molÁÅ 2 and 0.92 kcal/mol in PB, respectively [53][54][55]. The term (TΔS) in Eq 2 is the change in the conformational entropy upon ligand binding. Here, normal-mode analysis (NMA) was used for the calculation of the conformational entropy [56], which was calculated from the sum of translational, rotational, and vibrational components, and as there was high computational demand, only 125 snapshots extracted from the 10 ns of the MD trajectories after equilibration with 80 ps intervals were used. All the atoms in ADF1-actin complex were used for the normal mode calculation. Each snapshots were subjected to energy minimization for 10000 steps in the presence of a distance-dependent dielectric of 4r(i,j) (where r(i,j) is the distance between two atoms) until the root-mean-square of the elements of the energy gradient vector is less than 0.001 kcal mol -1 Å -1 . The mass-weighted Hessian matrix for each minimized snapshot was calculated and diagonalized by using normal mode analysis. The obtained frequency of the normal mode was used to calculate the entropy.
The calculation error bars are standard errors (SE) calculated using Eq 7; the STD is the standard deviation and N is the number of trajectory snapshots used in the calculation.
To further investigate the molecular determinants of ADF1 binding to actin, the effective binding energies calculated using the MM-PBSA method were decomposed into the contributions from individual residues [57].

Computational alanine scanning
The computational alanine scanning method involves replacing the side chain of a given residue (except glycine or proline) with a methyl group (alanine), then recalculating the absolute binding free energy of the mutated system [58,59]. In this work, the binding free energy of the alanine mutant was calculated using the MM-PBSA approach described above, with the snapshots for the wild type complex. All energy terms were calculated for 500 snapshots along the last 10 ns trajectory with 20 ps intervals. ΔΔG bind was defined by the following equation, where ΔG bind is the summation of the molecular mechanical energy term (ΔE MM ) and the solvation energy term (ΔG solv ).

Principal component analysis
Principal component analysis was carried out using the PTRAJ module of AmberTools. Five thousand snapshots were taken from the MD simulation trajectories. To obtain the proper trajectory matrix, overall translation or rotation motion were removed by fitting the coordinate data to the average structure. The trajectory data were then utilized to generate a covariance matrix between the Cα atoms i and j, defined as: Where x i and x j are Cartesian coordinates of the ith and jth Cα atom, N is the number of the Cα atoms considered, < x i > and < x j > represent the time average over all the configurations obtained in the MD simulation [60][61][62].

Results and Discussion
Dynamics behavior of the WT and mutant systems ADF1 is composed of four α-helices and six β-strands (Fig 1). ADF1 binds to actin with the long α-helix inserted into the hydrophobic cleft (groove) between subdomains 1 and 3 of the actin. The constructed complex was subjected to a 40-ns MD simulation to reach equilibrium. Subsequently, a snapshot was extracted from the equilibrium stage as the WT, and used as a starting point for the mutations. The protein stability of the four systems (ADF1-WT, ADF1-S6D, ADF1-S6 phos and ADF1-R98A/K100A in complex with actin) during the MD simulations was monitored by root mean square deviation (RMSD) of the backbone atoms (S1 Fig). The WT system reached equilibrium after 60 ns with average RMSDs of 1.80 Å. The ADF1-S6D and ADF1-S6 phos systems reached equilibrium after 30 ns with average RMSDs of 1.80 Å and 1.96 Å, respectively. However, the ADF1-R98A/K100A system stabilized after 70 ns with a larger average RMSD value of 2.56 Å. Thus, based on the RMSD results, our MD simulations are reliable enough for further investigation.
To identify which parts of the complex contributed most to protein mobility, the RMSF of the four systems versus the residue number was investigated (Fig 2). Three major sites were distinguished in ADF1 that interact with actin: the N-terminal loop, the α3-helix, and the β6-α4 loop [43]. It suggests that the R98A/K100A mutation in ADF1 induces a larger fluctuation on almost the whole ADF1 protein compared with the WT, especially in the N-terminal loop (residues 2-8) and the β6-α4 loop (residues 116-139) (Fig 2A). The protein mobility of ADF1 in the S6D mutation is slightly larger than in the WT. The residues with higher flexibility include residues 16-29 (belonging to the α1 helix) and Gln122 (belonging to the β6-strand). The protein mobility of ADF1 with phosphorylated Ser6 varies not significant from the WT, compared with the mutated systems. The flexibilities of the actins in the three systems are similar except for residues 36-50, which constitute a long loop belonging to subdomain 2 that is far away from the binding interface. Residues 349-354, which interact with the N-terminal loop of ADF1, became more flexible after ADF1-R98A/K100A mutation ( Fig 2B).
As can be seen from the RMSD and RMSF results, ADF1 was more flexible after mutation. To further explore the dynamic behavior of the complex, representative structures from the last 30 ns of the MD simulation trajectories of the mutation systems were superimposed on the WT according to the secondary structure of actin (Fig 3A). A noticeable rigid rotation of ADF1 versus actin was observed, which brought the N-terminal loop in ADF1 away from the actin subdomain 1, and brought the α4-helix in ADF1 close to the actin subdomain 3 in the ADF1-R98A/K100A system. The angle, represented by the Cα atoms of the three residues (Lys293 in actin, Ser102 and Ser136 in ADF1) was measured to represent the degree of this rotation (Fig 3A). The rotation angles of the four systems were monitored during the MD Fig 1. The overall structure of actin-depolymerizing factor 1 (ADF1) in complex with actin1 monomer. ADF1 is depicted with colored secondary structure (cyan for α-helix, magenta for β-sheet, and salmon for random coil). Actin is colored green. The mutated residues are shown as spheres.
Total binding free energy and per-residue contributions between ADF1 and actin predicted by MM-GB/PBSA methods Generally, a low dielectric constant of ε = 1 is used for solute in Molecular Mechanics Generalized Born Surface Area and Poisson-Boltzmann Surface Area (MM-GB/PBSA) methods [63]. Larger values ε = 2 or ε = 4 are also reported [50,64].  We compared the binding free energies calculated using three different solute dielectric constants, 1, 2 and 4. The results were listed in S1 Table. The predicted binding free energy for the WT system is always smaller than the ADF1-R98A/K100A system under three settings. The predicted binding free energy for the WT system is smaller than the ADF1-S6D system with the setting of ε = 1. When ε = 2 and 4, the predicted binding free energy for the WT system is larger than the ADF1-S6D system, which is not consistent with the conclusion in the experimental work [29]. The optimal value of the protein dielectric constant is still an problem need to be discussed in the literature [65]. In this work, ε = 1 made the best predictions. Based on the best results, the predicted binding free energy with the setting of ε = 1 is summarized in Table 1.
As shown in Table 1, both the MM-GBSA and MM-PBSA results suggest that the phosphorylation of Ser6, S6D and R98A/K100A mutations of ADF1 led to reduced binding affinity. The predicted results are in agreement with the experimental findings that after mutation the binding affinity of ADF1 and actin decreases [29,30]. Further analysis suggested that the main differences in binding free energy between the WT and the mutations arose from van der Waals interactions and the electrostatic energy (sum of the electrostatic solvation free energy and MM electrostatic energy, ΔG GB /ΔG PB + ΔE ele ) contribution. The calculated van der Waals contributions in the WT, ADF1-S6D, ADF1-S6 phos and ADF1-R98A/K100A systems were -107. 15  Therefore, to achieve a more precise quantitative interpretation of binding affinity, per-residue basis binding free energy decomposition was performed to determine the individual energy contribution to the interaction energy (S2 Table). A residue was reported only if the per-residue energy contribution difference between the WT and the mutation systems was larger than 1.00 kcal/mol (Fig 4). The ADF1 S6D mutation affected binding affinity mainly through residues Asp26, Asp27, Ser147, Arg149, Asp294, Lys328, Lys330, Arg337, Gln356, Lys375 and Phe377 in actin, and residues Ala2, Ser/Asp6, Asp93, Arg98, Lys100 and Thr124 in ADF1 ( Fig  4A). The summation of the energy contribution of Glu128 (ADF1) and Arg149 (actin) is -5.19 kcal/mol in the WT system versus 5.35 kcal/mol in the ADF1-S6D system (S2 Table). These two residues form hydrogen bond interactions in the WT system. The summation of the energy contribution of Ser6/Asp6 (ADF1) and Phe354 (actin) is -7.03 kcal/mol in the WT system versus -3.02 kcal/mol in the ADF1-S6D system (S2 Table). These two residues form OH-π interactions in the WT system. To calculate the OH-π interaction strength, 5000 snapshots were extracted to measure the angle of OH with the center of the benzene and the distance between the H atom of OH and the center of the benzene. Only an angle between 65°and 90°and a distance between 2.0 Å and 3.5 Å indicated a strong interaction. We observed that the hydroxyl of Ser6 in ADF1 formed a stable OH-π interaction with the benzene ring of Phe354 in actin in most of the MD trajectories in the WT system (S2A Fig). In the S6D mutation, the hydroxyl group (-OH) was substituted with a carboxy group (-COOH). The electrostatic energy contributions of Phe354 and Ser/Asp6 were 0.13 kcal/mol (WT) and 4.75 kcal/mol (ADF1-S6D), which means that a charged residue is unfavorable for the binding affinity. After S6D mutation, most of the interaction residues became unfavorable to the binding affinity except Asp26, Lys328 and Lys330 in actin, and Ala2 and Arg98 in ADF1. The phosphorylation of Ser 6 affected binding affinity mainly through many residues ( Fig  4B). The significant changes included the following interactions between ADF1 and actin. The summation of the energy contribution of Glu128 (ADF1) and Arg149 (actin) is -5.19 kcal/mol in the WT system versus 2.94 kcal/mol in the ADF1-S6 phos system (S2 Table). The summation of the energy contribution of Ser6 (ADF1) and Phe354 (actin) is -7.03 kcal/mol in the WT system versus 1.85 kcal/mol in the ADF1-S6 phos system (S2 Table). The electrostatic energy contributions of Phe354 and Ser6 were 0.13 kcal/mol (WT) and 9.34 kcal/mol (ADF1-S6 phos ).
We can see that with the R98A/K100A mutation, more residues changed the energy contribution (Fig 4C). The significant change was the increase of the summation of the energy contribution of Glu128 (ADF1) and Arg149 (actin), which is -5.19 kcal/mol in the WT system versus 4.36 kcal/mol in the ADF1-R98A/K100A system (S2 Table). Another significant change was the increase of the energy contribution of Ser6 (belonging to the ADF1 N-terminal loop), and T353, F354, and Q356 (belonging to actin subdomain 1). The significant reduction in the energy contribution was caused by the rigid rotation of ADF1, which brings the N-terminal loop away from the interaction surface and breaks the interactions between the residues mentioned above. In particular, the OH-π interaction between Ser6 in ADF1 and Phe354 in actin disappeared in most of the trajectories in the ADF1-R98A/K100A system (S2B Fig). The rigid rotation of ADF1 brings the α4-helix, the β6-strand, and some residues in the α3-helix close to subdomain 3 in actin, forming new interactions, including residues Gln122 (belonging to the β6-strand) and Arg137 (belonging to the α4-helix) and residues Asp290 and Asp294 in actin subdomain 3, and residues Lys107 and Lys111 (belonging to the α4-helix) and residue Glu169 in actin subdomain 3. The summation of the energy contribution of these residues changed considerably, which is also unfavorable for the binding of the ADF-actin complex. The overall energy contributions from the residues mentioned above were -3.56 kcal/mol (WT system) and 0.18 kcal/mol (R98A/K100A mutation system) (S2 Table).

Computational alanine scanning on ADF1 residues
Residues at 8 Å around the binding surface in ADF1 were calculated except for glycine and proline. The results of computational alanine scanning for ADF1 residues are shown in Fig 5. A positive ΔΔG bind value means that the residue was energetically more favorable than the corresponding alanine residue. The residues with a ΔΔG bind value larger than 2.00 kcal/mol were considered hot spots. As can be seen, there are twelve "hot spot" residues on ADF, which belong to the α3-helix (residues 97-111), namely Lys96, Val97, Arg98, Met101, Ile102, Lys107, and Lys111. The top five hot spots, with ΔΔG bind larger than 4.00 kcal/mol, were Lys98, Met101, Lys107, Thr124 and Glu128. Some of the predicted hot spots, including Arg98, Met101, and Ile102, are supported by a previous site mutagenesis study [66]. The other residues, including Lys96, Val97, Lys107, Lys111, Thr124 and Glu128, which were not identified in the former study [66], are worthy of verification by further mutagenesis experiments.
Larger magnitudes of the atomic fluctuations in the mutation systems than the WT system Principal component analysis identifies and quantifies which collective atomic motions contribute most to the overall motion of the molecule during simulation. PCA on the conformations of the WT, phosphorylated and mutation systems highlighted significant differences in the motion as a result of the mutations and phosphorylation. From the PCA plot, it is clear that eigenvectors computed from the MD trajectory for the WT/ADF1-S6D, WT/ADF1-S6 phos and WT/ADF1-R98A/K100A systems varied greatly, which clearly indicates the difference in protein motions as a result of mutations and phosphorylation (Fig 6). To get a direct observation of protein motion after equilibrium, we performed PCA and presented the results as porcupine plots. We can observe that the motions of the four systems were qualitatively different (Fig 7). The magnitudes of the atomic fluctuations in ADF1-S6D, ADF1-S6 phos and ADF1-R98A/ K100A systems were larger than that in the WT system. The mutation systems also induced large atomic fluctuations in actin. The ADF1-S6D had a more complicated motion direction that resulted from contributions of the N-terminal loop, the α1-helix, and the β4-β5 loop ( Fig  7A). The motion of ADF1-S6 phos varies from the WT and ADF1-S6D, which mainly resulted form the contributions of the α4-helix, and the β4-β5 loop (Fig 7B). It is clear that the motion of ADF1-R98A/K100A primarily resulted from contributions of the N-terminal loop, the β4-β5 loop, and the β6-β4 loop (Fig 7C).

Conclusion
In the present study, the binding mechanism of ADF1 with actin1 was explored using a combined computational protocol. Based on the snapshots from the MD simulations, the binding  free energy between WT, mutated and phophorylated ADF1 and actin was calculated by the MM-PB/GBSA method, which indicated that the binding between ADF1 and actin was tighter in the WT than in the mutated and phosphorylated forms (ADF1-S6D, ADF1-R98A/K100A and ADF1-S6 phos ). The van der Waals and electrostatic energy (the sum of the electrostatic solvation free energy and MM electrostatic energy) contribution differences were the main factors affecting the binding affinity. Further computational alanine scanning found critical residues, such as Arg98, Met101, Lys107, Thr124 and Glu128, and interactions that are important for ADF1 and actin binding affinity. The dynamic behavior studies using RMSF and PCA showed that ADF1-S6D and ADF1-R98A/K100A induced larger flexibility to the protein compared with the WT. Rigid rotation was triggered, which broke the interaction between the N-terminal loop and residues 353-356 in actin subdomain 1. In summary, the present study underlines the use of MD simulations in combination with MM-GB/PBSA free energy calculations to provide a detailed description of the binding mechanism of ADF1 and actin at the atomic level. Through per-residue decomposition and computational alanine scanning, a list of the important residues for ADF1 and actin binding were determined, which provides clues for further mutation experiments on the ADF/cofilin family.  Table. The contributions of the important residues for the binding of ADF1 with actin (kcal/mol). (DOCX)