Chemomechanical Simulation of Microtubule Dynamics With Explicit Lateral Bond Dynamics

We introduce and parameterize a chemomechanical model of microtubule dynamics on the dimer level, which is based on the allosteric tubulin model and includes attachment, detachment and hydrolysis of tubulin dimers as well as stretching of lateral bonds, bending at longitudinal junctions, and the possibility of lateral bond rupture and formation. The model is computationally efficient such that we reach sufficiently long simulation times to observe repeated catastrophe and rescue events at realistic tubulin concentrations and hydrolysis rates, which allows us to deduce catastrophe and rescue rates. The chemomechanical model also allows us to gain insight into microscopic features of the GTP-tubulin cap structure and microscopic structural features triggering microtubule catastrophes and rescues. Dilution simulations show qualitative agreement with experiments. We also explore the consequences of a possible feedback of mechanical forces onto the hydrolysis process and the GTP-tubulin cap structure.


Introduction
Microtubule (MT) dynamics is essential for many cellular processes, such as the positioning and separation of chromosomes in mitosis [1], or maintenance of cell polarity and cell shape [2]. An important feature, which enables MTs to exert pulling and pushing forces in these cellular processes, is their dynamic instability, which is the stochastic switching of MTs between states of growth by polymerization and states of fast shrinkage by depolymerization [3].
Switching from growth into shrinkage happens in catastrophe events, whose mechanism and triggers are not completely understood on the molecular level, but they are associated with a loss of the GTP-cap by hydrolysis within the MT [4,5] (see Refs. [6,7] for reviews). Hydrolysis is strongly coupled to mechanics of the MT, as is clearly seen in the curling of MT protofilaments into a "ram's horn" conformation after the catastrophe and during the shrinking phase [8]. The loss of the stabilizing GTP-cap triggers a release of binding energy and stored mechanical energy in the tubular MT structure. Therefore, shrinkage following a catastrophe is more than simple depolymerization of the MT; it is rather a rupture or crack propagation process between protofilaments, which releases chemical and mechanical energy while it propagates towards the minus end. The energy released during shrinking has biological functions and can be employed to exert pulling forces onto kinetochores during separation of MTs in mitosis [9].
The curling of hydrolyzed protofilaments into a ram's horn structure shows that GDP-tubulin dimers have a bent conformation [8,[10][11][12]. Tubulin dimers assembled within the MT body are in a straight conformation, on the other hand [13]. Hydrolysis of tubulin dimers embedded in a straight MT causes mechanical strains in the tubular structure because the surrounding MT lattice prevents these GDP-tubulin dimers from assuming their preferred bent conformation. This mechanical strain is released in a catastrophe via the rupture of lateral bonds.
There are different models explaining how the mechanical strain is increased by hydrolysis or how lateral bonds are weakened by hydrolysis such that the strained MT becomes more prone for catastrophes. The first cryo-electron microscopy (EM) studies showed blunt tips for growing MTs but curved tips for shrinking MTs [8] suggesting that GTP-protofilaments are straight while GDP-protofilaments are curved. Later evidence from cryo-EM showed that GTP-protofilaments are also curved, but significantly less than GDP-protofilaments [10]. The allosteric model is based on the assumption that hydrolysis of a tubulin dimer changes the dimer conformation from a rather straight GTP-conformation to a bent GDP-conformation. Hydrolysis of tubulin dimers embedded in a straight MT causes mechanical strain in the tubular structure because the surrounding MT lattice prevents these GDP-tubulin dimers from assuming their preferred bent conformation. This model was employed in almost all previous MT simulation models that consider MT mechanics [14][15][16][17][18][19] The lattice model, on the other hand, is based on evidence from X-ray and cryo-EM structures [20][21][22][23] and simulations [24,25] that also GTP-tubulin dimers assume a bent conformation and that hydrolysis rather affects the lateral and longitudinal dimer interaction energies. It is supported by recent experimental observations that both growing and shrinking MTs have bent protofilament ends [26]. Ref. [26] also presents first simulation results with a lattice model. But there is also recent evidence from molecular dynamics (MD) simulation pointing in a different direction and supporting an intermediate model, where hydrolysis affects interactions but also lowers GDP-tubulin flexibility [27]. If hydrolysis weakens lateral interaction energies, hydrolysis makes the structure more prone for a catastrophe. While in the allosteric model, the mechanical strain in the structure is increased by hydrolysis, in the lattice model, the mechanical strain that the MT structure can tolerate is reduced by hydrolysis. In both models, the result is an increased propensity for lateral bonds to rupture. Therefore, chemomechanical MT models with explicit bond rupture are a necessity to reproduce catastrophes. We build on existing modelling approaches based on the allosteric model [14][15][16][17][18][19] and include lateral bond rupture as explicit stochastic events with force-dependent rates, which can give important clues about how catastrophes are triggered in the MT structure.
The influence of tubulin dimer hydrolysis onto the mechanics of the MT lattice suggests that, vice versa, mechanical forces and torques acting on tubulin dimers via strains in the tubular structure could also affect hydrolysis rates, an effect which has been explored only in Ref. [17] previously. Although this interplay is plausible from a mechanochemistry point of view, experimental verification on the dimer level is extremely difficult and not possible yet, but we can employ chemomechanical MT models to explore and suggest possible implications for the dynamic instability.
The coupling between chemical events -namely polymerization events, dimer hydrolysis, bond rupture -and mechanical forces because of conformational changes due to these chemical events, is a characteristic of MTs and requires chemomechanical MT models on the dimer level in order to develop a microscopic understanding of their dynamic instability including catastrophe and rescue events [28]. In this respect, chemomechanical models go beyond a phenomenological description of MT dynamics in a four-parameter model based on growth and shrinking velocities and phenomenological catastrophe and rescue rates [29]. The challenge for microscopic chemomechanical models is to include all chemical events as stochastic processes, to perform conformational relaxation governed by MT mechanics following each chemical event, and, eventually, to also include the feedback of mechanical forces within the MT onto reaction rates of the chemical events.
We present a stochastic chemomechanical MT model on the dimer level. Our model includes (i) a mechanical model of the MT containing lateral elastic bonds between tubulin monomers in neighboring protofilaments and a harmonic bending energy between tubulin monomers with a nonzero equilibrium angle after hydrolysis (allosteric model), (ii) stochastic addition and removal of tubulin dimers, (iii) explicit stochastic lateral bond rupture and bond formation; the bond rupture rate is coupled to the mechanical stress state of the bond and thus via elastic interactions within the MT lattice also to the other bonds, (iv) stochastic hydrolysis of dimers with a rate that can also couple to the mechanical bending stress in the dimer. The stochastic kinetics (ii)-(iv) is handled by a Gillespie algorithm and after each stochastic event, a mechanical energy minimization mimicking the relaxational dynamics of the structure is applied to the MT.
In order to parameterize our model, we will focus on the simplified scenarios of a growing MT consisting of GTP-tubulin only and a shrinking MT consisting of GDP-tubulin only. In both cases, we can neglect hydrolysis (iv); in the growing GTP-MT, we can also neglect mechanics, which is generated by hydrolysis. In the presence of mechanics and hydrolysis, repeated catastrophe and rescue events are obtained and will be described and analyzed. One problem in chemomechanical MT models is the computational effort associated with the mechanical relaxation. We investigate in detail, which level of computational effort is necessary in our model to obtain a sufficient mechanical relaxation following each chemical event, on the one hand, and which simplifications can be taken to assure a finite simulation time for growing MTs, on the other hand. This will allow us to simulate arbitrarily long growing MTs at fixed computational speed.
Our chemomechanical model has to be compared to previous modelling approaches, which include the mechanics of the MT [14][15][16][17][18][19]: • Refs. [15,16] employ the allosteric model for dimer bending and include stochastic addition and removal of dimers. Hydrolysis is random. Mechanical energy minimization is performed only locally on randomly selected dimers. Lateral bond rupture is not implemented as explicit stochastic process but only included using a threshold energy criterion.
• The models in Refs. [14,19] focus on mechanics and do not include dimer addition and removal. They are also based on the allosteric model but consider fixed hydrolysis states. In Ref. [14], the lateral bond energy landscape is harmonic around a minimum but includes an energy barrier and a dissociated, i.e., ruptured state. Global energy minimization gives the final state of the static structure.
• In Ref. [18], the stochastic kinetics is added to a mechanical model similar to [14]. Here, the mechanical relaxation and lateral bond rupture is performed using Brownian dynamics (which include thermal fluctuations) with small time steps (equivalent to 2 × 10 7 minimization steps), which is only applied to 300 tubulin dimers at the plus end. Stochastic addition of dimers and removal by rupture of lateral and longitudinal bonds is included. The rupture of lateral bonds happens by activation over the bond energy barrier, the longitudinal rupture by a threshold criterion. Hydrolysis is random and stochastic with a rate that is independent of mechanics.
• Ref. [17] is also based on the allosteric model. Lateral bond rupture is possible using a threshold criterion. Mechanical energy minimization was performed globally. There is no addition or removal of dimers, but hydrolysis is included. In a first attempt to include a coupling of the hydrolysis rate to mechanical forces, the hydrolysis kinetics remained deterministic, however, with the most probable hydrolysis event determined by mechanical forces. In the present paper, we will add addition and removal of dimers and a fully stochastic hydrolysis kinetics.
Our chemomechanical model has also to be compared to previous purely chemical modelling approaches on the dimer level but without explicit mechanical model [30][31][32][33][34]. These models include attachment and detachment of tubulin dimers; some of these models [32][33][34] also include lateral bond rupture and are thus able to produce crack-like catastrophe events. Crack-like catastrophe events are, however, triggered by adjusting chemical rupture rates rather than including MT mechanics. The model by Margolin et al [33] has successfully reproduced features of the experimentally observed MT dynamic instability [35] but relies on a heuristic tuning of simulation parameters.

Microtubule structure and energy
Our MT model is formulated on the dimer level. The base units of the model are alpha-and beta-tubulin monomers. In our model, we represent each monomer as cylinder with radius r t = 2 nm and height t = 4 nm (see Table 1). Alpha-and beta-tubulin monomers form unbreakable tubulin dimers, which are arranged head-totail into protofilaments. 13 protofilaments form a 13_3 MT, i.e., a MT with a helical shift of 3 tubulin monomer lengths per turn. or the remainder of this paper, we will use triples (p, d, t) to address specific tubulin monomers within the MT with p ∈ {1, 2, . . . , 13} as protofilament number, d ∈ {1, 2, . . . , d(p)} as tubulin layer (with d = 1 denoting the minus end and d = d(p) denoting the plus end of the protofilament p), and t ∈ {1, 2} denoting the tubulin monomer within the dimer with t = 1 for the alpha-and t = 2 for the beta-tubulin monomers. For simplicity, we assume periodicity in p (i.e., p = 0 ≡ 13 and p = 14 ≡ 1) and combined periodicity in d and t (i.e., (p, d, 3) ≡ (p, d + 1, 1) and (p, d, 0) ≡ (p, d − 1, 2)). We will also generally refer to the lateral neighbors of tubulin monomer (p, d, t) using (p ± 1, d, t) even though at the seam, lateral neighbors differ in all three indices.
The MT is straight and oriented along the z-axis with the positive z-direction pointing to the plus end. Vectors m(p, d, t) and p(p, d, t) point to the to the lower (minus end) and upper (plus end) circular base of the tubulin monomer (p, d, t). The direction vector with length t = 4 nm points from m(p, d, t) to p(p, d, t) and is specified using spherical coordinates, i.e., azimuthal and polar angles, see Figure 1(A). The polar angle θ(p, d, t) is the only degree of freedom of each monomer, because we assume that monomers can only be displaced in radial direction, i.e., all azimuthal angles are fixed to φ(p) = 2π(p − 1)/13. As both alpha-and beta-tubulin have their polar angles as a degree of freedom, the model supports intra-and inter-dimer curling [36].  At the minus end of the MT each protofilament p starts with an alpha-tubulin arranged in a circle with mean MT radius R MT = 10.5 nm and with an offset z(p, 1, 1) = 3 t (p − 1)/13 in z-direction, such that the seam is between the 13th and the 1st protofilament. The protofilament length that will be used to calculate the growth and shrinkage velocities is the maximum z-coordinate max (p) of all tubulin monomers within the protofilament (see Supplementary Material for more details). The MT length is given by the average Every tubulin monomer has four interaction points: two in longitudinal direction and two in lateral direction. The longitudinal bond between alpha-and beta-tubulin monomers of the same dimer is considered unbreakable but the orientation of this junction can change via the beta-tubulin's polar angle θ(p, d, 2). In contrast, the longitudinal bond between adjacent tubulin monomers of different dimers can break and is modeled via the bond energy ∆G 0 * long (where the "0" refers to it being a standard energy [37] and the asterisk to the fact that it also includes the entropic cost of "immobilization" [30]). The lateral interaction points are located at the edge of the upper base (see Figure 1(A)). If there is a lateral bond between tubulin monomer (p, d, t) and its neighbor in the (p + 1)-th protofilament, the bond is modeled as a harmonic spring with base energy ∆G 0 lat : with the spring constant k lat of the bond and the vector s(p, d, t) connecting the lateral interaction points; s 0 1.47 nm is the rest length of the spring (see [17] and also consider the helical shift between two neighboring tubulin monomers of 3 t /13). Lateral bonds at the seam are assumed to have identical mechanical properties as other lateral bonds based on evidence that they do not constitute a weaker bond [22,38]. Additionally, there is a lateral repulsion term between neighboring tubulin monomers (regardless of whether they are bonded or not) to ensure a cylindrical form [17]: The bending of monomer junctions is described by a harmonic potential with bending constant κ: The bending angle ∆θ(p, d, t) = θ(p, d, t) − θ(p, d, t − 1) (see Figure 1(B)) is calculated with the neighboring monomer in the minus direction (using the periodicity convention in d and t, (p, d, 0) ≡ (p, d − 1, 2)), and ∆θ 0 (p, d, t) is its equilibrium value. For hydrolyzed beta-tubulin monomers and for alpha-tubulin monomers on top of a hydrolyzed beta-tubulin (and for the first alpha-tubulin monomers of each protofilament if the beta-tubulin in the same dimer is hydrolyzed), we use a rest angle ∆θ 0 (p, d, t) = 11°in order to reproduce the experimentally measured radius of curvature of 21 nm corresponding to an angle of 22°per dimer for a GDPprotofilament curling into the ram's horn configuration [10,39]. Otherwise (for an unhydrolyzed beta-tubulin monomer or an alpha-tubulin monomer on top of an unhydrolyzed beta-tubulin monomer), we assume a straight equilibrium configuration with ∆θ 0 (p, d, t) = 0°. This choice of rest angles implements the allosteric model, where GTP-hydrolysis leads to bending of tubulin dimers.
Our mechanical MT model is defined by the total energy where E lat (p, d, t) only contributes if there is a lateral bond between tubulin monomers (p, d, t) and (p + 1, d, t) and E rep (p, d, t) only contributes if tubulin monomer (p, d, t) has a lateral partner (p + 1, d, t).
There are four free parameters in our mechanical MT model (see Table 2): the longitudinal bond energy ∆G 0 * long , the lateral bond energy ∆G 0 lat , the lateral spring constant k lat , and the bending constant κ. For the repulsion constant k rep , we use the same value k rep = 10 −6 rad 2 nm 12 κ that has been found previously to ensure the overall cylindrical shape of the MT and only contributes a small portion to the MT energy [17]. In the simulation model, we do not use this mechanical energy to calculate forces for a microscopic dynamics such as Brownian dynamics on the dimer level (as opposed to [18]). We rather assume that mechanical relaxation dynamics is fast compared to chemical changes in the MT due to tubulin attachment and detachment, bond rupture and formation, or hydrolysis. The slowest mechanical process is relaxation of bending modes of protofilaments governed by small restoring bending moments. The basic time scale for this process can be estimated as [40], where η ∼ 10 −3 Pa s is the viscosity of water. This gives τ ∼ 10 −10 s, which is orders of magnitude smaller than typical time scales of seconds for chemical events. Therefore, even longer protofilaments relax fast compared to chemical changes. There is additional evidence from Brownian dynamics that bending mode relaxation is also much faster than immobilization in cryo-EM [41]. Therefore, we perform a quasiinstantaneous energy minimization of (6) between these chemical simulation steps. This is the computationally more efficient strategy to achieve mechanical relaxation. The rates of all chemical simulation events themselves determine the dynamics of the MT and are handled by a Gillespie algorithm as explained in more detail below.

Chemical simulation events
To simulate the dynamics of MTs, we include attachment of individual GTP-tubulin dimers and detachment of (laterally unbonded) tubulin dimers or whole (laterally unbonded) protofilament segments at the plus end, as well as lateral bond rupture and formation, and hydrolysis of tubulin dimers as stochastic chemical events into the simulation; Figure 2(A) summarizes the different possible events and the associated rates.

Attachment and detachment
At the plus end of each protofilament, a GTP-tubulin dimer can attach with an on-rate where k + is the pseudo-first-order polymerization rate and c tub is the concentration of free GTP-tubulin dimers. The on-rate is assumed to be independent of the hydrolysis state of the protofilament end.
For depolymerization, we assume that a tubulin dimer at the plus end can only detach if it has no lateral bonds. We also allow for detachment of whole protofilament segments starting from an interior dimer (d < d(p)) if the whole segment has no lateral bonds. Laterally unbounded dimers or segments can detach with a rate as given by Kramers theory with the longitudinal standard bond energy ∆G 0 * long (including the entropic cost of "immobilization") and the standard concentration c 0 = 1 M [37]. This approach differs from other models [15,30,42], where tubulin dimers can detach regardless of whether they have lateral bonds or not. In such models, if a tubulin dimer still has lateral bonds, its detachment rate decreases exponentially. In our model, we rather include lateral bond rupture and formation as separate stochastic events into the simulation (similarly to the purely chemical models in [32][33][34]); bond rupture can then be followed by detachment of laterally unbounded dimers or protofilament segments. Bond rupture enables dimer detachment and is necessary prior to a catastrophe; vice versa, bond reformation is necessary for a rescue event. Therefore, it is essential to also include the process of bond formation into the model. Moreover, it has been observed in MD simulations in [43] that lateral tubulin bonds can easily reform. The restriction that only laterally unbonded dimers can detach also causes an indirect increase of the effective off-rate if the last dimers of a protofilament are hydrolyzed because this tends to create stretched bonds, which rupture more easily.

Zipper-like lateral bond rupture and bond formation
We assume that bond rupture between protofilaments starts from the plus end and proceeds by a rupture front monomer by monomer towards the minus end; likewise, bonds can be reformed only monomer by monomer towards the plus end in a zipper-like fashion. As a result, we always have a rupture front between two neighboring protofilaments such that all monomers on top of the front toward the plus end are ruptured and all monomers below toward the minus end are intact. If tubulin monomer (p, d, t − 1) has a lateral bond with its neighbor in protofilament p + 1 but the tubulin monomer on top of it, (p, d, t), has no lateral bond with this neighbor, the rupture front can recede towards the plus end, and tubulin monomer (p, d, t) can form a bond with rate k form = k att (9) with the attempt rate k att . Vice versa, if the bond at (p, d, t) is intact and bond (p, d, t + 1) is broken, the rupture front can advance towards the minus end by rupturing this bond with a rate k rup = k att exp ∆G 0 lat + ∆G mech (10) which contains a chemical bond energy ∆G 0 lat and a mechanical energy ∆G mech , which accounts for the weakening of the lateral bond due to mechanical strain in the bond and enters according to Bell theory [44,45]. In our model, ∆G mech is due to the stretching of the Hookean springs representing the lateral bonds so that ∆G mech = F lat rup , where F lat = −∂E lat /∂| s(p, d, t)| is the force currently acting on the lateral bond and rup is the characteristic bond rupture length. We define rup as the length increase of the lateral bond from its rest length s 0 at which the stretching energy of the spring cancels the bond energy:

Hydrolysis without and with mechanical feedback
Lastly, GTP in beta-tubulin monomers can hydrolyze into GDP via a random (or scalar) hydrolysis rule meaning that almost every GTP-tubulin dimer in the MT can hydrolyze with a fixed rate k hydr regardless of the hydrolysis state of its longitudinal neighbor (which would be a vectorial hydrolysis rule). The "almost" in the previous sentence refers to the finding that the polymerization of the tubulin dimer (p, d) and thus the formation of a longitudinal bond between beta-tubulin (p, d − 1, 2) and alpha-tubulin (p, d, 1) catalyzes the hydrolysis reaction in beta-tubulin (p, d − 1, 2) [13]. As a consequence, only GTP-tubulin dimers that ever had another tubulin dimer on top of them can be hydrolyzed in our model.
We also consider the possibility that hydrolysis is mechanochemically coupled to the bending strain [17]. Then, the hydrolysis rate is modulated with a dimer-specific change ∆E hydr (p, d) in the energy barrier height of the hydrolysis reaction, which depends on the bending state of dimer (p, d). Because this bending state also depends via lateral bonds on the bending states in all neighboring dimers, and because the bending state of all neighboring dimers strongly depends on their hydrolysis state, the hydrolysis dynamics becomes effectively non-random but depends on the hydrolysis state of the neighbors.
The basis for our assumption of a tubulin dimer-specific mechanochemical hydrolysis rate is to view the equilibrium bending angle ∆θ 0 of a dimer as the reaction coordinate for hydrolysis which can be described by an energy profile F hydr (∆θ 0 ). F hydr (∆θ 0 ) has two local minima corresponding to the straight conformation with ∆θ 0 = 0°and the curved conformation with ∆θ 0 = 11°and a rate-limiting energy barrier of unknown height ∆F barrier hydr in between. We propose that hydrolysis of a tubulin dimer is eased if its actual bending angle ∆θ is closer to the equilibrium angle ∆θ 0 = 11°in the hydrolyzed state. We model this dependency by adding a dimer-specific bending energy contribution E hydr (∆θ 0 ) to F hydr (∆θ 0 ), which changes the energy barrier height from ∆F barrier hydr to ∆F barrier hydr + ∆E hydr (p, d), see Figure 3. ∆F barrier hydr can be absorbed into the constant rate k 0 hydr so that only ∆E hydr (p, d) remains in the Arrhenius factor in (12).
To calculate the change in the energy barrier height ∆E hydr (p, d), we now consider the total MT energy in (6) as a function of the hydrolysis reaction coordinate ∆θ 0 while keeping all polar angles {θ(p, d, t)} fixed. We simply assume that the energy barrier is centered between the minima at ∆θ barrier 0 = 5.5°resulting in Because hydrolysis of tubulin dimer (p, d) affects the rest bending angles of beta-tubulin monomer (p, d, 2) and alpha-tubulin monomer (p, d + 1, 1) and the rest bending angles only affect the bending energies (5), we finally  Figure 3: Schematic hydrolysis energy landscape with two local minima corresponding to the straight conformation (∆θ 0 = 0°) and the bent conformation (∆θ 0 = 11°) and an energy barrier ∆F barrier hydr at ∆θ 0 = 5.5°between them. ∆F release hydr is the energy released by hydrolysis. The dashed line represents the modified energy landscape due to the dimer-dependent contribution E hydr (∆θ 0 ).
so that only a local bending energy change has to be calculated. As a result, tubulin monomers in the MT lattice with larger bending angles ∆θ(p, d, t) tend to hydrolyze preferentially. For the terminal tubulin dimer of a protofilament (p, d(p)), the d + 1-term in (14) is missing because tubulin monomer (p, d(p) + 1, 1) does not exist. This results in an overall smaller energy barrier and, thus, a higher hydrolysis rate of the terminal tubulin dimer.
We also see that the base hydrolysis rate k 0 hydr in (12) is not the hydrolysis rate for a perfectly straight MT (∆θ(p, d, t) = 0°for all tubulin monomers) because there is still the constant contribution κ(5.5°) 2 to the energy barrier in (14) that reduces the hydrolysis rate. As these terms are proportional to the bending constant κ, we cannot simply absorb them into the constant factor k 0 hydr . We note that for almost all GTP-tubulin dimers in the GDP-body of the MT, we will typically find negative bending angles; these dimers bend inward in order to allow the longitudinal GDP-dimer neighbors to further bend outwards. For such negative bending angles the hydrolysis rate is reduced according to (14).
In addition to the previous four free parameters from the MT energy, the simulation events add three additional free parameters: the pseudo-first-order polymerization rate k + , the attempt rate k att , and the hydrolysis rate k hydr (or k 0 hydr ). In total, there are now seven free parameters, which are listed in Table 2.

Simulation and parameter determination
The actual MT simulation (implemented in C++) works as follows: 1. Initially, a MT with N GDP GDP-tubulin dimers followed by N GTP GTP tubulin dimers per protofilament is constructed with θ(p, d, t) = 0°for all (p, d, t).
2. Using the tubulin monomers' polar angles {θ(p, d, t)}, the MT's actual initial configuration is determined by minimizing its mechanical energy. Details on the minimization procedure will be discussed in the next section.
3. For all of the events described in the previous section, a list of possible events is determined and based on their rates k i , a "tentative" event time t i is calculated using Gillespie's first reaction method [46]: where r is a uniformly distributed random number from 0 to 1. The event i with the shortest event time t i is executed and the simulation time is increased by t i .
4. Assuming fast mechanical relaxation the MT's energy is minimized after any event.
5. The simulation terminates if a protofilament is shorter than two tubulin dimers. 1 Otherwise we go back to the third step to determine the next event.
There is a general agreement between different experiments [3,[47][48][49][50][51][52] that the MT growth velocity v gro increases linearly with the tubulin dimer concentration c tub and that the shrinkage velocity v shr is independent of c tub . We will use the results by Walker et al [47], which were measured for c tub ∈ [7.7 µM, 15 and lead to an individual critical concentration c tub,c 5 µM (below which v gro < 0).
To determine the values of the model parameters, we use a "divide and conquer" approach [15,30]. First, we consider MT growth, where mechanics are assumed not to play a significant role as protofilaments are not curling outward so that ∆G mech = 0. Thus, we use a GTP-only MT (N GDP = 0) and set k lat = 0 and κ = 0 so that the only free parameters left are k + , ∆G 0 * long , ∆G 0 lat , and k att . The goal of these simulations is to reproduce the measured growth velocity in (16) as function of the free tubulin dimer concentration c tub . Secondly, we consider MT shrinkage, where mechanics are now assumed to play a significant role, i.e., k lat > 0 and κ > 0. For a shrinking MT, we use N GTP = 0, N GDP > 0, and the parameter values already determined by the growth simulations to reproduce the shrinkage velocity in (17). In both cases, hydrolysis is ignored. A schematic overview of the entire parameter determination procedure can be found in Figure S7 in the Supplementary Material.
Comparing the number of free parameters and the amount of experimental data, we can already predict that we will not be able to determine one set of fixed parameter values but only restrict some parameter values to specific values if other parameter values are set to (arbitrarily but reasonably) chosen values. We will discuss this issue in more detail in the conclusion.

Energy minimization
In previous three-dimensional models, different energy minimization approaches have been used. VanBuren et al [15] used a local minimization approach in which they randomly selected individual tubulin dimers and then only locally minimized with respect to the parameters of this dimer. On average, each tubulin dimer was visited three times for minimization. Zakharov et al [18] employed a completely different approach by explicitly modelling the stochastic motion of tubulin monomers in space using Brownian dynamics (applied to the first 300 tubulin dimers at the plus end). They solve Langevin equations every 2 × 10 −10 s while using 10 −3 s as the time step for the events in their simulation resulting in O(10 7 ) dynamics steps between actual events. Using a parallel implementation run on a supercomputer, their simulation took more than a day to simulate 1 s of MT dynamics. There are drawbacks for both approaches: a local energy minimization scheme might not come close enough to a mechanically relaxed configuration, whereas a full Brownian dynamics simulation is computationally very costly. In this paper, we employ a systematic mechanical energy minimization between each stochastic chemical simulation event. We try to achieve a better mechanical energy relaxation than VanBuren et al [15] with significantly less computational steps than Zakharov et al [18].
In our simulation, we use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, a quasi-Newton method, provided by the GNU Scientific Library (GSL) [53] to minimize the total mechanical MT energy in (6) as a function of the polar angles {θ(p, d, t)}. If each protofilament in the simulated MT contains N GDP + N GTP tubulin dimers, there are a total of 26(N GDP + N GTP ) polar angles and thus the same number of minimization parameters. In realistic simulations, MTs can stay in the growing phase for a very long time resulting in an unbounded increase in the number of minimization parameters drastically slowing down the simulation. In essence, the average time for one minimization step increases with the MT length in this scenario making long-running simulations impossible.
To overcome this limitation, we will explore two possibilities to avoid having a MT length-dependent number of minimization parameters: 1. restricting the number of minimization steps per energy minimization to a small value but still considering all minimization parameters (this approach is similar to the strategy in [15]), 2. restricting the number of minimization parameters by only considering the tip of the MT but not restricting the number of minimization steps.
While the first strategy is easy to understand and implement, the second needs further specifications in terms of how we define the tip of the MT here. If a certain event is executed that affects tubulin dimer (p, d), we include all layers starting from min(0, d − d cutoff ) into mechanical energy minimization because mechanical interactions within the MT have a certain range, where d cutoff is a cutoff layer distance.
Below, we will compare these approaches of restricted minimization with respect to accuracy and speed and find that we obtain accurate energy minimization at a high simulation speed by using the second approach and restricting the number of minimization parameters with d cutoff = 10. We can compare with the approaches of Zakharov et al [18] and VanBuren et al [15] in terms of the average number of minimization steps between chemical events.
Zakharov et al [18] use O(10 7 ) Brownian dynamics steps between events and restrict the number of simulation parameters to 300 tubulin dimers at the plus end. With d cutoff = 10 we minimize on average with respect to a comparable number of 150 tubulin dimers at the plus end. To compare the efficiency, we consider a single quasi-Newton minimization step in our simulation to be equivalent to one time step of their Brownian dynamics (if we ignore the random thermal fluctuations in their Langevin equations, they are basically using a gradient descent method). We compare the event time t i divided by the number of minimization steps after the execution of that event to their Brownian dynamics time step of 2 × 10 −10 s. For shrinking MTs, one minimization step takes O(10 −5 s) after polymerization events, O(10 −4 s) after depolymerization events, and O(10 −7 s) after lateral bond events; all of these time steps are orders of magnitude larger than 2 × 10 −10 s and, thus, the simulation proceeds orders of magnitude faster, while we still achieve an accurate energy minimization. As a comparison with the 1 s of MT dynamics simulated in more than a day in a parallel computation in Ref. [18], we generally do not require more than a few hours for 1 min of MT dynamics (for a constant hydrolysis rate) using just a single CPU core.
VanBuren et al [15] apply a local minimization procedure and restrict minimization to, on average, three minimizations with respect to the parameters of each dimer. Because one step of their algorithm minimizes with respect to the parameters of a single tubulin dimer, a comparison to our quasi-Newton minimization steps which minimize the MT energy with respect to the parameters of, on average, O(150) tubulin dimers is not straightforward. In addition, VanBuren et al's model also contains longitudinal springs so that outward bending of single tubulin dimers as a consequence of local minimization can be compensated by stretching the next longitudinal spring. As our model does not contain such longitudinal springs, bending one tubulin dimer causes the whole protofilament part above the tubulin dimer to also bend outwards creating an effectively non-local, far-reaching interaction. Consequently, we are not able to also implement a local minimization procedure for comparison. To make a qualitative comparison between the two approaches, we assume that one minimization step of our BFGS algorithm, which acts on average on 300 parameters, i.e., 150 tubulin dimers, corresponds to 100 single tubulin dimer minimizations in the model of Ref. [15] as they consider three parameters per tubulin dimer. Between chemical events, we perform on average 150 BFGS minimization steps, which corresponds to 1.5 × 10 4 single tubulin dimer minimizations in Ref. [15]. Therefore, we apply the equivalent of 15000/150 = 100 single tubulin minimizations to each of the 150 tubulin dimers close to the plus tip on average as compared to three single tubulin dimer minimizations in the simulation model of Ref. [15]. Accordingly, we should achieve a more accurate mechanical energy relaxation.
We also compared our chosen minimization method, the BFGS algorithm, against the other multidimensional minimization algorithms using derivatives provided by GSL [53], including the conjugate gradient method, and found the BFGS algorithm to perform better. In particular, to fully minimize the initial configuration of a MT with N GDP = 20 and N GTP = 0, BFGS only required about a third of the time compared to the next best algorithm, a conjugate gradient method.

GTP-microtubule growth and model parameterization
MT growth mainly depends on the four parameters k + , ∆G 0 * long , ∆G 0 lat and k att , because the growing MT tip mainly consists of straight GTP-tubulin dimers. Therefore, we consider growth of a GTP-only MT (N GDP = 0) in the absence of hydrolysis and set k lat = 0 and κ = 0 so that the only free parameters left are k + , ∆G 0 * long , ∆G 0 lat , and k att . For k + = 2 µM −1 s −1 and k + = 4 µM −1 s −1 , we scanned the parameter space (∆G 0 * long , ∆G 0 lat , k att ) in steps ∆∆G 0 * long = 0.2 k B T . to find parameter values that reproduce the experimental growth velocity data of Walker et al in (16). The growth velocity v gro for each simulation was determined by fitting MT (t sim ) with a linear function. Experiments on MT growth show a linear dependence v gro (c tub ) = a gro c tub + b gro characterized by two parameters a gro and b gro from (16). If simulations reproduce a linear dependence of v gro as a function of c tub , we can determine two of the three model parameters (∆G 0 * long , ∆G 0 lat , k att ) by fitting to the experimental data (16) for a gro and b gro , i.e., two experimental constraints fix two model parameters as a function of the third parameter. This will allow us to parameterize a one-dimensional sub-manifold (a line) within the threedimensional parameter space (∆G 0 * long , ∆G 0 lat , k att ) where our model agrees with experimental growth data. This procedure is conceptually analogous to the approach of VanBuren et al [30], but we work in a higher-dimensional (three-dimensional) space of model parameters.
As a result, we obtain a line in the three-dimensional parameter space, which we parameterize by ∆G 0 * long , i.e., for a given value of ∆G 0 * long , a value of ∆G 0 lat (see Figure 4(A)) and a value of k att (see Figure 4(B)) is determined by the experimental growth data.
Afterwards, we will fix a particular value of ∆G 0 * long by the additional requirement that the simulation should exhibit an as linear as possible concentration dependence of the growth velocity v gro over a certain range of tubulin concentrations c tub (see Figure 4(D)) such that we arrive at parameter sets (∆G 0 * long , ∆G 0 lat , k att ) for k + = 2 µM −1 s −1 and k + = 4 µM −1 s −1 , see Table 3.
The results in Figure 4(A) show that the values of ∆G 0 * long and ∆G 0 lat depend only weakly on our chosen k + values. Figure 4(A) also shows that our data matches results obtained in [30] (this data was later re-used in [15,16,54]) but also differs from other results [31,55], which were all obtained by the same approach of fitting growth velocity data from Walker et al [47] (or their own growth data in [55]). Kononova et al [43] obtained bond energies from MD simulations of nano-indentation experiments; the values from Kononova et al [43] are much larger for both types of bonds (∆G 0 * long ∼ 2∆G 0 lat ∼ 25k B T ) and, thus, not shown in Figure 4(A). Qualitatively, the measured dependencies of ∆G 0 lat and k att on ∆G 0 * long can be understood as follows: the weaker longitudinal bonds are, the more likely it is that a tubulin dimer will depolymerize. To get the same growth velocity, this decrease in "longitudinal stability" has to be compensated by an increase in "lateral stability" by stronger lateral bonds (making it less likely that lateral bonds break and, thus, enabling depolymerization) or faster formation of lateral bonds (to stabilize newly polymerized tubulin dimers). Figure 4(C) shows the number of tubulin dimers d depoly that detach at once during depolymerization events. For increasingly stronger longitudinal bonds and, thus, weaker lateral bonds, multi-dimer depolymerization becomes more relevant. The data in the inset in Figure 4(C) is also compatible with results in Ref. [33] obtained with a purely chemical model.
Until now, we only considered free tubulin dimer concentrations c tub ∈ [7 µM, 16 µM] to use similar values as Walker et al [47], but there have also been other measurements with a larger range of c tub values [3,48,51,52]. In general, it is assumed that the growth velocity v gro increases linearly with c tub for the whole MT just as the polymerization rate in (7) increases linearly with c tub for individual protofilaments. Theoretically, it has been shown that, for multistranded polymers, lateral interactions give rise to a non-linear relation between growth velocity on monomer concentration [56]. For MT growth, a non-linear dependence on tubulin concentration was found in Ref. [31] using a two-dimensional model based on Ref. [30]. Over a larger range of c tub values, our simulations also exhibit a non-linear relation between v gro and c tub depending on the value of ∆G 0 * long , as shown in Figure 4 [47], see (16). To compare our lateral bond energies (per tubulin monomer) to other publications (lateral bond energy per tubulin dimer), the y-axis shows 2∆G 0 lat . (The numbers behind Ref. [30] refer to their value of k + .) (B) Lateral bond attempt rate k att as a function of the longitudinal bond energy ∆G 0 * long for our two values of k + from matching the concentration-dependent growth velocity data from Walker et al [47], see (16) large range of c tub values [3,48,51,52]. Therefore, we determined the remaining free parameter value of ∆G 0 * long for the two k + values from the condition that the concentration dependence of v gro is as linear as possible up to 50 µM. To determine these values of ∆G 0 * long , we ignored concentrations c tub below the individual critical concentration (for which v gro < 0) which violate our fundamental assumption of a growing MT.
In summary, we find a triple (∆G 0 * long , ∆G 0 lat , k att ) that fits the growth velocity data from Walker et al [47] and that gives a linear concentration dependence over a wide tubulin concentration range for two representative values of k + . Table 3 lists these parameter triples for k + = 2 µM −1 s −1 and k + = 4 µM −1 s −1 . For a given k + , these results fix four of the seven model parameters in Table 2 using experimental data on MT growth. To address the parameters κ and k lat , we now turn to MT shrinkage.

GDP-microtubule shrinkage and model parameterization
As opposed to MT growth, MT shrinkage also depends on the bending constant κ and spring constant k lat as protofilament curling and bond rupture become relevant processes for a shrinking MT. We consider a shrinking MT that initially only consists of GDP-tubulin dimers (N GTP = 0, N GDP > 0) with parameter values k + , ∆G 0 * long , ∆G 0 lat , and k att as already determined by the growth simulations and in the absence of hydrolysis (a shrinking, initially GDP-only MT acquires some GTP-dimers by attachment but remains GDP-dominated). To investigate shrinkage, MTs with N GDP = 20 and N GTP = 0 were used. For each parameter set, 20 simulations were run to get an average shrinkage velocity v shr . Experimental data on shrinking MTs show a shrinkage speed v shr that is independent of the tubulin dimer concentration. For each value of k + , we should be able to determine one of the two parameters (κ, k lat ) as a function of the other parameter by fitting such that the experimental value of the shrinkage velocity is reproduced in simulations (for parameters ∆G 0 * long , ∆G 0 lat , and k att fixed by the growth velocity data). We use the experimental shrinkage velocity of Walker et al, see (17), for this fitting procedure.
10 0 10 1 10 2 10 3 10 4    (17). All data points for each ∆G 0 * long fall on square root functions This functional dependence can be understood qualitatively by considering the mechanical contribution to the bond rupture rate (10), exp(F lat rup ), which, on average, should have the same value for all mechanical parameter combinations to produce the same shrinkage velocity. As the characteristic bond rupture length in (11) depends on k lat as rup ∼ √ k lat −1 , the average lateral bond force at rupture should depend on k lat like F rup ∼ k lat rup ∼ √ k lat . The lateral bond force F lat is a consequence of the lateral bonds stretching as the tubulin monomers curl outward to decrease the bending force  . Also, the resulting mechanical contribution F rup rup for the exponential function of the lateral bond rupture rate in Figure 5(D) is approximately constant as expected from our above argument.
As the experimentally measured shrinkage velocity v shr does not depend on the free tubulin dimer concentration c tub , we used constants to fit our v shr (c tub ) data. In reality, however, our data shows a linear dependence between v shr and c tub as shown in Figure 5(E) corresponding to a slowing down of depolymerization. This is caused by an increased probability for intermediate addition of tubulin dimers and lateral bond formation between them; these lateral bonds require additional time to rupture. While this dependency of v shr on c tub will have a small influence on the concrete value of the shrinkage velocity, we expect it to not have any qualitative effect on the overall MT dynamics. At higher tubulin concentrations, where the decrease of |v shr (c tub )| would become significant, the catastrophe rates decrease dramatically so that shrinking will rarely occur.
Comparing our results from figures 5(A) and (B) to other results is not always directly possible due to different modelling approaches but most find that k lat 1000 k B T /nm 2 and κ 1000 k B T /rad 2 [15,57,58], with some exceptions [43,59]. Previously, we used MD simulation data from Grafmüller et al [60] to calculate the bending constant κ [17]. Compared to Ref. [17], we have to adjust the calculation to consider both inter-dimer and intra-dimer bending resulting in κ 50 k B T /rad 2 . MD simulation in [43], on the other hand, give a persistence length of individual protofilaments of L p 6 µm, which corresponds to a significantly larger value of κ 1500 k B T /rad 2 for the bending constant. This discrepancy cannot be resolved at present. We use κ = 149 k B T /rad 2 in the following together with the corresponding value of k lat = 100 k B T /nm 2 according to Figure 5(B) which are values close to the ones used by [15].

Restricted energy minimization for efficient simulation
Until now, energy minimization was not restricted by either a maximum number of minimization steps or by only considering a subset of tubulin dimers at the MT tip so that we will consider this unrestricted minimization as the "gold standard" to which we will compare the two restricted energy minimization approaches described in Section 2.4. We use the shrinkage velocity v shr as the observable by which we judge the relevant cutoff values in the two approaches.
For restricting the number of quasi-Newton minimization steps, Figure 6  show that reducing the number of minimization steps by a factor of 10 can lead to deviating growth velocities. Therefore, the improved energy relaxation that we obtain in comparison to Ref. [15] by applying the equivalent of one order of magnitude more minimization steps should be relevant.
If minimization is restricted to a subset of minimization parameters at the tip of the simulated MT, this subset is defined by the cutoff distance d cutoff . To have a maximum improvement in simulation speed, d cutoff should be as small as possible. It is evident from the data shown in Figure 6(B) that values d cutoff < 5 have a detectable influence on the shrinkage velocity. We also ran some simulations with N GDP = 50 and also for k + = 2 µM −1 s −1 (see Figure S8 in the Supplementary Material) and based on all data, we choose d cutoff = 10 as a conservative value for the cutoff distance.
In summary, we are more confident in the second approach to only minimize the MT tip where actual conformational changes happen, because for this subset, the restricted energy is fully minimized. Additionally, the first approach still has the issue of slowing down with an increasing number of minimization parameters as all minimization parameters are considered. The second approach ensures that the number of minimization parameters does not scale with the MT length but remains bounded, which assures that we can simulate arbitrarily long growing MTs at a fixed minimal computational speed. In the first approach, the quality of the minimization will probably also decline because the number of minimization parameters increases while the number of minimization steps is kept constant. Lastly, the first approach, in contrast to the second approach, does not guarantee that the upper, i.e., the dynamic part of the MT is properly minimized.
We also note that in the presence of mechanical feedback onto hydrolysis, simulations take longer because minimizations after hydrolysis events need to consider more tubulin dimers if the hydrolyzed tubulin dimer is relatively deep in the MT lattice (see Supplementary Material for more details).

Full simulations exhibit repeated catastrophe and rescue events
Based on the previous section on energy minimization, we use d cutoff = 10 for full simulations in which the initial MTs have both a GDP body and a GTP cap, thus N GDP > 0 and N GTP > 0. We now aim for realistic MT dynamics with repeated phases of growth and shrinkage in the same simulation and catastrophe and rescue events in between. First, we only consider strictly random hydrolysis with a hydrolysis rate k hydr that is independent of tubulin dimers' position or mechanical forces and which is another unknown free parameter in our model. Hydrolysis coupled to mechanics via (12) will be considered later.
It poses a computational challenge for chemomechanical MT models to reach time scales of MT dynamics where repeated catastrophe events occur at realistic hydrolysis rates k hydr and tubulin dimer concentrations c tub . In Ref. [18], where mechanics was implemented via full Brownian dynamics, only short times scales could be reached (although the Brownian dynamics was applied to only 300 tubulin dimers at the plus end). Therefore, they increased the hydrolysis rate from their "normal" value of 0.5 s −1 (based on the 2 s delay between polymerization and phosphate release measured by [61], which is also used by [62]) into a range of 3 s −1 to 11 s −1 in order to trigger catstrophe events within computationally accessible time scales. They found a linear scaling of catastrophe rate with k hydr and employed a linear extrapolation to obtain catastrophe rates for realistic hydrolysis rates (see their Figure 3A). In our simulations, we observe that increasing k hydr beyond a certain (c tub -dependent) value leads to immediate MT shrinkage because the initial cap quickly hydrolyzes; this can be interpreted as an instantaneous catastrophe. In such cases (like in Figure 23 for c tub = 7 µM and k hydr = 0.5 s −1 ), there is no real growth phase based on which a catastrophe frequency could be determined. For these hydrolysis rates, the individual critical concentration c tub,c (where v gro = 0 is reached) has apparently increased above the given tubulin concentration.
The experimental data on the hydrolysis rate is limited, so that many publications determine the hydrolysis rate themselves by matching simulation results with experimental data [16,[30][31][32][33][63][64][65]. There are, however, more direct measurements in [61]. In most models and also in measurements from [61], the (random) hydrolysis rate is in the range of 0.1 s −1 to 0.5 s −1 (Ref. [30] use a relatively high value of 0.95 s −1 ). We explore exactly this range of hydrolysis rates, see Figure 23. Figure 23 shows MT growth curves (length vs. time) over simulation times up to t sim = 10 min for several representative tubulin concentrations and realistic hydrolysis rates. MT growth trajectories as in Figure 23 for other k lat values can be found in Figures S10, S11, S12, and S13 in the Supplementary Material. Simulations in Figure 23 were started with N GTP = 10 and N GDP = 20, but results are largely independent of the initial ratio N GTP /N GDP (see, for example, Figure S11 in the Supplementary Material). Our chemomechanical MT model is computationally efficient such that we can determine catastrophe and rescue rates as inverse average growth and shrinking times between repeated catastrophe and rescue events. In the Supplementary Material, we explain the algorithm that we used to identify catastrophe and rescue events and, thus, growth and shrinking times from MT simulation trajectories in detail. The results are shown in Figure 8. In comparison to typical experimental data [47,66], this decrease of the catastrophe rate with tubulin concentration seems too steep. Current phenomenological models for the MT catastrophe rate as a function of tubulin concentration can be found in [67,68], experimental data in [47,69]; the decrease of the catastrophe rate with GTP-tubulin concentration c tub appears steeper in the simulation for all hydrolysis rates k hydr = 0.1 s −1 to 0.5 s −1 .
In the following, we will discuss two aspects of MT growth and catastrophes in more detail, namely the dependence of growth velocity on hydrolysis rate and the detailed dynamics within single catastrophe events, which become accessible within a computational model and are impossible to address experimentally.

Growth velocity reduces linearly with hydrolysis rate because of cap structure
So far, we parameterized the model by fitting the growth velocity of GTP-only MTs, i.e., in the absence of hydrolysis to the experimentally measured velocity in (16). Hydrolysis reduces this growth velocity by increasing the probability of GDP-dimers dimers at the plus end. This increases the rate of bond rupture because hydrolyzed dimers tend to create stretched bonds which rupture more easily (there is no direct increase of the off-rate for hydrolyzed GDP-dimers in our model). As only laterally unbounded dimers can detach, hydrolyzed GDP-dimers at the plus end have an effectively higher detachment rate.
The last row of Figure 23 indicates and Figure 9(B) shows explicitly that increasing the hydrolysis rate decreases the growth velocity linearly although the growth reduction mechanism is indirect via the increased probability of bond rupture for hydrolyzed GDP-dimers. Our model parameterization was such that we obtain the experimentally measured growth velocities by Walker et al [47] at k hydr = 0 s −1 in Figure 9(B). Nevertheless, Figure 9(A) shows that there is still a linear relation between the free tubulin dimer concentration c tub and the growth velocity v gro so that it is possible to re-adjust parameters to reproduce the growth velocity in the presence of hydrolysis, once a particular hydrolysis rate can be reliably selected.
Because both the dependence on tubulin concentration in Figure 9(A) remains linear and the reduction by the hydrolysis rate in Figure 9   : Growth velocity v gro as a function of (A) the free tubulin dimer concentration c tub for different hydrolysis rates k hydr and as a function of (B) the hydrolysis rate k hydr for different free tubulin dimer concentrations c tub in comparison to the experimental data from [47]. (C) Average GTP-tubulin cap length N cap of protofilaments and (D) fraction of protofilaments without a GTP cap as a function of the hydrolysis rate k hydr . The standard set of parameters from Table 2 was used.
reached) increases linearly with the hydrolysis rate beyond the value c tub,c 5 µM of Walker et al [47]. Figure 23 clearly shows that increasing k hydr actually increases the individual critical concentration c tub,c . 2 The mechanism of growth velocity reduction by hydrolysis can be further elucidated by comparing the average GTP-tubulin cap length N cap of protofilaments (see Figure 9(C)), and the fraction of protofilaments without a GTP-cap (see Figure 9(D)): The higher the hydrolysis rate is, the smaller the GTP-cap and the higher the fraction of cap-less protofilaments is. 3 The increase in GDP-tubulin dimers depolymerizing from the protofilament tips for higher hydrolysis rates is due to an increase in the probability of uncapped protofilaments with the hydrolysis rate as shown in Figure 9(D). In Ref. [70], dependencies N cap ∝ c tub /k hydr and p(N cap = 0) ∝ k hydr /c tub have been predicted, which are in agreement with Figure 9(C) and (D).

Detailed dynamics within single catastrophe and rescue events
The chemomechanical model reproduces realistic MT dynamics including catastrophe and rescue events. Figure 10 shows typical MT growth paths featuring two catastrophe events and a rescue event in subfigure (C). Moreover, we observe "dips" in the growth path where a short phase of shrinking appears, which are similar to "stutter" events that have been observed in Ref. [35]. Videos of these two simulations with two-and three-dimensional representations of the MT structure can be found in the Supplementary Material. Using our computational model, we can systematically identify the point in a MT growth path, where a catastrophe becomes structurally unavoidable. This allows us to search for typical catastrophe-triggering features in MT growth. To analyze how probable it is at specific points in the simulation of MT dynamics that the MT continues a certain growth path, we chose two simulations with at least one significant event (meaning a catastrophe, rescue, or a "dip"/"stutter") and took configurations around such events as starting points for new simulations (similar to [33]). In these new simulations, MTs were allowed to grow (or shrink) for a maximum of 60 s, a sufficient amount of time to check if the new simulations show dynamics similar to the original simulation around the significant event.
The MT growth trajectory shown in Figure 10(A) has two significant events: a dip at t sim = 1.2 min and a catastrophe at t sim = 6.85 min; the trajectory in Figure 10(C) contains three significant events: a dip at the very beginning, a catastrophe at t sim = 9.15 min, and a rescue at t sim = 9.54 min. To determine whether newly run simulations with starting points from the initial simulation qualitatively follow the original simulation, we need criteria to identify dips, catastrophes, or rescue events. The exact criteria for these events in Figure 10(A) and (C) are stated in the Supplementary Material. In short, in order to identify whether a new simulation reproduces a catastrophe, we check after a time of 10 s to 15 s whether the MT is sufficiently short that a catastrophe must have happened; for a dip, we check whether the MT continued to grow without entering a catastrophe; for a rescue, we check that the MT did not completely vanish because it continued to shrink. For each initial configuration, we ran 20 new simulations and calculated the fraction of simulations that fulfilled these criteria. These fractions are the probabilities for the original growth path at different points in time, and they are shown color-coded in all the insets in Figure 10.
Both catastrophes and the rescue show that the transition from a high probability to stay in the current dynamic state to a high probability to switch into the other dynamic state occurs within a few seconds. In Figure 10(A) and (C), we first observe that catastrophes become practically unavoidable (red color code in (A.C) and (C.C)) after a phase of relatively slow shrinking by 50 nm to 100 nm; similar "transitional catastrophe" behavior has been observed in Ref. [35]. A dip, on the other hand, can only evade a catastrophe (yellow to red color code in (A.D) and (C.D)) if the MT length shrinks by significantly less than 50 nm. Because hydrolysis followed by straining and rupture of the lateral bonds is required before a laterally unbonded dimer can detach, MT shrinking by 50 nm suggests that roughly 6 dimer layers must hydrolyze in a row to trigger a catastrophe. This is, however, not sufficient to remove the entire GTP-cap. The GTP-cap length averaged over all protofilaments is still > 1 when the catastrophe becomes unavoidable (at points 3 in Figure 10(A) and 6 in Figure 10(C), see also Figure S1 in the Supplementary Material). As the corresponding MT snapshot insets 3 and 6 reveal, the reason for this discrepancy is the average over all protofilaments: it appears that typically only a "nucleus" of three neighboring protofilaments shrinks by more than 6 dimers, such that its GTP-cap is removed and its ends reach into the GDP-body of the MT, when a catastrophe is triggered. The MT snapshots in Figure 10(B) also suggest that rescue events require formation of a GTP-cap on almost all 13 protofilaments (with an average cap length ∼ 4) such that nuclei of three neighboring uncapped GDP-protofilaments are avoided. Further investigation of more catastrophe events will be necessary to definitely deduce catastrophe-and rescue-triggering structural MT features.

Hydrolysis coupled to mechanics changes the cap structure
Finally, we test how a mechanical feedback onto the hydrolysis rate as introduced in (12) and (14) changes the cap structure and dynamic behavior. In the presence of this mechanical feedback, tubulin dimers in the MT lattice with larger bending angles tend to hydrolyze preferentially.
Overall, we find a linear relation between k 0 hydr and the average hydrolysis rate k hydr (see Figure 11(A)) with k 0 hydr k hydr . When comparing MT growth with hydrolysis coupled to mechanics with average hydrolysis rate k hydr to MT growth with constant hydrolysis rate k hydr (for exmaple in Figure 11(D)-(F)), we use Figure 11(A) to choose the base hydrolysis rate k 0 hydr such that k hydr ≈ k hydr . Figure 11(B) shows the average hydrolysis rate k hydr as a function of the free tubulin dimer concentration c tub for k 0 hydr = 1.5 s −1 . Here, we observe a pronounced nonlinear concentration dependence with a decrease around the individual critical tubulin concentration c tub 10 µM. At the same concentration, also the porous cap length N pcap (see Figure 11(C)), which is defined as the difference between the number of tubulin dimers in a protofilament and the value of d of the first GTP-tubulin dimer counted from the minus end, starts to increase. As a result, the porous cap length for hydrolysis coupled to mechanics is much longer compared to a constant hydrolysis rate, even if the average effective hydrolysis rate is roughly the same. In the following, we argue that the reason for this increase in porous cap length is a decrease of the hydrolysis rate for GTP-dimers away from the tip. Mechanical feedback gives rise to preferential hydrolysis at the tip, i.e., the average hydrolysis rate k hydr (x) (over all actually executed hydrolysis events) is larger for small layer distances x ≡ d(p) − d from the tip, as can be seen in Figure 11(G). This is in line with previous results in Ref. [17] from a much simpler version of our model with a deterministic hydrolysis kinetics and without dimer attachment and detachment.
According to (14), GTP-tubulin dimers with larger bending angles tend to hydrolyze preferentially. If a straight GTP-dimer is bent inward (∆θ < 0°), its hydrolysis rate is reduced according to (14); if it is bent outwards (∆θ > 0°) the rate is increased. From the hydrolysis rates shown in Figure 11(G), it is possible to calculate the average bending angles using (12) and (14), The results for these bending angles as a function of the distance x from the top are shown in Figure 11(H). Surprisingly, almost all dimers are bent inwards (∆θ < 0°) on average apart from dimers close to the tip, We will try to interpret these results in the following.
An isolated GTP-dimer within the GDP-body can alleviate the bending stress of GDP-dimers by bending inward (∆θ < 0°), which allows longitudinally neighboring GDP-dimers to bend outwards (such that ∆θ > 0°) resulting  in an overall decrease of the MT energy (see Figures S14 and S15 in the Supplementary Material). Therefore, isolated GTP-dimers deep in the GDP-body hydrolyze with a reduced asymptotic rate k hydr ∞ k hydr .
We also find that, for several consecutive GTP-dimers in the same protofilament, GTP-dimers curl inward directly at the GDP/GTP interface resulting in a reduced hydrolysis rate (see Supplemental Figure S15), while GTP-dimers in the center of a GTP-island are straight so that they have a higher hydrolysis rate than at the GDP/GTP interfaces. Effectively, this hydrolysis rate distribution within a GTP-island results in a "anti-vectorial" hydrolysis mechanism with which GTP-islands are hydrolyzed from the interior in contrast to vectorial hydrolysis where hydrolysis happens at the GTP/GDP interfaces.
Also for GTP-dimers in layers closer to the MT tip, other longitudinally close-by GTP-dimers cooperate in alleviating bending stresses; then inward bending is still preferred, but the inward bending angle becomes smaller. This decrease in inward bending corresponds to an increase of the average hydrolysis rates k hydr (x) for GTP-dimers in these layers compared to GTP-dimers buried deeper in the MT body (see Figure 11(G) and (H)). For terminal tubulin dimers (x = 0), we observe a hydrolysis rate k hydr (x) higher than k hydr (while it is equal or lower than k hydr for all other layers x > 0). Hydrolysis in the first layer is enhanced because there are no tubulin dimers on top, such that hydrolysis has to overcome a smaller energy barrier as pointed out previously (the d + 1-term in (14) is missing corresponding to the δ d,d(p) -contribution in (19)).
As a result of the hydrolysis bias toward the tip, the spatial GTP-tubulin dimer distribution also differs. For concentrations where the MTs are growing only on time scales of several minutes (c tub ≥ 11 µM) for the chosen parameters, a constant hydrolysis rate leads to the expected exponential distribution of GTP-dimers shown in Figure 11(E) as observed in in vivo experiments [71]. Using an effective one-dimensional (or single protofilament) model similar to [63] to calculate the probability of tubulin dimers being GTP-tubulin dimers as a function of the polymerization rate k on , effective depolymerization ratek off , and hydrolysis rate k hydr matches the simulation results for concentrations at which the MTs can be considered in a steady state of growth (see Section 4 in the Supplementary Material). We use an effective depolymerization ratek off instead of k off , because we map onto the depolymerization process of a one-dimensional model so thatk off includes all effects from lateral bond formation and rupture and the actual depolymerization process in the full model.
If hydrolysis is coupled to mechanics, the spatial distribution is only exponential in its tail, has larger values at the MT tip, and GTP-tubulin dimers can be found much deeper in the GDP-body (see Figure 11(F)). These results reflect that the average hydrolysis rate k hydr (x) is decreasing towards the GDP-body and reaches a small limiting value k hydr ∞ k hydr for distances x = d(p) − d > 500 away from the tip, which governs the exponential tail (see Figure 11(G)). This can be rationalized by considering the probability p GTP (x) to find a GTP-dimer at distance x from the tip in a single protofilament and continuum approximation. The balance between attachment/detachment and hydrolysis leads to in the stationary state, which results in a sharp initial decrease of p GTP (x) because k hydr (0) is large at the tip but a much slower asymptotic exponential decrease when k hydr (x) ≈ k hydr ∞ k hydr , which explains the main features in Figure 11(F). In Section 4 in the Supplementary Material, we show that (31) describes simulations with a constant hydrolysis and with hydrolysis coupled to mechanics equally well. With p GTP (x), we can define an "average cap length" as¯ cap = ∞ 0 dx p GTP (x)x. This average cap length¯ cap is longer if hydrolysis is coupled to mechanics compared to a constant hydrolysis rate because p GTP (x) is much greater for larger x (see Figure 11(E) and (F)). As¯ cap < N pcap , this increase in average cap length also explains the increased porous cap length if hydrolysis is coupled to mechanics.
The relative increase of hydrolyzed GDP-dimers at the tip could make MTs more prone for catastrophes and give rise to an increased catastrophe rate and, eventually, a more realistic concentration dependence of catastrophe rates. Figure 12, however, shows that this is not the case. Instead, the same steep dependence on the (base) hydrolysis rate as in Figure 23 persists. In comparison to the MT growth trajectories with a constant hydrolysis rate shown in Figure 10, Figure 13 shows an example of a MT simulation in which the hydrolysis rate is coupled to mechanics. To calculate the probabilities shown in the insets, the same criteria as for Figure 10(A) were used. At first sight, these trajectories look similar to the corresponding trajectories for a constant hydrolysis rate Figure 10(A). There is, however, a significantly increased roughness of the trajectory during the growth phase, which could be interpreted as increased occurrence of "dips" or "stutter" events. A high probability of stutter events has also been observed in Ref. [35], which supports the existence of a mechanochemical coupling in hydrolysis. The catastrophe-triggering configuration of a "nucleus" of several neighboring protofilaments shrinking by more than 6 dimers is also similar as snapshots 4 and 5 in Figure 13(B) show.

Dilution experiments
In dilution experiments, the free tubulin dimer concentration c tub is reduced to c dil c tub at a certain point in time [5,72,73]. If the diluted concentration is sufficiently small or zero, the GTP-cap stops growing by polymerization (and depolymerizes) but continues to hydrolyze; after a characteristic delay time ∆t delay , the  We expect the delay time to be proportional to the GTP-cap length, ∆t delay ∝ N cap , as corroborated by Figure 14(C) and N cap ∝ c tub /k hydr according to Section 3.5 (see Figure 9(C) and (D)) [70]. This results in ∆t delay ∝ c tub /k hydr , which is in qualitative agreement with our simulation data in Figure 14(A) and (B). The comparison with the experimental dilution data from Ref. [73] in Figure 14(B) shows that delay times for a hydrolysis rate k hydr = 0.1 s −1 come close to the experimental data but appear to depend too steeply on c tub .

Discussion
We introduced, parameterized, and analyzed a chemomechanical model for MT dynamics in which, in addition to polymerization (attachment of dimers), depolymerization (detachment of dimers), and hydrolysis of dimers, the rupture of lateral bonds between monomers in neighboring protofilaments is explicitly modeled and coupled to the mechanics of the MT. The basis for this coupling is the allosteric model according to which a hydrolyzed dimer acquires a more bent configuration, which builds up mechanical stress in the MT tubular structue via lateral bonds between dimers.
As many model parameters as possible have been determined from the experimentally measured MT growth and shrinkage velocities measured by Walker et al [47]. To determine the values of the model parameters, we use a "divide and conquer" approach [15,30]. We used simulations of growing GTP-only MTs to parameterize longitudinal and lateral bond energies ∆G 0 * long and ∆G 0 lat and the attempt rate k att for lateral bond formation. By requiring a linear concentration dependence of growth velocity, we can fix all three parameter values for a given value of k + . We used simulations of shrinking GDP-only MTs to parameterize the bending constant κ and the spring constant k lat of the lateral bonds. Here, we can only fix one of the two parameters. Moreover, the hydrolysis rate k hydr is still a free parameter, for which we use values in the range 0.1 s −1 to 0.5 s −1 known from experiments [61].
The general philosophy of a divide-and-conquer approach is the successive fixation of simulation parameters by using first GTP-only growth, then GDP-only shrinkage and, eventually, catastrophe frequencies or dilution to fix the hydrolysis rate. This successive fixation is, however, problematic, as the corresponding experimental data is influenced by all simulation parameters in general. The problem becomes apparent when considering the hydrolysis rate: changes in the hydrolysis rate also affect the growth rate over a wide concentration range because hydrolyzed dimers have an effectively higher detachment rate, see Figure 9(B). Strictly speaking, all simulation parameters in Table 2 must be determined at once by fitting several experimental results simultaneously instead of the successive fixation in the divide-and-conquer approach or to apply the divide-and-conquer approach iteratively several times until a self-consistent parameter set is found. A simultaneous fixation of all parameters has been performed, for example, in Ref. [31] on a chemical model without bond rupture and, thus, with only four parameters (on-rate, bond energies, and hydrolysis rate). Future work on our model should include at least a re-adjustment of the parameters once a hydrolysis rate is selected such that the growth velocity of Walker et al [47] is reproduced in the presence of hydrolysis. If mechanical feedback onto hydrolysis is included, the model has to be re-parameterized again, in principle.
Our simulation model handles all chemical events, i.e., dimer attachment and detachment, bond rupture and formation, and hydrolysis using a Gillespie algorithm. After each chemical event we relax the resulting MT structure mechanically by mechanical energy minimization based on the assumption that the microscopic mechanical dynamics is much faster than the chemical steps. Therefore, mechanical energy minimization is the computationally most demanding step in the simulation. This is a common problem in all dimer-based chemomechanical MT models [15,16,18]. We address this problem by restricting the mechanical energy minimization to bounded number MT degrees of freedom near the plus end. We showed that restricting energy minization to a depth of d cutoff = 10 additional layers into the MT (in minus end direction) from the point of the last chemical event is an accurate and efficient choice. Computational efficiency of this procedure is better than performing a dedicated microscopic Brownian dynamics simulation [18] and better than random local energy minimization [15,16] (for the same accuracy in energy minimization). The restricted energy minimization strategy also ensures that the number of minimization parameters does not scale with the MT length but remains bounded, which assures that we can simulate arbitrarily long growing MTs at a fixed minimal computational speed using our approach.
Simulations do not require more than a few hours for 1 min of MT dynamics (for a constant hydrolysis rate) using just a single CPU core. Therefore, we can reach time scales of several minutes of MT dynamics which is the time scale for repeated catastrophe events for concentrations above the individual critical concentration, where the dynamic instability can occur. We performed a first systematic analysis of catastrophe and rescue rates in Figure 8, which indicates that the decrease of the catastrophe rate with tubulin concentration is too steep compared to experimental data [47,69]. It is also much steeper than simulation results of Ref. [18] but these results for the catastrophe rate relied on linear extrapolation from unrealistically high hydrolysis rates (3 s −1 to 11 s −1 ) down to realistic values (0.1 s −1 to 0.5 s −1 ). In the future, our computational model can also be used to measure the dependence of catastrophe rates on MT lifetime [66].
Within our model, we could also study single catastrophe and rescue events in detail, see Figure 10. The growth paths appear very similar to experimentally observed catastrophe and rescue events. Catastrophes typically feature an initial "transitional" phase of slow shrinking by 50 nm to 100 nm as also observed in Ref. [35]. Moreover, we observe "dips" in the growth paths resembling the "stutter" events from Ref. [35].
The most interesting results of chemomechanical models are possible statements about the typical catastrophetriggering configurations. In this respect, our simulations indicate that a catastrophe could be triggered by a "nucleus" of three neighboring protofilaments shrinking by more than 6 dimers, such that its GTP-cap is removed and its ends reach into the GDP-body of the MT. To rescue a shrinking MT the GTP-cap has to be re-established on almost all 13 protofilaments such that nuclei of three neighboring uncapped GDP-protofilaments are avoided. This shows that mechanical correlations in the dynamics of protofilaments are important in triggering catastrophe events. This is an aspect which is absent in the calculation of catastrophe frequencies based on simplified purely chemical models such as in Ref. [67], where protofilaments are regarded as effectively independent and uncorrelated.
Our model can achieve qualitative agreement with experimental data on dilution experiments (see Figure 14(C)) from Ref. [73] for relatively low hydrolysis rates of k hydr = 0.1 s −1 , which is an indication that the catastrophe mechanism is correctly captured by our chemomechanical model. This also constrains the hydrolysis rate, which is still a free parameter in our model, to lower values around k hydr = 0.1 s −1 .
Finally, we explored the consequences of a mechanochemical coupling in the hydrolysis of tubulin dimers. Because hydrolysis gives rise to bending of the GTP-dimers, we argue that mechanical forces on a dimer that increase its bending angle should also lead to higher hydrolysis rates, see (12) and (14). In the presence of mechanical feedback, hydrolysis gets a bias towards the MT plus end which, in turn, also causes an increase in porous cap length. At the same average hydrolysis rate, hydrolysis in the immediate tip of the GTP-cap is more likely while it is less likely in the remaining part of the cap such that GTP-tubulin dimers can be found much deeper in the GDP-body, see Figure 11. Individual catastrophe and rescue events (see Figure 13) look qualitatively similar in the presence of mechanical feedback but the probability of "dips" or "stutter" events is increased in agreement with Ref. [35]. The coupling of hydrolysis to mechanics does not increase catastrophe rates significantly such that the steep decrease of the catastrophe rate with tubulin concentration persists.
The main problem of our model appears to be the steep decrease of catastrophe rate with tubulin concentration, which could hint to a failure of basic assumptions. One possibility is that a direct effect of the hydrolysis state of the dimer onto the off-rate (as also suggested by atomistic simulations [74]) is relevant and not included in the model. Another possibility is a failure of allosteric models in general. The steep decline of catastrophe rates with tubulin concentration gives a hint that MTs are structurally too stable for GTP-rich caps. This might provide evidence for a shortcoming of the underlying allosteric model, which inserts GTP-dimers in a straight configuration that is more prone to form stable lateral bonds than a curved configuration. An alternative are so-called lattice models [20,21], according to which dimers are always bent but hydrolysis affects lateral and longitudinal dimer interaction energies. A systematic comparison of allosteric and lattice models towards the resulting concentration dependence of catastrophe rates within the framework provided here could help decide which class of models is more appropriate.
So far, almost all chemomechanical modelling approaches were based on the allosteric model [14][15][16][17][18][19] but recent experimental advancements in the analysis of the structure of MT tips [26] demonstrated that both growing and shrinking MTs have bent protofilament ends supporting similar earlier results [75][76][77][78]. Additionally, calculations using MT structures with different nucleotide content in the beta-tubulin [23] and all-atom MD simulations of GTP-and GDP-only MTs [24,25] revealed that hydrolysis weakens lateral bonds and strengthens longitudinal bonds. Both aspects support the lattice model for the influence of hydrolysis on MT mechanics. There is, however, also evidence from MD simulations for intermediate models, where hydrolysis affects interactions and also leads to a much lower GDP-tubulin flexibility [27]. Independent of these findings, our study based on the allosteric model is valuable for the following reasons: (i) In both the allosteric and the lattice model, catastrophes are cascades of lateral bond rupture and in both models, the bent shape of GDP-dimers is the dominating cause of mechanical strain in the MT structure. In the allosteric model, bending and mechanical strain is directly generated by the hydrolysis of GTP-dimers, whereas in the lattice model, the tubulin dimers are always bent but hydrolysis weakens lateral bonds. In both models, the result is an increased lateral bond rupture rate of mechanically strained bonds after hydrolysis. Therefore, an explicit modelling approach for lateral bond rupture as a stoachstic process under force generated by the bending of GDP-dimers will also be important in all future chemomechanical models based on the lattice model. So far explicit stochastic models of lateral bond rupture have only been included into two-dimensional models lacking explicit mechanics [32][33][34] or with heavy computational cost by explicitly simulating the Brownian dynamics of dimers and bonds [18]. (ii) The importance of lateral bond rupture becomes particularly clear for shrinking MTs or MTs entering a catastrophe. In these phases of the dynamic instability, GDP-tubulin dimers are significantly more relavant than GTP-tubulin dimers. As GDP-dimers are bent in both models and this bending gives rise to lateral bond stretching, we believe that both types of models will display a very similar behavior in these phases. The only difference in this scenario is that in the lattice model, the lateral bond energy ∆G 0 lat in the rupture rate (10) will depend on the nucleotide type of the bonded tubulin monomers, which also makes k rup an explicit function of the nucleotide state. Because the nucleotide state is predominantly GDP, the results for properly parameterized models will be very similar.
(iii) We also introduced a computationally efficient scheme to relax the mechanical energy between chemical events, which can also be employed in future chemomechanical lattice models. Within the allosteric model, we achieve a better mechanical energy relaxation than previous models [15] with significantly less computational steps than a full Brownian dynamics simulation requires [18]. (iv) The idea of a feedback of mechanical forces onto the hydrolysis rate can also be applied in future chemomechanical lattice models: if hydrolysis leads to a weakening of lateral bonds, one could expect mechanical strains that favor weakening of lateral bonds also to favor hydrolysis.
In the future, our model could be extended to also include regulating TIP+ proteins [79] for which different mechanisms of how they influence MTs could be implemented. Comparing the results of such simulations with experimental data could help to develop a mechanistic picture of the action of these proteins. Another future extension is MT polymerization under force [80,81]. So far, polymerization under force has been investigated using chemical models [56,[82][83][84][85]; the influence of an external force on the microscopic level, in particular the detailed dynamics of catastrophe events and the catastrophe-triggering configurations is unknown.

Supplementary Material
Microtubule structure and energy At the minus end of the MT, each protofilament p starts with an alpha-tubulin at with the mean MT radius R MT = 10.5 nm, such that the seam is located between the 13th and the 1st protofilament.
Using the direction vectors the plus end position p(p, d, t) of any tubulin monomer can be calculated by adding all direction vectors to the minus end vector, The protofilament length that will be used to calculate the growth and shrinkage velocities is the maximum z-coordinate of all tubulin monomers within the protofilament: For straight and slightly curved protofilaments, max (p) = p(p, d(p), 2) · e z is the position of the plus end of the protofilament. For strongly curved protofilaments exhibiting a ram's horn and curling backwards, the length can exceed the z-coordinate of the terminal beta-tubulin, max (p) > p(p, d(p), 2) · e z .
In order to define the lateral bond energies between tubulin dimers in neighboring protofilaments, we need to introduce interaction points, where the harmonic springs of the lateral bonds attach. The lateral interaction points are located at the edge of the upper base at p(p, d, t) + c(p, d, t) and p(p, d, t) − c(p, d, t) with the connection vector The connection vector points from the edge of the upper base of tubulin monomer (p, d, t) to its next neighbor (p + 1, d, t) and is used to define the harmonic spring energies of the lateral bonds.

Detailed dynamics of single catastrophes, rescues and dips
The MT growth trajectory shown in Fig. 10(A) in the main text has two significant events: a dip at t sim = 1.2 min and a catastrophe at t sim = 6.85 min. To determine whether newly run simulations with configurations from the initial simulation as starting points qualitatively follow the original simulation, we chose the following criteria for this particular simulation: • To identify whether a new simulation reproduces the dip, we checked that at the end of the simulation, the MT length is at most 400 nm shorter than the original simulation. As the relevant question for the dip is whether it could actually result in a catastrophe, we are only interested in the new simulations being shorter than the original simulation, hence there is no upper limit on the MT length difference.
• To identify whether a new simulation reproduces the catastrophe, we took the last entry in the MT length log and checked if its time t sim differs less than 10 s from the original simulation end and if MT < 200 nm, i.e., if the MT continued shrinking and depolymerized (almost) completely.
For each initial configuration, we ran 20 new simulations and calculated the fraction of simulations that fulfilled the criteria above. These fractions are the probabilities for the original growth path at different points in time and they are shown color-coded in the insets (A.D) and (A.C) in Fig. 10(A).
For the analysis of Fig. 13(A) in the main text, where the hydrolysis rate is coupled to mechanics, we used the same criteria. Figure 10(C) in the main text contains three significant events: a dip at the very beginning, a catastrophe at t sim = 9.15 min, and a rescue at t sim = 9.54 min. For this simulation, we used the following criteria: • For the dip, we used the same criterion as before.
• For the catastrophe, we checked that after 15 s (or 10 s for configurations after the catastrophe), the MT was at most 400 nm longer than the original MT at the same point in time.
• For the rescue, we took the last entry in the MT length log of the new simulations and checked if t sim > 59 s, i.e., if the new simulation finished due to the time constraint and not due to the MT having vanished meaning that the rescue actually happened.
Again, we ran 20 new simulations for each of the three initial configurations and obtained the probabilities for the simulation to follow the original growth path, which are shown as color code in the insets (C.D), (C.C), and (C.R) in Fig. 10(C). The color coding reveals that the dip in (C.D) is actually better characterized as a catastrophe immediately followed by a rescue. In addition, Fig. 10(B) shows snapshots of the MT configuration at characteristic points in the growth path, for example, before, in the middle, and after the catastrophe and rescue events. Figure 15 shows additional data on the catastrophe and rescue events from Fig. 10 in the main text (catastrophes from insets (A.C) and (C.C) and rescue from (C.R)). We show the mean protofilament length MT together with additional information on the length fluctuations σ of the individual protofilament lengths and the average GTP-cap length N cap . When catastrophes become unavoidable, the cap length has shrunken to around N cap ∼ 2, for rescues a cap length around N cap ∼ 4 seems necessary. Length fluctuations are also increased if catastrophe or rescues are triggered.  Fig. 10 in the main text is used and gray is used to show additional data before and after the highlighted parts in the insets of Fig. 10.

Determination of catastrophe and rescue rates
For the determination of catastrophe and rescue rates in Fig. 8 in the main text we employed the following algorithm. First a MT trajectory MT = MT (t sim ) is classified into growth and shrinkage intervals using a greedy threshold value ∆ 1 = 50 nm. We start at the beginning of an interval at t sim = t 0 and increase t sim searching for a suitable end t 1 of an interval. If MT (t 1 ) − MT (t 0 ) ≥ ∆ 1 , i.e., if a MT has grown more than ∆ 1 , the interval [t 0 , t 1 ] is classified as growth interval. Likewise, if MT (t 1 ) − MT (t 0 ) ≤ −∆ 1 , i.e., if a MT has shrunken by more than ∆ 1 , the interval [t 0 , t 1 ] is classified as growth interval. All plateaus, where the length changes by less than ∆ 1 are "absorbed" into surrounding growth or shrinkage intervals. This part of the procedure gives a complete classification into a (not necessarily alternating) succession of growth and shrinkage intervals.
If the (n-1)-th interval is a growth (shrinkage) interval and the n-th interval a shrinkage (growth) interval the n-th interval is marked as possibly containing a catastrophe (rescue).
Then we continue with the second part of the algorithm, where we employ a less greedy threshold value ∆ 2 = 300 nm. We search for a catastrophe time t c and an enclosing interval [t c− , t c+ ] in a potential catastrophe containing interval according to the following steps: 1. We find a t c− < t c with t c − t c− < 50 s and MT (t c ) − MT (t c− ) ≥ ∆ 2 , i.e., the MT grows by ∆ 2 within the previous 50 s or less. We select the largest t c− fulfilling these criteria.
2. We find a t c+ > t c with with t c+ − t c < 50 s and MT (t c ) − MT (t c+ ) ≥ ∆ 2 , i.e., the MT shrinks by ∆ 2 within the next 50 s or less. We select the smallest t c+ fulfilling these criteria. 3. Among all possible t c in the potential catastrophe containing interval, we choose the value producing the smallest enclosing interval [t c− , t c+ ] according to steps 1 and 2.
For rescues, we proceed analogously.
Two exemplary simulation trajectories that we analyzed with the algorithm are shown in Figure 16.
The algorithm identifies the points in time where catastrophes and rescues happen and, thus, also gives access to the times ∆t gr,k that a MT grows before the k-th catastrophe and times ∆t gr,k that a MT shrinks before the k-th rescue. Catastrophe and recue rates are obtained as inverse of the averaged average growth and shrinking times, Theoretical spatial GTP distribution We employ a similar approach as Ref. [63] and consider a one-dimensional MT (or single protofilament approximation) with polymerization rate k on , effective depolymerization ratek off , and hydrolysis rate k hydr (for all tubulin dimers, including the terminal one). In the steady state, the probability of the i-th tubulin dimer counted from the plus end (i = d(p) − d + 1) to be a GTP tubulin dimer is given by with i = 1 referring to the tubulin dimer directly at the plus end. As we are measuring the normalized probabilitỹ p i with ∞ i=1p i = 1, we have to compare our simulation results with While we have a constant polymerization rate k on and hydrolysis rate k hydr (with hydrolysis is not coupled to mechanics), there is no clear mapping of the effective one-dimensional depolymerization rate k off to our three-dimensional modelling because lateral bond formation and rupture results in an effective depolymerization rate. When comparing our simulation results with the theoretical prediction (30), we are using k off as a fitting parameter and using k on and k hydr from the simulation.  . Once c tub is sufficiently large so that the MTs can reach a steady state of growth, the prediction by the one-dimensional theory matches the simulation data.
In the main text, we give for the probability p GTP (x) to find a GTP-dimer at distance x = d(p) − d from the tip. (31) is a continuous version (with a continuous x ≈ i − 1) of the discrete master equation for p i that leads to the above result (29).
Again, a direct comparison with our data is not possible because of the unknown effective one-dimensional depolymerization ratek off . However, (31) can be rearranged for an explicit expression fork off so that we can calculatek off (x): It should be noted thatk off (x) is not the depolymerization rate of layer x = d(p) − d but the depolymerization rate of the last layer calculated using the data from layer x. If our data can be described by (31), we expect k off (x) to be independent of x. Figure 18 showsk off (x) for both a constant hydrolysis rate and hydrolysis coupled to mechanics. The derivative was calculated using the symmetric derivative, except for the first and last values, were a forward or backward derivative was used. Ignoring numerical issues due to the discrete derivative and insufficient data statistics for larger x values,k off (x) is sufficiently independent of x showing that (31) can also be used to describe the results of the simulations in which hydrolysis is coupled to mechanics.  Values withk off < 0 s −1 andk off > 100 s −1 are cut off here as they are due to insufficient data statistics.

Minimization time comparison between constant hydrolysis rates and hydrolysis rates coupled to mechanics
To gain further insight into the additional amount of execution time required by simulations in which hydrolysis is coupled to mechanics, Figure 19 shows a comparison of average times for minimization after each of the five possible chemical events polymerization, depolymerization, bond formation, bond rupture, and hydrolysis. While the average minimization time after polymerization, depolymerization, bond formation, and bond rupture is not affected by mechanical feedback, the average minimization time after hydrolysis increases in the presence of mechanical feedback. In both cases, at a certain point that roughly matches the point in time when the porous cap length N pcap reaches it steady state, the average minimization time after hydrolysis does not change significantly anymore either. For the two simulations shown in Figure 19, minimizations after hydrolysis take four times longer for hydrolysis coupled to mechanics than for a constant hydrolysis rate; the reason for this increase is that the porous cap length is on average more than six times longer such that minimization has to be executed for up to six times more layers if a dimer deep in GDP-body is hydrolyzed. In both cases, the simulation spends around 98 % of its time during minimization. Of that total minimization time, the simulation shown in Figure 19(A) uses roughly 30 % for minimizations after hydrolysis events, while simulation Figure 19(B) uses about 66 % for minimizations after hydrolysis events even though the difference between the percentage of hydrolysis events from all events only increased by little more than 0.1 percentage points.

Analysis of dilution simulations
For the determination of the delay time ∆t delay after GTP-tubulin dilution at time t dil , we employed the following algorithm.
First, we fit a linear growth law MT (t sim ) = v gro t sim + MT (0) to the MT length data for simulation times determines the catastrophe time τ cat > t dil as intersection point with the dilation plateau MT (t sim ) = MT (t dil ), which we fit for t dil < t sim < τ cat . The delay time is given by ∆t delay = τ cat − t dil .
Two exemplary simulation trajectories that we analyzed with the algorithm are shown in Figure 20.      Table 2 in the main text, we have created a sequence in which the following GDP-dimers in protofilament p = 3 were exchanged with GTP-dimers to measure the bending angles of the GTP-dimers: (B) (3, 31), (C) (3, 32), (D) (3,33), (E) (3,34), and (F) (3,30). The highlighted intervals show all tubulin monomers that bend inward due to all previous changes from GDP to GTP.