Computational Predictions of Volatile Anesthetic Interactions with the Microtubule Cytoskeleton: Implications for Side Effects of General Anesthesia

The cytoskeleton is essential to cell morphology, cargo trafficking, and cell division. As the neuronal cytoskeleton is extremely complex, it is no wonder that a startling number of neurodegenerative disorders (including but not limited to Alzheimer’s disease, Parkinson’s disease and Huntington’s disease) share the common feature of a dysfunctional neuronal cytoskeleton. Recently, concern has been raised about a possible link between anesthesia, post-operative cognitive dysfunction, and the exacerbation of neurodegenerative disorders. Experimental investigations suggest that anesthetics bind to and affect cytoskeletal microtubules, and that anesthesia-related cognitive dysfunction involves microtubule instability, hyper-phosphorylation of the microtubule-associated protein tau, and tau separation from microtubules. However, exact mechanisms are yet to be identified. In this paper the interaction of anesthetics with the microtubule subunit protein tubulin is investigated using computer-modeling methods. Homology modeling, molecular dynamics simulations and surface geometry techniques were used to determine putative binding sites for volatile anesthetics on tubulin. This was followed by free energy based docking calculations for halothane (2-bromo-2-chloro-1,1,1-trifluoroethane) on the tubulin body, and C-terminal regions for specific tubulin isotypes. Locations of the putative binding sites, halothane binding energies and the relation to cytoskeleton function are reported in this paper.

Tubulin proteins polymerize to form microtubules (MTs), nanoscale cylindrically shaped protein polymers that are part of the cellular cytoskeleton. The neuronal MT cytoskeleton, in particular, possesses a unique architecture [23], responsible for maintaining highly asymmetric neuron morphology and the intracellular transport of vesicles. Unlike MTs in all other cells, the MTs in dendrites are interrupted and oriented in local networks of mixed polarity. Moreover, the ionotropic GABA A receptor is anchored to the MT cytoskeleton via associated proteins [24], and an intact MT structure has been shown to be essential for the activity of GABA A receptors [25]. Likewise, MTs have been shown to modulate both sodium and calcium current within neurons [26][27][28][29][30].
A large number of neurological brain disorders bear the common feature of some kind of disruption to the neuronal cytoskeleton, particularly MTs, either directly or indirectly through associated proteins. These include neurodegenerative disorders such as Alzheimer's disease, Parkinson's disease, and Huntington's disease [31]. The relation between anesthesia, postoperative cognitive decline and dementia (in Alzheimer's disease and the other neurodegenerative disorders) is uncertain [32][33][34][35]. Anesthesia may exacerbate neurodegeneration, and/or have its own deleterious effect. Cytoskeletal MTs are a common link, affected by anesthetics and disrupted in neurodegeneration.
This suggests an important role for the MT cytoskeleton in anesthetic action, with potential side effects related to post-operative cognitive dysfunction. Recently it has been shown that genetic expression of tubulin is altered following exposure to desflurane [36], sevoflurane [37], halothane, and isoflurane [38]. LeFreche et al. [39,40] further showed that cognitive dysfunction following sevoflurane anesthesia was associated with hyperphosphorylation of the MT-associated protein (MAP) tau, separation of tau from MTs, and MT instability, the same signs associated with neurofibrillary tangles (NFTs) in Alzheimer's disease. Anesthetic binding to tubulin in MTs may be involved in anesthetic action and toxicity.
Tubulin, a peanut-shaped heterodimer with a and b monomers, has been identified as a direct binding target for halothane [22]. Additionally, experiment shows volatile anesthetics, particularly halothane, alter the tubulin self-assembly rates into MTs in a number of systems, both in vivo [41][42][43][44][45][46] and in vitro [47], albeit at extremely high concentrations. When polymerized into MT form, each tubulin subunit interacts with surrounding dimers, forming longitudinal contacts between dimers along the protofilament length, and lateral contacts between protofilaments (see Figure 1). There also exists an intradimer interaction between the aand bmonomers of a single tubulin dimer. In the more prevalent Blattice MT formation the lateral contacts are formed between like subunits (i.e. a-monomer to a-monomer, and b-monomer to bmonomer on adjacent protofilaments). In the less common Alattice, lateral interactions are between a and b monomers. Volatile anesthetics may inhibit MT assembly dynamics by a direct molecular interaction between the anesthetic molecule and the tubulin dimer hindering dimer-dimer, or intradimer interactions. It is also feasible that volatile anesthetics exert their action on MT dynamics through an alteration of the local environment affecting the highly flexible C-terminal tail regions of tubulin. Clearly, the site and mechanism of anesthetic action on MT assembly and stability remain to be determined.
Experimental methods to study anesthetic binding to proteins include NMR spectroscopy, photoaffinity labeling, and site directed mutagensis. NMR and photoaffinity labeling techniques can only be applied to purified proteins available in relatively large quantity. Currently, NMR methods are only capable of determining the structure of protein complexes with masses up to 20-30 kDa, which is well below the 110 kDa size of the tubulin heterodimer, rendering this method not directly applicable to this problem [48]. Without a priori knowledge of putative binding sites, site directed mutagenesis is also hindered by protein size. Direct photoaffinity labeling with halothane requires combination with methods such as protein digestion followed by mass spectroscopy to determine binding location, which also benefits from previously predicted locations to determine experimental protocols.
Surface geometry techniques used to predict anesthetic binding sites on proteins are based on static structural data [49,50]. These methods suffer a general weakness by not accounting for protein dynamics, neglecting rearrangements of local protein atoms and the resulting change in binding site availability. To address these issues we use a combination of molecular dynamics (MD) and surface geometry based binding site prediction to identify general putative volatile anesthetic binding sites on, or in, the tubulin protein. Blind docking followed site prediction to obtain halothane (2-bromo-2-chloro-1,1,1-trifluoroethane) binding energy estimates, as the majority of experiments between volatile anesthetics and tubulin investigate the interaction with halothane.

Putative Volatile Anesthetic Binding Sites on the Tubulin Body
A short 5 ns MD simulation was performed on two tubulin dimers in MT geometry with periodic boundary conditions, effectively modeling two infinite protofilaments. This 5 ns simulation is too brief to be representative of the entire conformational ensemble. However, it does serve to allow side chain motion while keeping the protein backbone of the structure relatively stable (see Figure 2).
Clustering of the 5 ns MD simulation trajectories resulted in 11 distinct protein conformations, with each cluster containing several snapshots of the protein at different time steps. Taking each snapshot in a cluster to be represented by the average conformation of the cluster resulted in the 11 dominant conformations existing for various portions of the simulation (see Table 1). Since each dominant conformation is represented by the average conformation of several snapshots at different timesteps, a given dominant conformation represents a certain percentage of the MD simulation (see Table 1 -% Simulation).
A modified PASS algorithm [51] was performed on each of the dominant conformations for both dimers in MT protofilament geometry, and for each dimer separately. PASS predicts putative binding sites (hydrophobic crevices and pockets) through an iterative coating of the protein surface with probe spheres. Potential sites are based on the burial depth of these spheres. Modification of this algorithm to additionally measure for hydrophobicity yields efficient prediction of volatile anesthetic binding sites [49]. This procedure yielded numerous putative anesthetic-binding sites on tubulin, which would be valid for any volatile anesthetic.
Due to the motion of side chains the predicted sites varied between the different protein conformations. The DBSCAN method [52] spatially grouped the predicted sites yielding 47 unique potential binding sites on the tubulin protein, however some sites were not found in all of the conformations. As such, each site was assigned a persistence value denoting the percentage of the MD simulation in which the potential binding site was found. The persistence of each site was calculated by taking the sum of the simulation percentages (see Table 1 -% Simulation) for each of the dominant conformations on which the site was predicted (e.g. If a site was found for dominant tubulin conformation 1, 2, and 3, the persistence would be 5.18% + 0.80% + 11.55% = 17.53% of the MD simulation). Persistence values varied greatly from 0.80% to 100% of the 5 ns simulation (see Table 2 and 3). Of the 47 predicted sites, 9 persisted for more than 70% of the simulation, and of these 5 persisted for the entire simulation (see Figure 3).

Halothane Docking
Tubulin body. Focused docking of a halothane molecule (see Figure 4) at each of the predicted putative volatile anesthetic binding sites on each of the tubulin clusters yielded varying binding energies between clusters for a given site. The binding mode with the lowest binding energy value for a given site was taken, generally resulting in one binding mode whose binding energy is listed in Table 2. In the case of multiple binding modes, the binding energy of the largest cluster was taken. Examination of the energy contributions yielded binding due mainly to van der Waals interactions with the Cl and Br atoms generally contributing the largest portion. Other energy contributions served to weaken this interaction.
C-terminal tails. For each tubulin isotype found in the brain, 50 distinct conformations of the C-terminal tail regions were generated sampling the tail conformational space (see    Table 4. The energy contributions yielded binding due mainly to van der Waals interactions again with the Cl and Br atoms contributing the largest portion. In general, binding energies increased with the number of available surrounding residues. Thus, tail conformations, which were compacted, forming loops or coils, provided more favorable binding conditions. Binding energies for these ideal binding-conditions were comparable to binding on the tubulin body.

Microtubule Polymerization Assay
Multiple polymerization assays were run for tubulin alone (control), tubulin with 10 mM paclitaxel, tubulin with 40 mM halothane, and tubulin with both 10 mM paclitaxel and 40 mM halothane (see Figure 6). A 40 mM halothane concentration was used to ensure interaction of halothane with tubulin. The control run showed optical density results consistent with normal tubulin polymerization. In the presence of 10 mM paclitaxel optical density increases sharply and falling off in the longer time limit. This is consistent with normal results obtained from mist experiment with quick polymerization showing the stabilizing effect of paclitaxel. In the presence of both 10 mM paclitaxel and 40 mM halothane the optical density curve is comparable to that of paclitaxel alone within standard error. This indicates that halothane has no effect on the interaction between paclitaxel and tubulin. In the presence of 40 mM halothane alone the optical density curve resembles the control situation, however the standard error indicates difference between the curves with a small increase in the optical density.

Interaction Mechanisms
The average brain volume is 1.4 L. Assuming approximately 10 9 tubulin dimers per neuron, with 10 11 neurons per brain, the average concentration of tubulin protein in the brain is 120 mM. However, in the neuron this may be 2-3 times higher since glial cells, possessing a MT density much less than neurons, are considered to comprise 50% of the brain. The remaining portion is composed of both neurons and other necessary structures such  as ventricles, blood vessels etc. Comparing an assumed total neuronal tubulin protein concentration of 240-360 mM, to a halothane concentration of 250-500 mM (,1 to 2 MAC) gives 1 to 2 halothane molecules per tubulin dimer at clinically relevant concentrations. However, occupancy depends not on equilibrated concentrations of protein and drug alone, but rather drug concentration and dissociation constant. The halothane binding energies suggest K d values between 6 and 16 mM. This yields a fractional occupancy of 1%-8% at 1 to 2 MAC, assuming the system is equilibrated. However, this does not account for the nonequilibrium state of biological systems, nor does it account for multiple affinity sites or sites of partial binding. The existence of many sites with similar binding energies made it difficult to assign binding to any particular site(s). In fact, it is likely that anesthetics bind non-specifically to many of the predicted binding sites. Low persistence of a binding site does not necessarily indicate that a potential site is invalid. Rather, it implies a lack of favorable conditions for binding, since these sites are associated with greater overall conformational free energies of the protein system. However, anesthetic molecules may bind to low persistence sites, potentially with a greater binding energy than to higher persistence sites. In light of this, it is expected that at a constant anesthetic concentration, the sites that are most occupied are determined by the sum of the conformational energy differences, as reflected in persistence, and binding free energy differences.
A total of 32 binding sites (16 on a-tubulin, and 16 on b-tubulin) were predicted on a single tubulin dimer, which were independent  of the dimer placement in MT geometry (see Table 2). Only one of the sites predicted on a-tubulin (site 23) persisted for more than 70% of the simulation, in fact lasting for the entirety of the simulation. Four sites predicted on b-tubulin (sites 1, 4, 7 and 21) lasted for more than 70% of the simulation, with two sites (4 and 21) lasting the entire simulation. An additional 8 sites were located at the intradimer a-b interface (see Table 3). These also did not depend on the placement of the dimer within the MT structure. Site 38 persisted in this region for 84.06% of the simulation, while sites 5 and 37 lasted for 100% of the time. Seven sites were dependent on the placement of the dimers within the MT lattice (see Table 3). One site was predicted at the dimer-dimer interface lasting for 88.05% of the simulation. The remaining 6 sites were located at the protofilament interfaces, either between aor btubulins, however all of these persisted for less than 50% of the simulation. Reasonable anesthetic concentrations may thus alter only longitudinal, or intradimer interactions, but only at sufficiently high concentrations would lateral attraction be affected.
Binding energy estimates for halothane with the tubulin Cterminal tails are comparable to binding on the main protein body suggesting another mode of interaction. Larger binding energies exist for more compact conformations of the C-terminal tail. As such, due to the flexibility of the C-terminal tails, interaction with halothane may sequester the tail region, holding them in more compact forms, and preventing normal tail movements. This is of importance to the function of MTs as evidence indicates the Ctermini play critical roles in regulating microtubule structure, function and interaction with MAPs [53][54][55][56][57][58][59][60]. Sequestration of the C-terminal tails by halothane into compact forms may indeed alter tubulin polymerization dynamics.
Of the key residues, anesthetic site 10 is within 5 Å of bL219, and site 42 is within 5 Å of bL217 and bP274. bP274 forms part of the hydrophobic pocket for the taxol 39-phenyl group, while bL217 makes hydrophobic contact with the 2-phenyl ring of taxol, a component of taxol shown to be absolutely required for its activity [62]. bL219 also makes hydrophobic contact with the 2phenyl ring.
Interaction of volatile anesthetics in this region may alter taxol binding resulting in loss of function. However, site 42 is negligibly persistent, only lasting for 0.80 % of the simulation, and site 10 persists minimally for ,27% of the simulation. Therefore, this behavior would only be expected in cases of high anesthetic concentrations.
The tubulin polymerization experiment suggests that this is indeed the case. Tubulin in the presence of 10 mM of paclitaxel produced similar optical density curves over time in both the presence and absence of 40 mM concentration of halothane. Thus, halothane does not affect the interaction of taxol with tubulin consistent with the computational prediction.
Colchicine. Colchicine inhibits MT polymerization by binding to free tubulin dimers at the intradimer interface and interacting with both aand b-tubulin residues. The colchicine site is mainly located within the b-tubulin, and is surrounded by residues bC241, the b-tubulin T7 loop, H8 helix (including bD251), S8 strand (including bV318), and S9 strand (including bK352), and interacts with residues aE71, aN101, and the atubulin T5 loop (including aV181) [63].
Steric hindrance between aV181 and colchicine prevents the atubulin subunit from occupying its normal position interfering with the straight conformation adopted by tubulin in protofilaments. The loss of this straight conformation prevents the tubulin M-loop from establishing lateral contacts between protofilament ends, leading to hindrance of MT dynamics and eventually MT depolymerization (or lack of MT polymerization from free tubulin). Additionally, binding of colchicine requires movement of the b-tubulin T7 loop and H8 helix to accommodate the drug molecule, and this movement also interferes with the a-b interaction.
Site 38 is within 5 Å of residue aV181, of the a-tubulin T5 loop, and the b-tubulin H8 helix and S8 strand, and persists for a considerable 84.06% of the simulation. Site 7 is located within 5 Å of residue bC241, and the b-tubulin T7 loop and S9 strand and persists for 78.88 %, while site 41 persists for 56.18% and is located within 5 Å of the b-tubulin H8 loop. Site 34 is located within 5 Å of the b-tubulin S8 and S9 strands including residue bK352, and site 46 is within 5 Å distance of the a-tubulin T5 loop, Table 4. Tubulin isotype, sequence and halothane binding energy range for the C-terminal tail regions found in the brain. including residue aV181, and bK352 of the b-tubulin S9 strand, however they persist minimally for less than 20% of the simulation.
The preferred binding mode of halothane at site 38 was shown to occupy a space surrounded by helix bH8 and strand bS8, while the binding mode at site 7 is within 3 Å of bC241 and surrounded by helix bH8 and strand bS9 (see Figure 7). Van der Waals interaction of halothane with the bH8 helix would reduce the mobility of these residues preventing the accommodation of the colchicine molecule in this pocket. This is consistent with findings that report a reduction in colchicine binding to tubulin in the presence of halothane [64].
Interaction of volatile anesthetics with the b-tubulin T7 loop or H8 helix can reduce the mobility of these regions preventing the movement required to accommodate the colchicine molecule within tubulin. Additionally, interaction of volatile anesthetics with residue aV181 can hinder the movement of the a-tubulin T5 loop produced by the steric clash with colchicine. As these sites are persistent it is thus expected that volatile anesthetics would hinder the binding of colchicine. The effect of volatile anesthetics on the binding of colchicine derivatives must also be considered as it is known that colchicine derivatives substituted at the methoxy positions of ring A can be cross-linked with bC241 and anesthetics at site 7 may reduce this interaction.
Vinblastine. Vinblastine is one of the vinca alkaloids and is an inhibitor of tubulin polymerization. The primary binding site of vinblastine is located at the MT plus (growing) end towards the inner lumen of the MT and works by inhibiting longitudinal (dimer-dimer) contacts in a protofilament. The vinblastine binding-pocket is boxed by the a-tubulin T7 loop (including aL248 and aN249), H10 helix and S9 strand (including aL352), and the b-tubulin T5 loop (including bV177 and bD179), H6 helix carboxy-terminal turn (including bY210), and the loop region between helices H6 and H7 (including bF214) [65]. Vinblastine affecting these regions alters the longitudinal interface between tubulin dimers constraining them in a curved assembly to avoid steric clashes with the drug molecule. The curvature induces displacement of the M-loop on subsequent dimers weakening lateral interactions leading to reduced MT dynamics and eventual MT depolymerization.
Site 21, which persists for the entirety of the simulation, is within 5 Å of bY210 in the carboxy-terminal turn of helix H6. Site 39, which persists for 88.05% of the simulation, and site 1, which persists for 70.12% of the simulation, neighbor residue bD179. Site 1 is also adjacent to bV177. Interaction of volatile anesthetics with these residues is expected to alter vinblastine binding either through competitive binding, or hindrance of the mobility required to avoid steric clashes. Sites 19 and 20 are in the region of the a-tubulin H10 helix, and sites 17 and 42 are within the btubulin H6 helix and H6-H7 loop regions. However, these sites all exhibit low persistence.
The effect of volatile anesthetics on derivatives of vinorelbine, a vinblastine analogue modified on the D' ring, should also be considered. Changes in these derivatives, which yield modified affinities and activities, occur in regions that interact with bY224 [66], a residue within distance of site 21.
Halothane binding in the vinblastine pocket is shown in Figure 8. The binding mode of halothane at site 21 is within 2 Å of bY210 and the molecule is surrounded by helices bH6 and bH7. Halothane binding at site 39 is within 5 Å of aK352 and bD179 and is cradled by the surrounding aT7 loop. Halothane at site 1 lies within 3 Å of bD177 and is surrounded on each side by loops bT5 and bT4. Due to the van der Waals interaction halothane can alter the mobility of these regions key to vinblastine binding potentially reducing vinblastine binding.
One point of interest is that vinblastine and colchicine are known to inhibit assembly at intermediate concentrations while inducing MT-to-macrotubule transformations at high concentrations [67][68][69]. Indeed it has been shown that halothane and other volatile anesthetics may induce the same aberrant forms of tubulin polymerization, albeit at extremely high concentrations [43,44,70]. Preferential binding of halothane to a less stable conformation of the vinblastine or colchicine cavity may result in destabilization of the longitudinal interaction resulting in an effect similar to that of these antimitotic drugs.
Peloruside A. The macrolide peloruside A is a MT-stabilizing agent that synergizes with taxoid drugs by acting at a unique site. The proposed peloruside A binding site is located on the exterior of the b-tubulin and is defined by the H9 helix (including bQ294), the H9-H9' loop (including bD297 and bR308), the H9' helix, the H9'-S9 loop, portions of strand S8, and the H10 helix (including bV335, bN339 and bY342) [71]. It is suggested that peloruside A stabilizes protofilament interactions by securing the a-tubulin T5 loop with the adjacent residues found in the btubulin H9-H9' loop, the H9' helix, and the H9'-S8 loop, although the exact mechanism is not known.
Site 4 contacts the H9 helix, the H9-H9' loop, the H9' helix and the H9'-S8 loop and persists for 100% of the simulation. Additionally site 37 and site 38 contact the H9' helix and H9'-S8 loop, and persist for 100% and 84.06% of the simulation, respectively. While the exact mechanism by which peloruside A acts is unknown, it can be expected that volatile anesthetics in this binding-region would alter both drug binding, and the resulting conformational changes which yield its stabilizing effect at reasonable anesthetic concentrations.
Sites 9, 17, 30, 31, 34 and 42 are all within this region as well, however their persistence is low. At high anesthetic concentration it is expected that the multitude of potential sites in this region would alter the effects of peloruside A significantly.
Laulimalide. Laulimalide is a potent, macrolide MT-stabilizing agent that binds to the exterior of the MT on b-tubulin near the C-terminal E-hook. The binding site, which encompasses the proposed peloruside A pocket, is surrounded by the b-tubulin H9-H9' loop (including bF296), H9' helix (including bR308), H10 helix (including bV335) and the H10-S8 loop (including bN339 and bY342) [72]. Laulimalide does not appear to bind to the atubulin.
Upon laulimalide binding both bR308 and bY342 reorganize to generate a cation-p interaction, increasing the stabilization of the loops in this pocket. In this configuration bR306 establishes polar contacts with the oxygen of the dihydropyran, in the laulimalide side chain. The hydrophobic bV335 and bF296 residues also reorganize to align in parallel with the laulimalide macrolactone ring, with the movement of bF296 creating an entrance into a deep cavity in b-tubulin. Finally, the bN337 residue, which also defines this cavity, forms a hydrogen bond with the hydroxyl group at laulimalide C15.
As with peloruside A, Site 4, 37 and 38 persist in this region for a considerable portion of the simulation suggesting that volatile anesthetics binding in this region would alter both drug binding, and the resulting conformational changes. Sites 9, 17, 30, 31, 34 and 42, again, would be expected to play a role only at high anesthetic concentration. One key point of interest is that site 4 is within 5 Å of bF296. Reduced mobility of this residue, by interaction with an anesthetic molecule, is expected to inhibit binding of the laulimalide molecule through unfavorable alignment of bF296 with the laulimalide macrolactone ring and subsequent shrinkage of the deep cavity entrance. Additionally, even though laulimalide may possess a greater binding affinity for tubulin than the volatile anesthetics, small anesthetic molecules occupying the deep cavity region may be prevented from escaping by an impinging laulimalide molecule. This would result in a binding ''stalemate'' preventing proper laulimalide effects.
Site 5, 23, 38 and 7 all persist for a significant portion of the simulation (100%, 100%, 74.06% and 78.88%, respectively), and all neighbor residues that directly interact with either GTP or the magnesium ion (Mg 2+ ) required for the intradimer ab stability. Sites 2, 24, 34, 35, and 46 are also in this same region, but persist for less than 21% of the simulation. Anesthetics at these sites may interfere with interactions between GTP, Mg 2+ , and tubulin. The effect of volatile anesthetics on GTP binding has not been studied for tubulin, but it has been found that volatile agents, at clinically relevant doses, have a direct effect on the conformation and stability of the GTP/ Mg 2+ bound state of some, but not all, GTP binding proteins [73]. This indicates the possibility of volatile anesthetics altering GTP dependent MT dynamics resulting in reduced polymerization, however further investigation is needed.
Also of note, site 5 borders residues aI171 and aT179, and site 38 is adjacent to residue aA180, which are all part of the sugar binding T5 loop. These residues are involved in longitudinal contacts between dimers, and the lateral, protofilament, interac-tions. This suggests a possible anesthetic induced MT depolymerization mechanism.
Exchangeable GTP. GTP bound to the b-tubulin, at the Esite, is exchangeable and can be hydrolyzed to GDP. GTP is required at the E-site in order for tubulin to polymerize, but is hydrolyzed to GDP upon polymerization. The result is a metastable MT structure stabilized by a so-called cap of remaining GTP-tubulin subunits. The loss of this cap results in rapid depolymerization.
Sites 1, 21 and 39 all persist for a significant portion of the simulation (70.12%, 100%, and 88.05%, respectively), and all neighbor residues that directly interact with either GTP or the magnesium ion (Mg 2+ ) in the non-hydrolyzed state, or with GDP in the hydrolyzed state. Sites 10 and 17 also make contact with key residues, but persist for less than 28% of the simulation. This indicates a potential for volatile anesthetics to affect GTP dependent polymerization dynamics as discussed above.
Similar to the N-site, sites 1, 21 and 39 neighbor residues in the T5 loop, potentially affecting lateral interactions between protofilaments. Additionally, site 21 encompasses bN206. Tubulin's preference for GTP is the result of hydrogen bonding of the 2exocyclic amino group in GTP to the hydroxyl groups of the N206 and N228 residues, and by hydrogen bonding of the 6-oxo group to the amino group of N206 [61]. Interactions of anesthetics with bN206 at the E-site could potentially reduce tubulin's specificity for GTP hampering the stabilizing nature of the GTP cap resulting in increased depolymerization events.

Comparison to MAP-binding Sites
MAPs bind to polymerized tubulin to regulate MT dynamics. The numerous MAPs identified carry out a wide range of functions dependent on the host tissue in which they are found. These functions include stabilizing and destabilizing MTs, guiding MT transport, cross-linking of MTs, and mediating interactions between MTs and other proteins or membranes. In the neuron MAP1A, MAP1B, MAP2 and tau are the most prominent MAPs, and all serve to stabilize MTs.
The exact binding location of these MAPs on tubulin is not known. However, it is known that MAP2 and tau share a sequence similarity in their MT binding repeats close to their carboxyl ends [74,75]. The proline-glycine-glycine-glycine (PGGG) motifs, believed to form tight turns upon association with MTs, are analogous to a sequence in the S9-S10 loop of a-tubulin (residues T361-L368), which fills the equivalent of the taxol site in the atubulin stabilizing the M-loop conformation [76,77]. It is believed that MAPs occupy the taxol-binding region in the b-tubulin securing the M-loop and stabilizing MTs in a manner similar to the action of taxol.
As discussed above, interactions of volatile anesthetics with MTs are not expected to alter the effect of taxol except at high anesthetic concentration. This can be extended to include MAPs acting by the same mechanism. The effect of volatile anesthetics on other MAP binding modes is unknown.

Comparison to Post-translational Modification Sites
In animals, tubulin may undergo numerous post-translational modifications involved in regulating MT stability and the interaction with MAPs. The a-tubulin may undergo acetylation, tyrosinolation/detyrosinolation, deglutamylation, palimotoylation, polyglutamylation polyglycylation, and phosphorylation whereas the b-tubulin may only undergo the latter three [78]. The majority of these modifications, including tyrosinolation/detyrosinolation, deglutamylation, polyglutamylation, and polyglycylation, occur in the C-terminal region of tubulin.
Post-translational modifications require enzymes to interact with the C-terminal region, and it is expected that this process would require the tails to exhibit normal flexibility. The interaction of halothane with the C-terminal tails indicates an increase in binding energy with the number of residues interacting with the molecule. Halothane may thus serve to reduce the motion of these tail regions interfering with normal post-translation modification processes.
Acetylation of a-tubulin takes place at the e-amino group of aK40. Palimotoylation occurs on a-tubulin C376. No predicted volatile anesthetic sites reside near these locations.
Phosphorylation takes place on both aand b-tubulin at tyrosine residues, or serine/ threonine residues on after residue 307 [79]. Several predicted sites reside near serine or threonine residues in this region, or near tyrosine residues, but the majority of these sites persist for less than 50% of the simulation. However, site 38, persisting for 84.06% of the time, is within 5 Å of bT314. Additionally, site 21, which lasts for the entire simulation, is within 5 Å of both bY210 and bY224. This indicates a potential interference of volatile anesthetics with tubulin phosphorylation, yet the exact location of this modification is not known. This is of interest as tubulin phosphorylation has been suggested as a molecular mechanism for memory encoding [80,81].

Comparison to Zinc Binding Sites
In regards to Alzheimer's Disease, zinc has been shown to promote the aggregation of b-amyloid, which then sequesters the metal making it largely unavailable, causing zinc dyshomeostasis in the vicinity of b-amyloid deposits [82]. While b-amyloid plaques are necessary to initiate the neurodegenerative process in AD, it is the NFTs that lead to neurodegeneration. Recently it has been argued that zinc sequestration by b-amyloid deposits deplete intraneuronal zinc which drives formation of NFT, MT destabilization and associated neuronal degeneration [83]. The study, which investigates zinc binding to MTs and their component tubulin proteins, predicts 6 zinc-binding sites under the highest stringency.
Cysteine, histidine, aspartic acid and glutamic acid amino acid residues are key to zinc binding and account for ,97% of all zincbinding amino acids [84]. The 6 predicted binding sites and key residues involved include the site associated with zinc-induced tubulin sheet formation (including aH192, aE420 and aD424) [61], a site near the non-exchangeable GTP site, an uncharacterized site on the a-tubulin (including aH266 and aD431), a site near the colchicine binding site (including bC243, bC356 and bD357), and two uncharacterized sites on the b-tubulin.
Site 33 neighbors the aH266 and aD431 residues of the uncharacterized a-tubulin site, but persists for a mere 1.59% of the simulation, indicating that volatile anesthetic interactions with zinc bound at this site are not expected to play a significant role unless high anesthetic concentrations are used. Site 8, which persists for 51% of the simulation, is within 5 Å of aD424 of the zinc-induced tubulin sheet formation site. This site is only expected to play a role in tubulin polymerization at zinc concentrations well above physiological conditions, however the zinc-binding site near the colchicine site is expected to play a significant role in tubulin polymerization at physiological zinc concentrations [83]. Site 7, and 44 border on the colchicine zinc-binding site. Site 7, which persists for 78.88%, is near bC241, while site 44, which persists for 22.71%, is adjacent to bC356 and bD357. Due to their large dipole moments, the presence of a volatile anesthetic in these regions is expected to interfere with van der Waals forces required for the residue coordination necessary for zinc binding. Should this occur volatile anesthetics would temporarily inhibit zinc binding, reducing the tubulin polymerization rate, resulting in a net loss of polymerized MTs, similar to the effect caused by b-amyloid zinc sequestration [83]. This could potentially cause memory impairment effects.

Comparison to Aromatic Amino Acids
While the mechanism of general anesthetic action responsible for erasing conscious awareness is still under debate, quantum mobility theory posits that anesthetics act by inhibiting signal propagation between the aromatic amino acids of tubulin [85]. Aromatic amino acids include tyrosine, tryptophan, phenylalanine and histidine, of which tubulin possesses 35, 8, 43 and 22, respectively. Numerous sites are near tyrosine residues (sites 11, 12, 16, 21, 22, 24, 32 and 33), phenylalanine residues (sites 4, 15, 16, 24, 25, 28, 32, 41, and 45), and histidine residues (sites 14, 22, 25, 26, 33 and 36) with some lasting for the entire simulation. Of these, tryptophan is posited to play a prominent role in quantum mobility theory. The closest predicted site to a tryptophan residue is site 37, persisting for 100% of the simulation, which is within 5.3 Å of bW346. This is consistent with the predictions of quantum mobility theory.

Experimental Validation
Numerous historical studies give evidence to the interaction between halothane and tubulin polymerization. The earliest studies showed concentrations of 2-2.8 % halothane in air by volume (,0.1 mM) yielded retraction of MT based axopods suggesting MT depolymerization [41,42]. Consistent results were found with larger concentrations of halothane. At 10 mM halothane, MTs were shown to decrease in length and density [44], and MT based flagellar structures were shortened [46], and 20 mM halothane likewise decreased MT density [45,86]. However, other investigations showed concentrations between 3 and 10 mM increased the number of MTs per area [44,86,87]. Other evidence indicates that halothane concentrations of 10 mM can produce aberrant MT formations, including ribbon structures and macrotubules [43,44,88].
This range of behaviors can be attributed to differing polymerization protocols. Tubulin polymerization is affected by many factors including pH, temperature and ionic concentrations. Similar variations in behavior were observed when first characterizing the effect of zinc on tubulin polymerization [89]. By varying buffer conditions and the ratio of tubulin to zinc, polymerization ranged from normal to aberrant, with various intermediate stages.
It is important to note that optical density measurements were not capable of quantifying polymerization of these zinc induced aberrant tubulin structures [89]. Optical density measurements also cannot detect the aberrant assembly of tubulin induced by the MT destabilizers cryptophycin 1 [90,91,92] and hemiasterlin [90]. The kinetics of MT assembly are often studied by optical density measurements since the geometry of long thin rods is such that absorption data is a function of the total mass of the assembled protein. As such, without knowledge of the shape of tubulin aggregates it would be dangerous to interpret absorption data quantitatively.
In the present study, a halothane concentration of 40 mM was used to ensure saturation of tubulin by halothane. The optical density measurements show that in the presence of both 10 mM paclitaxel and 40 mM halothane the optical density curve is comparable to that of paclitaxel alone within standard error. This suggests that halothane has no effect on the interaction between paclitaxel and tubulin, a result consistent with the computational predictions presented. In the presence of 40 mM halothane alone the optical density curve resembles the control, however the standard error indicates difference between the curves with a small increase in the optical density. While this seems to indicate that halothane does little to modulate tubulin assembly, previous evidence of altered tubulin structures and the inability of optical density to quantify such changes suggest that this may not be the case. Further investigation is clearly required.

Summary
Anesthetic binding to tubulin may be important for the mechanism of anesthetic action, and also for anesthetic side effects related to post-operative cognitive dysfunction and/or exacerbation of neurodegenerative diseases. The size of the tubulin macromolecule, its numerous preexisting non-polar, hydrophobic cavities, and the generally weak binding of volatile anesthetics hinder the experimental investigation of this molecular mechanism of interaction. For this reason, it is an important problem to be able to predict computationally the binding of anesthetics to tubulin, and of small molecules to proteins in general, although this is a challenge since binding is often nonspecific, with multiple binding sites being filled. Thus, we have used a combination of computational methods including molecular dynamics simulations, surface geometry based anesthetic binding site prediction, focused and blind docking to identify putative volatile anesthetic binding sites on, and in, the tubulin protein. We hope that this work will enable future experiments to resolve the necessary ambiguity in our computational results.
Multiple binding sites were found on the tubulin protein, but the availability of these sites for anesthetic binding was found to vary greatly. Since most binding energies were quite close in value, it is expected that those with lowest persistence will be filled only at large anesthetic concentrations.
Since the binding energy of volatile anesthetics in protein hydrophobic pockets is generally small, the potential for anesthetic molecules to bind by inducing fits through changes in protein conformation is extremely low. It is more likely that these molecules bind in preexisting non-polar, hydrophobic cavities. In this case destabilization of a protein, or protein structure, may result either from preferential binding of the anesthetic to a less stable conformation of the cavity, or the disruption of allosteric changes at protein interfaces.
We found that favorable thermodynamic conditions for the binding of anesthetics to tubulin result from van der Waals interactions. The nine sites predicted to persist for greater than 70% of the 5 ns simulation were located in the binding-pockets for colchicine, vinblastine, peloruside A, laulimalide, GTP, and GDP. These sites all reside in regions of either intradimer, or longitudinal interaction, indicating that at reasonable concentrations anesthetics do not alter lateral interactions. In fact, we do not predict strong binding in the taxol-binding site, which has been implicated in lateral interactions. This was also confirmed by experimental validation within the present investigation.
Our findings suggest that modification of intradimer and longitudinal interactions may be the general mode of MT destabilization by volatile anesthetics. This is consistent with observations that anesthetics are weak destabilizers of MTs under normal conditions. Steric hindrance caused by the antimitotics vinblastine and colchicine result in tubulin being constrained to a curved conformation preventing MT polymerization, and promoting the formation of macrotubules. At high concentrations anesthetics are shown to have the same effect, indicating the possibility of a similar mechanism, thus highlighting these putative binding sites. Indeed, our prediction of the interference of halothane with the binding of colchicine to tubulin has already been confirmed by an experimental study [64].
While anesthetic binding in protein hydrophobic pockets is the expected mode of interaction, we found comparable binding energy estimates between the C-tails and the tubulin body indicating an interaction between halothane and the tubulin Ctermini. Our results show an increase in binding energy with more surrounding residues of the tail suggesting that halothane may sequester the C-termini into compact forms. However, solvent effects are expected to play a prominent role in this interaction and docking alone cannot provide the detailed energy evaluations needed to investigate this mechanism. The role of C-termini in normal MT dynamics, including polymerization, post-translational modification, and interaction with MAPs, marks this potential interaction as one of interest.
The interaction of volatile anesthetics, including halothane, with tubulin and MTs is of great interest both for the mechanism of action of anesthetic gases, and for post-operative cognitive dysfunction and neurodegenerative diseases that present with dysfunctional neuronal cytoskeletons [31]. Neurons are unique as they are non-dividing cells, and thus their MTs are not required to repeatedly assemble and disassemble as in mitotic spindles. Part of this is due to the neuronal cytoskeleton, which is highly architectured and relatively stable. This stability both serves to prevent division and maintain neuron morphology, which is essential to overall neural function. Halothane has been shown to bind to tubulin, alter tubulin polymerization, and disrupt polymerized MTs [40][41][42][43][44][45][46][47]. This suggests a potential cause of anesthetic induced exacerbation of neurodegenerative disorders. In regards to postoperative cognitive dysfunction, our results suggest that anesthetics may not alter MAP binding directly to release tau and increase NFT formation. However, altered MT polymerization, either through the modification of intradimer/ longitudinal contacts, C-terminal tail dynamics, or through effects on cofactors relevant to polymerization, such as GTP/GDP or zinc, may result in either net loss of polymerized MTs, or MTs polymerized into aberrant forms. Due to the overall reduced number of normal MTs, an excess amount of free tau could arise. This free tau would be vulnerable to hyper-phosphorylation leading to an increase in NFT formation. Further investigation would be required to test these hypotheses.
The role of MTs in cell division marks tubulin as a prominent chemotherapeutic target. As such there exists the possibility of adverse drug reactions between volatile anesthetics, cancerous cell lines and antimitotic agents [93]. Thus, our results may not only shed light on the role of anesthetics in non-dividing neuron cells, but also on dividing cells, especially cancer cells, due to the prominent role of tubulin based MTs in these two classes of cells. Due to the potential involvement in post-operative cognitive dysfunction, and the potential adverse drug reactions between anesthesia and chemotherapy the volatile anesthetic-tubulin interaction warrants further investigation.

Molecular Dynamics (MD) Simulations
Tubulin body. Human tubulin isotypes TUBA4A and TUBB, being the most prevalent forms in human tissues [94] were modeled according to previous homology methods [95,96]. Initial set-up was performed in the leap and antechamber modules of AMBER9 [97][98][99] using the AMBER03 force field. To obtain globally minimized structures fifteen high temperature implicitlysolvated MT dimers were generated at a temperature of 5,000 K then slowly cooled to a target temperature of 300 K, using steps of 1,000 K over periods of 400 ps, then using a 0.8-power law cooling-schedule for temperatures below 2000 K. One of the resulting structures was chosen on which to perform MD simulations.
MD simulation included two adjacent protofilaments composed of two tubulin heterodimers. Tubulin dimers were placed in a configuration consistent with MT geometry [100]. Periodic boundary conditions were used to model the system as two periodic protofilaments. The protofilaments were aligned with their axis parallel to a periodic box of the length of one tubulin dimer, 81.2 Å . These were then relaxed for 300 steps without solvent, and another 300 steps with implicit solvent. Afterwards the structures were heated to 400 K, and then cooled to 300 K over a period of 400 ps to obtain further energy minimization.
Na + ions were added to neutralize the system, and following this 80 Na + and 80 Cl 2 ions were added to bring the ionic concentration to 100 mM. The system was placed in a periodic box with dimensions 81.2 Å by 107.5 Å by 152.5 Å , filled with pre-equilibrated explicit TIP3P water using the leap module of AMBER9 [97][98][99].
To relieve steric clashes with water the structure was further energy minimized in GROMACS 3.3.2 [101][102][103] using 5000 steps of steepest descent, followed by 1000 steps of conjugate gradient minimization. Following equilibration for 5 ns with a time step of 2 fs, MD simulation, with periodic boundary conditions at 300 K and constant pressure was performed over a period of another 5 ns with the same time step. The SHAKE algorithm [104] was used to constrain bonds involving hydrogen in all simulations and a non-bonded cutoff of 12 Å was used.
Tubulin C-terminal tails. The 1JFF [61] model of tubulin was repaired with the missing residues of 1TUB [105]. The repaired 1JFF preparation was performed in the leap module of AMBER9 [97][98][99] using the AMBER99SB force field. The structure was explicitly solvated with TIP3P water in a 25 Å box extending from the protein surface. Thirty-five Na + ions were added to neutralize the protein followed by the addition of 107 Na + and Clions to bring the salt concentration to 100 mM. The structure was then energy minimized with a conjugate gradient method using NAMD [106] over 40,000 time steps.

Volatile Anesthetic Binding Site Prediction
Trajectories from the 5 ns MD simulation were clustered into dominant conformations of the protein body using the g_cluster utility of the GROMACS [101][102][103] 3.3.2 program package with the single linkage method with 1 Å RMSD similarity cutoff comparing positions of all atoms, with C-terminal tails excluded. The average structures of the resultant clusters were subjected to further analysis. Percent simulation (% Simulation) of the average structures were calculated as the ratio between the number of timestep snapshots belonging a given cluster divided by the total number of timesteps in the MD simulation. Binding sites on the clusters were predicted using the surface geometry program PASS [51] modified specifically for volatile anesthetics [49]. The PASS [51] program was run on both the tubulin dimers in the protofilament conformation, as well as on each dimer separately to determine if MT geometry affected the prediction of a site.
Predicted binding sites between clusters were grouped via a Density Based Spatial Clustering of Applications with Noise (DBSCAN) [52] method with a minimum group size of 1 predicted site, and nearest neighbor distance of 5 Å . The center of each group was taken as the average position of all predicted sites within a given group. Not all predicted sites were found on all dominant conformations. Persistence of each site was determined by taking the sum of the simulation percentages for each of the dominant conformations on which the site was predicted.

Halothane Docking
Tubulin body. Focused docking runs were performed for halothane at all predicted sites for the middle structures of all clusters of tubulin. Halothane geometry was parameterized according to an ab initio structure calculation [108]. Docking was performed via AUTODOCK [109] 4.0 using a slow focused docking protocol [110], and box size of 10 Å per and grid spacing of 0.375 Å . Binding poses were clustered with a RMSD of 2 Å .
C-terminal tails. Blind docking runs were performed for halothane against each of the 50 tail conformations for each tubulin isotype. Halothane, as parameterized above, was docked against the C-terminal tails via AUTODOCK [109] 4.0 using a blind docking protocol [111,112]. Box sizes were adjusted to accommodate the C-terminal tail conformation with constant grid spacing of 0.375 Å . Binding poses were clustered with a RMSD of 2 Å .

Microtubule Polymerization Assay
All protein and reagents required for the assay were purchased from Cytoskeleton Inc. (Denver, CO). 97% pure lyophilized bovine brain tubulin reconstituted in de-ionized water (10 mg/ml) was used. 200 ml aliquots were snap frozen in liquid nitrogen to maintain maximum protein functionality. Assays were carried out in a 96 well plate. Plates were prewarmed at 37uC for 30 min prior to experiment. Using a 95% ethyl alcohol stock solution halothane (2-Bromo-2-chloro-1,1,1trifluoroethane) (Sigma Aldrich, St. Louis, MO) and paclitaxel (Hospira, Lake Forest, Ill) were prepared to 11 times final concentration in General Tubulin Buffer (80 mM PIPES, MgCl 2 , 0.5 mM EGTA, pH 6.9). The alcohol concentration of the control was adjusted accordingly. 10 ml of halothane concentration, General Tubulin Buffer (control), paclitaxel and paclitaxel with halothane were preheated in separate wells of 96 well plate for 2 min. Immediately before running an assay tubulin was thawed in 37uC water bath and put on ice. Tubulin protein was diluted in 450 ml of ice cold Tubulin Polymerization Buffer (80 mM PIPES, 2 mM MgCl 2 , 0.5 mM EGTA, pH 6.9, 1 mM GTP 10.2% glycerol). 100 ml of tubulin and polymerization buffer solution were added to each well containing 10 ml for a final volume of 110 ml. Assay were run for 1 hour at 37uC, with absorbance measurements taken every 30 s at 340 nm using a SPECTROmax 190 (Molecular Devices, Sunnyvale, CA) plate reader and collected using SOFTmax Pro version 4.0 software.