Allostery and molecular stripping mechanism in profilin regulated actin filament growth

Profilin is an actin-sequestering protein and plays key role in regulating the polarized growth of actin filament. Binding of profilin to monomeric actin (G-actin) allows continuous elongation at the barbed end (BE), but not the pointed end, of filament. How G-actin exchanges between the profilin-sequestered state and the filament state (F-actin) to support the BE elongation is not well understood. Here, we investigate the involved molecular mechanism by constructing a multi-basin energy landscape model and performing molecular simulations. We showed that the actin exchanging occurs by forming a ternary complex. The interactions arising from the BE binding drive the conformational change of the attached G-actin in the ternary complex from twist conformation to more flatten conformation without involving the change of nucleotide state, which in turn destabilizes the actin–profilin interface and promotes the profilin stripping event through allosteric coupling. We also showed that attachment of free profilin to the BE induces conformational change of the BE actin and facilitates its stripping from the filament. These results suggest a molecular stripping mechanism of the polarized actin filament growth dynamics controlled by the concentrations of the actin–profilin dimer and the free profilin, in which the allosteric feature of the monomeric actin plays crucial role.


Introduction
Directed growth of actin filaments plays crucial roles in cell morphology and many other fundamental cell processes, such as motility, division and endocytosis [1]. The monomeric actin (G-actin) attaches to the barbed end (BE) of the filament much easier than the pointed end, leading to asymmetric assembling dynamics, which is strictly controlled by the concentrations of available G-actin and the activity of regulatory proteins [1][2][3][4][5][6][7]. Profilin is a key protein regulating the polarized growth dynamics of actin filament [2,[8][9][10][11][12][13][14][15][16]. The monomeric G-actin available for filament growth is largely sequestered by profilin, which eliminates the non-productive filament polymerization [11]. It has been shown that binding of profilin allows continuous elongation at the BE and prevents the G-actin from attaching to pointed end of the filaments [17]. Revealing the molecular mechanism of the directed actin filament growth regulated by the profilin has been the long standing effort in the fields of biochemistry and biophysics [14,[18][19][20][21][22].
Since the G-actin available for polymerization mostly exists in the form of complex with the sequestering proteins, one key molecular event for the G-actin attachment to the BE involves the actin exchange between the profilin-sequestered state and the filament state (F-actin) ( figure 1(a)). In addition, it has been shown experimentally that high concentrations of profilin tend to accelerate the filament depolymerization event at the BE [22]. Consequently, the following questions arise and need to be addressed in order to reveal the molecular nature of the regulated growth dynamics of actin filament, including: (i) how the monomeric actin transfers from the profilin-sequestered state to the BE bound state to support the filament elongation, and whether the interactions arising from BE is enough to drive the conformational change from G-actin state to F-actin state; (ii) how the free profilin induces the depolymerization of actin filament. Experimentally addressing the above questions is extremely challenging due to the complex interplay of multiple molecular processes and the transient nature of the assembling dynamics.
Molecular dynamics (MD) simulations have shown great success in the studies of structural dynamics of actin filament due to the unprecedented spatial and temporal resolutions [23][24][25][26][27]. However, the molecular events involved in the filament growth typically occur at the time scales longer than millisecond, which are far beyond the capability of conventional all-atom MD simulations. On the other hand, the rich experimental data on the structure, thermodynamics, and kinetics of filament growth provide opportunity to construct a reduced computational model. The high computational efficiency of reduced model can significantly extend the accessible timescale [28][29][30][31], and therefore allows to simulate the whole molecular events involved in the filament elongation under the restraints from available experimental data, which will be useful to address the above key questions encountered in the studies of filament growth mechanism.
In this work, we constructed a multibasin energy landscape model to describe the allosteric feature of the actin monomer and the molecular events involved in the filament growth regulated by profilin. The atomic-interaction-based coarse-grained energy function (AICG2+) developed in our previous work [32,33] was used to describe the conformational fluctuations of proteins around the native structures. Some key parameters were further optimized according to available experimental data on the binding affinities between actin and its binding partners. By performing MD simulations and free energy calculations with the multibasin energy landscape model, we revealed an allosteric ternary-complex mechanism of the actin exchange. Particularly, we showed that the contacting interactions arising from the BE protomers are enough to drive the conformational change from G-actin to F-actin. We also revealed a molecular stripping mechanism of the profilin binding induced filament depolymerization. These results demonstrated the curial role of the allosteric feature of the actin monomer and provided new understanding to the molecular mechanism of the filament growth regulated by regulatory proteins.

Protein model and energy functions
In this work, we used the atomic-interaction based coarse-grained model with flexible local potential (AICG2+) to describe the protein motions around the native structures [32,33]. In the AICG2+ model, every residue was represented by a spherical bead located at the C α position. It was developed based on the perfect funnel approximation with the parameters optimized by using the atomistic simulations [34] and can well describe the conformational flexibility around the native structure. The AICG2+ energy function is given by where V bond describes the covalent connectivity of the neighboring coarse-grained beads along the protein chain. V flp loc is the local interactions describing the secondary structure propensity and flexibility of the protein chains [35]. V SB is the structure-based term describing the conformational fluctuations around the reference structure [29,34]. The V exv is the excluded volume interactions. In the above formula, R and R 0 are the coordinates of the coarse-grained beads in a given structure and in the reference structure, respectively. More details of the model can be found in appendix.
The conformational fluctuations of the G-actin can cover the G-actin state and the F-actin state, which requires an energy function with multibasin topography. Therefore, we employed the scheme developed in reference [36] to construct a multibasin energy function describing the allosteric motions of the actin (appendix). The inter-protein energy function for the proteins a and b is given by where λ ab was obtained by fitting the experimentally measured binding affinities of the actin and its partner proteins (appendix, table 1 and figure A.2). V ab SB and V ele are the structured-based term and the electrostatic term for the inter-protein interactions (More details were given in the appendix).
The reference structures of the profilin and G-actin used in the AICG2+ energy function were taken from protein data bank (PDB) with entry 2PAV [37]. The missing loop in the sub-domain 2 of the G-actin was added with the software modeller [38]. The structure-based interactions involving the residues of this loop were removed, which allows large fluctuations of the loop region at the G-actin state.
The reference structure of the actin filament was taken from the PDB entry 6DJM [39]. To simplify the molecular simulations, only two subunits from this structure were included to represent the filament BE. The AICG2+ energy function was applied for the two actin subunits, by which the two actin subunits remain bound together during the simulations.

Molecular simulations and free energy calculations
All the coarse-grained simulations reported in this work were carried out by Langevin dynamics with friction coefficient γ = 0.25 at 300 K using CafeMol package [40]. In the free energy calculations with two-dimensional umbrella sampling [41], two reaction coordinates were used, including the distance between the G-actin and profilin (R AP ) and the distance between the G-actin and BE of filament (R AF ) (appendix). The sampled conformations from the umbrella sampling were re-weighted to construct the free energy profile using MBAR [42]. In addition to the umbrella sampling simulations, we also conducted conventional MD simulations without applying biasing potential to simulate the profilin regulated actin attachment events. Around 400 independent simulations with different initial conditions were conducted. Each trajectory was lasted over 2 × 10 8 MD steps, with the step size of 0.2τ (here τ is the reduced time unit in CafeMol). Python libraries MDTraj [43] was used in the data analysis. The structures were visualized by using PyMOL [44]and VMD [45].

Multibasin energy landscape model capturing the allosteric feature of actin
Structural data showed that actin has a twist conformation when staying at free form or complexed with profilin [37,46], whereas it adopts a more flatten conformation in the filament [39,47,48], demonstrating allosteric behaviors [49,50] (figures 1(b)-(d)). Therefore, the computational model describing the filament growth should be able to capture the allosteric feature of the actin. In this work, we constructed a multibasin energy landscape model to describe the molecular events in the filament growth by combining the AICG2+ energy functions centered at the G-actin state and F-actin state [36]. The resulted double-basin energy function can well describe the large amplitude conformational motions between the F-actin conformation and G-actin conformation (figure A1 in appendix). By considering the allosteric feature of the actin monomer with the multibasin energy landscape model, we can then simulate the profilin coupled filament elongation and depolymerization with the consideration of the G-actin allostery. Such kind of multibasin energy landscape model has been successfully used in previous studies on the functional dynamics for various of protein machines, including molecular motors [51], transporters [52], and enzymatic catalysis [53][54][55][56][57].

Conformational fluctuations of the free and bound G-actin
The conformational flexibility of the free G-actin has been well recognized in previous studies. However quantitatively characterizing the conformational distribution of the free actin in experiment is difficult [58]. Previous all-atom MD simulations showed that although the free G-actin at ATP state mostly adopted a twist conformation similar to that observed in the profilin bound conformation, there is minor but significant probability that the twisted conformation switches to more flatten F-actin like conformation at room temperature [25]. The multibasin energy landscape model introduced above can well describe the conformational transitions between the G-actin and F-actin conformations (figures 1 and A.1). Following the all-atom MD results [25], the multibasin energy landscape model was designed such that the G-actin conformation is ∼2 k B T more stable than the F-actin conformation at free state (figures 1(e) and A.1), through which the nucleotide state can be implicitly included. Therefore the resulting energy function of the G-actin corresponds to the ATP binding state. Then we use this model to simulate the conformational fluctuations at the profilin bound state and BE bound state. The results showed that when the profilin was bound to the G-actin, the twist conformation becomes much more stable (Figures 1(f) and A.1), which is consistent with the all-atom simulations. Such effect of profilin binding is also consistent with the experimental observations that the profilin binding can significantly slow down the attachment of the G-actin to the filament [22]. In comparison, attachment of the free G-actin to the filament BE significantly stabilizes the flatten conformation (Figures 1(g) and A.1). Therefore the above results suggest that attachment of the G-actin to the BE of the filament during the filament growth is accompanied with the The root mean square distance of the actin-profilin interface relative to that in the native complex structure with (red) and without (green) binding to BE of filament. (c) Distribution of the fraction of formed native contacts (Q AF ) for the interface between the actin and filament with (red) and without (blue) the binding of profilin. (d) Root mean square fluctuation of the residues of the actin for the actin-profilin complex (green), actin-filament (blue), and ternary complex (red). The residues belonging to the actin-profilin interface (purple) and the actin-BE interface (black) were labeled by filled squares.
conformational motions of the actin, as featured by the flattening of the actin conformation, demonstrating a binding induced conformational motions. It is worth mentioning that whether the conformational switching from G-actin to F-actin is driven by the change of nucleotide state (involving ATP hydrolysis) or by the contacting interactions of the neighboring protomers is still in debate [39,59,60]. The above results demonstrate that the contacting interactions arising from the BE protomers can drive the conformational change of the attached actin, which is therefore consistent with recent discussions based on the cryo-EM data of the actin filament [39].

G-actin attachment to filament BE occurs via a transient ternary complex
With the success of the multibasin energy landscape model in describing the conformational distributions for the free and bound actin, we use this model to simulate the whole molecular event of the actin attachment to the BE of the filament. The simulation starts from the profilin bound state, which is the dominant state of the G-actin before attaching to the filament. As shown in figure 2(a), the profilin bound actin can successfully bind to the BE of the filament as indicated by the increasing of the Q AF (i.e. the fraction of formed native contacts between the F-actin and the BE of the filament; see appendix for details), leading to a actin-profilin-filament ternary complex. In the ternary complex, the interface between actin and the BE is quite flexible and samples wide range of conformations (figures 2(a) and 3(c)). At the same time, the BE binding induces moderate flattening of the actin conformation and deceasing of the distance between the sub-domains 2 and 4 (figures 2(a) and 3(a)). After forming such ternary complex, either the actin-profilin complex or the profilin alone may dissociate from the filament. Meanwhile, the BE binding induced conformational change of the actin alters the structure of the actin-profilin interface as shown by the slight increase of the interface RMSD and positional fluctuations (figures 3(b) and (d)). Such perturbation of the interface structure destabilized the actin-profilin interactions, facilitating the stripping of the bound profilin (figure 2(a) and movie 1) (https://stacks.iop.org/NJP/23/123010/mmedia), which is essential for the next round of binding of the profilin-sequestered G-actin. The above results clearly demonstrate the crucial role of allosteric coupling during the profilin stripping since the binding sites of the BE and profilin are not overlap. The two-dimensional free energy profile along the reaction coordinates R AF and R AP also showed a well defined free energy basin around the position of R AF = 32 Å and R AP = 53 Å ( figure 2(b)), corresponding to the ternary complex. In addition, the one-dimensional free energy profile showed that the barrier height for the profilin dissociation is much lower from the ternary complex than from the binary actin-profilin complex (figure 2(d)), which is consistent with the above allosteric coupling mechanism. We also conducted control simulations by switching off the electrostatic interactions. The results showed that the electrostatic interactions have significant contribution to the stability of the complex structure (figure A3).

Molecular striping mechanism of the profilin binding induced depolymerization of actin filament
Previous experimental data showed that at high profilin concentrations, the profilin tends to slow down and even depolymerize the filament. Although the steric effect plays crucial role as demonstrated previously, the molecular nature of such effect is not well understood. From the one-dimensional free energy profile shown in figure 2(e), we can also see that the dissociation barrier height of the actin monomer from BE is much lower when the profilin is bound compared to that without profilin binding. Such results suggest that the binding of the profilin to the BE not only prevents the further adding of the free monomers, but also tends to speed up the stripping of the attached monomers at the BE, leading to a molecular stripping mechanism. Apparently, such a molecular stripping mechanism is also mediated by the allosteric feature of the actin protein. Similar mechanism has been used in the genetic NF-κB mediated regulatory network [61].  shows a representative trajectory of a full molecular stripping event, in which the attachment of the profilin to the BE, forming a ternary complex, which is followed by the release of the actin-profilin complex (see also movie 2).
The above unbiased MD simulations of the profilin coupled filament elongation and depolymerization were independently conducted starting from different initial structures, by which we can observe successful elongation or depolymerization events within reasonable simulation time length. However, the results should be similar to those from one much longer MD simulation with continuous elongation and depolymerization events.
Notably, the profilin can either release without stripping the actin monomer or release as actin-profilin complex, with the proportion depending on the binding affinities of the two kinds of complexes. However, when the concentration of free profilin is high enough, the depolymerization of the actin by molecular stripping mechanism will become more significant, which can lead to the experimentally observed acceleration of the filament depolymerization at high profilin concentrations ( figure 5(b)) [22]. On the contrary, when the concentration of profilin-sequestered G-actin is high, the G-actin tends to continuously attach to the BE, with the bound profilin being stripped through allosteric coupling ( figure 5(a)).
To demonstrate the effects of the concentrations of free profilin and profilin-sequestered G-actin, we constructed a highly simplified kinetic model by assuming a steady state condition and neglecting the exchange between the free profilin and the actin-profilin dimer ( figure 5(c)). Based on the above unbiased MD simulations of the elongation and depolymerization, we extracted the transition rates for the individual steps in the kinetic model and extrapolated them to the conditions with different concentrations of free profilin and profilin-sequestered G-actin, with which we can calculate the flux of the elongation J e (appendix) at the steady state. Positive J e represents filament elongation, whereas negative J e represents filament depolymerization. As shown in figure 5(d), at low concentrations of free profilin and high concentrations of profilin-sequestered G-actin, the filament tends to elongate. On the contrary, at high concentrations of free profilin and low concentrations of profilin-sequestered G-actin, the filament tends to depolymerize. Such results support the above discussions on the concentration effect of the profilin regulated filament growth dynamics.

Conclusions
The assembly of actin filament is a delicately regulated dynamic process. In this work, we studied the underlying molecular events involved in the regulated growth of the filament by developing a multibasin energy landscape model with the restraints from available experimental data. We showed that the profilin sequestered G-actin attaches to the BE of filament through forming a ternary complex, in which the interactions arising from the BE protomers modify the interactions in the actin-profilin interface, promoting the profilin dissociation to support the continuous addition of other G-actin. We also showed that high concentrations of profilin can induce the depolymerization of the bound actin from the BE by using a molecular stripping mechanism. Both the G-actin attachment through ternary complex and the profilin induced depolymerization through molecular stripping mechanism rely on the allostetic feature of the G-actin. All these results provided new insights into the molecular mechanism of the filament growth controlled by regulatory proteins, which is essential for fully understanding many key biological processes involving the regulated actin filament growth dynamics.
It is worth noting that attachment of profilin-sequestered actin not always leads to successful elongation of the filament. After forming the ternary complex, when the profilin alone dissociates and the actin remains attached with the BE, an successful elongation step occurs. On the other hand, it was also observed in MD simulations that the profilin dissociates together with the actin, leading to failure of elongation. Therefore multiple trial attachment events may be needed for one successful elongation step. It is interesting to investigate whether other regulatory proteins are involved to increase the attachment efficiency. Undoubtedly, the partitioning of different pathways can be regulated by the concentrations of free G-actin and profilin.
Recently, all-atom level MD simulations were conducted to investigate the molecular mechanism of the polarized elongation of actin filament [23,62]. Such all-atom level simulations enable explicit description of the effect of nucleotide state on the actin conformations at the two ends of the filament. In our work, the nucleotide state was only implicitly considered. It is interesting in future studies to add the nucleotide molecules in the simulation system, by which we can explicitly investigate the effect of nucleotide state in the profilin regulated filament elongation. The current work mainly focuses on the events at the BE of the filament. In a recent work by Horan and co-workers [63], a coarse-grained model was used to simulate the actin polymerization along the BE and pointed end. It is also interesting in future studies to apply the multi-basin energy landscape model constructed in this work to investigate the actin polymerization at the pointed end.
In addition to the profilin discussed in this work, the assembly dynamics of actin filament can also be regulated by a number of other proteins, such as Tβ4, formin, Arp2/3, gelsolin, cofilin, and so on. Particularly, it has been shown that the attachment of the G-actin and its complex with profilin can also be mediated by formin, which makes the characterization of the involved molecule events in the assembly dynamics very difficult. It will be interesting in future studies to investigate the assembling mechanism using the computational model proposed in this work to consider the effect of other regulatory proteins.
In this work, the conformational dynamics of the proteins were described by the atomic interaction based coarse-grained model (AICG2+) [32,33], which is a structure-based model developed within the framework of energy landscape theory of protein folding utilizing a multiscale strategy and has been implemented in CafeMol software [29,34,40]. In this model, each residue was coarse-grained into one spherical bead locating at the C α position, and the energy function can be written as In the above energy function, V a (θ i ) and V dih (φ i ) are the flexible local potentials, which describe generic local propensity of given amino acids. These two terms were constructed based on a coiled library from PDB and were given in reference [35]. R and R 0 are the coordinates of coarse grained beads in a given snapshot and in the reference structure, respectively. r i , θ i , φ i represent the virtual bond length, virtual bond angle and dihedral angle respectively. r ij is the distance between the ith and the jth residues. r 0 i , θ 0 i , φ 0 i , r 0 ij are the corresponding values at the reference structure. ω and ω φ are the width of local potentials. 13 ij , 14 ij and nloc ij are the sequence-dependent interaction strength parameters. The default parameters in CafeMol were used except otherwise mentioned.
The conformational fluctuations of the G-actin can cover the G-actin state and the F-actin state, which requires a energy function with multibasin topography. Therefore, we employed the scheme developed in reference [36] to construct a multibasin energy function to describe the allosteric motions of the actin, which is given by where V R|R G 0 and V R|R F 0 are the AICG2+ potentials with G-actin and F-actin being used as reference structures, respectively. The parameter ΔV modulates the relative stability of the two conformational states. In this work, the value of ΔV was set as 21.0 kcal mol −1 , with which the free actin monomer is dominantly stayed at the G-actin conformational state. Δ was set as 168.0 kcal mol −1 , with which frequent conformational transitions can be observed within reasonable simulation time. The conformational motions of the actin can be monitored by the reaction coordinate χ, which was defined as with negative and positive values of χ representing the F-actin and G-actin states, respectively. In this work, the structure-based energy term between proteins a and b is given by Here the summation index i and j are the residue pairs forming inter-protein contacts in the reference PDB structures. The parameters nloc ij are the sequence-dependent interaction strength parameters as defined in equation (A.1). As the structure-based energy function V ab SB was applied to the native contacts formed in the reference PDB structure, it can reasonably describe the differences of the interactions for the long-pitch contacts and short-pitch contacts formed between the neighboring actin subunits in the filament.
Considering the important role of electrostatic interaction in the binding process for proteins, in this work, the electrostatic interactions between proteins was modeled by Debye-Hückel model [40], with the ionic strength being set to 150 mM (figures A1-A3).
In this work, the free energy profiles were calculated based on two-dimensional umbrella sampling simulations [41] and MBAR reweighting method [42]. The inter-protein distances were used as the reaction  coordinates in the umbrella sampling simulations. In calculating the inter-protein distances, we defined four residue groups and their centers of mass were used to calculate the distances. The groups for the actin and profilin were represented by the residues '106-111, 156-161' and '100-113', respectively (figure A4). The group for the BE of the filament is composed of the residues '191-202' of one subunit and '108-113, 176-181' of another subunit. In the umbrella sampling simulations, the spring constant and distance interval of the umbrella windows were set as 0.25 kcal mol −1 ·Å 2 and 2.0 Å, respectively for the distances less than 30 Å. For the windows with larger distances (larger than 30 Å), the spring constant and distance interval of the umbrella windows were set as 0.0625 kcal mol −1 ·Å 2 and 4.0 Å, respectively.
To more accurately reproduce the experimentally observed binding affinities of the actin and its partner proteins. The inter-protein energy function for the proteins a and b was tuned by a scaling factor λ ab , which were optimized to reproduce the experimental values of dissociation constants (table 1). In addition, the  default parameters for the strengths of the structure-based term of the actin monomer were scaled by 0.80 to give reasonable melting temperature.
The dissociation constant K D for a given protein pair is given by [64] where [Prot] is the protein concentration, which was estimated based on the size of the simulation box. p u and p b = 1 − p u are the fractions of unbound state and bound state, respectively. In order to calculate p u ,  Here Q ab represents the fraction of formed native contacts for the interface between the protein a and b. It was used to describe the similarity of the interface formed between protein a and b in the structures sampled by MD simulations compared to that in the reference PDB structure. It was defined by Here N nat is the number of contacts formed between two proteins a and b in the PDB structure. The summation index i and j represent the residue pairs forming the interface native contacts in the PDB structure. β and γ are parameters with the values of 10.0 and 1.2, respectively. In addition to the umbrella sampling simulations, we also performed long time scale MD simulations without introducing biasing potential to characterize the filament elongation dynamics containing the profilin, G-actin, and BE of the filament. Around 400 independent MD simulations with different initial conditions were conducted, with each simulation lasting over 2.0 × 10 8 MD steps.
To demonstrate the effects of the concentrations of the free profilin and profilin-sequestered G-actin, we constructed a highly simplified kinetic model by assuming a steady state condition and neglecting the exchange between the free profilin and the actin-profilin dimer ( figure 5(c)). The kinetic model contains four components, including the free BE of filament, the actin-profilin-BE ternary complex, the profilin-sequestered actin, and the free profilin, with the concentrations of [BE], [ternary], [AP], and [P], respectively. The total concentration of the free BE and the actin-profilin-BE ternary complex is fixed (C t ). In addition, the concentrations of the profilin-sequestered actin and the free profilin are the input quantities. Then the time evolving of the concentrations of the free BE and the actin-profilin-BE ternary complex can be given by:  1 [ternary] at the given concentrations of the profilin-sequestered actin and the free profilin. In the calculations, the C t was set as 1.0 C 0 , with C 0 corresponding to the effective concentration given by the box size of the MD simulations. k 1 , k 2 , k 3 , and k 4 are the transition rates, which were extracted from the unbiased MD simulations (table 2). The time unit t 0 used in the kinetic model was set as 1.0 × 10 9 MD steps.