Fascin − F-actin interaction studied by molecular dynamics simulation and protein network analysis

Abstract Actin bundles are an important component of cellular cytoskeleton and participate in the movement of cells. The formation of actin bundles requires the participation of many actin binding proteins (ABPs). Fascin is a member of ABPs, which plays a key role in bundling filamentous actin (F-actin) to bundles. However, the detailed interactions between fascin and F-actin are unclear. In this study, we construct an atomic-level structure of fascin − F-actin complex based on a rather poor cryo-EM data with resolution of 20 nm. We first optimized the geometries of the complex by molecular dynamics (MD) simulation and analyzed the binding site and pose of fascin which bundles two F-actin chains. Next, binding free energy of fascin was calculated by MM/GBSA method. Finally, protein structure network analysis (PSNs) was performed to analyze the key residues for fascin binding. Our results show that residues of K22, E27, E29, K41, K43, R110, R149, K358, R408 and K471 on fascin are important for its bundling, which are in good agreement with the experimental data. On the other hand, the consistent results indicate that the atomic-level model of fascin − F-actin complex is reliable. In short, this model can be used to understand the detailed interactions between fascin and F-actin, and to develop novel potential drugs targeting fascin. Communicated by Ramaswamy H. Sarma

For fascin (Figure1C), it was first found in sea urchin cells, and previous studies have proved its function of binding to F-actin (Kane, 1975;Otto et al., 1979).On the other hand, fascin is highly expressed in various cancers (Chen et al., 2022;Gupta et al., 2021;Lin et al., 2021).Fascin has four b-trefoil domains with residues of 1 À 139, 140 À 259, 260 À 382 and 383 À 493, respectively (Kureishy et al., 2002).Previous studies have made many efforts to explore how fascin binds to F-actin to form bundles.Although the structures of fascin and F-actin have been individually solved, the exact binding site of fascin is still unclear.
Site mutation experiments of fascin indicated that mutations on the phosphorylation sites (S39 and S274) and ubiquitination sites (K247, K250, K41, K43, K241, K353, K399, K464 and K471) affect the function of fascin (Cheng et al., 2021;Lin et al., 2016;Ono et al., 1997;Yamakita et al., 1996).It is speculated that the groove formed by domains 1 and 4 of fascin was a possible F-actin binding region, because the residue mutations in this region disable the function of fascin.In addition, various mutations for fascin were studied by Yang et al. (Yang et al., 2013) Based on their results, three actin binding regions were proposed, i.e. the region between domains 1 and 4, domain2 and domain 3.
As we have known that mutagenesis study can only be used for the speculate of binding region, and the golden standard is the analysis by fascin À F-actin structure.In the study of Aramaki et al., a complex structure of fascin and Factin was solved by cryo-electron tomography (Aramaki et al., 2016).According to their research, fascin straddles between the two F-actin chains in the form of monomer (Aramaki et al., 2016).However, due to the low resolution of 20 nm, the structure model was rather poor and atomic details were not clear.In addition, there were many conformational conflicts in this model and no metal ions reside in Factin.Previous studies have shown that Ca 2þ or Mg 2þ ions bind to actin, and Mg 2þ ion has higher binding affinity (Estes et al., 1987;Selden et al., 1986;Tobacman & Korn, 1983).
Metal ions play important role on the regulation the assembly and mechanical properties of F-actin (Hocky et al., 2016).To add the missing metal ions and obtain a relative reasonable structure, previous studies have built some models for F-actin, such as the well-known Oda model (Oda et al., 2009) and Namba model (Saunders & Voth, 2012).In the F-actin model constructed by Hocky et al., based on the rough structure of Kang et al., Namba-cation model was constructed (Hocky et al., 2016).In their model, authors not only replaced the low-resolution actin protein with the 1J6Z model (Otterbein et al., 2001) to obtain the structure with nucleotide (ADP), but also replaced the existing Ca 2þ with Mg 2þ ions and added another Mg 2þ ion in the 'stiffness' position (Hocky et al., 2016).However, no fascin À F-actin model was constructed.
To understand the detailed interaction between fascin and F-actin, a complex structure with high resolution is needed.However, it is very hard to solve it with experimental study.An alternative method is to improve the poor structure with theoretical methods, e.   Ye et al., 2021).Based on the improved structure, protein Àprotein/ligand interaction can be further investigated (Assaran Darban et al., 2017;Kalhori et al., 2022;Khashkhashi-Moghadam et al., 2022;Marjani et al., 2022;Sharifi-Rad et al., 2021).Thus, MD simulation is often used for the structure refinement and interaction between protein and protein/ligand (Kalhori et al., 2022;Maheri et al., 2022;Taheri et al., 2022;Zare-Feizabadi et al., 2021).
In this study, based on the poor cryo-ET structure, we first replaced the fascin and F-actin with high resolution structures to obtain a new complex, which is the same as previous studies (Aramaki et al., 2016;Saunders & Voth, 2011).Next, the structure was refined by MD simulation.To confirm the reliability of the model, we calculated the key residues on fascin (Figure S1) by the residue contribution to binding free energy and the network analysis, and our results are consistent with the experimental data.Thus, with this refined fascin À F-actin complex model, we can clearly see the interaction details between fascin and F-actin, which provides very important insights for the future study.

Model preparation
Based on the low-resolution structure of fascin À F-actin complex from Aramaki et al. (Aramaki et al., 2016), we replaced the two types of proteins with high-resolution structure (PDB ID:3A5M and PDB ID:3LLP for actin and fascin molecules) (Chen et al., 2010;Murakami et al., 2010).The detailed procedure can be found in supplementary documents.In addition, the missing residues in the structure of fascin were added by using Chimera program (Pettersen et al., 2004).For the two F-actin chains, there are an ATP molecule and a Mg 2þ ions in each G-actin.The magnesium ion is coordinate by two oxygen atoms of ATP and four oxygen atoms of surrounding water molecules, forming an octahedral Mg 2þ site.

Molecular dynamics (MD) simulations
Since the constructed structure has some conformational conflicts between fascin and F-actin chains interfaces, we run MD simulations to optimized the geometries.All MD simulations were performed with Amber20 programs for the fascin À F-actin complex system (D.A. Case et al., 2021).The complex was solvated in a 132 � 275 � 258 Å 3 water box.The ff14SB force field was used for protein complexes (Case et al., 2014).The protonation states of residues were determined by PROPKA3 (Olsson et al., 2011;Søndergaard et al., 2011).For histidine, His 65, 96, 193 and 310 on fascin and His 277 on actin were protonated on NE2 atom, while His 135 on fascin and His 89 and 163 on actin were protonated on ND1, and the remaining His residues were doubly protonated.Cl À and Na þ ions were added to neutralized the system and to maintain the ion concentration of 0.15 M. The system contains � 917,600 atoms.
To obtain a refined structure, 50000-step energy minimization (EM) was first performed to improve the conformational conflicts.In order to maintain the octahedrally coordinate of magnesium ions, harmonic restraints were used on magnesium ions and their coordinated atoms with a force constant of 1000 kJ/mol�nm À 2 .After EM, a 20-ps NVT MD simulation was run with the restraints.To equilibrate the system and obtain the geometries without conformational conflicts, the EM and NVT simulations were performed and iterated for 5 cycles, and followed by a 100ps NPT simulation.Finally, a 100-ns NPT production simulation was carried out with the restraints (2000 kcal/mol�Å À 2 ) for the same atoms.The coordinates were saved every 2.5 ps, and total 4000 frames were obtained.

Binding free energy calculation
The absolute binding free energy of fascin is calculated with the end-point method of MM/GBSA (molecular mechanics combined with generalized Born and solvent-accessible surface area solvation).In this method, the binding free energy DG bind is calculated by (1) The binding free energy (DG bind ) of the complex is obtained by complex energy (G complex ) subtracting the protein energy (G protein ) and the ligand energy (G ligand ).DG bind is contributed by the internal energy (DE internal ), the electrostatic energy (DE ele ), van der Waals energy (DE vdW ), the electrostatic solvation energy (DG GB ), the nonpolar contributions (DG SA ), and the entropy.In this paper, the entropy was not considered.In addition, the residue decomposition of free energy was performed to explore the important residues on fascin (Gohlke et al., 2003).

Network analysis of the MD ensembles
Protein structure network (PSN) is important for understanding protein À protein interactions (Zhu et al., 2018).In our study, PyInteraph 2 was used to analyze the PSN of fascin (Sora et al., 2020;Tiberti et al., 2014).The network was calculated with salt bridge (J� onsd� ottir et al., 2014) and hydrogen bond used as molecular interactions.For the salt-bridge interaction, it was considered as a connection when the distance between two residues with opposite charges was less than 4.5 Å.In addition, the duration of this distance should account for more than 20% of the MD trajectory.For the hydrogen-bond interaction, the distance between the acceptor atom and the hydrogen atom on donor should be less than 3.5 Å and the donor-hydrogen-acceptor angle must be greater than 120 � .Likewise, the duration of this interaction was required to be more than 20%.Finally, we combined the data from salt-bridge calculation and hydrogen-bond calculation.The detailed steps can be found from the tutorial of Sora et al (Sora et al., 2020).All PSN data were analyzed by the module in PyInteraph suite, and 3D structure was visualized with PyMOL program (Pasi et al., 2012;Sora et al., 2020).

Construction the model of fascin 2 F-actin complex
To construct the model of fascin À F-actin complex, the position of fascin on F-actin must be identified.In the structure from Aramaki et al., the accuracy of the position of fascin was quite low; many conformation conflicts can be found and no metal ions and nucleotide reside in F-actin (Aramaki et al., 2016).For the latter problem, we construct a new Factin chain with a high-resolution structure of G-actin (PDB ID: 3A5M) (Murakami et al., 2010).The G-actin was used to fit into the F-actin structure from Aramaki et al.This structure of G-actin is similar to the structure of 1J6Z (Figure S2) which has been used by previous study (Oda et al., 2009;Saunders & Voth, 2012).The reason we used 3A5M was that Mg 2þ ion in this structure had a standard six coordination structure (Murakami et al., 2010).In the structure of Aramaki et al., a single F-actin contains 13 actin monomers, that is, the conventional helical cycle (Castaneda et al., 2018).For the later all-atom MD simulation, we reduced the number of actin monomers from 13 to 7. Since the polymerization of three actin monomers is called nucleation (Pollard, 2016), seven actin molecules in each chain have contained more than two actin trimers, which is enough for our study on the binding of fascin.
Next, the structure of fascin from a high-resolution crystal structure 3LLP (Chen et al., 2010) was fitted in the original structure from Aramaki et al. to get an improved complex structure, although a few conformational conflicts can be still found on the interface in this conformation.We obtained a stable complex conformation without conformation conflict (Figure 2) by several geometry optimizations of energy minimization and MD simulations with different ensembles, e.g.short-time NVT and NPT, and a finally 100-ns MD simulation with NPT ensemble was carried out.Compared with the initial conformation, no conformational conflict is found.Moreover, in our model, magnesium ion binds to each actin with ATP and forms a standard six coordination.

Molecular dynamics simulation
Based on the trajectory from 100-ns MD simulation, we calculated the RMSD for Ca atoms of the system.It can be seen from Figure 3A that the geometries of complex have converged.Likewise, Radius of gyration (Rg) of the whole system shows that the overall structure changes only in the range of 0.5 Å (Figure S3).For each actin protein, we also calculated the relative root-mean-square fluctuation (RMSF), and found that the actin proteins contact with fascin have different fluctuations (Figure S4).The residues close to fascin have lower RMSF values (Figure S4), indicating the fascin has significant effect on F-actin.
In addition, we examined the coordination of Mg 2þ ion.In actin, Mg 2þ ion is coordinated with two oxygen atoms  (O1B and O2G) of ATP and four water molecules.In our MD simulation, as shown in Figure 3B, ATP and the four water molecules are stable around the Mg 2þ ion; the distance between Mg 2þ ion and the O1B/O2G of ATP is �1.9 Å, and the distances of Mg 2þ À O(water) are �2.0À 2.2 Å.These parameters are consistent with the previous study (Hocky et al., 2016).All those data indicates that the structure of the fascin À F-actin complex is reliable anc can be used for subsequent calculations, viz., MM/GBSA and Protein structure network (PSN).

MM/GBSA and energy decomposition
To identify the key residues on fascin, we run the MM/GBSA calculation and binding free energy decomposition by residue.The binding affinity of fascin is À 198 kcal/mol, indicating that fascin has a strong interaction with F-actin and this is consistent with the fact that fascin bundling F-actin to from microfilament (Otto et al., 1979).Then, decomposition of binding free energy of fascin by residue was calculated to further understand the interaction between fascin and F-actin.In this work, the residue with a contribution � À 1 kcal/mol was considered as a key site for binding (Huang et al., 2020).
In the mutagenesis study by Yang et al, the effects of 51 single-site mutations were examined (Yang et al., 2013).As is shown in Figure 4 and Table S1, 13 residues labeled in red had effect for the binding of fascin, and 38 residues labeled in blue had no effect.The decomposition of energies is shown in Figure 4.Among these sites, the experimental strong response sites K22, S39, K41, K43, R149, K358, K379, H392, R398, R408 and K471 are calculated by residual energy decomposition of MM/GBSA, and their contributions to binding free energy are À 0.21, 0.02, À 0.01, À 0.09, À 0.09, 0.34, 0.07, À 0.03, À 0.03, À 0.02 and À 0.02 kcal/mol respectively (Figure 4, Table S1).Our results show that MM/GBSA underestimates the contributions of these 11 residues to fascin binding F-actin.For those residues, all the them except K358 do not locate on the interface between fascin and F-actin (Figure S5 red balls).On the contrary, residues F31, R68, R109 and R344 which locate on the contact surface are overestimated as important sites for the binding with large contributions to binding free energy, i.e.À 3.47, À 1.17, À 1.62 and À 5.72 kcal/mol (Figure 4, Table S1, Figure S5).Therefore, for the residues which are far from the interface, MM/GBSA tends to give small contribution for the binding energy, whereas the contributions are overestimated for the residues which locate in the interface (Figure S5 blue balls).MM/GBSA fail to give correct results for those sites may be caused by the limitation of this method.

Network analysis
Since MM/GBSA gives poor results for the residues those are far from the interface, we calculated the interaction network of residues to construct the relationship among all residues.It is well known that noncovalent interactions in protein are important for its structure and function.The relation between function and noncovalent interaction can be expressed by the network of protein structure (Tiberti et al., 2014).In this study, the network was described by salt-bridge and hydrogen-bond interactions (Fas et al., 2021;J� onsd� ottir et al., 2014).
To identify the key residues, the hubs degree is the only parameter which should be set carefully.In this work, several values of hubs degrees (HD) were preliminary tried and it was finally set to 10 (Figure 5).Comparing to experimental data, residue with the HD > 10 is predicted as a key site, whereas residue is a noneffective site with HD � 10 (Figure 5).(Yang et al., 2013).The effective sites are red whereas the no effective sites are blue.
Overall, for the 51 single-point mutations, residual network results give 5 overestimated sites (K32, R63, D88, R90 and K250 are experimentally no effective sites) and 6 underestimated sites (F29, S39, Q291, K379, H392 and R398 are experimentally effective sites).In order to better understand the results, we analyzed the calculation of each site and the interaction between residues.Our results show that strong associations among strong effect sites are found and mutation of these sites will change the ability of fascin to bind Factin (Table S2).For the interaction of noneffective residues, since they are not involved in those associations, mutation of them do not affect the ability of fascin to bind F-actin (Table S2).For example, E27 and K43 are both experimentally strong-effect sites of fascin (Yang et al., 2013), and they form a hydrogen bond with a duration of 63.18% in the simulation; K42 and are experimentally noneffective sites and they form salt-bridge interactions with a duration of 78.23% (Table S2).The detailed residual network of fascin is shown in Table S2.
However, there are 11 residues those cannot be predicted accurately by this method.It should be noted that the key molecular interaction must maintain for a certain period of time.In our study, salt-bridge or hydrogen-bond interactions with the duration >50% of MD simulation was deemed to be effective.First, we analyzed the 6 underestimated sites, i.e.F29, S39, Q291, K379, H392 and R398.For F29, it forms hydrogen bonds with the backbone of A349 of actin with a duration of 90% in the simulation (Figure S7A), which means that the hydrogen bond is quite strong and the connect to actin is should be very important.However, no other interaction is found for F29, which leads to the underestimation of the role of F29.For S39, the mutation of S39D had strong effect to cause the loss of function experimentally (Yang et al., 2013).However, our previous results showed that the mutation of S39A did not affect the formation of filopodia (Zeng et al., 2017).Our calculation result failed to be consistent with the experimental data of S39D, because this mutation would lead to a new formation of salt bridge (D39 À K41) (Wu et al., 2021).For Q291 and H392, our results show that they do not form any interactions with other residues, which leads to no hub nodes as noneffective sites.For K379, it forms a salt bridge with E288, but the duration of the salt bridge is 40% which lower than the interaction duration threshold we set (Figure S7B).In addition, the backbone of K379 forms a hydrogen bond with V264.E288 and V264 is predicted to be a noneffective site with HD ¼ 8 and 0 (Figure 5).Thus, this region is calculated to be noneffective.For R398, it forms salt bridge with D405 and the duration is 70% in the simulation (Figure S7C).However, D405 is predicted as a non-effect residue with HD ¼ 5.In addition, our experimental results show that D405A is a non-effective mutation (Figure S6).In addition, R398 also forms a hydrogen bond with T403 (Figure S7D).Since the hubs degree of T403 is only 9, T403 is also considered to be an ineffective site.It can be found that the residue would be deemed as noneffective site by network if it strongly connects to a noneffective site, e.g.K379 and R398.
Next, we analyzed the five overestimated residues, i.e.K32, R63, D88, R90 and K250.For K32, its backbone forms hydrogen bond with the backbone of G30 of fascin and S350 of actin (Figure S8AÀ C).However, these hydrogenbond interactions are formed by the backbone of K32, which means that these interactions do not vanish even if K32 is mutated into alanine.For R63, it forms salt bridge with E49 (Figure S8D), but E49 is shown to be an ineffective site in the binding experiment (Yang et al., 2013).Therefore, even if the mutation of R63A destroys the salt bridge, it should not affect the function of fascin.The residue D88 and R82 form a salt bridge, and the hubs degree of R82 reaches 20 (Figure 5).Thus, R82 is deemed as key residue by our calculation, although no experimental data for it.In addition, D88 forms hydrogen-bond interaction with A72 (Figure S8F).Again, A72 is predicted as key residue since its HD is 10 (Figure 5).Therefore, D88 is calculated as very important residues incorrectly because its interacting residues have high HD value.For residue R90, it forms salt-bridge interactions with both D53 and E106, and latter interaction in our MD simulation persists for a longer time (Figure S8G, H).According to the experimental results, E106 is a no-effect site (Yang et al., 2013).Therefore, the mutation of R90 should not affect the fascin binding function.For K250, only one hydrogen bond is formed between its backbone and L253 (Figure S8I) and HD of L253 is only 6 (Figure 5).Therefore, the hub degree of K250 should be wrong.
To clearly show the key residues with high values of HD, we fitted them to the structure of fascin (Figure 6).The important sites are shown as stick model (Figure 6B).Based on our model, the regions where fascin interacts with F-actin have larger values of HD.The central region of fascin is mainly with small value of HD black circle in Figure 6B).Our network results for the distribution of key residues are in good agreement with the previous experimental data, which shows the reliability of the model and calculation results.In addition, it should be noted that hubs degrees and interactions of each residue should be checked carefully and the duration of these interactions must be considered (for this system, duration should be larger than 50%).

Discussion
The invasion and metastasis of malignant tumor cells in vivo has always been an important reason for the difficulty of treatment of tumors.Actin cytoskeleton plays an important role in the movement and invasion of tumor cells (Fife et al., 2014;Kelley et al., 2008).Filopodia and lamellar pseudopodia composed of actin promote the metastasis of cancer cells (Machesky, 2008;Yamaguchi et al., 2005).The polymerization of actin requires the participation of a certain concentration of divalent oxygen ions (Angelini et al., 2005;Castaneda et al., 2018;Tang & Janmey, 1996).The presence of cations has an effect on the hardness and elongation of F-actin (Castaneda et al., 2018).In vivo, there are also many ABPs involved in the polymerization of actin (Zhang et al., 2021).Understanding the interaction between these ABPs and Factin can help us better to understand the metastasis process of cancer cells and to design inhibitors to block the metastasis and diffusion of cancer cells in this process more specifically.
As an ABPs, fascin directly participates in the binding of F-actin in the form of monomer, helping to form a stable and robust bundle structure (Maier et al., 2015).However, the exact actin binding region on fascin has not been determined.Mutagenesis studies have also been carried out for fascin to explore the effect of mutations on the binding function of fascin (Yang et al., 2013).Three actin binding regions on fascin were inferred based on the mutagenesis studies.On the other hand, cryo-electron microscopy study gave a possible pose for fascin on F-actin (Aramaki et al., 2016;Yang et al., 2013).
In this study, based on a poor cryo-electron microscope structure of Aramaki et al., we refined this structure by MD simulation (Figure 2).The combined MM/GBSA and residual network results show that this model is rational and can be used for the further study (Yang et al., 2013).
First, we obtained the converged geometries of fascin À Factin complex.Whole complex system becomes stable in 100 ns (Figure 3A).For the metal center of actin, the Mg 2þ ion maintains octahedral coordination condition with four water and ATP molecules (Figure 3B).
Second, based on the trajectory of MD simulation, we run MM/GBSA calculation.Our results show that fascin has a strong binding affinity of À 198 kcal/mol.To compare with the experimental results of a large number of residual mutations (Yang et al., 2013), binding free energy decomposition was carried out.Our results show that key residues on fascin have large contribution to binding free energy only if the residues on the interaction surface between fascin and Factin (Figure 4).However, MM/GBSA failed when the key residues are far from the contact surface (Figure S5).
Third, the effect of long-range interactions on protein Àprotein recognition and binding was described by the residual network (Sengupta & Kundu, 2012).In this method, the key residues are identified by the value of hubs degrees (HD).In this work, the residues with HD > 10 are deemed to be key residues.Comparing to the experimental data of 51 single-point mutations, 40 of all mutations are consistent with the experimental results directly based on the value of HD.For the other 11 sites, the detailed interactions must be analyzed to give reasonable results.For example, some interactions formed by backbones of residues lead to overestimation of the importance of residues since mutation does not change these interactions, such as K32.Therefore, the analysis of these interactions in detail for the residual network is very necessary.
In addition, as is shown in Figure 5, the key residues far from the interface can be identified by residual network.For the key residues K22, K41, K43, R149, K358, K408 and K471, MM/GBSA gives bad result that all of them are predicted to be noneffective sites.With residual network, all of them are calculated as key residues, which is consistent with experimental data.
In summary, the fascin À F-actin model we constructed can reliably show the interaction between fascin and F-actin.Based on this model, the key residues can be identified by MM/GBSA and residual network methods, which are consistent with experimental data.On the other hand, the methods we used in this paper are transferable to other systems for studying and explaining the important residues in protein Àprotein interactions.It can provide prior knowledge for subsequent experiments and the development and optimization of inhibitors.

Figure 2 .
Figure 2. The complex of F-actin with fascin.Each F-actin chain contains seven monomers.The two F-actin chains are crosslinked by fascin.

Figure 3 .
Figure 3. (A)The convergence of the whole system simulation process.(B) The distances of six coordination bonds formed by magnesium ion are calculated over time.The coordination bond formed by magnesium ion and the oxygen atom of ATP is maintained at about 1.9 Å, and the other four with water molecules are maintained in 2.0 À 2.2 Å.

Figure 4 .
Figure 4. Energy decomposition of 51 single site mutation residues with experimental data(Yang et al., 2013).The effective sites are red whereas the no effective sites are blue.

Figure 5 .
Figure5.Protein structure network analysis of Fascin/F-actin system.We use 10 as the threshold to distinguish effective residues from ineffective residues.

Figure 6 .
Figure 6.The tube model of fascin with values of HD setting to b-factor (plotted by pymol program).The thick tube means the large value of HD.