Conformational Changes in Acetylcholine Binding Protein Investigated by Temperature Accelerated Molecular Dynamics

Despite the large number of studies available on nicotinic acetylcholine receptors, a complete account of the mechanistic aspects of their gating transition in response to ligand binding still remains elusive. As a first step toward dissecting the transition mechanism by accelerated sampling techniques, we study the ligand-induced conformational changes of the acetylcholine binding protein (AChBP), a widely accepted model for the full receptor extracellular domain. Using unbiased Molecular Dynamics (MD) and Temperature Accelerated Molecular Dynamics (TAMD) simulations we investigate the AChBP transition between the apo and the agonist-bound state. In long standard MD simulations, both conformations of the native protein are stable, while the agonist-bound structure evolves toward the apo one if the orientation of few key sidechains in the orthosteric cavity is modified. Conversely, TAMD simulations initiated from the native conformations are able to produce the spontaneous transition. With respect to the modified conformations, TAMD accelerates the transition by at least a factor 10. The analysis of some specific residue-residue interactions points out that the transition mechanism is based on the disruption/formation of few key hydrogen bonds. Finally, while early events of ligand dissociation are observed already in standard MD, TAMD accelerates the ligand detachment and, at the highest TAMD effective temperature, it is able to produce a complete dissociation path in one AChBP subunit.


Introduction
The transmembrane nicotinic acetylcholine receptors (nAChRs), belonging to the so-called 'Cys-loop' super-family of ligand-gated ion channels (LGICs) [1,2], are involved in a variety of biological processes [3][4][5] and have been implicated in the onset of Alzheimer's disease [4,6] and nicotine addiction [7]. The nAChR channel pore, located in the transmembrane domain (TMD), opens following the binding of agonist ligands in the orthosteric site located in the extracellular domain (ECD) of the protein. Conversely, the channel is mainly closed in the resting state, when either no ligand is present or an antagonist is bound, and in the desensitized, agonist-bound state. The main endogenous agonist for nAChRs is acetylcholine (ACh). Nicotine also binds nAChR as an agonist while lobeline (Fig. S1) has a full or a partial agonist effect [8,9].
Currently, only limited information is available on the atomic structure of nAChRs. Low resolution structures of Torpedo acetylcholine receptor with closed [10] and open [11] pore were obtained from electron microscopy data. Structures of distant prokaryotic homologues of nAChRs, present in Gloeobacter violaceus and in Erwinia chrysamthemi, GLIC and ELIC, were solved during the last years by high-resolution crystallography [12][13][14][15][16][17][17][18][19]. Interestingly, the GLIC structures put in evidence a binding site for anesthetics in the TMD [20], ketamine binding to the orthosteric site [21] as well as a new semi-closed state of the channel [22]. The ELIC structures showed the binding of anesthetics both in the TMD and in the ECD [23] and the binding of acetylcholine to the ECD [21]. However, despite these new pieces of information, the nAChR gating conformational transition is not yet fully understood. Valuable insight on the ligand binding mechanism came from studies of a water-soluble homologue of the nAChRs ECD, the pentameric acetylcholine binding protein (AChBP). Indeed, AChBPs have been crystallized bound to different ligands displaying agonist or antagonist effects on nAChRs [24,25]. The structures of unliganded (apo) [26,27] and liganded (holo) state of AChBP [27][28][29][30][31][32][33][34][35][36][37][38][39][40] revealed how the ligands bind to the orthosteric pocket. This pocket, present in each of the subunits, is lined by the so-called loop C and extends at the interface between subunits (Fig. 1). The AChBP structures revealed the influence of the ligand type on the degree of Cloop closure against the protein core [32], as the loop arrangements cluster into three groups: i) the antagonist-bound ''open'' conformation; ii) the apo ''intermediate'' conformations; iii) the agonist-bound ''closed'' conformations. It is worth noting however that a recent AChBP structure [41] questions this model, in the case of the small antagonist ligand dihydro-beta-erythroidine (DHbE).
Although the ECD/TMD interface region [31], involved in discriminating between agonist versus antagonist binding and in the long-range communication between the extracellular binding site and the transmembrane pore, is not present in AChBP, similar structural patterns and pharmacological responses have made AChBP a valuable proxy for investigating ligand recognition by nAChRs.
Due to the lack of high-resolution structural information on the nAChRs, computational methods have been extensively used to probe the conformational transition. Nanosecond time scale Molecular Dynamics (MD) simulations have been performed on homology models of nAChRs [42][43][44]. Along with normal mode analysis, performed by using an elastic-network representation for the protein [45,46] or on a full atomistic protein model [44,47,48], they provided first useful insights into the gating mechanism. However, because the binding-to-gating process in nAChRs takes place in the millisecond timescale in physiological conditions [49], standard MD approaches are not useful. Targeted MD has been used, in which the C-loops have been forced to move from the outward to the inward, agonist bound, position [50,51]. 10 ns MD simulations have been performed on AChBP in complex with acetylcholine [52], providing a picture on the ligand interaction with binding site residues; with nicotine and carbamylcholine [53], focusing mostly on the dynamics of residues and water in the binding pocket rather than on large global motions, and with partial agonists [54]. Spontaneous, although preliminary, conformational changes have been observed in unbiased MD simulation of a homology model of the a7 nAChR ECD both in the apo state and bound to an antagonist toxin [55]. Steered MD was used to study the unbinding of nicotine from the AChBP protein [56]; the unbinding of acetylcholine from a homology model of the ECD of the nicotinic receptor has been investigated, by using as reaction coordinate the distance between the ligand and the binding site [57].
In the present work we use full-atom standard MD and the enhanced sampling method temperature-accelerated MD (TAMD) [58] to study the AChBP holo-to-apo transition. In particular, we investigate the conformational changes involved in the opening and closing of the C-loops. In TAMD, a set of extra variables are introduced, coupled to the original system via collective variables (CVs). The original and new variables are then evolved concurrently in condition of effective adiabatic separation and at two different temperatures, the physiological temperature for the physical system and a higher artificial temperature for the new variables. In this way, the system navigates the free energy landscape associated to the extra variables overcoming barriers that are even much higher than the energy at the physiological temperature. TAMD can be used to reconstruct the free energy landscape from direct sampling via reweighting [58,59], or combined with a mean-force interpolation method [60]. TAMD has already been applied to a variety of rare events sampling studies [61][62][63][64][65][66][67] and has proved to be particularly useful in protein conformational searches [68][69][70][71][72][73][74][75]. Here, we use it to search for AChBP conformations by exploring the space of a set of suitable variables.
The conformational transition between the lobeline-bound (holo) and the apo states of Aplysia californica AChBP is investigated. First, unbiased MD simulations on the hundred of nanoseconds time scale allowed to characterize the relative mobilities of the different AChBP regions and to give an insight on which CVs are suitable to describe the transition between the holo and apo conformations. MD simulations starting from ad hoc modified holo conformations put in evidence the influence of few key residues on the nature of the metastability, or in other words, let us to locate the transition bottlenecks at molecular level. Then, TAMD simulations were performed using as CVs the centers of mass of the C-loop and the cys-loop, to accelerate sensibly the holo-to-apo transition of AChBP.

MD Simulations
MD simulations of AChBP were carried out for several systems in explicit water, each 150 ns in length: the apo pentamer (P), the lobeline-bound pentamer in the presence (P1+L) or in the absence (P1) of lobeline. Details on the starting structures and simulation methods are given in Section ''Materials and Methods'' and in Text S1.
As first, an analysis of the protein secondary structure has been performed along the trajectories. The content in a helices and b strands, calculated using STRIDE [76], mostly superimpose with the corresponding contents in the X-ray crystallographic structures (data not shown) and is quite similar in all trajectories and subunits except for the two small 3-10 helix (residues 61-63) and a-helix (residues [67][68][69][70]. While the shortest disappears along the simulations, the residue stretch at 67-70 still shows a significant percentage of 3-10 helix content along the whole trajectory in each simulation done. We attribute these secondary structure changes to the absence of crystal packing. To further assess the stability of the model systems along the trajectories, the Root Mean Square Deviation (RMSD) of the individual subunits with respect to the starting conformations after equilibration was calculated after removing the roto-translational motions of each individual subunit (Fig. S2). All profiles show a plateau after 20 ns, at about 1.2 Å for P1+L and 1.7 Å for P and P1. In P1+L, the RMSD curves of all subunits are superimposed within 0.5 Å , whereas they are superimposed within 1.0 Å and 1.5 Å in P1 and P, respectively: this means that the absence of ligand allows larger conformational changes and more variability among subunits. This is in agreement with results from previous simulations of AChBP [53].
The atomic Root Mean Square Fluctuations (RMSFs) profiles ( Fig. 2) are similar for the three trajectories P1+L, P1 and P, with the exception of few regions: the C-loop, the F-loop, the a1-b1 loop (residues 9-27) connecting the helix a1 and the strand b1, and the b1-b2 loop. With respect to P, in P1+L the RMSFs averaged among subunits are smaller in these regions. In particular, the loop C displays the largest decrease, followed by the loop F. Indeed, because of the interaction with lobeline the loop C is locked in a specific conformation with limited fluctuations [27]. On the other hand, the fluctuations in cys-loop region displays quite similar values in all simulations. Despite the overall similar RMSFs profiles in trajectories P1+L, P1 and P, the differences between P1 and P show that removing ligand does not produce the same RMSFs profile observed for P in the time window considered. The large difference in the C-loop fluctuations is due to the opening of C-loop in P (the distance d intra C , defined in Section ''Trajectory Analysis'', is 15.1 Å , averaged over the trajectory and over the five subunits) to a larger extent than in P1+ L (averaged d intra C = 13.5 Å ) and P1 (averaged d intra C = 13.8 Å ). Differences in opening correspond to those observed in the starting X-ray crystallographic structures: 2BYS (d intra C = 12.7 Å , averaged over the five subunits) and 2W8E (averaged d intra C = 14.7 Å ). The cys-loop does not display significant displacement with respect to the remaining part of the pentamer, as the distances d inner cys and d outer cys defined in Materials and Methods (Section ''Trajectory analysis'') display variations smaller than their corresponding standard deviations (data not shown). This is most likely related to the structure of AChBP, which encompasses only the ECD of nicotinic receptors.
Distributions of d intra C were calculated along the trajectories, by joining together data from all five subunits. Results are shown in Fig. 3. The P1+L distribution is unimodal, centered on the value 13.4 Å , whereas the P distribution is bimodal with two peaks located at 13.2 and 15.6 Å (Fig. 3, first row). The P1 distribution is also uni-modal, although broader than P1+L and displaying a small shoulder at 12.4 Å . The distributions calculated along the three time intervals: 0-50, 50-100 and 100-150 ns (Fig. 3, second row) are almost superimposed for P1+L and P whereas they display more heterogeneity for P1. This is the sign of a reached stable condition of the P1+L and P trajectories. P1 is not yet in stable conditions and this is probably due to the destabilization produced by removing lobeline. This destabilization is not overcome in 150 ns, as is also visible in RMSD and RMSFs plots (see Fig. S2 and Fig. 2).
Starting from biased conformations of AChBP (conf1 and conf2, described in Text S1, section ''New AChBP models'') removes this destabilization and makes the d intra C distribution evolve along 100 ns toward the bimodal distribution observed in P (Fig. 3, third and fourth rows). The drift is even more pronounced in P1conf2, where the LYS143 sidechain was modified. At variance, in the presence of lobeline (P1conf1+L, P1conf2+L, see Table 1), the d intra C distribution stays unimodal as in the simulation P1+L. The presence of lobeline is thus able to counterbalance the introduced perturbation, even if the orientation of LYS143 is changed.
To summarize, no major change in the average pentamer architecture is observed along the MD trajectories of P1+L, P1 and P. Larger conformational drifts and heterogeneity in atomic fluctuations appear in the absence of lobeline and in the apo system P. The distribution of C-loop openings in P1+L/P1 on one side, and in P on the other side, are respectively unimodal and bimodal, revealing the existence of two metastable states which cannot easily inter-convert in the MD timescale here considered. Biasing the initial sidechain conformations of P1 induces the Cloop opening, in particular if LYS143 is perturbed.

TAMD Simulations
The data discussed so far confirm that the closed to open transition of the C-loop occurs on very long time scales. To accelerate this transition, three sets of TAMD simulations at four different effective temperatures were carried out using systems all with closed C-loop conformations (see Table 1). The first two sets are for unliganded AChBP and were started from the conformations corresponding to the 40 ns snapshot of the P1 and P1+L free MD simulations, removing the ligand in the second case. We indicate them as (P1) i and (P1-40 ns) i , where i indicates the effective energies b b {1 equal to 3, 5, 7 and 10 kcal/mol. The third set is for liganded AChBP, starting from the 40 ns snapshot of the P1+L free MD and is indicated as (P1+L) i .
As a first check on the TAMD simulations we calculated the protein secondary structure content along the trajectories. Similarly to what was observed along the standard MD, the a helix and b strand content along all TAMD trajectories mostly superimpose with the corresponding contents in X-ray crystallographic structures (data not shown), with the exception of the two small 3-10 helix (residues 61-63) and a-helix (residues 67-70), as already observed along the MD trajectories. This points out that during TAMD the protein secondary structure is not altered, with respect to the MD.
In order to check if TAMD has been successful in driving the Cloop opening, the d intra C distributions along the TAMD trajectories were calculated and compared with the bimodal distribution obtained in the standard P simulation (see Fig. 3). Results are shown in Fig. 4  An almost wide unimodal distribution is observed over 0-2 ns at 10 kcal/mol in the unliganded case (see Fig. S3, upper panel); in the presence of lobeline, bimodal distribution is achieved over the first 2 ns (Fig. S3, lower panel). Note, however, that, already at 7 kcal/mol, the d intra C distributions quickly evolve toward a profile Simulation of AChBP PLOS ONE | www.plosone.org much more different than the one observed for P, with d intra C values up to 36 Å . The presence of lobeline is no more an obstacle to the C-loop opening, but TAMD seems to evolve the system to states distinct from the natural apo state P.

Hydrogen Bonds between AChBP Residues
Let us discuss now the role of specific residue-residue interactions in the stabilization of the different C-loop conformations. The comparison of the X-ray crystallographic structures 2W8E (apo AChBP) and 2BYS (holo AChBP) points out that the C-loop closed conformation in the presence of lobeline is stabilized by a network of hydrogen bonds, among which is the one between LYS143 in b strand 7 and TYR188 in the C-loop (see Fig. S4F). In Fig. 5 we report the persistence of this hydrogen bond along our standard MD simulations, defined as the percentage of time along a full trajectory that the bond is formed. For each simulation, data are plotted for the five different protein subunits. As shown in Fig. 5A, the LYS143/TYR188 bond is very rarely formed in P, while it is more formed in most of the subunits in the presence of lobeline. Modifying the initial P1 protein structure makes the hydrogen bond percentage evolve toward the one in the native apo P state, particularly for P1conf2. These structures relaxes plausibly toward the target state, i.e. the unbound one, on the time scale of the unbiased MD simulation, pointing out where the transition bottlenecks could be located, at molecular level. Along the TAMD simulations started from the P1 and P1-40 ns initial conditions (Fig. 5B) the behavior is much more similar to the one observed in the target P protein with respect to the one along the standard MD P1. TAMD produces the evolution toward P spontaneously already at b b {1 equal to 3 kcal/mol, and even more at to 5 kcal/mol. We would like to stress that TAMD is not used here to validate the method against the simulation results on the ad hoc structures, as it could appear, but to explore the free energy underlying the C-loop transition, with the aim to overcome the known bottlenecks. A posteriori, the results obtained remark the usefulness and powerfulness of the TAMD method per se.
The relation between the LYS143-TYR188 hydrogen bond and the degree of C-loop closure is better evidenced by plotting the two dimensional histograms of d intra C and the distance between LYS143 and TYR188 donor/acceptor atoms. In Fig. 6 we show such maps calculated over all unbiased MD trajectories (panels A-C). When the bond is formed or partially formed, as in P1 and P1+L (Fig. 6, panels B and C), the C-loop is in the closed state; in correspondence to larger LYS143-TYR188 distances as in P (panel A), the C-loop is mainly in the open form.
The strength of the LYS143-TYR188 bond and the possibility that its formation may contribute to the C-loop closure was already investigated for the apo AChBP in [50], where the potential of mean force for the interaction of the two above mentioned residues was calculated using umbrella sampling on the LYS143NZ-TYR188OH distance. In that work, the bound-like conformation was stabilized by about 2.0 kcal/mol at a distance of 3Å . The strength of the hydrogen bond was found somewhat weak; according to the authors, this could arise from the flexibility of the C-loop and/or concurrent interactions of both residues with water molecules.
At variance with monodimensional PMF calculation, we study the transition in the full protein, i.e. we accelerate all C-loops together. According to our results, the formation and persistence in time of the bond between LYS143 and TYR188 is strictly anticorrelated to the presence of the lobeline. The hydrogen bond is weak and transiently formed in P1+L, as distance values larger than 3Å are observed (see Fig. 6, panel C), while the distribution in P1 show a larger population in the range 2-4 Å (Fig. 6, panel B). The weakness in P1+L is brought about by the transient hydrogen bond that the lobeline could establish with TYR188 (see below). Therefore, in AChBP bound to lobeline, the hydrogen bond between LYS143 and TYR188 does not have a dominant structural role in maintaining the C-loop closed.
The relation between the LYS143/TYR188 donor/acceptor distance and d intra C was followed also along the TAMD trajectories. Results are shown in Fig. 6, for the TAMD starting from P1 and P1-40ns at b b {1 equal to 3 and 5 kcal/mol (panels D to G), confirming i) how TAMD was able to efficiently push the P1 protein structure out of the basin in which it was trapped along the 150ns MD trajectory (see Fig. 6, panel B), already at b b {1 equal to 3 kcal/mol (see Fig. 6, panel D); ii) the higher energy cost for opening the C-loops and disrupting the LYS143/TYR188 bond when the lobeline is present in the binding pocket, as pointed out by the comparison of the maps in panels H, I.
Another hydrogen bond that has been long discussed for AChBP is the one between SER146 and TYR93 (see Fig. S4A). It has been suggested that TYR93 has a role as ''gatekeeper'' in AChBP [77]: depending on its position, it makes the lobeline pocket inaccessible (as in 2W8E) or open (as in 2BYS). Furthermore, experimental results on full nicotinic receptor [78] suggest that mutating TYR93 (in particular aY93W) causes a change in the structure of the binding pocket that compromises both the mobility of the ligand into and out of its docking site, and the speed of channel opening.
In Fig. 7 we report the persistence in time, along our simulations, of the SER146/TYR93 bond: in P the bond is always present (Fig. 7A), while in P1+L it is never formed; this is due to the presence of the lobeline in the binding pocket. As it will be shown below, the carbonyl group of SER146 interacts with the lobeline hydroxyl group along the entire P1+L trajectory, in four out of five subunits. In the unbiased P1 trajectory, an intermediate behavior between P and P1+L is observed. Once again, biasing the initial P1 protein structure makes the hydrogen bond percentage to evolve toward the one in the natural P conformation, particularly for P1conf2.
Along the TAMD simulations started from the P1 and P1-40ns initial conditions (Fig. 7B) the persistence in time becomes much more similar to the one observed in the target P protein with  respect to the one along the standard MD P1. TAMD produces the evolution toward P spontaneously already at b b {1 equal to 3 kcal/mol. In analogy with the LYS143/TYR188 bond, the two dimensional histograms shown in Fig. 8 point out a relation between close/open C-loops and the formation/disruption of the SER146/TYR93 bond. It is clearly shown how TAMD was efficient in driving the system toward the natural P state already at b b {1 3 kcal/mol; this was better achieved by starting from the P1-40ns initial condition.
The features of the SER146/TYR93 bond in the TAMD exploration are further discussed in Text S1 (see section '' TAMD exploration and the role of SER146/TYR93'', Fig. S8, Fig. S9 and Fig. S10).

AChBP-lobeline Interaction
We now turn to discuss the lobeline/AChBP interactions and conformations along our MD and TAMD simulations of the complex systems. At first, we monitored [27]: (i) four hydrogen bonds connecting lobeline carbonyl and the hydrogen atom He1 of TRP147, lobeline amide and carbonyl of TRP147, lobeline hydroxyl and carbonyl of SER146, and lobeline hydroxyl and carbonyl of TRP147; (ii) three van der Waals interactions: between lobeline piperidine ring and TRP147 indole, between lobeline methyl and TYR188 aromatic ring, and between lobeline methyl and TYR195 aromatic ring. The hydrogen bonds are monitored through the distances between hydrogen donor and acceptor atoms. The van der Waals interactions are monitored through the distances between the centers of mass of the two atom groups. Among the hydrogen bonds observed in the crystallographic structure 2BYS [27], only one, between the lobeline hydroxyl group (LOB-OH) and the SER146 carbonyl (SER146-CO) is formed during the MD simulations (Table S1). Between the rings of TYR188 and TYR195, the lobeline methyl (LOB-CH3) seems to prefer to interact with TYR195 ring, as the distances LOB-CH3/  (Table S1). This marked difference from the X-ray crystallographic structure can be explained by the fact that only the global electronic envelop of the ligand was visible in the Xray crystallographic structure 2BYS, thus hampering a precise definition of the ligand conformation. Among the protein-lobeline distances monitored on the trajectory P1+L, in the subunit P3, the distance between the lobeline hydroxyl (LOB-OH) and the carbonyl groups of SER146 and TRP147, displays a strong transition around 64 ns (Fig. 9A, green and red curve, respectively). This transition is accompanied by the formation of a hydrogen bond between the lobeline amide (LOB-NH) and the carbonyl group of TRP147 (Fig. 9A, black curve). An examination of frames extracted in the 65-100 ns range of P1+L reveals that one of the aromatic  Simulation of AChBP PLOS ONE | www.plosone.org rings of lobeline changed its rotameric state around 64 ns, and the LOB-OH then pointed toward the complementary subunit, as represented in Fig. 9D, establishing transient hydrogen bonds with the hydroxyls of the TYR55 (Fig. 9A, purple curve), TYR93 (Fig. 9B, black curve) and TYR188 (Fig. 9B, green curve). This movement of lobeline is an isolated event, which might describe the initial steps of lobeline dissociation.
The ligand lobeline has been shown to have a partial agonist effect [79], on a7 [8] and a4b2 [9] nAChRs. Although it is not clear what partial agonism means at the molecular level [80], the observation of the lobeline transition would support a model where a partial agonist is an agonist ligand that interacts weakly with the receptor and is thus more prone to dissociation. On the other hand, it has been observed that other partial agonists cocrystallized with AChBP cause an incomplete closure of loop C against the binding site [32]; at variance, in 2BYS the C-loop is capped around the lobeline, similar to what observed with full agonists. This has been also observed in the recent Xray crystallographic structure of AChBP from Capitella teleta, bound to lobeline [81]. Our results (see the d intra C histogram in P1+L in Fig. 3) show that the lobeline-induced tight closing is kept during the simulations. The lobeline OH transition would reconcile these two points of view as the lobeline weakened interaction would arise from isolated dissociation events rather than from different, induced, conformations of C-loop.
In the biased MD simulations P1conf1+L and P1conf2+L, the lobeline transition was not observed. This could be explained by considering that the starting points were produced by modifying, among others, the orientation of GLN186, which strongly influences the position of TYR188 (data not shown).
The lobeline transition is however observed in TAMD trajectories (Table S2), and significantly accelerated in at least one subunit (see Fig. 10 for a representative subunit in each TAMD trajectory). In particular, it is observed progressively in more than one subunit as the effective temperature increases. At  finally in P5 about 27 ns. At b b {1 equal to 7 kcal/mol, it occurs in four out of five subunits: P4 first at 2.5 ns, then P3 at about 5 ns, followed by P5 at about 10 ns and finally in P1 at about 18 ns. At the highest b b {1 , 10 kcal/mol, the event is observed in the five subunits at the following time frames along the trajectory: P4 at about 0.6 ns, P3 at about 3 ns, P2 at about 3.3 ns, P5 at about 6.4 ns, and P1 at about 6.6 ns. Remarkably, we observe that the same chain of events characterizes the lobeline transition in MD and TAMD, i.e. we observed the same key residues bonds rupture/formation events (see Figs. 9 and 10).
In the TAMD at b b {1 equal to 10 kcal/mol, the dissociation of the lobeline from the binding pocket is complete after few nanoseconds (Fig. 11). The ligand exits from the binding site and migrates downward to the solvent; its motion is assisted by a corresponding downward motion of the C-loop. A similar path was obtained for the unbinding of acetylcholine from the ECD of a7 nicotinic receptor using Steered MD in [57]. In that work, the free energy difference between the bound and unbound states was estimated to be about 10-15 kcal/mol, consistent with the effective energy adopted by us. At variance with the paths proposed by Zhang et al. [57] however, the dissociation path observed in the present work was not obtained by forced pulling but occurred spontaneously in the TAMD trajectory. A downward motion like the one we observed would be not supported by the analysis of several crystal structures of AChBP both in the holo and the apo state, in which it is shown that the loop C assumes different intermediate positions mainly in the direction outwards perpendicular to the main axis of the protein. We would like to stress that the ligand dissociation trajectory observed here is a possible one as it is provided by TAMD. Indeed, TAMD trajectories are likely to pass close the most probable reaction path [69] although there is no explicit requirement for this. In particular, the first step of lobeline dissociation is additionally supported as the most plausible part of the dissociation because the same event is spontaneously

Discussion
In the present work, the transition of the Aplysia californica AChBP from the holo to apo state was explored using MD and TAMD simulations. The transition has been qualitatively described using: (i) opening distances of C-loop, (ii) individual hydrogen bonds. The TAMD simulations have been shown to dramatically accelerate the transition observed in some MD simulations only after some ad hoc modifications of the structure. Indeed, to get the holo/apo transition of the loop C in unbiased MD simulations it was needed to modify the native structure of the apo protein (P1), by disrupting the pattern of some key hydrogen bonds. Possible first steps of spontaneous lobeline dissociation by TAMD have been observed. At the highest effective temperature, a path of complete ligand dissociation has been observed in one subunit.
The analysis of MD simulations has shown that the AChBP states bound to lobeline and apo are stable within an interval of 150 ns. This agrees with the timescale of channel opening in nAChRs, which is on the order of 1 ms [49]. AChBP bound to lobeline shows along the MD trajectories limited oscillations around a well defined mean structure. On the contrary, the apo AChBP obtained from 2W8E has a much less defined conformational state due to the large C-loop internal mobility, although nAChRs in the corresponding apo state display a well-defined physiological state with a closed channel.
The initial change of orientation of few residues (ARG57, ARG95, GLN103, GLN184, LYS143) is sufficient to reduce, at least in part, the long metastability of AChBP in standard MD trajectories. These residues have been selected as they display a variability of orientations in the AChBP crystal 2BYS. This variability in the crystal is probably the sign of an equilibrium of AChBP in solution among several conformations. It is interesting to see that pushing the conformations of these mobile residues toward the ones observed in the absence of lobeline (PDB entry: 2W8E), induce a whole transition of AChBP toward the apo state, as seen by monitoring the C-loop behavior. Thus, the relevance of these residues modifications is supported by the plausible effects the modifications have on the ACHBP. One should notice that, among the modified residues, only GLN184 is located in C-loop, and only LYS143 directly interacts with C-loop. Most of the residues are thus located apart from the region displaying the major drift during the transition. This is the sign of a domino effect relating the conformation of these residues to the C-loop opening [82].
The TAMD simulation produces spontaneously the holo-to-apo transition, accelerating by a factor of 10 the transition obtained as described above by altering the initial configuration in standard MD. The TAMD results give very qualitative estimation of the energy barriers which are crossed during the transition, through the observation of the effective energy b b {1 values for which the conformational transition is observed. The value is in the range 3-5 kcal/mol for the AChBP transition, which is found to be in qualitative agreement with the free energy barrier of about 2 kcal/ mol between the open and closed C-loop states of the ECD of a7 nAChR, calculated by umbrella sampling using the LYS143/ TYR188 bond as collective variable [50]. In the presence of lobeline, the acceleration of the transition by TAMD is reduced, and the transition, not seen at b b {1 = 3 kcal/mol, starts to be observed from b b {1 equal to 5 kcal/mol. The increase of about 2 kcal/mol induced by the presence of lobeline is in qualitative agreement with the affinity of 0.3 nM [27] experimentally measured for this ligand on Aplysia californica AChBP. Interestingly, the TAMD acceleration also influences the early steps of lobeline dissociation. Indeed, the complete dissociation of lobeline from the binding pocket is achieved, for b b {1 equal to 10 kcal/mol, in the nanosecond timescale.
The analysis of the hydrogen bonds along standard MD and the TAMD trajectories permitted to describe the transition from the holo to the apo state in terms of the evolution of a network of interactions. Among the hydrogen bonds, the one between LYS143 and TRP188 has been shown previously in the literature Figure 11. Snapshots along the TAMD trajectory. Snapshots taken at 1 ns, 3 ns, 4 ns, 5 ns, 6 ns and 6.3 ns along the (P1+L) 10 trajectory for S3-S2 interface. S3 is in green; S2 is in red, cartoon representation; lobeline in VDW representation. Lobeline is fully in the solvent at 6.3 ns. doi:10.1371/journal.pone.0088555.g011 [50] to be important. But, also other hydrogen bonds like SER146/TYR93 have been found even more relevant to the mechanism of the transition, as pointed out by the TAMD simulations. Most of the residues involved in hydrogen bonds are conserved in nAChRs and we can thus expect similar behavior for them in nAChRs. Their importance on the binding/gating events in nAChRs has been already tested by site direct mutagenesis [83,84], validating the model proposed here.
All the observed conformational changes are quite local, but nevertheless relevant according to what it has proposed in the literature, in particular by Auerbach and coworkers, in the socalled conformational wave model [85,86] where the ensemble of variations in local interactions is thought to produce the overall allosteric mechanism. The results obtained here on AChBP could be then fully relevant for the nAChRs. Furthermore, the effects of lobeline on the TAMD transition agrees with the lobeline function on nAChRs as full agonist, as the presence of lobeline decrease or suppress the acceleration of the C-loop opening.

Starting Structures
It is known that Acetylcholine (ACh) is the natural ligand of AChBP, but, up to now, there is no X-ray crystallographic structure of AChBP bound to ACh. As it is quite important to start from a ligand position assessed by precise information from structural biology, a high-resolution structure of AChBP bound to lobeline was used, rather than a model obtained by docking ACh to an empty AChBP structure. Moreover, the choice of lobeline as ligand is based on its therapeutic importance [8,9]. As starting conditions for the apo and holo AChBP conformation we used the X-ray crystallographic structures PDB entry names: 2W8E [33] and 2BYS [27] of Aplysia californica AChBP. They were the highest resolution AChBP structures available (2.05 Å for 2BYS and 1.9 Å for 2W8E), at the time this work started (a new AChBP apo structure has been deposited later, PDB entry: 2Y7Y at 1.89 Å ). Missing residues in structures 2BYS (first loop after the helix a1) and 2W8E (in loop C) were reconstructed using homology modeling with Modeller v9.8 [87,88]. In what follows, we refer to the original residue numbering in 2BYS [27]. Note that residue numbering is shifted by +2 with respect to the original numbering in 2W8E [33].
In order to investigate the role of selected residues on the stability of the AChBP conformations, we modified their arrangement in the initial protein conformations thus creating alternative starting points for the lobeline-bound AChBP structure. Details are fully given in Text S1 (see Section ''New AChBP models'').

MD and TAMD Simulations
MD simulations of AChBP were carried out for several systems in explicit water: the apo pentamer (P), the lobeline-bound pentamer in the presence (P1+L) or in the absence (P1) of lobeline. The starting points for P and P1+L/P1 simulations were respectively the first five chains, called A to E, of the PDB entries 2W8E and 2BYS, completed as explained. Starting from the new AChBP structures, additional MD simulations were carried out both in the presence and in the absence of lobeline. Set up of the systems and the standard MD simulation protocol are fully described in Text S1 (see section ''Molecular Dynamics simulations'').
As for TAMD, the theoretical basis was originally presented by Maragliano and Vanden-Eijnden [58]. The coupled system of equations describing TAMD, for a system whose configuration is specified by x[R N , is: where i = 1,…3N and j = 1,…m. h(x)~(h 1 (x),:::,h m (x)) are collective variables (CVs) that are functions of the atomic Cartesian coordinates and z = (z 1 ,…,z m ) is a set of particular values of these variables; V (x) is the inter-atomic MD potential; k is the spring constant coupling the CVs with their target values; c and c c are friction coefficient; g is a white-noise process, i.e. a Gaussian process with mean zero and covariance the Boltzmann constant and the T is the system temperature; b b~1 k B T T where T T is an artificial temperature, T TwwT. The aforementioned set of equations describe the motion of x(t) and h(t) over the following extended potential: The advantage of TAMD [58], is that if (i) c c is chosen sufficiently large so as to guarantee that z(t) evolves much more slowly relative to the fundamental variables, and (ii) k is sufficiently large such that z j (t)&h j (x(t)) at any given time, then the z variables effectively evolve according to the equation: where F (z) is the free energy associated to the CVs h(x). Therefore the points in the CVs space are sampled according to the Boltzmann factor exp({ b bF) which is tuned to be more uniform than the associated physical Boltzmann factor exp({bF ) by taking T TwwT. TAMD therefore accelerates the sampling of the conformational landscape through the CVs space by increasing the probability of visiting points with relatively low physical Boltzmann factors.
The procedure for determining suitable values of c c and spring constant k is described in Text S1 (see Section ''Determining the value of fictitious friction c c and spring constant k'', and Fig. S6). TAMD trajectories were initiated from: (i) the 40 ns frame of MD trajectories P1+L and P1, (ii) the 40 ns frame of the P1+L trajectory, to which the lobeline has been removed.
As collective variables we used the Cartesian coordinates of the center of mass of the loop C (residues 183-194) and cys-loop (residues 128-140) of each subunit of AChBP, for a total of 30 variables. We could do this because, when compared with other available enhanced sampling methods, TAMD allows the use of a larger number of variables. Since we used a large number of collective variables and we observed the opening mechanism in the TAMD simulations, we did not compare alternative choices. A posteriori, the analysis of two different distances to assess C-loop opening, including the distance between the loop and the facing monomer surface (see below, ''Trajectory analysis''), shows that although these distances are not directly evolved in TAMD, they sample the range of values from the closed to the open state. In particular, we have considered ''open'' those conformations that showed a value of the distance similar to the one in the crystal structure 2W8E.
TAMD simulations were run from 10 up to 30 ns, with temperature T T of about 1500 K, 2500 K, 3500 K and 5000 K, corresponding to b b {1~3 , 5, 7 and 10 kcal/mol. In TAMD runs, the atomic force coming from the coupling potential and the evolution of the CVs are implemented via a Tcl script linked to the NAMD code.
A list of the MD simulations started from the original conformations and of the TAMD simulations is given in Table 1, together with other MD simulations started from modified conformations. Systems set-up details are given in Table S3.

Trajectory Analysis
In order to characterize the protein conformational changes, we used, following [89], three distances between the centers of mass of Ca atom groups. For the loop C, the distance d intra C (Fig. S5A) internal to each subunit, is calculated between the center of mass of residues 183-194 of loop C in the principal subunit and the center of mass of the residues 143 and 144. For the cys-loop, two distances d inner cys (Fig. S5B) and d outer cys (Fig. S5C) are defined, between the cys-loop center of mass (residues 128-140) and the center of mass of residues 179-181,198-204 and 120-127,49-60,30-43, found respectively in the inner and outer b-sheet of the same subunit. These b-sheet residues were picked up as the ones displaying average Mean Square Fluctuations smaller than 0.8 Å in the P simulation (see Results section), to guarantee a ''rigid'' back-wall against which the cys-loop displaces. To study the role of specific residue-residue interactions, an intra-protein hydrogen bond analysis has been performed along all MD trajectories (see Fig. 5, Fig. 7 and Fig. S7). Details are fully described in Text S1, see Section on ''Hydrogen bond analysis''.    Text S1.

Supporting Information
(PDF)