Atomistic Basis of Microtubule Dynamic Instability Assessed Via Multiscale Modeling

Microtubule “dynamic instability,” the abrupt switching from assembly to disassembly caused by the hydrolysis of GTP to GDP within the β subunit of the αβ-tubulin heterodimer, is necessary for vital cellular processes such as mitosis and migration. Despite existing high-resolution structural data, the key mechanochemical differences between the GTP and GDP states that mediate dynamic instability behavior remain unclear. Starting with a published atomic-level structure as an input, we used multiscale modeling to find that GTP hydrolysis results in both longitudinal bond weakening (~ 4 kBT) and an outward bending preference (~ 1.5 kBT) to both drive dynamic instability and give rise to the microtubule tip structures previously observed by light and electron microscopy. More generally, our study provides an example where atomic level structural information is used as the sole input to predict cellular level dynamics without parameter adjustment.


INTRODUCTION
Biological processes occur over a wide range of spatial and temporal scales, spanning from angstroms and picoseconds at the atomistic level to micrometers and minutes at the cellular level. While rapid advances in structural biology are yielding unprecedented insight into biological structures at the angstrom scale, converting this information into quantitative prediction of cellular level dynamics on the time scale of minutes remains a major challenge. Fortunately, several multiscale modeling approaches have been developed to address this challenge. 5,6,11,16,46,52,56,62,71,73 For example, association kinetics of small biochemical ligandreceptor systems have been explored computationally by a multiscale approach that matches experimental values, connecting molecular dynamics (MD) simulations to Brownian dynamics (BD) simulations with milestoning theory. 70 Another example of a multi-scale framework initiating from atomistic details studied actin filaments by integrating all-atom MD simulations and coarse-grained (CG) techniques to study the impact of the nucleotide state on the filament's conformation, and the hydrolysis rate constant. 71,75 These methods, although bridging information across two or three scales temporally and spatially, do not yet provide a single framework to connect across length-time scales from atoms to cells.
Microtubules, are a prime example of a complex biological system with spatial and temporal scales ranging from atomistic phenomena such as nucleotide hydrolysis to cellular level organization and behavior that mediate cellular functions such as cell division and migration. 29 Microtubules self-assemble from abtubulin heterodimers that form both longitudinal and lateral bonds to build a hollow cylinder. 14,48,72 Tubulin heterodimers have conserved protein structures with defined binding zones which makes them ideal to study at the atomistic-molecular level. 24,41,51 A key feature of microtubule assembly is the phenomenon of ''dynamic instability,'' where the highly dynamic plus-end switches abruptly and stochastically between extended phases of net growth and net shortening, which depends on the guanosine triphosphate/diphosphate (GTP/GDP) nucleotide state in the b-tubulin subunit. 4,32 Whereas microtubules grow via net addition of GTP-tubulin subunits, the GTP-tubulin soon hydro-lyzes within the microtubule lattice to GDP-tubulin resulting in a so-called ''GTP cap'' that stabilizes the growing microtubule. 18,19,55,57 When the GTP cap is lost through a combination of GTP hydrolysis and stochastic loss of GTP subunits, the labile GDPtubulin core of the microtubule is exposed and the microtubule undergoes ''catastrophe'' followed by rapid shortening. 48 Because dynamic instability is essential for mitosis and cell migration, a key question is what are the fundamental thermodynamic and mechanical differences between the GTP and GDP states of the tubulin subunit that lead to net growth and shortening, respectively.
Previous studies of microtubule dynamic instability posited that the energetic difference between the two nucleotide states of tubulin in models of microtubule assembly 13,21,25,34,44,45,47,49,59,60,67,68,76,77 is due to differences in the tubulin-tubulin lateral bond strength, 25,44,45,47,60,68 tubulin-tubulin longitudinal bond strength, 25,44,77 the bending preference or flexibility, 13,21,34,49,59,67,76 or indirectly dependent on the lateral bond as a result of nucleotide-dependent bending strain energy. 7 At the level of atomic structure, a recent cryo-electron microscopy (cryo-EM) study showed that GDP, GDP-Pi and GTPcS all have compacted lattices (~82 Å ) compared to the GMPCPP extended lattice (~83.7 Å ), 44 consistent with conclusions from earlier cryo-EM studies. 1,77 This study also suggested that the compaction in GDP-lattice causes stronger longitudinal interactions and perturbs the lateral bond between the protofilaments, thus resulting in splayed protofilament ends and microtubule depolymerization. This view was recently challenged by a study 20 showing that MT lattice expansion is not related to the GTP-state, but rather is induced by the GMPCPP-state and that the GDP-and GTP-lattices have inherently similar dimer compactions. However, in this, and the other structural studies cited above, the conclusions drawn from cryo-EM studies regarding the mechanism of dynamic instability are based on the assumption that the contacts of the residues remain stationary as captured by the cryo-EM structure.
To gain insight into the atomistic dynamics, MD simulations have been performed using the protein data bank (PDB) structures as the sole input. Using the structure of tubulin protofilaments as a function of nucleotide state, a multiscale approach was taken to integrate the results of 100 ns equilibrium MD simulation with a CG-MD analysis, 25 representing a molecular group of atoms as one particle. 2 This study concluded that GDP-microtubules are less stable because GDP-tubulin has weaker longitudinal and lateral contacts in the lattice compared to GTP-tubulin. However, this conclusion was drawn from MD anal-ysis of contact maps and CG analysis of bond strength and lengths, while more recent work 30 using 10 allatom MD replicates per nucleotide state combined with BD and thermokinetic modeling revealed no evidence of nucleotide dependence of the lateral bond. More generally, as stated by the authors, 25 analysis of residue contact map and CG bond strength and length is not a measure of the actual binding strength of tubulin dimers; rather, an energy landscape of the lateral and longitudinal interaction of the dimers needs to be calculated. In addition, while a number of other MD studies have contributed toward an atomistic understanding of intradimer bending mechanics, 34,59 to our knowledge, previous studies have not examined the free energy of interdimer bending as a function of its nucleotide besides equilibrium trajectory analysis. 21,42,49 Finally, previous studies have not estimated the nucleotide dependence of the strength of the longitudinal bond, generally regarded as stronger than the lateral bond. Thus, an integrated multiscale dynamics approach is needed to understand the atomistic basis of microtubule dynamic instability, and, more generally, to connect atomic structures to cellular level behavior.
To address this challenge, we performed simulations at multiple scales where the output of one scale was used as the input to the scale above it. Using published crystal structures as inputs, we performed MD simulations of tubulin-tubulin interactions to obtain an estimated longitudinal potential of mean force (PMF), with which we then simulated tubulin addition and loss from the microtubule lattice via BD to obtain estimated k on , k off , and DG 0 . Using these k on and k off values as inputs, we used thermokinetic and mechanochemical modeling to predict microtubule dynamic instability at the scale of micrometers and minutes. Thus, we were able to simulate the cellular dynamics using the published crystal structures as the only inputs with no adjustable parameters. We find that the GTP-tubulin longitudinal bond has a stronger PMF than GDP-tubulin by DU long % 6.6 ± 2.8 k B T, which translates into a standard Gibbs free energy difference of % 4 ± 0.5 k B T. While this difference largely explains the atomistic basis of dynamic instability, we further found via thermokinetic and mechanochemical modeling that an outward bending preference of GDP-tubulin relative to GTP-tubulin contributes to dynamic instability and is necessary for realistic tip structures (~1.5 k B T). Overall, our multiscale approach (Fig. 1) presents a methodology, using microtubule dynamic instability as a specific example, for using atomistic structural information to predict cellular level behaviors 63 without parameter adjustment.

Multi-scale Approach
Our multi-scale framework starts from atomistic structural information of tubulins to obtain the ensemble-averaged potentials of lateral 30 and longitudinal interactions of tubulin, as the two heterodimers mainly interact with each other laterally and longitudinally.
The potentials obtained from MD simulations as well as dimer dimensions from tubulin structures were then used as inputs to our BD simulations to model Brownian motion of these heterodimers in solution as well as their assembly into a preassembled MT lattice. Since the time-averaged potentials obtained from MD simulations are equivalent to the canonical ensembleaveraged potentials of tubulin dimers if there is ergodic convergence, 9,12,54 we used umbrella sampling to ensure sufficient sampling of the ensemble. In our multiscale framework, each model is separated by different time scales, with the MD simulation ranging from~30 to 100 ns (PMF vs. unconstrained simulations), the typical Brownian dynamics simulation ranging from 10 ls to 0.5 s, and the thermokinetic modeling ranging from 100 ls to 20 min.
The MD potentials (U) along with dimer dimensions from tubulin structures were then used as inputs to our BD simulations to simulate Brownian motion of the heterodimers and their assembly into a preassembled MT lattice, based on the following equations (Eqs. 1-3). 8 First, since the interaction zones were not perfectly aligned at all simulation times, the MD potentials were time-averaged in a BD simulation while a subunit was within the binding zones to obtain an intrinsic bond energy DG 0 B À Á . In Eq. 1, m is the total number of unbinding events simulated and n is the number of steps taken before unbinding for an unbinding event.
Next, the standard Gibbs free energy change of tubulin binding DG 0 À Á was defined based on the sum of the intrinsic bond energy and a rigid-body entropic penalty of binding DG 0 S À Á (Eq. 2). Last, the standard Gibbs free energy change was related to the kinetic rates of assembly (Eq. 3), where k on is the association rate constant (M 21 s 21 ) and k off is the dissociation rate constant (s 21 ). Finally, the BD-calculated kinetic information of tubulin assembly as well as binding potentials from MD simulations were used in our previously developed pseudo-mechanical thermokinetic 68 and mechanochemical model 67 to capture MT dynamic behavior and detailed MT tip structures. In the pseudo-mechanical model, the hydrolysis of GTP-tubulin to GDP-tubulin was incorporated to model transitions as well as the mechanical effects of hydrolysis, i.e. addition of a dimer was still necessary to stimulate hydrolysis of the dimer below, without the instantaneous coupling between the two. In the mechanochemical model, the bending of the subunits is explicitly modeled in 3-dimensional coordinates, shown in Eq. 4, where k curl is the outward curling spring constant, and u is the angle between the preferred and actual orientation of the dimer.

Molecular Dynamics Simulations
MD Simulations of all systems were run using NAMD 2.10 software package 37 using the CHARMM 36 force field 69 for parametrization. VMD 1.9 33 was used for visualization and trajectory analysis. Longitudinally-paired tubulins, both dimers with GDP-or GMPCPP-nucleotides, were extracted from the published cryo-EM structures of microtubules by Zhang et al. 77 (PDB ID 3JAS, 3JAT). GTP-tubulin structure was built based on GMPCPP-tubulin structure by swapping GMPCPP out for a GTP, which was necessitated by lack of a true GTP-tubulin structure in a microtubule lattice. The protein complex along with the nucleotides were all parametrized using the CHARMM-GUI interface. 35 Each simulation system was initially energy minimized for 12,000 steps using the conjugate gradient algorithm, and then they were solvated in TIP3P water, 36 using a 10 Å margin from each side. The simulation systems were neutralized with MgCl 2 ions at 2 mM concentration based on physiologically-relevant salt concentrations.
The solvated systems were heated to 310 K for 1 ns using a Langevin thermostat, 26 and then run in an NPT ensemble (T = 310 K and P = 1 atm). The simulations were followed by a total production run of 350 ns for each system (after 50 ns equilibration). All simulations were run with 2 fs time step and a cutoff radius of 12 Å for van der Waals interactions, using Particle Mesh Ewald (PME) for long range non-bonded interactions. 15 The equilibrium run trajectories were stored every 3000-time steps (6 ps). Root-meansquare deviation (RMSD), root mean squared fluctu-ations (RMSF), hydrogen bonds (with the criteria of donor-acceptor distance of less than 3 Å and the angle cutoff of 20°) and salt bridges were calculated using manually written tcl scripts and plugins available in VMD. The buried solvent-accessible surface area was calculated as previously described. 30 To simulate the effect of lateral neighbors in the microtubule lattice without increasing the number of atoms, we identified lateral bond residues from Ref. 30 and longitudinal bond residues from our unconstrained simulations and applied a harmonic constraint to those atoms. The stiffness of the constraints of lateral and longitudinal neighbors were chosen according to the stiffness of the lateral PMF, 30 j lat = 1 kcal/mol Å 2 and unconstrained longitudinal PMF, j long = 0.6 kcal/mol Å 2 .
For free energy calculations, we employed the umbrella sampling method 65 combined with weighted histogram analysis method (WHAM) 27,40 to be able to sample the ensemble sufficiently and have independent simulations that each can be run for longer sampling time in parallel, considering the large number of atoms.
A PMF, a free energy landscape as a function of a specified reaction coordinate, was obtained for each nucleotide case. The reaction coordinate was defined as the longitudinal center of mass (COM) to COM distance of the dimers, without any rotations, since that is the most probable path of longitudinal unbinding of the dimers according to our BD simulations (Fig. S5). In addition, choosing the COM as the metric for measuring inter-dimer distances allowed us to consistently mesh our MD simulation potentials to our coarse-grained BD simulations. We believe that the calculated effective COM-to-COM PMFs are reasonable representatives of many-dimensional free energy surfaces for each protein-pair, which is very challenging to calculate. 62 However, due to the constraints of the parameter space of our multi-scale model and how the model outputs should eventually be consistent with the in vitro/vivo MT dynamics, we believe that the MD-derived lower-dimensional CG potentials are not far from the actual interaction potentials of tubulins needed for calculating the thermodynamic properties via BD. 22,52 In addition, the entropic information is obtained via BD simulations and the potentials are further revised. Our approach relies on targeting macroscopic behavior (thermodynamic data) to justify our CG path selection (COM-to-COM reaction coordinate) when compared to other ''bottom-up'' CG methods. 62 The bias potential stiffness was tuned to be 10 kcal mol 21 Å 22 to give sufficient overlap of the histograms of the windows. Fifteen windows were created, each being 1 Å separated from their nearest window. The reaction coordinate in free energy simulations was recorded every 200-time steps (0.4 ps).
To ensure the time convergence of the PMFs to the ensemble-average PMFs, we employed similar methodology in Ref. 30 to run multiple replicates with different initial conditions and to gradually increase window sampling time until the PMF change did not exceed a threshold of 1.5 k B T, determined by the Monte Carlo bootstrap error of the PMFs. A window sampling time of 30 ns was sufficient to produce a convergence in each PMF. The last 200 ns of the equilibrium run was used to choose equilibrated initial structures for creating replicates of umbrella windows for each system.
NVIDIA Tesla K40 GPUs were used to accelerate the simulations on the Mesabi cluster at the Minnesota Supercomputing Institute (MSI), University of Minnesota, and NVIDIA Kepler K80 GPUs were used on Comet and Bridges, Extreme Science and Engineering Discovery Environment (XSEDE) 66 dedicated clusters at the San Diego Supercomputing Center (SDSC) and Pittsburgh Supercomputing center (PSC).

Brownian Dynamics Simulations
Tubulin dimers' association and dissociation from a protofilament in microtubule lattice was simulated using the BD model of Castle et al. 8 The entropy corrected PMFs, with entropy calculated based on the rigid body Shannon entropy method as described in Ref., 30 were used as the input to the BD simulations. Simulations were run for two nucleotide cases and each was run for a total of 500,000 iterations of binding simulations and 20 to 20,000 iterations of unbinding simulations (total time of 0.1 to 1 s). For binding simulations, half-force radius was calculated for the PMF energy profiles and used as the binding radius in BD simulations. 30 For unbinding simulations, we used a separation distance criterion of R U = 11 nm, according to, Ref. 8 where the probability of rebinding is very low (p <0.01). BD of dimer's dissociation from a protofilament with two lateral neighbors was not simulated due to extremely high stability and long unbinding time (no unbinding event up to 0.1 s).

Thermokinetic and Mechanochemical Modeling
Microtubule assembly dynamics were simulated using a pseudo-mechanical thermo-kinetic modeling, as previously developed 68 and modified. 7 For all simulations, an in vitro parameter set was used (Table S4) with variable energy penalty for hydrolysis. Seed length and starting GTP layers were set to 2 lm. Onrate penalties of 2 and 10 were added for one and two lateral neighbor cases, respectively. 8 Microtubule tip structures were obtained using the mechanochemical model as previously described. 67 All the shortening microtubules were simulated starting from uncapped protofilaments with blunt tips, ran for 500 events. Model parameters were similar to Ref. 67 with modification of preferred bending angles and flexibilities. Higher flexibility in GTP-dimer was chosen according to bending angle data variance (Table S1A) as double the flexibility of GDP-dimer. As for bending preference, the preferred radial bending angle was selected as 22°, 67 and the preferred tangential bending angle was set to 11°, according to Table S1A, where the tangential mean value is almost half of the radial mean value.

Tubulin Dimers Reach Equilibrium in the Absence and Presence of Microtubule Lattice Constraints
To investigate the longitudinal interaction between tubulin heterodimers, an interaction essential to protofilament formation, we modeled a pair of tubulin dimers stacked longitudinally via MD simulations. Since tubulin structures were initially obtained from straight protofilaments 77 (Fig. 2a), we were interested to see how removing the lattice constraints would influence oligomer (dimer of heterodimers) bending and conformation. During the simulation, we monitored the RMSD of the backbone atoms to ensure global equilibrium (Fig. 2b). We also considered the possibility that microtubule lattice constraints affect the conformations of the dimers, mostly by confining their bending motions. These constraints are believed to be stored as strain energy that keeps the dimer's conformation straight in the lattice. To test this hypothesis, we simulated the effects of lateral neighbors and the bottom longitudinal neighbor as harmonic constraints on the atoms mainly involved in the lateral 30 and longitudinal bond (highlighted in red in Fig. 2c). As for harmonic stiffness, we first calculated the longitudinal PMF for unconstrained dimers and then used that along with our previously calculated lateral PMF 30 to estimate a stiffness for the bonds by fitting a harmonic potential around the potential well. We obtained stiffnesses of j = 1 and 0.6 kcal/mol/Å 2 (j = 0.7 and 0.4 nN/nm) for lateral and longitudinal bonds, respectively. In this way, we simulated the extreme case of the lattice constraints for a dimer at a protofilament tip, which is having two lateral neighbors and one longitudinal neighbor.
We found that these lattice constraints do not deform the protein structure, as evidenced by the plateau of the RMSD of the backbone atoms after~70 ns for both nucleotide states (Fig. 2d). Our simulations confirmed the two major bending modes for the oligomers, tangential and radial outward bending, regardless of the nucleotide or lattice constraints (Fig. S1A, Supplemental Movies S1 and S2), as previously shown for single tubulin dimers. 30,34 To compare our data to previously published tubulin MD simulations, 21,34 the last 300 ns of our trajectories were analyzed further in terms of RMSF and bending angle dynamics (Figs. S1-S3, Tables S1A, S1B, S2A, S2B), which confirmed behavior consistent with these previous studies. We avoid drawing strong quantitative conclusions about nucleotide dependence of bending an-gles based solely on these results since the angle data is highly variable and is not likely to be converged due to longer relaxation times. Nonetheless, we can conclude that our protein backbone structure is globally equilibrated, and the presence of the lattice constrains the fluctuations of the dimers, resulting in a straighter oligomer configuration, as evidenced by the reduction of the average backbone RMSD (Figs. 2b, 2d) and average structures obtained from the trajectories (Fig. S1B). The global equilibration allowed us to proceed to probe the non-covalent binding nature of the equilibrated longitudinal interface and the longitudinal bond strength, i.e. the longitudinal PMF.

GTP-and GDP-Tubulin Have Two Longitudinal Interaction Zones in Common But Differ on a Third Zone
High-resolution cryo-EM microtubule structures identified three main longitudinal contact zones between the dimers in a protofilament. 50 Three zones also reproduced the best estimates of experimental onoff rates of protein assembly when modeled via BD simulation. 8,53 To investigate whether tubulin oligomers would maintain their longitudinal contacts despite the initial relaxation of the lattice compaction observed in GDP-lattice upon equilibration (Fig. S4A-F), 77 we searched for residues involved (> 25% of the trajectory time) in H-bond and ionic interactions at the interdimer interface. Our computational results, shown in Table 1 and Figs. 3a and 3b, identified three nonsymmetrical longitudinal interaction zones in each equilibrated nucleotide state in the case where the bottom dimer is lattice-constrained: (1) S9 and H10-S9 loop with T5 loop, (2) H11-H11¢, H11¢, H11¢-H12 with H8-S7 and H4-S5 loops, and (3T) H6 with H10 in GTP-tubulin only and (3D) residues 2 to 4 in the asubunit's N-terminus with T2 loop in GDP-tubulin only. In contrast to the results of Nogales et al. (1999), the interaction of T7 loop with T2, T1, and the nucleotides (termed ''zone C'' by Nogales et al.) was less than 10% and 18% of the total trajectories in GDP and GTP states, respectively. This is due to the relaxation of the simulated dimers in water without the constraints of the lattice around all the dimers, as clearly observed in the compaction release of GDP dimers, which increases the distance between the nucleotide and interacting residues. Based on our results of constrained dimers (Fig. S4E), although the GTPand GDP-dimers had initial interdimer distances of 42.1 Å and~40.6 Å respectively, they both plateau to the same expanded interdimer distance at the end of 400 ns of our simulation (~42 Å ) (Fig. S4E). The comparison of our unconstrained dimer distances (Fig. S4F) to our constrained distances (Fig. S4E) showed that in the unconstrained simulations, both GDP-and GTP-dimers reached a similar higher interdimer distance at the end (~43.5 Å ) due to lack of constraints in solution (Fig. S4F). Our results agree with a recent study 20 in that both GDP and GTP-lattices have similar compaction, and differ in the final state of those lattice repeat distances (our results showed a final expanded lattice vs. a compacted lattice for both nucleotides). This discrepancy might be due to the fact that our simulation system of two dimers differs from a complete lattice structure in that it lacks the long-range constraints of the neighbors that leads to lattice-stored mechanical stress and has expanded due to thermal fluctuations.
These results also suggest that GTP-and GDP-dimers differ in one of their major longitudinal interaction zones, potentially causing different longitudinal bond strengths as a function of nucleotide state. To identify the relative contribution of various non-covalent interaction types at the longitudinal interface, we calculated the number of H-bonds, salt-bridges and extent of hydrophobic interactions quantified as buried solvent accessible surface area (BSASA) involved (persistence > 25% of the trajectory time) in the interdimer interaction, as shown in Figs. 3c and 3d. We found that the longitudinal interface is dominated by hydrophobic interactions (~25-30 nm 2 BSASA), which is qualitatively distinct from the lateral interface which is dominated by H-bonds and ionic interactions. 30 As expected, lattice constraints make the number of longitudinal contacts higher by preventing the dimers from bending considerably. However, in contrast to previous studies which argue that higher contact numbers equals a stronger bond, 44 we note that this metric is not necessarily representative of the total longitudinal bond strength and is more an indicator of the quantity of longitudinal contacts without knowing their relative strengths. While these analyses inform our qualitative understanding of the nature of the bond and the relative shifts that occur upon nucleotide hydrolysis, assessing the quantitative strength of a bond requires calculating bond potentials (PMFs) and comparing the potential well-depths.  The Longitudinal Bond is Significantly Stronger for GTP-Tubulin Than it is for GDP-Tubulin The differential strength of the longitudinal interactions between GTP-and GDP-tubulin heterodimers could contribute significantly to the phenomenon of dynamic instability. Therefore, we estimated the PMF of the longitudinal COM to COM distance for GDPand GTP-tubulins, as free oligomers or lattice constrained oligomers. Our distance-based reaction coordinate choice, although limited to low dimer rotational degrees, was found to have the highest efficiency of longitudinal binding (Fig. S5) and is more common in the literature. 28 We used umbrella sampling with WHAM 27,40 (see ''Methods'') to efficiently sample the ensemble, which has the advantage over traditional time-averaging techniques of overcoming energy barriers inaccessible to traditional MD simulations. To make conclusions about bond strengths on the scale of k B T, we found it necessary to run multiple replicates (Table S3). As shown in Figs. 4a and 4b, our longitudinal PMFs clearly show a stronger GTP tubulintubulin interaction compared to GDP, regardless of the presence of lattice constraints. We found that the potential well-depths for the two nucleotide states differ by~6.6 k B T (p < 0.02; Table 2).
To investigate the stiffnesses of the potentials, we calculated the half-force radius (described in Ref. 30) in Table 2 and, as predicted by constraints confining the atom movements, all lattice-constrained potentials were slightly stiffer. In addition, the GDP-tubulin potentials were stiffer than the GTP-tubulin potentials. Since the potential well-depths were not statistically different between the constrained and unconstrained cases (p > 0.9), we calculated a mean longitudinal PMF for GDP-and GTP-tubulin, averaging all 15 replicates of constrained and unconstrained simulations together (Fig. S6). To test whether different nucleotide states had different potential minima after the averaging, we performed a statistical analysis of the minima locations. In the same simulation configuration, i.e., constrained or unconstrained, none of the potential minima locations for either nucleotide state was significantly different (p > 0.8), indicating that the GDP lattice compaction had relaxed into the GTP extended lattice. Overall, we conclude that the stronger longitudinal bond strength of GTP-tubulin compared to GDP-tubulin can explain why GTP-capped protofilaments would favor polymerization while GDP protofilaments would dissociate faster from a protofilament end due to weaker longitudinal bond strength.
Brownian Dynamics Simulations Imply a Free Energy Difference, DDG 0 long , of % 4

k B T Between Nucleotide States
To predict the influence of the PMF on the subunit addition-loss kinetics and thermodynamics as quantified by the standard Gibbs free energy difference between the GTP and GDP nucleotide states, we performed BD simulations to model dimers' diffusion and assembly kinetics on time scales up to 0.5 s, instead of only~1 ls via MD. Following our multi-scale approach, we calculated the full entropic penalty of longitudinal binding, which allows us to account for the fact that the stronger longitudinal bond (~25-32 k B T well depth) will have a higher entropic penalty compared to the weaker lateral bond (~11-12 k B T well depth; Hemmat et al. 2019). MD entropy-corrected potentials (Fig. 5a) were used as input to our BD simulations of microtubule lattice assembly (Figs. 5b, 5c). Binding and unbinding of an incoming tubulin dimer to the tip of a protofilament with zero or one lateral neighbor were then simulated. Our zeroneighbor simulation results, as depicted in Table 3 show that a longitudinal energetic difference DG 0 long of~3.5 ± 0.5 k B T exists between GDP-and GTPstates using the different input potentials obtained from the MD simulations for unconstrained dimers. To test the hypothesis that the slightly stiffer potentials from constrained MD simulations would change the BD model outputs, we ran another set of BD simulations using lattice constrained MD inputs. While the energetic difference mean value is higher for the constrained case, DG 0 long % 4.7 ± 0.6 k B T, it is not statistically different from the unconstrained value (p > 0.6). We conclude that GTP tubulins have stronger longitudinal bonds as measured by DG 0 long than GDP tubulins by 3.5 to 4.7 (% 4 on average) k B T.  p-values < 0.02 constrained dimers compared to unconstrained dimers of the same nucleotide state.

BIOMEDICAL ENGINEERING SOCIETY
To account for the stabilizing effect of a lateral neighbor, we ran BD simulations of dimer incorporation into a ''cozy corner,'' i.e. a protofilament tip with a longitudinal and a single lateral bond, using our lateral 30 and longitudinal potentials as inputs. The results, summarized in Table 4, yield the estimated lateral bond value for GDP-and GTP-states in the microtubule lattice, from comparing the total energy DG 0 À Á value to the zero-neighbor case. The lateral bond strength, as expected from the lateral potential inputs from our previous study, 30 is found to be nucleotide independent, based on the values for unconstrained dimers DG 0 GDPÀlat ¼ À 4:2k B T and DG 0 GTPÀlat ¼ À 4:9k B T, and DG 0 GDPÀlat ¼ À 4:0k B T and DG 0 GTPÀlat ¼ À 4:2k B T for constrained dimers, which are not statistically different (p > 0.7). In addition, the association rate constants, k on , and dissociation rate constants, k off , of GDP-tubulin are generally slower for the constrained input potentials as compared to the unconstrained input potentials. This is because as the potential becomes stiffer, it becomes harder for the dimer to diffuse in and out of the potential well. 8 However, for GTP-tubulin, a significant change in the kinetics is not expected as the stiffness of the longitudinal potential is not significantly different in constrained vs. unconstrained simulations. Based on our previous BD simulations of bimolecular association in crowded environments, we do not expect the kinetics and thermodynamics to be significantly altered in in vivo environments due solely to macromolecular crowding. 61 Altogether, these BD results give us insight into the possible mechanisms by which GDP-tubulin differs from GTP-tubulin in its stabilizing behavior in microtubule assembly that can be used in larger scale thermo-kinetic modeling of microtubule dynamic behavior. Specifically, these results are consistent with a model where a potential well-depth difference DU long À Á of~6.6 k B T in the input longitudinal potential creates a 3.5-to 5-fold decrease in the on-rate constant, and an 11-to 22-fold increase in the off-rate constant of GDP-tubulin compared to GTP-tubulin, depending on their lateral neighbor case.
Since a stronger longitudinal bond will favor oligomer formation, we tested the hypothesis that GTPtubulin oligomers are potential intermediate structures.

BIOMEDICAL ENGINEERING SOCIETY
We calculated the average length of oligomers (in number of subunits, including single subunits) that is resulting of the stronger GTP-tubulin longitudinal bond 7,31 (Eq. 5), and showed that even with the stronger bond DG 0 long ¼ À 9:5 and at [Tub] of 5.6 lM, GTP oligomers are not significant in solution (average length of 1.3 subunits).
Thermokinetic Modeling Identifies a Preferred Bending Angle Difference to Reproduce Experimental Microtubule Tip Structures Collectively, the MD and BD simulations show that GDP-tubulin pays an energetic penalty through its weaker longitudinal bond while the lateral bond strength remains nucleotide-independent. These results led us to hypothesize that the energetic penalty due to GTP hydrolysis in the microtubule lattice must either exist on the longitudinal bond only or on both longitudinal bond and bending flexibility and/or angle preference. Unfortunately, we do not yet have detailed structural information regarding the bending flexibility or angle preference of the dimers, and previous studies have not reported a potential for various bending modes of tubulin. To eliminate one of our possible hypotheses for dynamic instability, we ran our previously described thermokinetic model 7,30,68 (Fig. 6a) with DDG 0 implemented on the longitudinal bond alone or in combination with lateral bond energetic differences due to possible nucleotide-dependent bending flexibility/mechanics, which allowed us to predict microtubule net assembly rates, and tip structures during net polymerization and depolymerization. We first ran a single state thermokinetic model of microtubule assembly using our BD estimated on rate constants and in vitro parameter set (Table S4). The heatmap of the resultant microtubule net assembly rates (Fig. 6b) indicates that our MD-BD estimated free energies are in reasonable agreement with published in vitro assembly rates. 23 Based on our previous analyses, microtubules undergo dynamic instability with zero net rate overall (averaged over multiple rounds of dynamic instability) with an apparent microtubule tip net assembly diffusion coefficient of~0.8-16 9 10 4 nm 2 /s, extracted from experimental microtubule length variance measurements. 7,23 As shown in Figs. 6c and 6d, the calculated microtubule net assembly rates and apparent diffusion coefficients as a function of DDG 0 lat and DDG 0 long are consistent with those reproducing dynamic instability (highlighted inside the dashed line borders). Dynamic instability is found over a range of lateral energetic penalties of less than 1.5 k B T combined with the MD-BD estimated DDG 0 long of 3.5 ± 0.5 k B T. Additionally, using the parameters estimated, our model shows catastrophe frequency over a range of tubulin concentration (Fig. S7), consistent with previous in vitro studies. 10,72 Interestingly, if the energetic penalty of hydrolysis that destabilizes the lattice is only paid through the longitudinal bond, microtubules would grow and shorten with protofilaments having very low protofilament length standard deviation, r tip (Fig. 6e), i.e. ''blunt'' tips (while a larger r tip reflects ''tapered'' tips). 17,23 However, if the penalty is applied to both longitudinal and bending angle preference and/ or flexibility (i.e. the lateral bond in our model), tapered tips with variable protofilament lengths and assembly dynamics appear in the simulated microtubules, manifested as a larger r tip (Fig. 6f). This finding is consistent with previous in vitro 23 and in vivo 17 experiments that reported highly variable tapered tips for microtubules by high resolution imaging and microtubule tip tracking. It is also in line with previous work that suggested there is a flexibility difference between the two nucleotide-states, with GTP being softer at its intra-and inter-dimer interface, 21,34,42 or a bending preference 67,76 with GTPprotofilaments growing less curved or nearly straight compared to the extensive outward peeling observed for GDP-protofilaments.  (Table S4). The netrates are obtained using a pseudo-mechanical model without hydrolysis. Blue and orange circles show our predicted areas where GDP-and GTP-tubulin are located, respectively. Black circle shows the reference point for dynamic instability.  Since both bending flexibility and angle preference are directly related to mechanics of a microtubule lattice, using a more detailed model where mechanics are directly modeled can help to pinpoint the possible mechanisms of dynamic assembly. The previously developed mechanochemical model of microtubule assembly 67 was used in our study to investigate whether a difference in bending angle preference (in both radial and tangential directions) and/or in bending flexibility combined with a DDG 0 long will recreate experimentally observed assembly rates as well as tip structures. We first tested the possibility of GDPtubulin having higher radially outward bending preference (h GDP x ¼ 22 , chosen according to electron micrographs evidence of a curling outward GDP vs. straight GTP protofilaments, 10,39,43,67 Fig. 7a). The results show microtubule dynamic instability (Fig. S8A), blunt growing microtubules (Fig. S9A) and three major types of shortening tips: blunt, splayed and tapered (Fig. 7a), consistent with experimentally observed EM tip structures. 10,43 Tangential bending preference was then investigated as a possible mechanism by itself (Fig. 7b) or combined with radially bending kink (Fig. 7c). Our model showed that dynamic instability is disrupted, microtubules rarely undergo catastrophe, and mostly grow in the case of a tangential bending kink only, with h y ranging from 11 to 22° (Figs. S8B, S9B). However, combined with radial outward bending, microtubules show dynamic behavior (Figs. S8C, S9C) and tip structures similar to Fig. 7a, meaning that a combination of bending preference both radially and tangentially is plausible for microtubule assembly.
Next, as argued in previous studies, 21,34 we assessed the possibility of a flexibility difference between different nucleotide states as the main cause of dynamic instability (Figs. 7d, 7e). Interestingly, assuming that GTP-dimer is more flexible than GDP-dimer (~2-fold based on Table S1A) and both dimers prefer a radially kinked or straight conformation (h x = 22°, 0°) did not result in normal microtubule dynamics, i.e. abrupt microtubule length changes with extended periods of dynamic pause for both dimers kinked, and single state assembly with microtubules growing consistently without any shrinkage for both dimers preferring to be straight (Figs. S8D, S8E, S9D, S9E). Lastly, a combination of flexibility difference along with a differential radial bending preference was applied to our model. Dynamic instability was observed (Figs. S8F, S9F) and shortening microtubules showed either splayed, blunt or tapered tip structures (similar to Figs. 7a, 7c).
These results highlight the importance of a radial outward bending angle preference for GDP-tubulin that exceeds that of GTP-tubulin in mediating microtubule dynamic assembly and tip structures. Although we cannot rule out that GDP-tubulin's radial bending preference might be accompanied by other mechanisms such as flexibility differences or tangential bending preference as well, we conclude that radially outward kinking of GDP-tubulin is necessary for the model to agree with the experimental dynamics and tip structures.

DISCUSSION
Here we investigated the fundamental atomistic and molecular mechanics underlying a complex biological phenomenon, microtubule dynamic instability, by using a multiscale approach, integrating structural, mechanochemical, and kinetic perspectives that span from atoms to cellular scales. Our MD results, using previously published structures of tubulin, show that the longitudinal bond, in contrast to the lateral bond, is nucleotide dependent and is~3 times stronger than lateral bond. However, we find that a longitudinal  bond difference is insufficient by itself to produce the experimentally observed tapered growing tip structures, and a nucleotide-dependent radial preferred angle is essential to recreate curling protofilaments commonly found at the tips of shortening microtubules and blunt tip structures in growing microtubules. Thus, by using this multiscale approach without parameter adjustment, we conclude that dynamic instability occurs primarily by weakening of the longitudinal bond (~4 k B T) and secondarily by outward curling between dimers (~1.5 k B T) upon GTP hydrolysis in the microtubule lattice. More generally, these results show how dynamic simulations can be used to leverage atomistic structural data to identify the mechanistic origins of cell-level dynamics that are important to cell behavior.
Although the possibility of a nucleotide-regulated bending flexibility has received much attention recently, 3,21,34,74 the suggested mechanism did not reproduce microtubule dynamics and predicted tip structures within our experimentally-constrained approach. An implication of this finding is that there is another important factor playing a role in maintaining dynamic instability, as observed in experiments. Nucleotide-dependent lateral and longitudinal bonds, with GDP-tubulin having a stronger longitudinal and a weaker lateral bond, suggested by a recent cryo-EM study 44 were also ruled out based on the results of multiscale model and our previous MD studies of the lateral bond PMF. 30,77 It was only the addition of a radial bending preference to our nucleotide-dependent longitudinal bond in our model that captured both predicted microtubule dynamics and tip structures. Our results agree with the findings of Rice et al. (2008) where soluble dimers have similar intradimer conformations in solution, and further we show that the intradimer bending angles remain small and nucleotideindependent upon lattice constraints. In addition, in the presence of the lattice constraints, an outward radial interdimer bending preference (% 22°) along the protofilament is modulated by the nucleotide state. As an alternate mechanism, a GTP-dimer with a stronger longitudinal bond, a weaker radially outward bending angle and higher dimer flexibility cannot be ruled out within our methodology, consistent with experimental boundaries.
Despite the fact that our results are in good agreement with experimental microtubule dynamic assembly measurements and previously estimated microtubule computational models, 67,68 it could be argued that our MD estimated potential energies can be affected both by the initial crystal structures, and by the selected reaction coordinate. The initial cryo-EM structures used in our study were obtained in the presence of kinesin, a tubulin dimer marker, 77 which has been speculated to affect GDP-tubulin longitudinal compaction, 58 or GMPCPP-lattice spacing. 78 We note that although longitudinal compaction in the lattice might be regulated by kinesin to a degree, the instant compaction release and plateauing after 100 ns in our MD simulations in the absence of kinesin justifies that any external effect is quickly relaxed (Figs. 2, S4). The fact that we obtain consistent PMFs in both constrained and unconstrained simulations further shows that our PMF calculation is reliable within the replicate-toreplicate variability (± 1.8 to 3.5 k B T SEM) and is not dependent on the initial conditions. Since we were unable to investigate all possible reaction coordinates, we chose the most probable distance-based path to investigate the bond energies between the dimers as determined by BD simulation with compliant longitudinal zones (Fig. S5), and found that a COM reaction coordinate can best mesh together the MD to BD scales. Additionally, it will be important in the future for our thermokinetic model to include minus end dynamics as another way to validate our modeling results as well as long-range effect of neighbors 38 (without ad hoc parametrization) and the seam interactions 64 to investigate their potential effects on the interaction energies.
The predictions made by our multiscale modeling approach show that although microtubule assembly dynamics have multiple parameters governing their behavior in the thermokinetic model, it is tightly regulated and only a limited range of parameters would agree well with the experimental results. Having a multiscale approach is significantly beneficial in narrowing down the predictions made from atomic/molecular level to cell-level behavior. This multiscale methodology can be a framework to predict the dynamics of other self-assembled polymers, especially those whose subunits are relatively rigid and ordered, such as F-actin, deoxygenated sickle hemoglobin fibers, amyloid fibrils, and virus capsids to make physiologically relevant predictions from atomic changes, such as the effects of hydrolysis, mutations, post-translational modifications, microtubule-associated proteins, drugs, and therapeutic agents.

ACKNOWLEDGMENTS
The authors thank Dr. Jonathan Sachs for advice and helpful discussions. This study was supported by National Institutes of Health under Award Number RF1-AG053951 and the Institute for Engineering in Medicine (IEM) award at the University of Minnesota to DJO. The authors acknowledge the Extreme Science and Engineering Discovery Environment (XSEDE), Comet system at the San Diego Supercomputing Center (SDSC) and Bridges system at the Pittsburgh Supercomputing Center (PSC) through allocation MCB160060, and the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for providing resources that contributed to the research results reported within this paper.
This work is based on the preprint submission on bioRxiv.org: ''Atomistic basis of microtubule dynamic instability assessed via multiscale modeling'', https://d oi.org/10.1101/2020.01.07.897439Copyright holder for this preprint is the author/funder and it is made available under a CC-BY-NC-ND 4.0 International license.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://crea tivecommons.org/licenses/by/4.0/.