Electronic Energy Migration in Microtubules

The repeating arrangement of tubulin dimers confers great mechanical strength to microtubules, which are used as scaffolds for intracellular macromolecular transport in cells and exploited in biohybrid devices. The crystalline order in a microtubule, with lattice constants short enough to allow energy transfer between amino acid chromophores, is similar to synthetic structures designed for light harvesting. After photoexcitation, can these amino acid chromophores transfer excitation energy along the microtubule like a natural or artificial light-harvesting system? Here, we use tryptophan autofluorescence lifetimes to probe energy hopping between aromatic residues in tubulin and microtubules. By studying how the quencher concentration alters tryptophan autofluorescence lifetimes, we demonstrate that electronic energy can diffuse over 6.6 nm in microtubules. We discover that while diffusion lengths are influenced by tubulin polymerization state (free tubulin versus tubulin in the microtubule lattice), they are not significantly altered by the average number of protofilaments (13 versus 14). We also demonstrate that the presence of the anesthetics etomidate and isoflurane reduce exciton diffusion. Energy transport as explained by conventional Förster theory (accommodating for interactions between tryptophan and tyrosine residues) does not sufficiently explain our observations. Our studies indicate that microtubules are, unexpectedly, effective light harvesters.

Site identification and docking. After determining the dominant conformations of the dimer from the MD ensemble, we performed a scan for potential binding sites on the dimer using the SiteFinder tool provided in MOE. To focus the analysis not only on the most likely binding sites, but also on pockets with a potentially lower affinity, we investigated all predicted binding sites with a positive PLB score 13 , and docked the four different ligands within these regions. In addition, to further complement the analysis, the ligands were also docked to three experimentally known and clinically relevant binding sites on the tubulin dimer, namely the (a) Colchicine binding site, located near the inter-monomer interface and defined by loop βT7, helix βH8, strands βS8 and βS9 and the αT5 loop; (b) the taxane binding site, located on the β subunit surrounded by helix H1, the H6-H7 loop, the H7 helix, the M loop and the S9-S10 loop; this site is partially overlapping with site number 5 found by the SiteFinder tool (c) the vinca alkaloid binding site, at the longitudinal dimer-dimer interface within protofilaments, surrounded by loop αT7, helix αH10 and strand αS9. Docking runs were carried out in MOE, with the following methodology: we performed the initial placement of the ligands using the triangle matcher algorithm, which describes the active site using α spheres. The generated poses were then scored using the London dG scoring function 3 . Subsequently, for each ligand, the 30 best poses found using the triangle matcher algorithm were further energy-minimized starting from their initial placement, using a Molecular Mechanics approach with the protein kept rigid. The final poses were finally re-scored using the GBVI/WSA dG scoring function 14 , and the top 5 poses were considered for analysis.

S3
(A) Quality of the Homology Models. The human alpha 1A and beta 1 tubulin chains shared a sequence identity of 99.3% and 97.9% with the corresponding chains of the 3J6F template, respectively, thus ensuring the reliability of the homology modelling approach. Indeed, on the phipsi plot, the generated human alpha/beta dimer model featured 94.18% of residues in core regions, 5.35% in allowed regions and 0.47% outliers. In comparison, the same figures for the 3J6F template were 96.81%, 2.83% and 0.36% respectively. The QMeanDisCo global score of the model was 0.76 ± 0.05, which is comparable within error to the corresponding score of the crystallographic template (0.79 ± 0.05) (see Figure S12: comparison of QMean scores). Both quality assessments confirm the comparable structural quality of the homology model with respect to its 3J6F template.
(B) Docking. Clustering of the last 100 ns of the MD ensemble at a 0.15 nm RMSD cutoff yielded a single cluster, suggesting overall conformational stability and no major rearrangements in the protein structure. The centroid of this cluster, representing the dominant conformational state of the tubulin dimer throughout the simulation, was chosen as the conformation to perform docking against. First, we scanned for possible binding sites on the tubulin dimer using the SiteFinder tool in MOE, and found a total of 43 putative binding sites, with PLB scores ranging from 3.66 (best) to -0.79 (worst). All binding sites with positive PLB values were kept as docking sites, resulting in a total of 13 analyzed binding sites. A summary of these binding sites is provided in Table S6 together with the experimentally known binding sites on tubulin for colchicine, taxanes and vinca alkaloids. In the following, the binding site numbering of Table S6 will be used to identify the sites.
Etomidate Predicted binding energies for Etomidate range from -6.77 to -4.36 kcal/mol (mean±std:-5.71±0.51 kcal/mol). The best binding energy was obtained on site 5 (-6.77 kcal/mol), and the interaction between Etomidate and the protein in this site is highlighted in Figure S9A below. The average results in all 16 binding sites are summarized in Table S7.
Isoflurane Predicted binding energies for Isoflurane range from -4.64 to -3.66 kcal/mol (mean±std:-4.11±0.22 kcal/mol), indicative of a slightly weaker predicted interaction with respect to etomidate. The best binding energy was obtained on site 3 (-4.64 kcal/mol), and the corresponding interaction between isoflurane and the protein in this site is highlighted in Figure  S9b, although it is to be underlined how the differences between the different binding poses and binding sites are realistically all within the error of the docking methodology. The average results in all 16 binding sites are summarized in Table S8.
Picrotoxinin Predicted binding energies for Picrotoxinin range from -5.88 to 0.56 kcal/mol (mean±std: 4.22±1.79 kcal/mol). The best binding energy was obtained on the taxol binding site (-5.88/ kcal/mol), and the corresponding interaction between the ligand and the protein in this site is highlighted in Figure S9C. The average results in all 16 binding sites are summarized in Table  S9.
Picrotin Predicted binding energies for Picrotin range from -5.72 to -2.82 kcal/mol (mean±std:-4.89±0.64 kacl/mol). The best binding energy was obtained in the taxol binding site (-5.72 kcal/mol), and the corresponding interaction between the ligand and the protein in this site is highlighted in Figure S9D. The average results in all 16 binding sites are summarized in Table  S10.
Summary of docking studies. Overall, for all compounds except isoflurane, the highest predicted affinity is consistently found to be in or near (site 5) the taxol binding site, which is a known and relevant binding site on the tubulin dimer. Predicted energies in the other putative binding locations are nevertheless similar and most are within the thermal noise level of 0.6 kcal/mol. Hence, no hard conclusions can be drawn regarding a distinctive preferred binding site for either of the analyzed ligands on human tubulin.
Picrotoxin is a 1:1 molar mixture of picrotin and picrotoxinin, two highly similar compounds. Hence, docking was performed for both ligands, and no substantial difference in predicted binding affinities emerged, coherently with their high chemical similarity. No conclusions can be drawn from the docking simulations to distinguish a more active compound on tubulin among the two; in fact, both compounds showed the best predicted dG value in the taxol binding site, where they are predicted to form at least two hydrogen bonds with residues ARG318, ARG359 and ARG276 with differences well within the error of docking predictions.

KINETIC MONTE CARLO SIMULATIONS Preparation of microtubule crystal structure
We constructed the all-atom structure of the full microtubule using the following protocol. Firstly, the tubulin sheet of PDB structure 7SJ7 15 was preprocessed using MOE 3 to fix missing residues in the structure, assign the correct protonation states and perform energy minimization to relief atomic clashes. We then extracted the central tubulin dimer of this structure and used it to construct an initial template of a single tubulin ring, composed of 13 dimers, by fitting 13 copies of the dimer onto the original electron density map deposited in the Electron Microscopy Data Bank (accession code EMD-25156), using the "Fit in Map" tool of UCSF ChimeraX 16 . To construct the final microtubule, we first carried out 200-ns Molecular Dynamics simulations of the tubulin sheet to equilibrate the structure and obtain a diversified structural ensemble. In detail, simulations were carried out in GROMACS 2021.4 9 in the NPT ensemble, using the AMBER 99SB-ILDN force field 6 , PME electrostatics with a cutoff of 1.2 nm, the velocity rescale thermostat 10 at 300K and the Parrinello-Rahman barostat 12 . We subsequently clustered the last 150 ns of this MD simulation with a cutoff of 0.10 nm using the gmx cluster tool, and extracted the centroids of each cluster, focusing only on the central dimer of the tubulin sheet. These centroids represent a set of dominant conformations of the tubulin dimer within a microtubule, and were used to build the final microtubule. In greater detail, we generated each of the final tubulin rings by randomly choosing centroids from the MD clusters and RMSD-fitting them onto the template tubulin ring described above. This was repeated for each of the 31 tubulin rings composing the final microtubule, which was assembled by spacing the rings axially by 8.15 nm, as described in earlier literature for the 13:3B MT lattice 17 . In summary, we eventually obtained a final B-lattice microtubule with a length of 252.65 nm and composed of an ensemble of different tubulin dimer conformations, extracted from all-atom MD simulations and representative of the conformational heterogeneity of tubulin within a microtubule.
Simulation protocol for homo-and hetero-transfer among tryptophan and tyrosine residues. The monomers of TYR and TRP were optimized in the CAM-B3LYP/AUG-CC-PVDZ level of theory and the optimized geometry was used to calculate the transition dipole moments using TDDFT (CAM-B3LYP/AUG-CC-PVDZ). The calculations were carried out in Gaussian 16 program package 18 .The obtained transition dipole moments were reoriented at the respective molecular coordinates for calculating the long-range coulombic coupling strengths in the 31-long dimer microtubule crystal structure. We took the dissipation rate (kDis) as 3.06 x 10 8 s -1 19 .
Charge transfer Integral calculations were carried out on 5 representative TYR dimers selected from the microtubule crystal structure with Inter-aromatic (RIA) distances below 6.3 Å. A constrained optimization (B3LYP/6-311+G(d)) were performed on the selected TYR dimers by freezing the carbon skeleton prior to charge transfer integral calculations. The charge transfers Integral calculations (B3LYP/6-311+G(d)) were carried out using the dimer projection (DIPRO) method 20 using the CATNIP program package 21 available from https://github.com/JoshuaSBrown/QC_Tools).        Table S1. Nearest neighbor distances (Å) between tryptophan residues in different tubulin polymers. CD2 (carbon at the delta-2 position) atoms in pairs of tryptophan residues were used to measure distances between residues.