Mechanics and kinetics of dynamic instability

  1. Thomas CT Michaels
  2. Shuo Feng
  3. Haiyi Liang  Is a corresponding author
  4. L Mahadevan  Is a corresponding author
  1. Paulson School of Engineering and Applied Sciences, Harvard University, United States
  2. Department of Modern Mechanics, University of Science and Technology of China, China
  3. IAT Chungu Joint Laboratory for Additive Manufacturing, Anhui Chungu 3D Institute of Intelligent Equipment and Industrial Technology, China
  4. Department of Physics, Harvard University, United States
  5. Department of Organismic and Evolutionary Biology, Harvard University, United States

Abstract

During dynamic instability, self-assembling microtubules (MTs) stochastically alternate between phases of growth and shrinkage. This process is driven by the presence of two distinct states of MT subunits, GTP- and GDP-bound tubulin dimers, that have different structural properties. Here, we use a combination of analysis and computer simulations to study the mechanical and kinetic regulation of dynamic instability in three-dimensional (3D) self-assembling MTs. Our model quantifies how the 3D structure and kinetics of the distinct states of tubulin dimers determine the mechanical stability of MTs. We further show that dynamic instability is influenced by the presence of quenched disorder in the state of the tubulin subunit as reflected in the fraction of non-hydrolysed tubulin. Our results connect the 3D geometry, kinetics and statistical mechanics of these tubular assemblies within a single framework, and may be applicable to other self-assembled systems where these same processes are at play.

Introduction

Microtubules (MTs) are polar tubular polymers formed by the self-assembly of the protein tubulin. MTs are ubiquitous in eukaryotic cells, where they are a major component of the cellular cytoskeleton, and participate in a number of essential cellular functions, such as cell migration, morphogenesis, transport within cells and cell division (Desai, 1997; Gupta et al., 2015; Huber et al., 2013; Fletcher and Mullins, 2010). They are also involved in regulating the shape and dynamics of axons, cilia and flagella (Mimori-Kiyosue, 2011).

The basic building blocks of MTs are tubulin heterodimers. These are formed by α-tubulin and β-tubulin, two structurally similar globular proteins with mass of about 55 kDa. αβ-tubulin dimers are arranged longitudinally into flexible tubulin filaments called protofilaments (PFs). A number (between 9 and 16, typically 13) of such PFs then assembles by lateral interactions to form the MT lattice (Mandelkow et al., 1986; Chrétien and Wade, 1991; Mitchison, 1993; Tilney et al., 1973). MT growth occurs by the addition of tubulin dimers mainly at the plus end, where β-tubulin is exposed. Upon hydrolysis of guanosine-tri-phosphate (GTP), tubulin subunits undergo a structural conversion that weakens lateral bonds, destabilises the subunit in the MT lattice and converts the relatively straight tubulin state into a state that is bound to guanosine-di-phosphate (GDP) and is characterised by an increased longitudinal curvature (Wang and Nogales, 2005; Alushin et al., 2014).

MTs are not static assemblies. They can repeatedly and stochastically vary their length by undergoing alternating phases of assembly and disassembly both in vivo and in vitro. This phenomenon is termed ‘dynamic instability’ and it is essential to a number of cellular functions, such as chromosome separation, the remodelling of spatial organisation of the cytoskeleton during mitosis or the exploration of extracellular environment (Mitchison and Kirschner, 1984; Gildersleeve et al., 1992; Cassimeris et al., 1988; Sammak and Borisy, 1988). Understanding the factors that regulate MT dynamic instability is central to cell physiology and disease. Yet a detailed understanding of dynamic instability still remains elusive (Aher and Akhmanova, 2018; Hemmat et al., 2018). This difficulty originates in part from the fact that dynamic instability is the result of several mechanical and kinetic aspects operating at multiple time and length scales (Figure 1).

Connecting mechanical and kinetic aspects of MT dynamic instability across multiple time and length scales.

MTs are composed of tubulin heterodimers formed of α-tubulin and β-tubulin. Both α-tubulin and β-tubulin are bound to GTP, but only the GTP that is bound to a β-tubulin is hydrolysable (Nogales et al., 1998). When the tubulin dimer is part of the MT lattice, GTP hydrolysis increases the spontaneous longitudinal curvature along the dimer axis. This causes GDP-tubulin dimers to be less tightly bound in the MT lattice (Wang and Nogales, 2005; Alushin et al., 2014). In the scheme, GTP-tubulin dimers are shown in green, while GDP-tubulin dimers are shown in blue; dark color indicates β-tubulin and light colour indicates α-tubulin. Tubulin dimers are connected head-to-tail into PFs. Typically, 13 such PFs align laterally to form the MT lattice, which is a long and hollow cylindrical shell with an outer diameter of approximately 25 nm and a thickness of about 5 nm (Mandelkow et al., 1986; Chrétien and Wade, 1991; Mitchison, 1993; Tilney et al., 1973). The mechanical stability/instability of the MT tube structure results from a competition between lateral and longitudinal curvatures. While MTs made of GTP-tubulin are relatively straight (about 5° per subunit), MTs made of GDP-tubulin tend to curve outward longitudinally at the plus end due to the longitudinal spontaneous curvature of GDP-tubulin dimers (about 12° per subunit) (Chrétien et al., 1995). Consequently, MTs consisting of GDP-tubulin (hydrolyzed MTs) tend to be mechanically less stable that their non-hydrolysed counterparts. The images show: (a) Schematic structure of an unhydrolyzed αβ-tubulin dimer with bound nucleotides highlighted in orange (Alushin et al., 2014). (b) EM image showing the characteristic shape of a depolymerising MT plus end, resembling a ram’s horn (VanBuren et al., 2005). (c) Sequence of TIRF microscopy images of a MT end illustrating the switching between phases of polymerisation and depolymerisation during dynamic instability; catastrophe is associated with the loss of the GTP-caps (green fluorescence is from a protein that is believed to associate with the GTP cap). (a) is adapted from Alushin et al. (2014); (b) is reproduced from Austin et al. (2005). (c) is reproduced from Duellberg et al. (2016) under the Creative Commons Attribution License CC-BY-4.0 (https://creativecommons.org/licenses/by/4.0/).

© 2005 Company of Biologists Ltd. All rights reserved. Panel B is reproduced from Austin et al. (2005) with permission. It is not covered by the CC-BY 4.0 licence and further reproduction of this panel would need permission from the copyright holder.

As such, dynamic instability of MTs has been the focus of extensive experimental and theoretical work (Mandelkow et al., 1986; Chrétien and Wade, 1991; Mitchison, 1993; Tilney et al., 1973; Mitchison and Kirschner, 1984; Fygenson et al., 1994; Gildersleeve et al., 1992; Cassimeris et al., 1988; Sammak and Borisy, 1988; Hemmat et al., 2018; Aher and Akhmanova, 2018; Dogterom and Leibler, 1993; Brun et al., 2009; Antal et al., 2007; Mahadevan and Mitchison, 2005; Nogales et al., 1998; Wang and Nogales, 2005; Alushin et al., 2014; Zovko et al., 2008; Hoog et al., 2011; Mandelkow et al., 1991; Chrétien et al., 1995; Austin et al., 2005; Gardner et al., 2011; Hunyadi and Jánosi, 2007; Rice et al., 2008; Kueh and Mitchison, 2009; Jánosi et al., 1998; Molodtsov et al., 2005; VanBuren et al., 2005; Zapperi and Mahadevan, 2011; Bertalan et al., 2014b; Wu et al., 2009; Duellberg et al., 2016; Seetapun et al., 2012; Cheng et al., 2012; Cheng and Stevens, 2014; Jain et al., 2015; Aparna et al., 2017; Zakharov et al., 2015). Studies initially invoked a kinetic view capturing MT dynamic instability phenomenologically by different rates of polymerisation and depolymerisation depending on the state of the tubulin-phosphate complex (Dogterom and Leibler, 1993; Brun et al., 2009; Antal et al., 2007). Evidence of changes in the structural properties of MT subunits during hydrolysis later showed how lattice-bound tubulin dimers undergo a structural conformational change that increases their curvature (Wang and Nogales, 2005; Alushin et al., 2014; Mahadevan and Mitchison, 2005), suggesting that outward curving tips of hydrolysed MTs can cause such structures to be mechanically unstable (Zovko et al., 2008; Hoog et al., 2011; Austin et al., 2005; Mandelkow et al., 1991; Chrétien et al., 1995; Gardner et al., 2011; Hunyadi and Jánosi, 2007; Rice et al., 2008; Kueh and Mitchison, 2009). Several coarse-grained computer simulations have since then adopted this structural-mechanical view, considering MT elasticity explicitly, to understand different aspects of dynamic instability, including hydrolysis-driven mechanical deformations near the cap (Jánosi et al., 1998), force generation by shrinking microtubules (Molodtsov et al., 2005), or 3D sheet-like/blunt tips (VanBuren et al., 2005), or stochastic microtubule tip configurations and their relation to catastrophe (Zakharov et al., 2015). Few studies, however, have attempted to capture this mechanical view with the aim of emphasizing the qualitative features necessary for dynamic instability and providing phase diagrams that delineate the zones where dynamic instability is seen. Those that exist, for example (Zapperi and Mahadevan, 2011; Bertalan et al., 2014b), focus on the 1D limit, where MTs are modelled as adsorbed chains: the predictions from these models can be qualitatively different from those that correctly account for the 3D geometry of MTs (see ‘Mechanical stability of 3D MTs in the presence of quenched disorder’).

Here, we use a combination of theory and simulations to establish how the kinetics of polymerization and the mechanics associated with the 3D geometry of MTs act together across multiple scales to regulate dynamic instability. We also establish the role of disordered remnants of GDP-tubulin in determining the statistics of MT rescue. All together, our study provides a set of qualitative phase diagrams that delineate the regions of parameter space where dynamic instability is seen, consistent with previous observations while providing experimentally testable predictions.

Methods

Computational model

To complement the theory (see ‘A phase diagram for mechanical stability of 3D MTs’), we developed a minimal coarse-grained computational model of MT mechanics and dynamics. To characterize the tubulin heterodimers, we use two patchy spheres linked together by a flexible hinge. We derive the patchy particles from coarse-grained representations of colloidal particles that are decorated by patches on their surface; the patches represent specific anisotropic interactions that promote binding with patches on other particles. In our model, interactions between dimers are described by patches carrying two types of interactions: longitudinal and lateral contacts (Huisman et al., 2008). Longitudinal contacts link dimers head-to-tail, arranging them into PFs. Lateral contacts connect parallel PFs to form the cylindrical shell of the MT. Each dimer has two longitudinal contacts and four lateral ones, corresponding to three patches per monomer (Figure 2a). Interactions between two patches on monomers i and j are described by the following potential (Feng and Liang, 2012):

(1) Vpatchy(r,θi,θj)=VS(r,θi,θj)+VB(r,θi,θj)+VT(r,θi,θj).

We see that there are three distinct contributions associated with the stretching (VS), bending (VB) and twisting (VT) modes, and we choose the following forms for these potentials:

(2) VS(r,θi,θj)={ϵ[(1ea(rr0))21]r<rmϵ[(1ea(rr0))21]Dθ(θi)Dθ(θj)rmr<rc0rrc
(3) VB(r,θi,θj)=b(θi2+θj2)Dr(r)Dθ(θi)Dθ(θj)
(4) VT(r,θi,θj)=c(ϕ/r)2Dr(r)Dθ(θi)Dθ(θj)

where

(5) Dr(r)={1sin4(π4rr0rdr0)r<rd0rrd
(6) Dθ(θ)={1sin4(π4θθd)θ<θd0θθd

Here, r denotes the center-to-center distance between monomers, while the angles θi and θj describe the spatial directions of the patches (Figure 2b). The Morse potential term VS, defined in Equation 2, describes the non-covalent interaction between patches, with ϵ being the depth of the potential well, r0 the equilibrium distance between monomers, and a is a parameter that controls the curvature of the potential well and, hence, determines the stretching modulus (see Equation 7). When r<rm, where rm=r0-log(2)/a, the potential VS behaves as an isotropic repulsive interaction. In the range rmr<rc, where rc=5r0 is the cutoff for VS (VS is set to zero for rrc), VS is modified by multipliers Dθ(θi) and Dθ(θj). This yields an anisotropic attractive potential that exists only when the patches are aligned. Indeed, the multipliers Dθ(θi) and Dθ(θj), which are defined in Equation 6, weaken the attraction between the patches when these are not aligned (Figure 2c): Dθ reaches its maximum value when θ=0. The cutoff of Dθ is θd=π/3 and limits the influence of the patches within particular range of spatial directions. The potential terms VB and VT, defined in Equations 3 and 4, characterise bending and twisting deformations respectively. They are described as classical harmonic potentials with curvature b, respectively, c, and are modified by the multipliers Dr(r), Dθ(θi) and Dθ(θj), which are defined in Equations 5 and 6. These multipliers limit the range of VB and VT to specific spatial locations and directions. The cutoff of Dr is set as rd=2.7r0, which is smaller than rc (the cutoff of VS). This choice makes VB and VT shorter-range interactions compared to VS.

Definition of computer simulation and mechanical model.

(a) The GTP-tubulin subunit (αβ-heterodimer) in our coarse-grained computer model. The blue patches (1 and 2) are longitudinal connecting points, while the pink patches (3 to 6) are lateral ones. This setup allows the dimer to curve inward along the MT longitudinal direction and outward laterally. The positions of the patches allow us to control the longitudinal and lateral curvatures of the subunit, giving a strong curvature (12°) to hydrolysed dimers and a weaker curvature (5°) to non-hydrolysed dimers. (b) The interaction potential Vpatchy(r;θi,θj) between two patchy particles as a function of their centre-to-centre distance r for θi=0 and increasing θj=θ. Inset: the bending angles θi and θj are defined as the angles between the centre-to-patch line and the centre-to-centre line; the twisting angle φij is defined by the projection of the normals 𝒏i and 𝒏j on the plane perpendicular to the centre-to-centre line. (c) Our mechanical model describes a MT as a continuous, thin elastic sheet. This is parameterised as a surface of revolution obtained by rotating the function R(x)=R0+u(x) along the MT long axis (x-axis), where R0 is the natural radius of the MT and u(x) describes the local deformation of the surface away from R0. We approximate lateral interaction potentials between tubulin subunits in neighbouring PFs by a set of spring potentials (springs with stiffness S). Bending of the PFs in the MT causes the connecting springs to stretch. This competition between bending energy and stretching energy determines the mechanical stability of MTs. Inset: hydrolysed tubulin dimers curve naturally in two directions as described by the longitudinal and lateral spontaneous curvatures κ and ω, respectively.

The parameters in our coarse-grained computational model are linked to the mesoscopic mechanical properties of MTs (see section ’Mechanics’) including the interfilament spring stiffness S, the filament bending stiffness B and the filament torsional rigidity K as

(7) S=2a2lϵ,B=bl,K=8cl,

where l is the length scale of tubulin dimers. l takes different values depending on whether we are calculating longitudinal or lateral properties. In particular, we set l=2r0=8 nm when calculating longitudinal properties and set l=r0=4 nm for lateral ones (see Appendix 1 and Table 1 for further details on the computer simulation model and a summary of parameter choices). We choose the potential-well parameters in our simulations such that the resulting mesoscopic mechanical parameters Equation 7 are consistent with the experimentally measured values of the mechanical properties of typical MTs (Gittes et al., 1993; Mickey and Howard, 1995; Felgner et al., 1996; Tolomeo and Holley, 1997; de Pablo et al., 2003; Sept and MacKintosh (2010); Deriu et al., 2010). Simulations of this coarse-grained model of MTs were performed using Molecular Dynamics (MD), as described in Appendix 1.

Table 1
Values of parameters used in our coarse-grained simulations.
ParameterValueDescription
r04 nmEquilibrium distance between tubulin monomers
ϵ07.665 × 10-20JPotential energy depth for longitudinal interactions
ξ5° or 12°Angle between α-tubulin and β-tubulin within a GTP- or GDP-dimer
χchainSee Equations 24-29
χside1103.45° or 90°See Equations 24-29
χside213.8°See Equations 24-29
llong2r0Longitudinal size of tubulin dimers
llatr0Lateral size of tubulin dimers
along,alat20/r0Parameters that control longitudinal and lateral stretching stiffness
blong20ϵ0Parameter that controls longitudinal bending stiffness
blat10ϵ0Parameter that controls lateral bending stiffness
clong500ϵ0r02Parameter controls longitudinal twisting stiffness
clat0.5ϵ0r02Parameter that controls lateral twisting stiffness
ϵlongϵ0Potential energy depth for longitudinal interactions
ϵlat0.0906ε0Potential energy depth for lateral interactions

Results

A phase diagram for mechanical stability of 3D MTs

Our analytical model to study the mechanical stability of 3D MTs is a function of the underlying parameters describing MT mechanics, MT growth kinetics and subunit hydrolysis, which we first consider in the deterministic limit (see ’A phase diagram for mechanical stability of 3D MTs’). This forms the basis for studying the mechanical stability of MTs when there is heterogeneity in the state of tubulin, that is it could be either GTP or GDP bound (see ‘Mechanical stability of 3D MTs in the presence of quenched disorder’), and allows us to investigate the role of GTP-remnants (containing random fractions of non-hydrolysed subunits) on rescue (see ‘Role of GTP-remnants in rescue’).

The starting point of our mechanical model is that MTs exist as individual polymers with persistence lengths in the O(mm) range, which is much larger than the typical length of MTs (µm range) (Fletcher and Mullins, 2010; Huber et al., 2013). This observation suggests that MTs can be considered to have a well-defined shape that is not affected significantly by thermal fluctuations. Moreover, previous studies indicate that lateral bonds between tubulin dimers are considerably weaker than longitudinal ones (VanBuren et al., 2002; VanBuren et al., 2005; Molodtsov et al., 2005). We model MTs as a set of adherent PFs that have bending stiffness B; we approximate lateral interactions between PFs by a series of spring potentials (springs of stiffness S) and assume that extending these springs beyond a critical displacement causes MTs to become mechanically unstable and break, potentially leading to dynamic instability.

Kinetics

Previous studies (Wang and Nogales, 2005; Alushin et al., 2014) suggest that tubulin dimers that are part of the MT lattice have different mechanical properties depending on their hydrolysis state. In particular, upon hydrolysis the tubulin dimer undergoes a structural transformation from a relatively straight state to a state with finite curvature. In a 3D setting, the tubulin dimer can be curved both in the longitudinal direction and in the lateral direction. Let κ(0) denote the longitudinal curvature of tubulin dimers in their GTP-state and let κ() be the longitudinal curvature in the hydrolysed state. For simplicity, we assume that the hydrolysis reaction affects primarily the longitudinal curvature of the tubulin dimers, such that their curvature ω in the azimuthal direction can be considered to be constant. This assumption can be relaxed, see Appendix 2. Thus, as a result of GTP-hydrolysis, the longitudinal curvature of tubulin dimers, κ, changes with time, which we assume follows first order kinetics so that

(8) dκdt=kH(κ()-κ),

where kH is the rate of hydrolysis. While bound tubulin changes its structure, unbound (bound) tubulin can attach (detach) to (from) the free end of the MTs, which we also describe using a minimal first order kinetic law for the evolution of the length of the MT n(t) (expressed in number of subunits) so that

(9) dndt=k+[m]-k-=kG.

Here kG is the net growth rate, k+ is the elongation rate constant, k- is the dissociation rate constant and, for simplicity, we have assumed a constant subunit concentration [m] in solution.

Mechanics

In addition to MT growth and subunit hydrolysis, we also need to account for the elastic deformation of the MTs since the geometric state of the assembly is linked to its mechanical state. Individual PFs can bend, but are also constrained by inter-filament interactions, so that there are two contributions to elastic energy: (i) curvature energy associated with the bending of PFs, (ii) stretching energy of the springs connecting neighbouring PFs. We capture these energy contributions in a continuum picture that describes a MT as a thin elastic surface of revolution obtained by rotating the function R(x)=R0+u(x) along the long MT axis (x-axis), where R(x) is the local radius of MT and R0 is the natural radius (Figure 2c). In the small gradient approximation, corresponding to u1, the total elastic energy can be written as (see Appendix 2 for details):

(10) tot=0(B2[u′′-κ]2+Σ2[u]2+S2u2)𝑑x,

where Σ=B(1+ωR0)2/R02 and =/x denotes derivative with respect to x. The first term in Equation 10 is the energy of MT that penalises deviations from its natural curvature κ(t) which itself reflects its state of hydrolysis. The second term is a surface energy term that penalises area increase due to the outward curving of the surface. The third term is the stretching energy of the springs. The minimum energy configuration results from a competition between bending energy, which favours a natural curved MT state, and elastic spring energy, which favours a straight cylindrical MT configuration. The overall shape of the axisymmetric tubule is then obtained by solving the Euler-Lagrange equation associated with Equation 10 (see Appendix 2 for details):

(11) Bu′′′′-Σu′′+Su=0

subject to the boundary conditions u(,t)=u(,t)=0 (fixed minus end), u′′(0,t)=κ and u′′′(0,t)=0 (free plus end), and is coupled to the kinetic Equations 8 and 9.

Condition for mechanical stability

There are three natural dimensionless parameters (two mechanical parameters and one kinetic parameter) in our model that read:

(12) α=(Bκ2Suc2)1/2,β=(ΣucBκ)1/2,γ=kGkH.

The first parameter α describes the effect of longitudinal curvature. The second parameter β pertains to lateral curvature. Coupling these mechanical parameters to the kinetics of subunit hydrolysis and MT growth introduces an additional relevant dimensionless parameter γ, which is the ratio of the rate of hydrolysis of GTP-tubulin dimers to the net rate of addition of GTP-subunits to the MT plus end (see Appendix 2 for details).

These parameters serve as the basis for a phase diagram for the mechanical stability of a 3D MT (Figure 3a). Assuming that a mechanical instability arises when the elastic MT is deformed so that the radial displacement crosses a critical value (u(0)>uc), we can solve Equation 11 in terms of the maximal deformation u(0) to yield a condition for when the MT is mechanically unstable. Rewriting this in terms of the scaled longitudinal curvature yields a critical value (see the Appendix 2 for details):

(13) α=1+4β21-e-n/γ-12β2.

above which MTs are mechanically unstable (the transition curve in the αβ-plane in Figure 3a). This critical value depends on the lateral curvature parameter β, which is related to MT radius through β1/R0. In particular, the critical value for α is maximal (α=1) when β=0, that is R0. This situation corresponds to the limit of a one-dimensional MT (Zapperi and Mahadevan, 2011). The critical value for α then decreases with increasing lateral curvature β, that is decreasing MT radius. Overall, these results suggest that MTs with smaller radius are mechanically less stable than MTs with larger radius, consistent with the intuition that increasing azimuthal curvature increases the mechanical strain on the MT and thus makes it more likely to fracture. Furthermore, increasing the rate of MT growth over hydrolysis acts to stabilise MTs mechanically. In the mechanically stable phase, hydrolysis is slower than the addition of GTP-tubulin at the plus end, leaving a stabilising GTP-cap of size n (see Appendix 2, Video 1). In the mechanically unstable phase, hydrolysis is faster than subunit addition at the plus end. Consequently, the ‘hydrolysis front’ takes over the ‘growth front’, which destabilises MTs. In this case, the PFs curve outward near the plus end, leading to a characteristic morphology of depolymerising MTs that resembles rams’ horns (VanBuren et al., 2005) (see Appendix 2, Video 2). The transition curve separating these fast and slow hydrolysis regimes depends on both the longitudinal and lateral curvatures of the MT.

Role of lateral and longitudinal curvatures in mechanical stability of a 3D MT in the absence and presence of quenched disorder.

(a) Schematic phase diagram for mechanical stability/instability with three axes (two mechanical axes and one kinetic axis): longitudinal curvature parameter α=(Bκ2/Suc2)1/2, lateral curvature parameter β=(Σuc/Bκ)1/2, and ratio of growth rate to hydrolysis rate γ=τκ/τG. (b) 2D phase diagram separating regions of mechanically stable and unstable MTs as a function of α and β. Data points are from computer simulations and the solid line is the prediction of Equation 13. (c)-(d) Implementation of lateral (c) and longitudinal (d) disorder at the level of the spring stiffness S in our coarse-grained computer model. The data points are simulation results and the solid lines are fits to Equation 16 and Equation 17. The data show that lateral disorder acts to destabilise MTs mechanically (c), while longitudinal disorder strengthens MTs (d). In both cases, disorder was generated using the distribution of spring constants in Equation 14 with σ=0.02S, that is the disorder parameter was k=50. The amount disorder of each tubulin subunit represents the average value of the relative stiffness S/S of its two lateral interactions. See Videos 14 for movies illustrating the mechanical stability of a MT with and without quenched disorder (see Videos 14).

Video 1
MT growth and hydrolysis kinetics.

This movie shows the interplay between MT growth kinetics (we focus here only on the addition of GTP-tubulin subunits at the plus end) and the subsequent hydrolysis of the incorporated subunits in the older parts of the MT. The interplay between growth and hydrolysis can be seen in the emergence of a growing front and a hydrolysis front that move with different speed. In this case, the hydrolysis front does not catch up with the growth front at the plus end, yielding a stabilising GTP-cap and a mechanically stable MT configuration for the entire duration of the simulation. The growth rate is rG=10-4 steps−1, while the probability of hydrolysis for each dimer is given by pH=e-rHt, where rH=10-10 steps−1 is the rate of hydrolysis. Note the separation of timescales between MT growth and subunit hydrolysis (rGrH). To keep the helical structure stable, we use the following parameters blong=8ϵ0, blat=4ϵ0 and clong=5ϵ0r02, while the other mechanical parameters are the same as listed in Table 1

Video 2
MT catastrophe.

This movie shows the interplay between MT growth and subunit hydrolysis. The rate of hydrolysis is rH=5×10-6 steps−1 while the growth rate is decreasing from rG=5×10-3 steps−1 to rG=2×10-6 steps−1. Hydrolysis destroys the GTP-cap of MT and causes catastrophe. The mechanical parameters are the same as in Video 1

In Figure 3b, we show that crossing the transition curve given by Equation 13 causes a switch from the mechanically stable phase into the mechanically unstable phase. This could result from variations in either α, β or γ. For MTs in the mechanically stable phase, catastrophic failure can still occur via thermal activation. In this regime, the rate of catastrophe follows Arrhenius’ law rcexp(ΔE/kBT), where kBT is the thermal energy and ΔE is the energy barrier given by ΔEα1+β2 is a measure of the ‘distance’ from the transition curve Equation 13 in the phase diagram of Figure 3b; the further away a MT is from this transition curve, the less likely it is to undergo catastrophe. The dependence of the rate of catastrophe on MT radius is through the parameter β, such that lnrc1/R02. Thus, at constant temperature and at fixed values of the mechanical parameters, the rate of catastrophe increases with decreasing MT radius or, equivalently, decreasing PF number.

Comparison with computer simulations

We used our computer simulations to test the prediction from Equation 13 for how the critical value of α varies with β, which is a function of MT radius R0 that is controlled by changing the number Nf of PFs in the MT. The results (Figure 3b) show that the critical value for α is maximal when β=0, and decreases with increasing β, in agreement with the theoretical prediction of Equation 13 (solid line).

Mechanical stability of 3D MTs in the presence of quenched disorder

Having considered the deterministic limit where growth kinetics, hydrolysis and longitudinal/lateral curvatures characterise the mechanical stability of a 3D MT, we now consider the role of randomness by including a random fraction of GTP-tubulin dimers in their lattice. We can model this situation by introducing quenched disorder in the state of the tubulin subunit. Quenched disorder describes the general situation when certain parameters in the system become random variables; disorder can be considered to be ‘quenched’ when the probability distribution of parameter values either does not vary with time or it varies with time slowly compared to some underlying fast dynamics, and thus cannot be described solely using equilibrium statistical mechanics. In the context of MTs, such a separation of timescales emerges very naturally when comparing fast polymerization/depolymerization kinetics and the comparatively slower GTP turn-over.

Since the primary mode of MT instability is due to the breaking of lateral bonds, a natural parameter for discussing the role of disorder is the spring stiffness S. The underlying motivation for this choice is that mechanical forces can influence the rates of chemical reactions. Indeed, mechanical work contributes to the free energy, which in turn determines the rates of a chemical reaction (Howard, 2001). In our context, GDP-tubulins within the MT lattice experience mechanical stresses due to the strong curvature. These mechanical stresses can shift the polymerisation-depolymerisation equilibrium and favour free monomers. This is consistent with the idea that GDP-tubulins are less tightly bound to MTs than GTP-tubulins (Wang and Nogales, 2005; Alushin et al., 2014), even if the chemical bonds are identical.

Thus, instead of having a well-defined spring constant S throughout the MT, we consider a MT with varying S. Each lateral interaction is characterised in principle by a different spring constant S, which is drawn from a time–independent probability distribution p(S) of spring constants. For convenience, we choose the Gamma-distribution

(14) p(S;k,S)=kkSk-1exp(-kS/S)SkΓ(k),

where Γ(x) is the Gamma function, S is the average spring stiffness, and the parameter 1/k=σ/S, with σ being the standard deviation of the distribution, is the coefficient of variation that describes the degree of disorder in the system. The choice of the Gamma-distribution admits a simple parameterisation in terms of the coefficient of variation that allows us to explore a range of different extreme value statistics. For instance, for k=1 the Gamma distribution p(S;1,S) yields the exponential distribution with intensity λ=1/S, while for k1 it yields a normal distribution with mean μ=kS and variance σ2=kS2. We distinguish two limiting modes for disorder: lateral (Figure 3c) and longitudinal (Figure 3d). Any realisation of disorder can then be decomposed into a combination of these two limiting modes.

Lateral disorder

In the presence of disorder in S, lateral interactions are characterised by variations in S azimuthally but not longitudinally. In this case, whether the MT undergoes catastrophe depends on the breaking of the weakest lateral bond along the circumference of the MT, resulting in a MT that is ‘cut-open’ along the longitudinal direction (Figure 3c and Appendix 2, Video 3). This situation is fully analogous to what happens when pulling a one-dimensional chain by its ends: the chain will break as soon as its weakest link breaks. The mechanical stability of a MT with lateral disorder is thus equivalent to the mechanical stability of a MT with uniform spring stiffness Smin, where Smin denotes the average value of the weakest spring stiffness along the MT circumference. This replacement maps the study of the mechanical stability of a MT with lateral disorder onto a problem of extreme value statistics (Zapperi and Mahadevan, 2011; Bertalan et al., 2014b): the determination of Smin for a system of N independent and identically distributed links with spring constants S1,,SN. In Appendix 2, we show that Smin can be calculated from Equation 14 using extreme-value statistics, yielding:

(15) SminS=1kΓ(k+1k)[Γ(k+1)N]1/k.

The condition for mechanical instability of a MT with lateral disorder in S is thus obtained by replacing S by Smin in Equation 13, yielding:

(16) α=1+4β21-e-n/γ-12β2SminS.

Since Smin<S, Equation 16 predicts that the transition curve between mechanically stable and unstable MT regions shifts towards the instability region. Thus, lateral disorder weakens MTs.

Video 3
Role of lateral quenched disorder in MT mechanical stability.

This movie illustrates the stability of a MT in the presence of lateral disorder in the spring stiffness S, which describes the strength of lateral contacts between PFs. Note that lateral disorder causes the MT to be mechanically unstable along directions with weakest lateral bonds (the MT ‘cut opens’ along these weak directions). The color code indicates the local value of S for each subunit, which is obtained by averaging over the relative stiffness of its two lateral bonds (drawn from the Gamma-distribution (14) with k=50). The simulation parameters are the same as listed in Table 1, except that we set ξ=5.8 for all tubulin dimers.

Longitudinal disorder

A different situation arises when quenched disorder is distributed longitudinally. Here, the strongest lateral bond determines MT stability. In fact, longitudinal disorder leads to the presence of ‘rings’ of particularly strong bonds that prevent the MT from depolymerising completely (Figure 3d and Appendix 2, Video 4). The mechanical stability of a MT with longitudinal disorder is thus equivalent to that of a MT with uniform spring stiffness Smax, where Smax is the expected value of S associated with the strongest lateral bond. Using extreme-value statistics, one finds (see Appendix 2) Smax/S=(γe+logN)/k, where γe0.5772 is the Euler-Mascheroni constant (Taloni et al., 2018; Bertalan et al., 2014b; Zapperi and Mahadevan, 2011). Hence, the curve separating mechanically stable and unstable regions in the presence of longitudinal disorder is Zapperi and Mahadevan (2011):

(17) α=1+4β21-e-n/γ-12β2SmaxS.

Since Smax>S, longitudinal disorder in S reinforces MTs.

Video 4
Role of longitudinal quenched disorder in MT mechanical stability.

This movie illustrates the stability of a MT in the presence of longitudinal disorder in the spring stiffness S. Note that the MT is in this case mechanically more stable than in the presence of lateral disorder. This is because longitudinal disorder leads to the presence of ‘rings’ of strong bonds that prevent the MT from peeling off completely. The simulation parameters are the same as in Video 3

It is important to note that longitudinal disorder is the only mode of disorder present in a 1D MT (Zapperi and Mahadevan, 2011). Lateral disorder is thus a defining feature of the 3D geometry of MTs. Our results thus reveal a fundamental role of MT dimensionality: while in a 1D setting quenched disorder stabilises MTs mechanically, in a 3D setting it can destabilise MTs.

Comparison with computer simulations

We have tested the theoretical predictions of Equation 16 and Equation 17 using our coarse-grained simulations (Figure 3c,d). Quenched disorder was realised using Equation 14 with disorder parameter k=50. These simulations confirm that longitudinal quenched disorder increases the mechanical stability of MTs (Figure 3c), whereas the effect of quenched disorder in the lateral direction is to destabilise MTs mechanically (Figure 3d).

Role of GTP-remnants in rescue

Using our theoretical model of mechanical stability of 3D MTs in the presence of quenched disorder, we are now in the position to investigate the role of remnants of GTP-tubulin in rescue. Rescue refers to the transition from depolymerisation to polymerisation during MT dynamic instability but is still poorly understood (Brouhard, 2015). Experiments indicate that GTP-tubulin addition at the plus end is not critical for rescue (Gardner et al., 2013; Walker et al., 1988), but that the presence of remnants of GTP-tubulin along the MT lattice in so-called ‘GTP-islands’ can lead to MT rescue (Tropini et al., 2012; Dimitrov et al., 2008; Aumeier et al., 2016; Gardner et al., 2013; Vemu et al., 2018). In particular, experiments in vivo have revealed a strong correlation between rescue probability and the presence of remnants of GTP-tubulin in older parts of the MTs, suggesting that these ’GTP-remnants’ can function as rescue sites (Dimitrov et al., 2008). This view was further supported by the observation that the presence of a slowly hydrolyzable analogue of GTP bound to tubulin subunits contributes to MT rescue (Tropini et al., 2012. Aumeier et al., 2016) also demonstrated the possibility to generate GTP-islands along the MT lattice in a controlled manner by means of laser damaging and subsequent repair of the damaged site by incorporation of GTP-tubulin from solution (Schaedel et al., 2015): rescue occurred at laser-damaged sites in the presence of free GTP-tubulin (Aumeier et al., 2016). Separately, recent studies (Vemu et al., 2018; Vemu et al., 2019) reported of a damage-repair mechanism that stabilises MTs mediated by the enzymes spastin and katanin. Overall, these studies suggest that disordered GTP-islands in an otherwise structurally periodic lattice are involved in rescue regulation. Since these GTP-remnants are characterised by a random mixture of different states of tubulin, we ask if our framework might help to quantify these observations.

Computer simulations

We first used our coarse-grained simulations to study the role of disordered GTP-remnants in MT rescue. We generated reinforcing islands by inserting, in the middle of a fully hydrolysed, depolymerising MT, a ring consisting of several layers of GTP-tubulin dimers (Figure 4a). We then observed whether the reinforcing GTP-islands were able rescue the depolymerising MTs as a function of two parameters: 1) the length of the GTP-island Nrf (defined here as the number of layers in the island) and 2) the fraction ϕ of GTP-tubulin in the island. The results of these simulations (Video 5) are shown in Figure 4b. Note that the parameter ϕ controls the amount of disorder present in the island at the level of GTP-hydrolysis. This mimics both the scenario when rescue islands are formed because not all GTP-tubulin is able to hydrolyse, as suggested in Dimitrov et al. (2008), or when rescue islands result from the incorporation of GTP-tubulin during the repair process of a damaged site, as suggested in Aumeier et al. (2016). If disorder varies slowly over time compared to the characteristic timescale of polymerisation/depolymerisation, we can model slow changes of MT mechanical properties by making the relevant mechanical parameters explicit functions of time. For the parameters in our simulation, we find that when the reinforcing island is one layer long (Nrf=1), the probability of rescue is close to zero, irrespective of the GTP-fraction in the reinforcing island. Interestingly, when Nrf>1, we observe that rescue probability prescue increases with ϕ in a highly nonlinear manner. Specifically, prescue is either close to zero or close to one for most values of ϕ, with a sharp increase in the transition region.

Role of disordered GTP-islands in rescue.

(a) Schematics showing the rescue of a depolymerising MT reinforced by a disordered GTP-tubulin island. In this example, the length of the reinforcing island is Nrf=3 and the GTP-island consists of 20 GTP-tubulin dimers, which corresponds to a GTP-fraction of ϕ=20/(3×13)0.51 (see Video 5). (b) Simulated rescue probability by GTP-islands as a function of GTP-fraction ϕ and reinforcing island length Nrf. The rescue probability prescue was estimated by repeating simulations six to ten times for each pair ϕ and Nrf. (c) Percolation model for the role of disordered reinforcing GTP-islands on MT rescue. (d) Simulated site percolation on a square lattice with dimensions Nrf×13 as a function of GTP-fraction ϕ and island length Nrf. The probability of percolation was obtained by averaging over 103 realisations for each pair ϕ and Nrf.

Video 5
MT rescue by disordered GTP-island.

This movie illustrates the successful rescue of a depolymerising MT by a disordered GTP-island placed along the MT lattice. The mechanical parameters in this case are the same as in Video 1.

Percolation model of rescue

To qualitatively understand the observed nonlinear behaviour of rescue probability with GTP-fraction in the reinforcing island ϕ, we propose a site percolation model of rescue (Figure 4c). Site percolation is concerned with the following question: given a random graph, in which each site is (independently) occupied with probability q or empty with probability 1-q, what is the probability that a connected path of occupied sites exists between the boundaries of the graph? In our percolation model of rescue, each site of the reinforcing island is occupied by a GTP-tubulin dimer with probability ϕ, while it is occupied by a GDP-tubulin dimer with probability 1-ϕ. In Sec. ’Role of quenched disorder in mechanical stability of MTs’, we have shown that the presence of randomly distributed weak lateral bonds (mediated by GDP-tubulin) can destabilise MTs mechanically when disorder is longitudinal. As such, a MT will be mechanically unstable when a connecting path of GDP-tubulins runs longitudinally through the reinforcing island (Figure 4c). The question of whether a reinforcing island with GTP-fraction ϕ is able to rescue a depolymerising MT is thus analogous to site percolation with q=ϕ. The rescue probability prescue thus relates to 1-pperc, where pperc is the probability of percolation of a longitudinal path of GDP-tubulin subunits through the length of the reinforcing island. Figure 4d shows that the results of site percolation on a square lattice of dimensions 13×Nrf with varying ϕ are in qualitative agreement with simulated rescue probabilities (Figure 4b).

Discussion

Our multi-scale approach to dynamic instability incorporates the mechanics and 3D geometry of MTs, the kinetics of tubulin addition and GTP-hydrolysis, and quenched disorder in the state of the tubulin subunit. Our results provide a series of phase diagrams for the presence of dynamic instability, revealing the dimensionless mechanical and kinetic parameters controlling the problem. Compared to previous analytic studies of dynamic instability, our results reveal the key role of the 3D geometry of MTs. In particular, we find that the mechanical stability of a MT is strongly affected by its radius (Mahadevan and Mitchison, 2005); a MT with a smaller radius is mechanically less stable than an identical MT with a larger radius.

Since MT radius has been shown to vary because the number Nf of PFs composing MTs typically ranges between 9 to 16 (Chrétien et al., 1992; Chalfie and Thomson, 1982; Chaaban and Brouhard, 2017; Cueva et al., 2012; Díaz et al., 1998), our study suggests quantitative experimental tests via measurements of catastrophe rates rc as a function of PF number Nf to verify the theoretical prediction for the rate of catastrophe log(rc)1/Nf2. Another key prediction from our study is that the rescuing power of GTP-islands displays a sharp drop at intermediate values of GTP-fraction. In particular, the percolation model predicts that there is a critical point for the GTP-fraction, ϕ=ϕc, below which reinforcing islands lose their ability to rescue MT disassembly. The numerical value of the threshold ϕc depends on the thickness of the reinforcing island as well as on the MT lattice structure. This critical GTP-fraction could be determined experimentally and compared to theory by using non-hydrolyzable analogs of tubulin (Tropini et al., 2012), to control the amount of disorder at the level of the state of tubulin subunit in the island, in combination with super-resolution microscopy (Huang et al., 2009) to establish island length. Finally, we note that our assumption of the form of the quenched disorder in terms of the Gamma distribution is just that - an assumption. Further experimental work will be required to solve the inverse problem of estimating the average GTP-fraction in GTP remnants from rescue probabilities determined experimentally (Dimitrov et al., 2008; Aumeier et al., 2016; Vemu et al., 2018) to see if we might determine both the form of the disorder and the intrinsic parameters characterising it, thus allowing future research to address the question of how to control dynamic instability.

Appendix 1

Computer model

We assign a local coordinate system to each dimer to describe its overall rotation, as well as the relative position of the patches within the dimer. To describe the whole dimer, we need three basis vectors t,s and n is directed along the long dimer axis and points from the 𝜶-tubulin to the 𝜷-tubulin. 𝒔 and 𝒏 lie in the plane perpendicular to 𝒕. This coordinate system is particularly convenient for describing the rotation of a dimer as a whole. However, it cannot capture the relative rotation of the 𝜶- and 𝜷-tubulin monomers in the dimer. This inner rotation can only be along the axis 𝒏, and thus two patches coordinate systems (one for each tubulin in the dimer) can be obtained by rotating the local coordinate system around 𝒏 as:

(18) tβ=tcos(ξ/2)ssin(ξ/2),
(19) sβ=tsin(ξ/2)+scos(ξ/2),
(20) nβ=n,
(21) tα=tcos(ξ/2)+ssin(ξ/2),
(22) sα=tcos(ξ/2)+ssin(ξ/2),
(23) nα=n,

where ξ is the angle between the vectors 𝒕𝜶 and 𝒕𝜷 (the axes of the two monomers), as shown in Appendix 1—figure 1a.

Appendix 1—figure 1
Computer model.

(a) The local coordinate systems of a dimer. (b) Top, left and front views of two spheres forming a dimer and the six patchy points.

There are two types of connections between dimers: longitudinal and lateral contacts. Longitudinal contacts link dimers head to tail and arrange them into PFs. Lateral contacts connect parallel PFs arranging them into a cylindrical MT shell. Each dimer has two longitudinal contacts and four lateral ones, that is there are six patchy points on each dimer. These are numbered as in Appendix 1—figure 1b or Figure 2a of the main text. The coordinates of the six patchy points are described by the following equations:

(24) R1=tβcos(χchain)sβsin(χchain),
(25) R2=tαcos(χchain)sαsin(χchain),
(26) R3=(tβcos(χside1)nβsin(χside1))cos(χside2)+sβsin(χside2),
(27) R4=(tβcos(χside1)+nβsin(χside1))cos(χside2)+sβsin(χside2),
(28) R5=(tαcos(χside1)nαsin(χside1))cos(χside2)+sαsin(χside2),
(29) R6=(tαcos(χside1)+nαsin(χside1))cos(χside2)+sαsin(χside2).

The values of the angles χchain,χside1,χside2 are given in Table 1.

Appendix 1—figure 2
Mesoscopic mechanical properties of MTs as functions of the parameters in the coarse-grained simulations model.

(a) Young’s modulus E vs a, (b) torsional rigidity K vs c and (c) bending stiffness B vs b. Solid lines correspond to the theoretical predictions (see Equation 7 in Section ’Computational model’ of main text), and are in good agreement with the computer simulations.

Appendix 2

Theoretical model

Here we provide additional mathematical details pertaining to the theoretical model discussed in the main text, in terms of the elastostatics of the structure and kinetics of polymerisation/depolymerisation and hydrolysis. Finally, we study this model in the presence of quenched disorder.

Mechanical stability of 3D MTs (statics)

Our model considers two contributions to the total elastic energy of MTs: 1) a curvature energy term describing the bending of PFs in the MT and 2) the stretching energy of the springs that connect neighboring PFs. We now consider these two energy contributions in detail.

Curvature energy

We describe the curvature energy term by means of the Helfrich (1973):

(30) curvature=𝑑A[k12(2H)2+k2K],

where k1 and k2 are bending rigidities, H=(κ1+κ2)/2 is the mean curvature, K=κ1κ2 is the Gaussian curvature, κ1 and κ2 are the principal curvatures, and dA is the surface area element. We describe the MT as a continuous elastic sheet modelled as a solid of rotation obtained by rotating the function R(x) along the x-axis, which describes the microtubule growth axis (Appendix 2—figure 1a). We assume that the MT possesses a natural radius R0 (when flat) and, for convenience, we write R(x)=R0+u(x), where u(x) describes the deviation of local radius R(x) from R0. The MT surface is thus parameterised as (x,R(x)cosφ,R(x)sinφ) where x[0,) and φ[0,2π). In this parametrisation, x=0 corresponds to the MT plus end, while x= corresponds to the minus end. We then use the following expressions for the mean and Gaussian curvatures of a surface of revolution Gray, 1997:

(31) H(x)=R(x)R′′(x)-(1+R(x)2)2R(x)[1+R(x)2]3/2,K(x)=-R′′(x)R(x)[1+R(x)2]2,

where =d/dx denotes derivative with respect to x, and the surface area element is dA=2πR(x)1+R(x)2dx. From these expressions, we note that the Gaussian curvature term, K(x)dA=-R′′(x)/[1+R(x)2]3/2dx, is fully integrable, that is it gives rise to boundary terms only. We thus focus on the mean curvature term only. To this end, we consider a small gradient approximation, which corresponds to u(x)1. We then write R(x)=R0+u(x) in Equation 31 and after keeping only the leading order terms, we find

(32) (2H)2dA2π(R0u′′(x)2+u(x)22R0)dx+,

where stands for higher order terms in u(x) or for boundary terms (i.e. terms that are fully integrable and thus, after integration, simply shift the curvature energy by a constant value).

Appendix 2—figure 1
Continuum elastic shell model of microtubule.

(a) To describe the curvature energy, MTs are modelled as thin elastic surfaces of revolution obtained by rotating the function R(x)=R0+u(x) along the MT long axis (x axis). R0 is the natural radius of the MT. (b) Neighbouring PFs within the MT are connected by Nf Hookean springs with stiffness ks.

Our calculation so far has accounted for bending energy relative to a flat cylindrical MT. To account for the natural longitudinal and lateral curvatures of subunits (κ and ω, respectively), we extend the Helfrich Hamiltonian as follows:

(33) curvature=𝑑A[k12(2H-2H0)2+k2K],

where H0=(κ+ω)/2 is the natural (or spontaneous) mean curvature. In the small gradient approximation (u1), the leading contribution to the curvature energy in the presence of natural curvatures κ and ω is found to be

(34) (2H-2H0)2dA2π(R0[u′′(x)-κ]2+u(x)22(1+ωR0)2R0)dx+.

Thus, the curvature energy is

(35) curvature0(B2[u′′(x)-κ]2+Σ2u(x)2)𝑑x,

where B=2πR0k1 and Σ=B(1+ωR0)2/R02.

Elastic spring energy

The second contribution to the total elastic energy of MTs is due to the stretching energy of the springs connecting neighbouring PFs. To construct this energy contribution, we model the cross section of the MT as a set of PFs connected by Nf harmonic springs with stiffness ks (Appendix 2—figure 1b). The stretching energy is proportional to the square of the deviations of the distance d between PFs from their equilibrium position d0. By symmetry in the φ direction, springs are stretched by d(x)θNR(x)=2πR(x)/Nf (θN is defined in Appendix 2—figure 1b). The rest length of springs is d0=θNR0=2πR0/Nf. Thus, the stretching energy per unit length of the system of Nf springs can be estimated as:

(36) Nf×ks2[d(x)-d0]2Nfks2(2πNf)2[R(x)-R0]2.

Thus, writing R(x)=R0+u(x), the stretching energy is found to be:

(37) stretchS20u(x)2𝑑x,

where S=4π2ks/Nf.

Total elastic energy

In summary, the total elastic energy of a 3D MT can be written in terms of u(x) by combining the contributions from curvature energy and spring energy:

(38) tot=curvature+stretch0(B2[u′′(x)-κ]2+Σ2u(x)2+S2u(x)2)𝑑x.

The first term describes the longitudinal bending energy of the MT away from the natural curvature κ. The second term comes from the lateral curvature. At leading order, this term is a surface energy term that penalises the increase in surface area when the MT curves out; the parameter Σ plays the role of a surface tension. Finally, the third term is the stretching energy of the springs, which gives rise to an energy density contribution proportional to u(x)2.

Euler-Lagrange equation and minimum energy configuration

Having defined the elastic energy functional for the problem, tot[u]=0(u,u,u′′)𝑑x (Equation 38), we now consider the associated Euler-Lagrange equation that describes the minimal energy configuration of MTs:

(39) δtot[u]δu=u-ddx(u)+ddx2(u′′)=0.

Using Equation 38, Equation 39 is found to be:

(40) Bu′′′′-Σu′′+Su=0.

Subject to the following boundary conditions

(41) u()=u()=0 (MTflatatx=,i.e.atMTminusend),
(42) u(0)=κ (naturalcurvatureatx=0,i.e.atMTplusend),
(43) u(0)=0 (noshearforceintheMT)

the analytical solution to Equation 40 reads:

(44) u(x)=κλ2λ12(λ2-λ1)e-λ1x+κλ1λ22(λ1-λ2)e-λ2x,

where

(45) λ1,2=1α2±1α4-4β42,

and the two relevant length scales in the problem are

(46) α=(BΣ)1/2,β=(BS)1/4.

The first length scale is associated with the lateral curvature of the MT. In fact, from Σ=B(1+ωR0)2/R02 it follows α=(B/Σ)1/2R0.

Condition for mechanical stability/instability of 3D MTs

To study the mechanical stability/instability of MTs as a function of its mechanical parameters, we assume that the springs connecting PFs break if their extension exceeds a critical value uc. From Equation 44 we thus obtain the following condition for mechanical instability:

(47) u(0)=κβ2(1+β2α2)>uc.

To rewrite this in a more transparent way in terms of the dimensionless parameters, we first define the scaled longitudinal and lateral curvatures:

(48) α=Bκ2Suc2,β=ΣucBκ.

so that the parameters α and β are related to the characteristic length scales defined in Equation 46 through

(49) α=(βc)2,β=cα,

where c=uc/κ is a length scale corresponding to the geometric average of the longitudinal radius of curvature and the critical extension uc.

The condition Equation (47) can be reformulated most conveniently as

(50) β>1-ααα=1+4β2-12β2,

The resulting curve is shown in a phase diagram in Appendix 2—figure 4a.

Appendix 2—figure 2
Role of kinetics on mechanical stabillity of MTs.

(a) Phase diagram for mechanical stability of 3D MTs, as described by Equation 50. (b) Effect of the relative rates of tubulin hydrolysis and MT growth (γ=kG/kH) on mechanical stability of MTs, as described by Equation 62. See: link to Videos 1 and 2.

Scaling behaviour of critical curvature with mechanical parameters

Appendix 2—figure 3
Scaling behaviour of critical curvature.

(a) Scaling behavior of the critical longitudinal curvature κc with bending stiffness B and spring stiffness S (inset). The scaling relationship over the range of B and S values investigated in the simulations is found to be κc(B/S)-0.53, in agreement with the predictions of Equation 55. (b) Scaling exponent σ (Equation 55) as a function of β/α interpolates between the limits σ=-0.5 for β/α0 and σ=-1 for β/α.

Equation (50) predicts that there is a critical value for the longitudinal curvature, κc, above which MTs are mechanically unstable. We now study the scaling behaviour of κc with system parameters such as the bending stiffness B and spring stiffness S. To this end, we performed a series of computer simulations by varying the longitudinal bending stiffness B of PFs and the lateral connecting strength ϵ, which in turn determines the spring stiffness S. The results are shown in Appendix 2—figure 3a. We find that κc follows a scaling law with the longitudinal bending stiffness B and S, κc(B/S)σ. This scaling behaviour is shown in the double-logarithmic plots in Appendix 2—figure 3a, where the slope corresponds to the scaling exponent σ-0.53. The observed scaling behaviour can be rationalised using (50), which can be rewritten as

(51) κc(B/S)σ,

where σ is the scaling exponent, given by

(52) σ=log(κc)log(B/S).

After rewriting Equation 47 as

(53) κc=ucβ2(1+β2α2)

we find

(54) σ=log(B/S)[log(uc)-log(β2)-log(1+β2α2)]=-12-1(1+β2α2)log(B/S)(β2α2),

where in the second step we have used β=(B/S)1/4 (Equation 46). Finally, using xa/logx=axa, we arrive at:

(55) σ=-12(1+11+α2β2)={-0.5 when βα-1 when βα

The scaling exponent σ interpolates between −0.5 for βα (limit of a single-filament model of MT Zapperi and Mahadevan, 2011) and −1 for βα (limit of strong lateral curvature), see Appendix 2—figure 4b. This scaling behaviour for κc has been verified using our coarse grained simulations of depolymerising MTs in Appendix 2—figure 3a of the main text. Representative values for the mechanical parameters in the simulations are B=1.23×10-26 Nm2, S=91 MPa and R=10 nm. These values give β=(B/S)1/43.4 nm, α=10 nm and, therefore, σ-0.55, in close agreement with the simulations in Appendix 2—figure 4a.

Coupling mechanical stability with MT kinetics

We now combine our static calculation of mechanical stability of MTs with the dynamics of subunit hydrolysis and MT polymerisation/depolymerisation. This allows us study how kinetics affect the phase diagram of Appendix 2—figure 3a. Let kG be the rate of addition of subunits to the MT end and let kH be the rate of hydrolysis of GTP-tubulin dimers. With these parameters we can construct a dynamic dimensionless parameter as the ratio between the times for hydrolysis and growth

(56) γ=kGkH.

To understand how the dynamic parameter γ modifies the static stability of the MTs, our starting point are the kinetic equations in the main text:

(57) dκ(t)dt=kH(κ()-κ(t)),dω(t)dt=kH(ω()-ω(t)),
(58) dn(t)dt=k+[m]-k-=kG.

Equations 57 and 58 can be combined together by expressing time as t=n/kG (follows directly from the solution to Equation 58). The solution to Equation 57 can thus be expressed as

(59) κ(n)=κ(1-e-kHn/kG)+κ(0)e-kHn/kG=κ()(1-e-n/γ)+κ(0)e-n/γ,
(60) ω(n)=ω()(1-e-n/γ)+ω(0)e-n/γ.

Using this parametrisation, we can express the condition for mechanical instability, Equation 47, as

(61) (1-e-n/γ)α[1+αβ2(1+ω0R1-ω0Re-n/γ)2]>1.

This condition can be solved with respect to α to yield the following expression for the phase boundary:

(62) α=1+4β2(1+ω0R1-ω0Re-n/γ)21-e-n/γ-12β2(1+ω0R1-ω0Re-n/γ)2.

The resulting phase diagram as a function of γ is shown in Appendix 2—figure 4b. We see that increasing γ increases the region of mechanical stability of MTs. The physical interpretation of this result is that MTs remain mechanically stable as long as the transition curve, Equation 47, is reached after a stabilising cap of length n is added onto the MT end. Increasing the rate of MT growth, kG, over subunit hydrolysis, kH, favours the formation of the stabilising cap. In the limit when ω does not vary with time, Equation 61 reduces to

(63) (1-e-n/γ)α(1+αβ2)>1,

which can be solved to yield Equation 13 of the main text.

Mechanical stability of MTs in the presence of laterally-distributed quenched disorder

As argued in the main text, in the presence of lateral disorder, the mechanical stability of a MT is controlled by its weakest link. To understand this quantitatively, we replace the MT with lateral disorder with an equivalent one with uniform spring stiffness Smin, where Smin denotes the average value of the weakest spring stiffness. This replacement maps the study of the mechanical stability of a MT with lateral disorder onto a problem of extreme value statistics (Zapperi and Mahadevan, 2011; Bertalan et al., 2014a): the determination of Smin. To this end, consider N independent and identically distributed links with spring constants S1,,SN. We assume that the values of spring constants are random and distributed according to a Gamma-distribution

(64) p(S)=kkSk-1exp(-kS/S)SkΓ(k),

where S is the average spring stiffness (Appendix 2—figure 4a). The ratio between the standard deviation σ and the average S

(65) Cv=σS=1k

is the coefficient of variation, a key parameter which we use to describe the degree of disorder in the system; small values of k correspond to nearly ordered system while large values correspond to a strongly disordered system.

With Smin=miniSi being the smallest value of spring constants and letting

(66) PN(S)=Pr[Smin<S]

be the cumulative probability distribution for Smin over the N links, we can calculate PN(S) directly from the definition of PN(S), yielding

(67) 1PN(S)=Pr[SminS]=i=1NPr[SiS]=(1P(S))N,

where

(68) P(S)=0Sp(S)dS

is the cumulative probability distribution of S. For large N, we can approximate the exact expression in Equation 67 as

(69) PN(S)=1(1P(S))N1exp[NP(S)].

The interesting behaviour of PN(S) happens for small values of S,when PN(S) is controlled by the low-value tail of P(S). The cumulative distribution function P(S) can then be calculated explicitly from Equation 64 as

(70) P(S)=γ(k,kS/S)Γ(k),

where γ(k,x) is the incomplete gamma function. Since we are interested in the low-value tail of P(S), we can use the small x expansion of γ(k,x) 

(71) γ(k,x)xkk+,x0

to arrive at the following expression

(72) P(S)=γ(k,kS/S)Γ(k)(kS/S)kΓ(k+1),

which is valid for S0 (Appendix 2—figure 4b). The cumulative probability distribution for Smin then converges to a Weibull distribution

(73) PN(Smin)=1-exp[-N(Smin/S0)k], with S0=Γ(k+1)1/kSk.

The average value for the weakest spring constant Smin is therefore

(74) SminS=1kΓ(k+1k)[Γ(k+1)N]1/k.

Having solved the extreme-value statistics problem of determining Smin allows us to estimate the transition curve separating regions of mechanical stability and instability for a MT exhibiting lateral disorder. The curve separating these zones is obtained by replacing S by Smin in Equation 50, yielding:

(75) α=1+4β2-12β2SminS,

where we note that Smin<S. Thus, by comparing Equation 75 with the deterministic result Equation 50, we see that, in the presence of lateral disorder, the curve separating mechanically stable from mechanically unstable MTs shifts in such a way as to increase the region of mechanical instability. Lateral disorder thus weakens MTs, making them more susceptible to catastrophic failure.

Appendix 2—figure 4
Gamma distribution.

(a) Gamma distribution Equation 64 with k=5. (b) Cumulative probability distribution P(S) for the Gamma distribution, Equation 70, and low value-tail, Equation 72 (dashed line).

Data availability

All data generated or analysed during this study are included in the manuscript and supporting files.

References

    1. Felgner H
    2. Frank R
    3. Schliwa M
    (1996)
    Flexural rigidity of microtubules measured with the use of optical tweezers
    Journal of Cell Science 109 ( Pt 2:509–516.
    1. Gildersleeve RF
    2. Cross AR
    3. Cullen KE
    4. Fagen AP
    5. Williams RC
    (1992)
    Microtubules grow and shorten at intrinsically variable rates
    The Journal of Biological Chemistry 267:7995–8006.
  1. Book
    1. Gray A
    (1997)
    Modern Differential Geometry of Curves and Surfaces with Mathematica
    CRC Press.
  2. Book
    1. Howard J
    (2001)
    Mechanics of Motor Proteins and the Cytoskeleton
    Sinauer Associates, Inc.

Article and author information

Author details

  1. Thomas CT Michaels

    Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, United States
    Contribution
    Conceptualization, Formal analysis, Investigation, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6931-5041
  2. Shuo Feng

    1. Department of Modern Mechanics, University of Science and Technology of China, Hefei, China
    2. IAT Chungu Joint Laboratory for Additive Manufacturing, Anhui Chungu 3D Institute of Intelligent Equipment and Industrial Technology, Wuhu, China
    Contribution
    Software, Validation, Visualization, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0530-8381
  3. Haiyi Liang

    1. Department of Modern Mechanics, University of Science and Technology of China, Hefei, China
    2. IAT Chungu Joint Laboratory for Additive Manufacturing, Anhui Chungu 3D Institute of Intelligent Equipment and Industrial Technology, Wuhu, China
    Contribution
    Software, Supervision, Validation, Investigation, Methodology
    For correspondence
    hyliang@ustc.edu.cn
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7458-8036
  4. L Mahadevan

    1. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, United States
    2. Department of Physics, Harvard University, Cambridge, United States
    3. Department of Organismic and Evolutionary Biology, Harvard University, Cambridge, United States
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    lmahadev@g.harvard.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5114-0519

Funding

Swiss National Science Foundation

  • Thomas CT Michaels

National Natural Science Foundation of China (11272303)

  • Shuo Feng
  • Haiyi Liang

National Natural Science Foundation of China (11072230)

  • Shuo Feng
  • Haiyi Liang

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We acknowledge financial support from the Swiss National Science foundation (TCTM), the National Natural Science Foundation of China (11272303,11072230) (SF, LH).

Copyright

© 2020, Michaels 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.

Metrics

  • 4,098
    views
  • 538
    downloads
  • 25
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Thomas CT Michaels
  2. Shuo Feng
  3. Haiyi Liang
  4. L Mahadevan
(2020)
Mechanics and kinetics of dynamic instability
eLife 9:e54077.
https://doi.org/10.7554/eLife.54077

Share this article

https://doi.org/10.7554/eLife.54077

Further reading

    1. Computational and Systems Biology
    2. Physics of Living Systems
    Ju Kang, Shijie Zhang ... Xin Wang
    Research Article

    Explaining biodiversity is a fundamental issue in ecology. A long-standing puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individual-based modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.

    1. Computational and Systems Biology
    2. Physics of Living Systems
    Divyoj Singh, Sriram Ramaswamy ... Mohd Suhail Rizvi
    Research Article

    Planar cell polarity (PCP) – tissue-scale alignment of the direction of asymmetric localization of proteins at the cell-cell interface – is essential for embryonic development and physiological functions. Abnormalities in PCP can result in developmental imperfections, including neural tube closure defects and misaligned hair follicles. Decoding the mechanisms responsible for PCP establishment and maintenance remains a fundamental open question. While the roles of various molecules – broadly classified into “global” and “local” modules – have been well-studied, their necessity and sufficiency in explaining PCP and connecting their perturbations to experimentally observed patterns have not been examined. Here, we develop a minimal model that captures the proposed features of PCP establishment – a global tissue-level gradient and local asymmetric distribution of protein complexes. The proposed model suggests that while polarity can emerge without a gradient, the gradient not only acts as a global cue but also increases the robustness of PCP against stochastic perturbations. We also recapitulated and quantified the experimentally observed features of swirling patterns and domineering non-autonomy, using only three free model parameters - the rate of protein binding to membrane, the concentration of PCP proteins, and the gradient steepness. We explain how self-stabilizing asymmetric protein localizations in the presence of tissue-level gradient can lead to robust PCP patterns and reveal minimal design principles for a polarized system.