Microtubules (MTs) are key components of the cytoskeleton and play a central role in cell division and development. MT assembly is known to be associated with a structural change in -tubulin dimers from kinked to straight conformations. How GTP binding renders individual dimers polymerization-competent, however, is still unclear. Here, we have characterized the conformational dynamics and energetics of unassembled tubulin using atomistic molecular dynamics and free energy calculations. Contrary to existing allosteric and lattice models, we find that GTP-tubulin favors a broad range of almost isoenergetic curvatures, whereas GDP-tubulin has a much lower bending flexibility. Moreover, irrespective of the bound nucleotide and curvature, two conformational states exist differing in location of the anchor point connecting the monomers that affects tubulin bending, with one state being strongly favored in solution. Our findings suggest a new combined model in which MTs incorporate and stabilize flexible GTP-dimers with a specific anchor point state.https://doi.org/10.7554/eLife.34353.001
Microtubules (MTs) are dynamic cytoskeletal filaments formed from -tubulin heterodimers that are abundant in eukaryotic cells and involved in many processes that are essential for cell physiology, for example, cell division and neural growth (Hyams and Lloyd, 1994). Unlike other cytoskeletal components in the cell, for example, actin filaments which continuously grow as long as enough G-actin is present, MTs stochastically switch between growing and shrinking phases even under sufficient free tubulin concentrations. This dynamic instability is observed both in vitro and in living cells and allows the MT cytoskeleton to be rapidly remodelled in response to internal and external signals (Mitchison and Kirschner, 1984; Gardner et al., 2011).
Free tubulin requires guanosine triphosphate (GTP) to initiate the formation of new MTs (nucleation) and elongate already existing ones (growth) (Weisenberg et al., 1968; Weisenberg, 1972). In the latter case, GTP-tubulin stacks head-to-tail at the tips of growing MTs (primarily those capped by -tubulin) and forms one-dimensional protofilaments that associate laterally to create a hollow cylinder typically comprising 13–14 protofilaments. Every dimer contains two GTP-binding sites: one nonexchangeable site that is buried in the intradimer contact and one exchangeable site exposed on the outer surface of the -subunit. Upon addition of a tubulin dimer to a growing MT end, -tubulin in the incoming dimer completes the GTP-binding pocket of -tubulin in the preceding dimer, which enables the hydrolysis of GTP in the pocket, i.e. cleavage of GTP into a guanosine diphosphate (GDP) and a -phosphate. MTs grow steadily as long as GTP hydrolysis in the lattice lags behind the arrival of new dimers, which creates the so-called GTP cap at the MT tip (Mitchison and Kirschner, 1984). If a MT occasionally grows at a slower rate than GTP hydrolysis, the MT lattice depolymerizes rapidly. This catastrophe is opposed by a rescue mechanism, possibly due to ‘islands’ of GTP-dimers embedded in the depolymerizing MT lattice, which serve as nucleation checkpoints (Dimitrov et al., 2008). GDP-tubulin released during depolymerization exchanges GDP for GTP and becomes again competent for polymerization and nucleation.
Despite decades of intensive research, the precise mechanism of how the conformational dynamics of tubulin contributes to MT assembly remains elusive (Brouhard, 2015). Analysis of electron microscopy images of MTs and structural data have revealed that tubulin dimers are straight when locked in the MT lattice (Nogales et al., 1998; Li et al., 2002), slightly kinked ( per dimer) at the tips of growing MTs (Chrétien et al., 1995), and strongly kinked outwards ( per dimer) at the tips of shrinking MTs (Mandelkow et al., 1991; Arnal et al., 2000). X-ray crystallography and cryo-electron microscopy (cryo-EM) structures support this evidence and can be tentatively assigned to each of the above polymerization states (Nogales et al., 1998; Löwe et al., 2001; Gigant et al., 2000; Wang and Nogales, 2005). It has hence become clear that tubulin straightening takes place during polymerization, whereas strongly kinked conformations are characteristic of depolymerized GDP-tubulin. Nevertheless, the role of GTP binding in the complex process of MT assembly still remains unclear. One of the central questions being controversially discussed is: does unassembled tubulin exist in different nucleotide-dependent conformational states which may modulate its polymerization kinetics?
The debate of whether tubulin conformation is or is not determined by the nucleotide state of its -subunit originated from the studies of small tubulin oligomers and rings formed during cold depolymerization of pre-assembled MTs (Melki et al., 1989; Müller-Reichert et al., 1998). According to the allosteric model (Figure 1, right branch), there is an equilibrium between two conformations of free tubulin, straight and kinked, under the allosteric control of the nucleotide state. Here, free tubulin dimers bind GTP prior to assembly, which triggers a kinked-to-straight conformational change such that dimer integration into the MT lattice is sterically compatible. This model has been supported by the observation of an intermediate curvature state of tubulin at low temperatures using a nonhydrolyzable GTP analog (Müller-Reichert et al., 1998; Wang and Nogales, 2005). Alternatively, the lattice model (Figure 1, left branch) postulates that free tubulin adopts a kinked conformation in solution irrespectively of the nucleotide state, and only upon integration into the MT lattice, tubulin dimers are forced into a straight conformation. Here, the role of GTP binding is not to control the intrinsic tubulin conformation but to strengthen the dimer-dimer bonds in the lattice. This model relies on observations from MT drug binding assays and small-angle X-ray scattering (SAXS) experiments (Manuel Andreu et al., 1989; Buey et al., 2006; Rice et al., 2008). Thus, the main difference between the two models is that the allosteric model treats the kinked-to-straight transition in tubulin as a cause and the lattice model as a consequence of polymerization. This difference originates from the fact that the unconstrained (MT-free) dynamics of GTP-tubulin are largely unclear. For the lattice or allosteric model to hold, different properties of unassembled dimers are required. Therefore, depending on the allosteric response of GTP-tubulin in solution compared to GDP-tubulin, it should be possible to rule out one or both canonical models, at least in their original form.
Currently, indirect evidence for both models exists. Crystal structures of unassembled tubulin in all relevant nucleotide states, albeit bound to MT-depolymerizing proteins or drugs, have been resolved by several research groups (Gigant et al., 2005; Nawrotek et al., 2011; Ayaz et al., 2012; Prota et al., 2013). These structures did not show marked differences between the different nucleotide states, thus supporting the lattice model. The intrinsic conformational dynamics of both GTP- and GDP-tubulin in solution have also been simulated for several tens of nanoseconds using all-atom molecular dynamics (MD) and both were found to be kinked (Gebremichael et al., 2008). Longer MD simulations (up to 100 ns) have not revealed GTP-induced tubulin straightening either (Bennett et al., 2009; Grafmüller and Voth, 2011; André et al., 2012; Grafmüller et al., 2013). In another computational study, lateral binding of two GDP-tubulin dimers has been shown to shift the conformations of individual dimers toward the straight state; the GTP-state was not analyzed though (Peng et al., 2014). The available structural data combined with the simulation evidence therefore indirectly favor the lattice model. However, external factors (e.g., MT-affecting drugs or proteins) are still critical for the crystallization of tubulin, whereas nanosecond time scales of tubulin dynamics so far covered by MD simulations might be too short for large-scale nucleotide-induced effects to occur.
Recently, the allosteric model has been supported by tubulin polymerization assays and structural experiments which showed that two -tubulin mutations (D417H and R262H) related to ocular motility disorder in humans substantially affect MT dynamics (Ti et al., 2016). Interestingly, while these mutations were not seen to cause any substantial changes within the MT lattice, they significantly promoted MT polymerization. Ti et al., 2016 have concluded that the mutations might exploit an allosteric mechanism that reduces the fraction of kinked vs. straight tubulin dimers relative to wild-type tubulin. Another -tubulin mutation (T238A) known to yield hyperstable MTs in yeast has been recently shown to uncouple the conformational dynamics of tubulin from its GTP-hydrolyzing activity in MTs (Machin et al., 1995; Geyer et al., 2015). The authors demonstrated that this mutation also reduces the propensity of tubulin dimers bound to the MT cap to be kinked. Although the observations of Ti et al., 2016 and Geyer et al., 2015 are indirect and their relation to GTP hydrolysis is not yet clear, they suggest that at least a certain degree of allosteric regulation through GTP binding may still be beneficial for tubulin to enter the MT lattice.
In light of these new findings and the absence of strong evidence for or against either of the current models, the dynamics and flexibility of unassembled tubulin in solution need closer investigation. Here, we address the effect of the nucleotide state on the intrinsic dynamics and energetics of tubulin in the absence of external factors (e.g., MT-affecting drugs or proteins) by all-atom, explicit-solvent conformation and free energy calculations. Currently, high-resolution GTP- and GDP-tubulin structures in various curvature states are available, allowing direct simulations of the corresponding conformational basins without the need for manual nucleotide cleavage or docking. We have performed multiple microsecond long simulations of a tubulin dimer in solution in two different incarnations, GTP- and GDP-tubulin, and compared them to a set of all tubulin structures currently available in the Protein Data Bank (PDB). Analysis of these simulations has identified two states associated with nucleotide-dependent modes of tubulin bending and separated by a high free energy barrier, and has enabled us to quantify their energetics and relative occupancies under equilibrium conditions. Our results resolve previous seemingly contradictory findings and suggest a new combined mechanism of MT assembly.
To study how the nucleotide state affects the dynamics and intrinsic bending of tubulin, we performed multiple microsecond long MD simulations starting from tubulin structures in two different nucleotide states: GTP and GDP (see Materials and methods). For each nucleotide state, two independent -long simulations were started from a straight and kinked conformation, respectively, yielding a total of of tubulin dynamics.
In order to compare the simulated tubulin ensembles to experimentally known structures, we performed a principal component analysis (PCA) on a set of all tubulin structures currently deposited in the PDB (Figure 2, Figure 2—video 1; see also Materials and methods). The PCA revealed the major conformational motion in this set and served as a measure of how intrinsic bending of the tubulin dimer in our simulations compared with that of experimental structures. Figure 2A shows all considered experimental structures on a two-dimensional (2D) plane spanned by the first two PCA conformational modes which describe the largest conformational variations among these structures (see also Figure 2—figure supplement 1 for a simplified 2D sketch demonstrating the basic idea of this analysis). Each data point in this 2D plot represents one experimental structure, meaning that similarly kinked dimers are co-located along the x-axis. The PCA and visual inspection of all investigated structures indicate that these clearly fall into two subpopulations: straight tubulin or TUB and kinked tubulin or TUB (Figure 2B, green and pink points, respectively). Interestingly, no intermediate structures were observed in the known experimental data, suggesting that the population of such intermediates is probably very low and may not be resolved by cryo-EM and X-ray crystallography.
The first conformational mode contributes % to the total structural variation in the set and describes a dimer bending motion (Figure 2C, Figure 2—video 2). Bending occurs around an 'anchor point' located at the intradimer interface that does not move during the transition. This anchor point involves hydrophobic interactions between the H8 helix of -tubulin and the surface of -tubulin. Interestingly, a similar unaffected contact, but between adjacent dimers in MT protofilaments, has been previously shown to be involved in MT lattice compaction upon GTP hydrolysis (Alushin et al., 2014; Zhang et al., 2015). The second conformational mode contributes % to the total structural variation and describes a ‘breathing’ motion of the dimer (Figure 2—video 3), which might be indicative of scaling differences in the experimental structures and/or different protein environment, particularly in those derived from cryo-EM reconstructions of MTs and 2D sheets of crystallized tubulin (Figure 2A, 3JAT/3JAS vs. 1TUB/1JFF). Notably, this difference is not present in the TUB subpopulation which is exclusively constituted by X-ray structures co-crystallized with small ligands, stathmin-like factors or darpin.
Figure 3 compares the four simulated tubulin ensembles to the experimental structures in terms of the above PCA results (green and pink trajectories vs. black dots). Most importantly, the nucleotide state seemed to strongly affect tubulin bending flexibility. While the GDP-tubulin ensembles were on average more 'confined' during the simulation time, GTP-tubulin was seen to explore a broader range of curvatures (Figure 3, motion along conformational mode 1), with many conformations being more kinked than is evident from the experimental structures (Figure 3, left panel, pink ensemble). Interestingly, all simulated structures started from straight TUB conformations (i.e. instantaneously removed from the MT lattice where they are stabilized by neighboring dimers) quickly adopted similarly kinked conformations within several nanoseconds (Figure 3—figure supplement 1, top inserts), which suggests that these straight conformations may be highly unrelaxed in solution. Moreover, none of the simulations started from TUB fully converged onto the curved conformation of unpolymerized tubulin; and conversely, none of the simulations started from TUB reached MT-like conformations within several microseconds. There are, we think, two possible explanations for the observed lack of connection between the TUB and TUB simulations: (a) insufficient sampling of one global state corresponding to the kinked dimer (Figure 3, hypothetically located at nm (GTP) and nm (GDP) along the x-axis), or (b) the presence of two global states of tubulin dynamics separated by a high free energy barrier (Figure 3, minima hypothetically placed at nm and nm (GTP) and nm and nm (GDP) along the x-axis). As shown in Figure 3—figure supplement 1 (main four panels), independent equivalent simulations started from different experimental structures are locally converged, which motivated us to have a closer look at whether there was a hidden energy barrier separating the two sets of simulations in Figure 3.
To test the possibility of two global states and an associated free energy barrier, we systematically derived a reaction coordinate (RC) that simultaneously takes into account both the experimental structures (Figure 2A) and our simulated ensembles (Figure 3), and reveals the highest barrier the dimer has to cross in order to transit from one state to the other (Materials and methods; see also Voß, 2014). This coordinate is termed the ensemble separation RC, and the main idea behind its derivation is sketched in Figure 4—figure supplement 1 for a simple 2D case. We note that, strictly speaking, the nature of the hypothesized barrier in the simulated ensembles may not necessarily be related to tubulin bending because conformational modes 1 and 2 in Figure 3 were derived using only a limited set of experimentally known conformations.
A surprising result of this analysis is that the conformational change associated with the identified ensemble separation RC, indeed, was unrelated to tubulin bending or nucleotide state but rather involved rearrangements in the individual monomers and at the intradimer interface (Figure 4—video 1; discussed in detail in the last Results section). This fact is shown in Figure 4 as 2D free energy profiles along the ensemble separation RC and a common bending RC. This bending coordinate was defined by performing a PCA on the sum of both simulated ensembles as was done for the experimental structures in Figure 2. First, Figure 4 clearly shows no overlap between the TUB and TUB simulated ensembles along the ensemble separation RC, confirming the lack of overlap in Figure 3. Second, the motion along the ensemble separation RC is uncoupled from the bending motion both for GTP- and GDP-tubulin, confirming that the identified free energy barrier is unrelated to both tubulin curvature and nucleotide state. Finally, it can be inferred from these 2D profiles that, in our simulations, tubulin easily probes a wide range of curvatures (motion along the y-axis), independently of the starting conformation and with GDP-tubulin sampling a much more restricted range of curvatures (Figure 4, left vs. right panel).
We assume the reason for the poor ensemble overlap in Figure 3 is likely the exclusive use of experimental structures to deduce the bending transition (Figure 2), which is less informative and in part masks the presence of the high free energy barrier identified by our ensemble separation search. Hence, the combination of experimental data and simulation yields a more complete view on the bending transition and reveals an additional conformational change in unassembled tubulin (Figure 4—video 1) – to our knowledge so far unknown – that is observed both in experimental and simulated structures.
We next focused on a quantitative characterization of the nucleotide-dependent bending motions (Figure 4, vertical axis). The fact that we observed two global free energy minima for tubulin in both nucleotide states suggests that the bending dynamics may be different in either of these. To test this possibility, we treated both minima separately. Figure 5A shows molecular representations of the largest-amplitude bending motions of GTP-tubulin in both free energy minima (see also Figure 5—video 1 and Materials and methods for details on the derivation of these motions). Although projections of these motions onto the corresponding 2D free energy profile look very similar (Figure 5B, green and pink lines), both considerably differ in the way the tubulin dimer bends. While the motion derived from the TUB ensemble (pink mode) mainly describes bending 'orthogonal' to the MT wall, resembling protofilament splaying at the MT end, the motion derived from the TUB ensemble (green mode) includes a ‘tangential’ component, i.e. twisting of the -subunit relative to the -subunit around the imaginary MT axis. To distinguish between these two different motions later on, we refer to them as the splay-bend (SB) and twist-bend (TB) mode, respectively.
Having identified the SB and TB modes of tubulin bending enabled us to directly test the lattice and allosteric models by calculating free energy profiles along these modes (Figure 5B) and to assess which of the two scenarios is more feasible from the energetic point of view. According to the lattice model (Figure 1, left branch), one would expect the profiles for GTP- and GDP-tubulin to be very similar, with kinked conformations having lower free energy values in either case. In contrast, the allosteric model (Figure 1, right branch) predicts that straight conformations of GTP-tubulin would have lower free energy values than kinked conformations, with the opposite being true for GDP-tubulin. Both lattice and allosteric models are unsupported by our free energy calculations, which suggests a different mechanism of allosteric control through GTP binding by tubulin.
Comparison of the free energy profiles for GTP- and GDP-tubulin revealed that GTP-tubulin favors a much broader range of intrinsic curvatures, all of which are almost isoenergetic (Figure 6 and Figure 6—figure supplement 1; free energy differences of ). This observation is valid both for TB and SB motion. Contrary to what is postulated by the allosteric model, GTP does not shift the bending preference of tubulin toward low-curvature states according to our calculations. Rather, the range of curvatures energetically accessible to GDP-tubulin becomes strongly restrained (steeper free energy curves in Figure 6), while intermediate dimer kinking is preserved in either nucleotide state (approximate coincidence of the free energy minima in Figure 6). This suggests that tubulin might gain considerable flexibility through binding GTP without changing its average, intermediately kinked conformation.
We also estimated the free energy of kinking () stored in GTP- and GDP-tubulin that are held in the straight, MT-like conformation. For this aim, the free energy values for the GTP- and GDP-tubulin structures derived from cryo-EM reconstructions of MTs and used in our free simulations (Zhang et al., 2015) were calculated (Figure 6, left panel, black and gray squares). Our estimations yield and , suggesting that it costs less free energy for GTP-tubulin to adopt the straight MT conformation. Importantly, the calculated is in agreement with the most recent estimate of curvature strain in straight GDP-protofilaments of per dimer obtained by laser tweezer measurements (Driver et al., 2017). We also emphasize that these energy costs do not include energetic contributions from neighboring dimers in the MT lattice.
The free energy barrier identified by the ensemble separation search does not depend on tubulin curvature or nucleotide state (Figure 4). It separates the free energy basins with different bending dynamics in either of them (Figure 5; Figure 6). It is hence likely that this barrier, in addition to the nucleotide, controls tubulin’s activity by making one or the other mode more favorable without interfering with the dimer curvature. We therefore analyzed the conformational change associated with the barrier and calculated the energetics of this transition.
Visual inspection of the respective collective motion (Figure 7A and Figure 4—video 1) revealed several large-amplitude internal rearrangements in the - and -tubulin monomers that might contribute to the observed barrier. Those include an upward-downward movement of the central H7 helix in the intermediate domain (residues 206–371) accompanied by a translocation of the H7-H8 loop, which leads to a shift in the position of the anchor point (motions of non-interacting flexible loops on the outer surface were not considered). Analysis of the same transitions in GDP-tubulin yielded a qualitatively similar picture of the intradimer rearrangement (not shown). Hence, crossing the barrier along the ensemble separation RC does not require the dimer to bend but rather induces changes in the monomer structure and at the intradimer interface.
It was not possible to estimate individual basin depths ( and ) as well as the free energy difference between the two states () based on the free simulations (Figure 4), because the transition state was only sparsely populated. We therefore performed free energy calculations to overcome sampling issues in the transition region and to reconstruct the full free energy landscape along the ensemble separation RC (Figure 7B and Figure 7—figure supplement 1). The individual basin depths (relative to the transition state at nm) were found to be and , thus confirming the absence of spontaneous transitions in the -long free simulations. Consequently, the free energy difference between the states was estimated to be , suggesting that the basin corresponding to the SB motion is energetically much more favorable (% of the free GTP-tubulin population will exhibit the SB bending motion).
Our results suggest a new mechanism by which GTP binding by unassembled tubulin is linked to its conformational dynamics. This mechanism reconciles the previously proposed lattice and allosteric models (Figure 1) into a new unified model, hence resolving the discrepancies between previous experimental observations. We propose that tubulin enters the MT lattice via a combination of both mechanisms (Figure 8). In solution, tubulin can exist in two anchor point states, each capable of nucleotide-dependent bending but differing in the way the dimer bends, with the SB bending mode being strongly energetically favored over the TB mode (Figure 4; Figure 7B; Figure 8, horizontal transition in the bottom cycle). The average dimer conformation in either state is intermediately kinked, irrespective of the bound nucleotide (Figure 6), which might explain experimental observations that colchicine binding and SAXS profiles of soluble GTP- and GDP-tubulin are almost identical (Manuel Andreu et al., 1989; Rice et al., 2008). In this respect, our findings are also consistent with the early dimer (Gebremichael et al., 2008) and follow-up protofilament simulations (Grafmüller and Voth, 2011; Grafmüller et al., 2013) where both dimers and protofilaments (GTP and GDP) attained largely kinked conformations on the several tens of nanoseconds timescale. We have now shown that bending flexibility at the intradimer interface is controlled by the nucleotide state and gets strongly enhanced upon GTP binding (Figure 6; Figure 8, upward transition in the bottom cycle). Although allosteric control through GTP binding is still present in our model, it differs conceptually from that proposed by the canonical allosteric model (Figure 1, right branch), wherein GTP binding forces tubulin to undergo a kinked-to-straight conformational change.
In our model, the GTP-induced increase in flexibility does not force tubulin to adopt the straight MT-like conformation, which is in part compatible with the lattice model (Figure 1, left branch). Straightening of kinked tubulin at the MT tip and within the lattice is accomplished through induced fit, i.e. the energy of kinking stress (Figure 6, left) is compensated by lateral and longitudinal bond energies. This compensation is merely a lattice effect and does not involve allosteric regulation via GTP. Notably, this scenario does not contradict the experimental data where MTs polymerized with a non-hydrolyzable GTP-analog, GMPCPP, were seen to depolymerize into significantly less curved protofilaments than normal MTs (Müller-Reichert et al., 1998) and intermediate-curvature states of polymers made of GMPCPP-tubulin were observed by cryo-EM (Wang and Nogales, 2005; Nogales and Wang, 2006). Formation of straight or less kinked GTP-tubulin conformations at MT tips, we think, might not strictly require that those conformations be energetically favored in solution. Holding GTP-tubulin in the straight lattice conformation is less challenging for the MT lattice, as evidenced by calculated for GTP-tubulin () and GDP-tubulin () (Figure 6, left). Less kinked, cap-bound GTP-dimers would be more likely to be observed, given at least one lateral neighbor. Finally, the loaded-spring-like stress caused by being non-zero both for GTP- and GDP-tubulin, together with the axial compaction in the lattice induced by GTP hydrolysis (Alushin et al., 2014; Zhang et al., 2015), could generate destabilizing strain that would be released during MT depolymerization.
The coexistence of the two curvature- and nucleotide-independent anchor point states separated by a large free energy barrier is striking (Figure 4; Figure 7). It also raises the question how such a high barrier of in unassembled tubulin is compatible with its polymerization dynamics. We think there are two possible scenarios. The first possibility is sketched in Figure 8 (top cycle) and assumes that the activation energy to cross the SB-TB barrier for a newly arrived dimer at the MT tip comes from favorable lateral and longitudinal bond formation. The second possibility is that SB and TB dimers are equally well incorporated into the MT lattice, and that the transition may be accomplished much later in the MT body. While we cannot discriminate between the two possibilities at present, we favor the first because: (a) no tubulin structure derived from the MT body has (yet) been deposited showing a dimer in the SB state; and (b) the values for longitudinal and lateral bond energies found in literature seem to be sufficient to antagonize the high SB-TB barrier (VanBuren et al., 2002; VanBuren et al., 2005; Gardner et al., 2011; Castle and Odde, 2013; Kononova et al., 2014). More specifically, the longitudinal bond energy ranges between and (not including the entropic cost of binding from the bulk to the MT tip which is ([Castle and Odde, 2013)), while the lateral bond energy is estimated to be between and . Taking the lower bounds of these estimates, one arrives at a combined (one lateral and one longitudinal) energy of −20.4 kT, which is already enough to counterbalance the high SB-TB activation energy.
Although we cannot provide strong evidence for the functional role of the barrier separating the two anchor point states (Figure 7), it is tempting to speculate that this barrier, in addition to the nucleotide state, may control the kinetics of MT assembly and disassembly, i.e. MT catastrophe rates. In fact, this barrier might resolve the so far unexplained inability of modern mechanochemical computational models of MT assembly (VanBuren et al., 2002; VanBuren et al., 2005; Margolin et al., 2011; Margolin et al., 2012; Castle and Odde, 2013; Zakharov et al., 2015) to reproduce typical MT lifetimes (several minutes), lengths (several microns), as well as their moderate dependence on free tubulin concentration observed in vitro (Walker et al., 1988; Walker et al., 1991; Gardner et al., 2011), which are the key features of MT dynamic instability. This dependence of MT lifetimes and lengths on free tubulin concentration is predicted to be too steep (VanBuren et al., 2002), implying a characteristic time for a catastrophe event under physiologically relevant concentrations being on the order of years (Zakharov et al., 2015). There must be therefore a lack of important mechanical features which causes a ‘hyperstabilization’ of the MT in these models. Our results obtained with atomistic MD simulations may now provide some of these missing features.
In conclusion, our results, combined with previous structural knowledge, extend our understanding of MT dynamic instability. In particular, the new combined model of MT assembly cross-bridges previous contradictory experimental observations and may address one of the longstanding questions in the MT field, namely: why does GTP-tubulin polymerize and GDP-tubulin does not? As shown by our work, assuming that GTP-induced tubulin flexibility, and not the dimer conformation as such, is the driving force for MT assembly resolves this discrepancy and reconciles contradictory experimental data reported previously. Thus, we believe that this new mechanism is an important step toward revisiting our understanding of the MT life cycle and accounts for unexplained complexities of MT growth and catastrophe.
The set of tubulin structures was extracted from the PDB using the BLAST method (Altschul et al., 1990). The template sequences of - and -tubulin were extracted from the structure of a straight GMPCPP-tubulin dimer (PDB ID: 3JAT [Zhang et al., 2015]). At the time of performing the search (April, 2017), the PDB contained tubulin structures sharing at least 80% sequence similarity with the template sequences (sus scrofa). This preliminary set was then aligned by filtering out structures whose chains shared less than 96% sequence identity and whose aligned parts of the sequences covered less than 85% of the template sequence using ProDy (Bakan et al., 2014), i.e. structures that deviated strongly from the reference (3JAT) in terms of sequence or that had too many missing residues were excluded. This additional filtering yielded a set with a total of 91 structures. This final set covered a broad range of tubulin structures: 62 structure in the kinked conformation, 29 structures in the straight conformation as well as various nucleotide states (GTP, GDP, GTP-analogs like GMPCPP etc.). Most of the kinked structures were in the so-called T2S complex with the MT depolymerizing factor Op18/stathmin (Belmont and Mitchison, 1996), whereas all straight structures originated from cryo-EM reconstructions of MTs or 2D sheets of crystallized tubulin. Residue numbering mismatches in the -subunits were fixed. All residues of -tubulin are referred to in the main text according to the corrected numbering.
Principal component analysis (PCA) (Amadei et al., 1993; de Groot et al., 1996) was performed on structural sets/ensembles using only the backbone atoms and excluding the flexible C-termini. To clarify the nomenclature, we briefly describe the essence of the PCA below. Given a set of atomic configurations (e.g., a set of PDB structures or a MD trajectory), each consisting of atoms , the main goal of a PCA is to reduce the dimensionality of the conformational space by finding unit vectors in () which describe most of the ensemble variance. For this aim, after removing translational and rotational motions by alignment of the ensemble with a reference structure (here, PDB ID: 3JAT), the covariance matrix of atomic positions is constructed,
where represents the average configuration. Diagonalization of yields a set of orthogonal unit eigenvectors and eigenvalues (in descending order),
The projection of the ensemble onto the eigenvector, , is termed the principal component or conformational mode and denotes the collective motion along this eigenvector with the variance given by . For molecular simulations, the first 1–20 conformational modes usually account for 80–90% of the ensemble variance (Amadei et al., 1993).
Six tubulin structures were selected for subsequent MD simulations: two straight dimers (PDB IDs: 3JAT (GMPCPP), 3JAS (GDP)) and four kinked dimers (PDB IDs: 5JQG (GTP), 5FNV (GTP), 4ZHQ (GDP), 4ZOL(GDP)). For convenience, we refer to the straight structures as TUB and TUB, and to the kinked structures as TUB and TUB. Missing atoms and residues as well as the C-termini (:437–451, :426–445; presumably unstructured) were added using MODELLER version 9.17 (Fiser et al., 2000). Protonation states of the histidines were assigned using the GMCT package (Ullmann and Ullmann, 2012). The straight tubulin structures were extracted from lattice patch models (3JAT, 3JAS) of which one (3JAT) contained the non-hydrolyzable GTP analog (GMPCPP) in the E-site of -tubulin. GMPCPP was converted into GTP by replacing the carbon atom between - and -phosphate with an oxygen atom, and the new bond lengths and angle relaxed to their equilibrium values during minimization. The kinked structures were extracted from T2S complexes coassembled with multiple small ligands (calcium and chloride ions, ethanesulfonic acid, glycerol, etc.). Stathmin as well as the small ligands were stripped out while preserving the GDP and Mg-coordinated GTP molecules.
GROMACS version 4.6 (Szilárd et al., 2015) and CHARMM22* force field (Piana et al., 2011) were used for all simulations. The simulation setup was similar to that described in (Rauscher et al., 2015). Briefly, for every nucleotide/curvature state TUB, the simulated system consisted of a tubulin dimer centered in a rhombic dodecahedral box filled with CHARMM-modified TIP3P water (MacKerell et al., 1998) and 0.15M KCl ( atoms in total). Center-of-mass subtraction was applied to the solute’s atoms (including GTP, GDP, and Mg in GTP) in order to prevent the dimer from drifting away from the periodic box, which permitted more convenient trajectory analysis.
Prior to the production runs, each system was subject to an initial equilibration phase involving steepest-descent energy minimization followed by 1 ns of MD simulation in the NVT ensemble at 100K with position restraints applied to the solute’s non-hydrogen atoms (only protein) and 6 ns of simulation in the NPT ensemble at 1 atm during which the temperature was gradually increased from 100K to 300K. The LINCS algorithm (Hess et al., 2008) was used to constrain the lengths of bonds with hydrogen atoms, allowing a 2 fs integration time step. The cutoff radius for the Lennard-Jones and short-range electrostatic interactions was set to 0.95 nm, and long-range electrostatics were calculated using the particle-mesh Ewald (PME) method (Essmann et al., 1995) with a 0.12 nm grid spacing. The velocity rescaling thermostat (Bussi et al., 2007) was used for all simulations. Berendsen pressure coupling (Berendsen et al., 1984) was used to maintain the atmospheric pressure during equilibration.
For each TUB, the last frame of the equilibration phase was extracted, and virtual sites (Hess et al., 2008) for the solute were introduced to remove the fastest degrees of freedom. Together with constraining all bond lengths with the LINCS algorithm, this allowed a 4 fs integration time step. The simulation was then continued in the NPT ensemble at 300K using the Parrinello-Rahman algorithm (Parrinello and Rahman, 1981) for pressure control for a total duration of . For TUB and TUB, two independent production runs per system were conducted and combined. For TUB and TUB, production runs started from two different crystal structures were combined. This yielded an accumulated trajectory of for each TUB, covering a total of of free tubulin dynamics. Source files and a step-by-step guide necessary to reproduce the simulations as described above are freely available (https://github.com/moozzz/simulation-protocols/tree/master/free-tubulin-simulation [Igaev, 2018a]; copy archived at https://github.com/elifesciences-publications/simulation-protocols). Unless differently specified, VMD (Humphrey et al., 1996) and UCSF Chimera (Pettersen et al., 2004) were used for visualization and native GROMACS tools for structure analysis.
To study the energetics of and transitions between different conformational states, the ensemble separation search was employed (Voß, 2014), yielding a unit vector that best separates two minima in the multidimensional free energy landscape. Essentially, given two simulated ensembles projected onto a PCA-based vector space, the method seeks the best linear reaction coordinate such that the projections of the probability densities of the two ensembles onto this vector have the smallest overlap. The optimization problem consists in finding the minimum of the following overlap integral:
where are projections of the multidimensional probability densities sampled by the two ensembles onto :
and represents the reaction coordinate along the axis defined by .
In practice, a search in the -dimensional conformational space is unfeasible given a molecular system with thousands of atoms due to the curse of dimensionality. Using the basic property of the PCA, that is, dimensionality reduction, it is however possible to accelerate the search by approximating with a linear combination of only the first PCA eigenvectors. Here, serves as a regularization parameter which is incremented until the solution to Equation 4 becomes independent of .
The ensemble separation RC search algorithm was implemented in a custom-made script (https://github.com/moozzz/orc_search_gauss [Igaev, 2018b]; copy archived at https://github.com/elifesciences-publications/orc_search_gauss). For a given , the search consisted of two major steps. First, 10 random unit vectors were generated (500 scans 2000 vectors) with the coordinates being drawn from independent, one-dimensional normal distributions . In each scan, the vector that corresponded to the minimal overlap in Equation 4 was selected, and a local, gradient-free minimization was carried out by the downhill simplex method (Nelder and Mead, 1965). The vectors resulting from the independent scans were then averaged to yield the final ensemble separation vector . Convergence of the search was assessed by monitoring the variances of individual coordinates of as functions of the number of scans.
A new vector basis was first derived by a PCA on the combined TUB and TUB trajectory for each nucleotide state. The first PCA eigenvectors were then used as a basis set to find the ensemble separation vector that minimized the overlap between the two ensembles, where was varied to account for all relevant motions in the ensembles. In both cases, less than 20 PCA eigenvectors were sufficient to separate the ensembles (Figure 4—figure supplement 2). Finally, the ensembles were projected onto a 2D plane constituted by the ensemble separation vector and an orthogonal vector with the largest variance in the combined ensemble (i.e. the first PCA eigenvector in the -dimensional subspace; Figure 4).
Functional mode analysis (FMA) (Hub and de Groot, 2009; Krivobokova et al., 2012) was applied to derive the collective motion of protein dynamics related to a particular functional quantity – here, tubulin bending. In its partial least-squares (PLS) form (Krivobokova et al., 2012), the FMA yields a vector basis with the lowest possible dimensionality that guarantees optimal correlation between the fluctuations in the functional quantity of interest and free protein dynamics. Briefly, in the PLS-based FMA, a regression problem of the form
is solved, where denotes the functional quantity of interest represented as a vector, is the matrix of atomic positions for an ensemble consisting of structures, is a matrix composed of basis column vectors, and is a vector of coefficients to be optimized. The basis is defined iteratively according to both variance in (similarly to a PCA; see above) and linear correlation between and . Like in the ensemble separation search, serves as a regularization parameter which is incremented until the maximal predictive power of the regression model in Equation 5 is reached. This is achieved by simultaneous cross-validation against an independent set of data. Finally, an ensemble-weighted mode of motion is constructed that correlates best with and has sufficient variance to contribute to . Details on the GROMACS implementation of the PLS-based FMA are given elsewhere (Krivobokova et al., 2012).
We chose the RMSD relative to the straight tubulin structure derived from the MT lattice (3JAT) as the functional quantity of interest (Figure 7—figure supplement 1A). Any deviation from the straight conformation increases the RMSD, irrespective of the bending direction. However, since the RMSD is a nonlinear function of the atomic coordinates, nonlinear correlation contributions cannot be assessed. Nevertheless, an acceptable quality model was obtained for this nonlinear case as shown in Figure 7—figure supplement 1B,C. Correlation coefficients between regression model and data are shown in Figure 7—figure supplement 1B for model building () and cross-validation () as a function of the FMA basis dimension for both TUB and TUB half-ensembles (from to ). The cross-validation showed that the FMA converges at and components for the TUB and TUB ensembles, respectively. Figure 7—figure supplement 1C demonstrates the overlayed RMSD data half-set and the FMA fit shown for both model building and cross-validation parts, using the optimal number of components in each case. High correlation coefficients observed for model building and cross-validation in both cases () confirmed that the obtained FMA models are adequate and reflect the general features of RMSD fluctuations. The second TUB and TUB half-ensembles (from to ) were used as additional and independent cross-validation sets, which yielded and , respectively. Hence, the found FMA bending modes are a robust representation of the relation between the RMSD to the straight reference structure and the ensemble dynamics.
To calculate the free energy profiles along the ensemble separation RC and the reaction coordinates defined by the FMA bending modes, umbrella sampling simulations were carried out (Torrie and Valleau, 1977; Kästner, 2011). In each case, starting configurations (seeds) were derived by projecting a free simulation trajectory onto a FMA bending mode or the ensemble separation RC and selecting structures placed along the bending coordinate with a step of nm. For those simulations, where the projected trajectory did not cover the desired range along the bending coordinate, additional seeds were generated from preceding umbrella windows using the structures closest to the desired value along the reaction coordinate. Harmonic potentials with spring constants were used to restrain the seed configurations, where is the vector defining the reaction coordinate. For technical reasons, inverted flooding potentials were employed to approximate (Lange et al., 2006). Each of the restrained structures was simulated for 120–200 ns, and the first five ns in each simulation were discarded as equilibration. The free energy profile was reconstructed using the WHAM method (Kumar et al., 1992) implemented in the GROMACS g_wham tool, and its uncertainty was assessed by Bayesian bootstrapping of complete probability histograms (Hub et al., 2010).
As a further control, a conceptually unrelated method proposed by (Zhu and Hummer, 2012) was used to crosscheck the uncertainty of the free energy profiles. Unlike in Bayesian bootstrapping, where new sets of histograms are generated based on the sampling in each of the umbrella windows, this method relies on the statistical error of the mean force in every individual window. The mean force errors, in turn, are obtained from block averages, with the optimal block size calculated as , where is the total duration of the trajectory and is the autocorrelation time of in window .
The state of the guanosine nucleotide allosterically affects the interfaces of tubulin in protofilamentJournal of Computer-Aided Molecular Design 26:397–407.https://doi.org/10.1007/s10822-012-9566-x
Structural transitions at microtubule ends correlate with their dynamic properties in Xenopus egg extractsThe Journal of Cell Biology 149:767–774.https://doi.org/10.1083/jcb.149.4.767
Molecular dynamics with coupling to an external bathThe Journal of Chemical Physics 81:3684–3690.https://doi.org/10.1063/1.448118
Dynamic instability 30 years later: complexities in microtubule growth and catastropheMolecular Biology of the Cell 26:1207–1210.https://doi.org/10.1091/mbc.E13-10-0594
Structure of growing microtubule ends: two-dimensional sheets close into tubes at variable ratesThe Journal of Cell Biology 129:1311–1328.https://doi.org/10.1083/jcb.129.5.1311
An extended sampling of the configurational space of HPr from E. coliProteins: Structure, Function, and Genetics 26:314–322.https://doi.org/10.1002/(SICI)1097-0134(199611)26:3<314::AID-PROT7>3.0.CO;2-D
Nucleotide-dependent lateral and longitudinal interactions in microtubulesJournal of Molecular Biology 425:2232–2246.https://doi.org/10.1016/j.jmb.2013.03.029
GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular SimulationJournal of Chemical Theory and Computation 4:435–447.https://doi.org/10.1021/ct700301q
g_wham—A Free Weighted Histogram Analysis Implementation Including Robust Error and Autocorrelation EstimatesJournal of Chemical Theory and Computation 6:3713–3720.https://doi.org/10.1021/ct100494z
MicrotubulesNew York: Wiley-Liss.
Tubulin bond energies and microtubule biomechanics determined from nanoindentation in silicoJournal of the American Chemical Society 136:17036–17045.https://doi.org/10.1021/ja506385p
The weighted histogram analysis method for free-energy calculations on biomolecules. I. The methodJournal of Computational Chemistry 13:1011–1021.https://doi.org/10.1002/jcc.540130812
Flooding in GROMACS: accelerated barrier crossings in molecular dynamicsJournal of Computational Chemistry 27:1693–1702.https://doi.org/10.1002/jcc.20473
Refined structure of alpha beta-tubulin at 3.5 A resolutionJournal of Molecular Biology 313:1045–1057.https://doi.org/10.1006/jmbi.2001.5077
Microtubule stability in budding yeast: characterization and dosage suppression of a benomyl-dependent tubulin mutantMolecular Biology of the Cell 6:1241–1259.https://doi.org/10.1091/mbc.6.9.1241
All-atom empirical potential for molecular modeling and dynamics studies of proteinsThe Journal of Physical Chemistry B 102:3586–3616.https://doi.org/10.1021/jp973084f
Microtubule dynamics and microtubule caps: a time-resolved cryo-electron microscopy studyThe Journal of Cell Biology 114:977–991.https://doi.org/10.1083/jcb.114.5.977
The mechanisms of microtubule catastrophe and rescue: implications from analysis of a dimer-scale computational modelMolecular Biology of the Cell 23:642–656.https://doi.org/10.1091/mbc.E11-08-0688
The determinants that govern microtubule assembly from the atomic structure of GTP-tubulinJournal of Molecular Biology 412:35–42.https://doi.org/10.1016/j.jmb.2011.07.029
Structural intermediates in microtubule assembly and disassembly: how and why?Current Opinion in Cell Biology 18:179–184.https://doi.org/10.1016/j.ceb.2006.02.009
Polymorphic transitions in single crystals: A new molecular dynamics methodJournal of Applied Physics 52:7182–7190.https://doi.org/10.1063/1.328693
UCSF Chimera--a visualization system for exploratory research and analysisJournal of Computational Chemistry 25:1605–1612.https://doi.org/10.1002/jcc.20084
Structural Ensembles of Intrinsically Disordered Proteins Depend Strongly on Force Field: A Comparison to ExperimentJournal of Chemical Theory and Computation 11:5513–5524.https://doi.org/10.1021/acs.jctc.5b00736
Tackling Exascale Software Challenges in MolecularDynamics Simulations with GROMACSEASC 2014 Conference Proceeding. pp. 3–27.
Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella samplingJournal of Computational Physics 23:187–199.https://doi.org/10.1016/0021-9991(77)90121-8
GMCT : a Monte Carlo simulation package for macromolecular receptorsJournal of Computational Chemistry 33:887–900.https://doi.org/10.1002/jcc.22919
Mechanochemical model of microtubule structure and self-assembly kineticsBiophysical Journal 89:2911–2926.https://doi.org/10.1529/biophysj.105.060913
Conformational Changes in Ligand Binding ProcessesUniversity of Göttingen.
Dynamic instability of individual microtubules analyzed by video light microscopy: rate constants and transition frequenciesThe Journal of Cell Biology 107:1437–1448.https://doi.org/10.1083/jcb.107.4.1437
Molecular and Mechanical Causes of Microtubule Catastrophe and AgingBiophysical Journal 109:2574–2591.https://doi.org/10.1016/j.bpj.2015.10.048
Convergence and error estimation in free energy calculations using the weighted histogram analysis methodJournal of Computational Chemistry 33:453–465.https://doi.org/10.1002/jcc.21989
Luke M RiceReviewing Editor; UT Southwestern Medical Cetner, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Microtubule assembly governed by tubulin allosteric gain in flexibility and lattice induced fit" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Guest Reviewing Editor and Anna Akhmanova as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: William O Hancock (Reviewer #2); Gregory A Voth (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
The manuscript from Igaev and Grubmuller uses molecular dynamics (MD) simulations to investigate whether/how the identity of the bound nucleotide (GTP vs GDP) affects the conformation of ab-tubulin. This is an interesting and challenging question that has not been definitively settled, and where new insights can impact the understanding of microtubule dynamics and regulation. The authors perform atomistic molecular dynamic (MD) simulations of GTP/GDP-tubulin dimers starting from different conformations to probe the conformational dynamics and free energy profiles of unassembled tubulin dimers. They find that both GTP- and GDP-tubulin have an ensemble average of intermediate bending conformation, and that GTP-tubulin has much larger bending flexibility compared to GDP-tubulin. Igaev and Grubmuller conclude that increased bending flexibility and reduced free energy of kinking make GTP-tubulin more favorable for incorporation into the MT lattice. Finally, the authors also use functional mode analysis to identify two distinct bending modes for tubulin – twist-bending (TB) and splay-bending (SB) – are identified through functional mode analysis in both nucleotide states. Igaev and Grubmuller speculate that the transition from SB state to TB state when tubulin is bound at the assembling MT tip generates destabilizing strain that could contribute to the kinetics of the MT disassembly process. The ideas that nucleotide state controls tubulin flexibility, and that nucleotide-dependent differences in flexibility may be important for microtubule dynamics, are interesting and new. All reviewers felt that the work was technically performed at a very high level, and that a suitably revised version of the paper that satisfactorily addresses the concerns articulated below has the potential to make a nice contribution to the field.
1) While the work was very well presented, there was some concern that as written the manuscript might not be accessible to a broad audience. Additional efforts at simplification and clarification would be extremely helpful. 'Conformational mode' or 'reaction coordinate' axes are not always easy to relate to the known structures even given the black dots in Figure 3, Figure 4, and Figure 5 that show where the experimental structures fall. Figure 4—figure supplement 1 was a good example of helpful simplification of what are otherwise complicated concepts explained with too much technical jargon.
2) The results from the simulations, upon which all the analysis rests, are hard to assess: whether convergence has been reached and/or whether different equivalent simulations yield overlapping results, and how some of the conformations observed relate to the known structures is not easily discerned from what was presented. This should be addressable through additional panels (possibly supplemental) made from existing data: convergence or lack thereof could be illustrated using color coding or shading overlaid onto the pink and green distributions of Figure 3 and Figure 4 to show simulation time and/or to separate the equivalent simulations, and by annotating the rotation angle between a- and b-tubulin in Figure 3, Figure 4, and Figure 5, and possibly also by proving reference structures in the accompanying supplemental videos (which are helpful). Alternately, or in addition, it would be useful to show the extent of bending transition (conformational mode 1) along the simulation time as the evidence of the simulations are converged. Would the differences between Video 2 and Video 3 be better perceived in a side-by-side presentation?
3) The manuscript is in part framed around the idea that it is testing competing models for how nucleotide affects tubulin conformation: an 'allosteric' model positing that GTP-binding causes tubulin to become straighter, and a 'lattice' model positing that tubulin binding to the GTP lattice is what causes tubulin to become straighter. Strictly speaking, because the authors did not perform simulations of tubulin assemblies, their data may not really have much to say about the lattice model. Furthermore, some of the simulations start from 'straight' conformations that were taken from recent cryo-EM structures of microtubules in different nucleotide states. It seems that these simulations do not converge onto the curved conformation of unpolymerized tubulin (and conversely, that simulations starting from curved tubulin do not reach straight conformations). It is interesting that there appears to be such a large difference in free energy between curved and straight conformations, but at the same time the lack of connection between the two sets of simulations also leaves room for doubt about whether something is missing (or being missed). To address these issues surrounding the 'straight' simulations, the authors are requested to more explicitly acknowledge (i) the potential limitations of starting simulations from straight conformations of tubulin 'plucked' from the lattice (i.e. removed from its natural context where the conformation is stabilized by longitudinal and lateral interactions in the lattice), and (ii) the potential limitations of not having simulated assemblies of tubulin. Along these lines, the reviewers recommend removing or at least substantially condensing the penultimate paragraph (and associated discussion in Materials and methods section). It seemed overly detailed for how speculative it was, given the limitations articulated here along with a lack of simulations of polymerization dynamics.
4) In Figure 7B, the uncertainty of the free energy profile is large at the two local minimums and nearly zero at the transition state around 2.5 nm. This seems counterintuitive. The authors should also evaluate the uncertainty using the block average method and compare with the results presented here. In Figure 8, the transition from SB state to TB state when GTP dimer has arrived at the MT tip requires an activation energy of around 20 kT, assuming the free energy profile in Figure 7B is correct. This seems to be much too high considering the MT assembly speed to be around 2 μm/min, or 50 tubulin dimers/sec. Could the authors justify the model using existing experimental data?https://doi.org/10.7554/eLife.34353.040
- Maxim Igaev
- Helmut Grubmüller
- Maxim Igaev
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
This work was supported by the Max Planck Society and the German Research Foundation (DFG) via the research grant No. IG 109/1–1. The authors gratefully acknowledge computer time provided by the North-German Supercomputing Alliance (HLRN). The authors also thank Gregory Bubnis, Nicholas Leioatts, Thomas Ullmann, and Andrea Vaiana for helpful discussions.
- Luke M Rice, Reviewing Editor, UT Southwestern Medical Cetner, United States
© 2018, Igaev et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.