Tubulin response to intense nanosecond-scale electric field in molecular dynamics simulation

Intense pulsed electric fields are known to act at the cell membrane level and are already being exploited in biomedical and biotechnological applications. However, it is not clear if electric pulses within biomedically-attainable parameters could directly influence intra-cellular components such as cytoskeletal proteins. If so, a molecular mechanism of action could be uncovered for therapeutic applications of such electric fields. To help clarify this question, we first identified that a tubulin heterodimer is a natural biological target for intense electric fields due to its exceptional electric properties and crucial roles played in cell division. Using molecular dynamics simulations, we then demonstrated that an intense - yet experimentally attainable - electric field of nanosecond duration can affect the bβ-tubulin’s C-terminus conformations and also influence local electrostatic properties at the GTPase as well as the binding sites of major tubulin drugs site. Our results suggest that intense nanosecond electric pulses could be used for physical modulation of microtubule dynamics. Since a nanosecond pulsed electric field can penetrate the tissues and cellular membranes due to its broadband spectrum, our results are also potentially significant for the development of new therapeutic protocols.

Being able to control protein-based cellular functions with an electromagnetic field could open an exciting spectrum of possibilities for advancing biotechnological processes. Besides, it paves the way for the development of new biomedical theranostic approaches to treat various diseases where specific proteins are known targets. Electrostatic interactions in proteins are of paramount significance for protein function 1 . Strong molecular electric fields are known to play an essential role in protein folding 2 , protein-ligand and protein-solvent interactions 3 as well as protein-protein interactions. Furthermore, the enzymatic activity of proteins 4 exploits local molecular electric fields to affect the potential energy surface of the reaction involved. A molecular electric field can also become crucial for protein stability since even a small imbalance in electrostatic interactions can cause malfunction of the protein 5 . Therefore, it is reasonable to assume that external electric fields (EFs) with appropriately chosen parameters of strength, frequency, and duration could modulate protein function. To overcome thermal noise effects and avoid heating side-effects, short (<100 ns) intense (>MV/m) electric pulses lend themselves as a suitable form of the electromagnetic field that can be utilized to modulate protein function 6,7 . Indeed, it has been demonstrated through molecular dynamics simulations that EFs can affect conformation of pancreatic trypsin inhibitor 8 , insulin [9][10][11] , lysozyme [12][13][14][15] , β-amyloid and amyloid forming peptides 16,17 , and soybean hydrophobic protein 18 . Further, EFs also unfolded myoglobin 19,20 , induced transition of peptides from a β-sheet to a helix-like conformation 21 , and caused structural destabilization of (a short peptide) chignolin 22,23 in molecular dynamics simulations. Moreover, recent studies demonstrated that EFs can affect water diffusivity and ion transport across transmembrane proteins such as aquaporins [24][25][26][27][28] , ion channels 29 and voltage sensors 30 . Experiments showed that EF applied directly by electrodes can catalyze the reaction 31 in a similar manner as a molecular electric field at enzyme sites 32 and switch protein conformational states 33 . Proteins from the tubulin family seem to possess an unusually high magnitude of charge and dipole electric moment 34,35 , thus possibly being a suitably sensitive target for the action of EFs. αand β-tubulin monomers, which form stable heterodimers, are also crucial components of the cytoskeleton structures called microtubules, which are essential for cell division 36 among many other roles they play in living cells.
Several experimental works have been reported, which demonstrate the effects of EFs on microtubule structures. It has been shown that EFs in the intermediate frequency range (100-300 kHz) had a profound inhibitory effect on the growth rate of a variety of human and rodent tumor cell lines 37 and also in vivo in the case of human brain tumors 38 presumably by interfering with the polymerization of mitotic spindle microtubules which are composed of tubulin dimers. In recent work, ultra-short intense pulses of EFs were shown to cause Ca 2+ independent disruption of dynamic microtubules in glioblastoma 39 . However, it is not clear whether the observed effect of EFs was direct or indirect through the action on the membrane channels first and then transmitted downstream into the cell interior. Hence, an exact molecular-level mechanism of this electric field action on microtubules and tubulin remains unknown.
Microtubules, due to their expected special electric 40 and vibrational [41][42][43][44] properties, were proposed to be involved in endogenous electrodynamic processes in cells [45][46][47] . However, all-atom molecular simulations of external EF effect on tubulin have been carried out only recently 48,49 . They specifically investigated the EF effects on protein mechanics but did not include the C-terminal tail. The C-terminal tail is a highly flexible unstructured domain of tubulin, which (i) is essential for tubulin-protein interactions 50,51 , (ii) is the main site of protein mutation and post-translational modifications 51 , and (iii) greatly contributes to the overall electric properties of tubulin accounting for approximately 30-40% of the total charge 35 . A very recent study 52 presented results of the molecular dynamics simulation analysis of the electric field effect on tubulin. This recent study investigated electric field strengths in the range between 50 and 750 MV/m, which overlaps with the values used in our simulations, and carried out simulations over 10 ns while our study reports on simulations that ran up to 30 ns. The previous study only examined displacement effects on key secondary structure motifs such as tubulin's C-termini and alpha helices. In contrast, in the present study, we aim to unravel the mechanisms of interaction of intense EFs on the tubulin protein family, given their essential and attractive role as drug targets for cancer therapy, due to their involvement as key-players in cellular self-organization. It is already known that various stabilizing/ destabilizing tubulin-binding drugs such as taxanes, colchicines, and vinca alkaloids, bind to different sites on the tubulin dimer, modulating microtubule-based processes 53 . This type of binding is assumed to be based on purely electrostatic interactions like those employed in the recognition between proteins 2 . Conversely, what is presently unexplored, is the possibility to modulate such electrostatic environment through EF. Therefore, we first employ bioinformatics tools to systematically compare electric charges and dipole moments of the various isoforms of tubulin to the so-called PISCES set, which represents all unique chains in the whole protein database. Then, we quantitatively evaluate the effects of intense nanosecond EFs on the overall shape of the tubulin dimer as well as the time-dependent evolution of its dipole moment. Finally, we closely analyze the dipole moment of individual residues with particular attention to those forming the binding pockets for the most important tubulin-modulating pharmacological agents (e.g., paclitaxel or vinca alkaloids) and GTP/GDP molecules. Thus, in this paper, we cast light on the mechanism of direct EF action on tubulin at the molecular level.

Results
proteins from tubulin family have exceptional electric properties. In this study, we first compared the structural charges and dipole moments of proteins from the tubulin family to those of the PISCES set of proteins (all unique protein chains available in the RSCB protein database). It can be seen from Fig. 1a that tubulin proteins possess a generally much higher structural electric charge (with a mean value of −22 e per monomer) than most PISCES proteins (mean value of −4 e). This is likely due to an excess of acidic over basic amino acid residues present among the tubulin proteins compared to the PISCES set of proteins, see SI 1. Taking as an example the 1TUB tubulin structure, we found that while a relative content of acidic residues is similar in tubulin versus the average from the PISCES set, tubulin contains fewer basic residues than an average protein. It is worth noting that a large fraction (on average around 34 % across the tubulin dataset) of electric charge of tubulin is due to the unstructured and highly flexible C-terminal tail (CTT). Such high electric charge is remarkable since the www.nature.com/scientificreports www.nature.com/scientificreports/ CTT of tubulins is rather short -having on average around 20 residues. A rather high value of the electric charge of tubulin CTT is due to a high content of acidic residues. Figure 1b shows that tubulins possess a much higher dipole moment than PISCES proteins: with a mean of 2,166 debye versus a mean of 555 debye, respectively. The high dipole moment arises from an asymmetric electric charge distribution, i.e., an asymmetric distribution of acidic versus basic residues. This asymmetric charge distribution is partially due to the fact that a significant fraction of charge is located on the C-terminal tail. However, even the tubulin structures where the CTT is not resolved, possess a rather high dipole moment, e.g., both 1TUB and 1JFF crystal structures result in a dipole moment above 2,000 debye. The results in Fig. 1 demonstrate that the tubulin proteins have exceptional electric properties. Hence it is reasonable to analyze in detail if these properties can be exploited to manipulate tubulin's structure and hence its function using an external EF with specifically designed characteristics.
Rounding effect on tubulin shape induced by electric field. We have employed molecular dynamics (MD) simulations, see Methods for details, to analyze the effects of an EF acting on tubulin. The tubulin structure used in these simulations consists of αand β-monomers, which form a stable noncovalently bonded heterodimer, see Methods for details.
In the MD simulation performed, we have added an EF as an external Coulomb force acting on every atom in the system. We analyzed the effect of the EF strength starting from 20 MV/m up to 300 MV/m. The rationale for the selection of this range of field strength is both theoretical and experimental. Theoretically, the interaction , where μ is the tubulin dipole moment vector and E the electric field vector, exceeds thermal energy for field strengths in the range >MV/m. The range <70 MV/m is also experimentally attainable since it is below the field strength of dielectric breakdown of water-like media with an exact value depending on the ionic strength and the electric pulse duration 54,55 . The length of the simulation was selected to be up to 30 ns. This time scale is short enough to be computationally tractable and corresponds accurately to the duration of the electric pulse, which can be applied in experiments. Furthermore, electric pulses of this duration typically do not cause any appreciable heating in experiments conducted even for the experimentally attainable field strengths mentioned above.
The most basic integral parameter of protein is its shape. Hence, we first analyzed how the EF affects the shape of the tubulin dimer, see Fig. 2 (note that no water molecules and counterions are displayed). To this end, we approximated the shape of tubulin by an ellipsoid and obtained a new coordinate system within the frame of this ellipsoid, see Methods. This approach allows a sound analysis of protein polarization process, excluding the trivial roto-translational effects taking place due to the external EF. At first, we observed that the effective shape of the equivalent tubulin ellipsoid is affected by a 100 MV/m EF in a way that the medium axis is elongated by 11% and the long axis is shortened by 2%, see Table 1.
We also observed that tubulin undergoes rotation -this will be discussed in more detail in the next section as the rotation is related to the dipole properties of the tubulin analyzed there.
Orientation effect of electric field on tubulin dipole moment. As the next step in this investigation, we focused on the effect on the electric dipole moment of the tubulin dimer, see Fig. 3. It is readily seen that the tubulin dipole moment under zero-field conditions has an average value of around 2,500 debye (Fig. 3a). However, in the presence of a 100 MV/m EF, the dipole moment is increased to more than 6,500 debye, i.e., it more than doubles in the process. In the Fig. 3b the polarized C-terminal tails extending from the structure are also visible. To get more in-depth insight into time evolution and field strength dependence, we further analyzed the tubulin y0 dipole component, the most affected one (i.e., with the highest polarization change) among the three dipole components shown in Fig. 3. The external electric field vector was oriented in the z-direction (Cartesian reference system, see Fig. 2a) during the whole course of the MD simulation while the initial orientation of the tubulin dipole had an anti-parallel component to it. At the zero electric field strength the effective y′ component of the tubulin dipole moment fluctuates around 2,000 debye, with a somewhat steady profile for the whole simulation time (see Fig. 4). This stable dipole value evolution is due to the coordinate system transformation adopted (i.e., from the external Cartesian system to the internal protein system), which removes the dipole roto-translational effects, highlighting internal polarization effects. The EF tends to act by exerting a torque on the tubulin's dipole moment so that the dipole component in the direction of the field vector is maximized. Data reported in Fig. 4 suggest that EF in the range of 20 MV/m induce a slow, but significant, polarization in the last 5 ns of simulation, so that the y′ dipole component increases from 2,000 to 3,000 debye. Higher external fields amplify and speed up this polarization effect up to a value above 5,000 debye after 7 ns for a field strength of 100 MV/m. Moreover, a second field-effect appears on the β-tubulin CTTs. While in the zero field simulation and within approximately 20 ns of 20 MV/m condition, the CTT remains more or less close to the surface of tubulin, for the field strengths ≥50 MV/m the CTT is pulled away from the tubulin surface and becomes outstretched (see Fig. 3).
No unfolding effect on tubulin up to 100 MV/m. The fundamental and functionally-crucial characteristics of any protein structure are the type and location of its secondary structure motifs. Therefore, we have focused next on the analysis of how an EF affects the count of tubulin residues being part of a certain secondary structure motif. The results of this analysis can be seen in Fig. 5, which displays the count of residues in alpha helices, beta-sheets, coils and turns as averaged over the last 5 ns of the corresponding trajectories. The field strength of 100 MV/m seems not to cause any effect on the structure -just slightly increasing the content of the tubulin residues present in coils at the expense of beta-sheets, alpha helices and turns. Only 300 MV/m induces a substantial effect on tubulin secondary structures, where more than 70 residues are converted from alpha-helices to coils. Such conversion is a sign of unfolding, where charged CTTs are experiencing an electrostatic force, which is strong enough to undergo pulling out from the tubulin body. Since there are two alpha helices (H12 and H11) www.nature.com/scientificreports www.nature.com/scientificreports/ close to CTTs, those are the first ones which are unfolded. Furthermore, we selected the CTT of β-tubulin since it is longer and carries more charge than the α-tubulin CTT for a more detailed view. In Fig. 6a, we quantified the count of CTT residues being in coil structures within the 30 ns time scale for 0, 20, 50, and 100 MV/m field strength conditions. We can see that the 20 MV/m condition does not tend to change the secondary structure of CTT compared to the no field condition. However, both 50 and 100 MV/m field strengths tend to increase the  www.nature.com/scientificreports www.nature.com/scientificreports/ count of residues being in a coil structure at the expense of α-helices and turns. In these cases, the higher the field strength, the shorter the time needed to achieve the maximum amount of residues in a coil structure.
Electric field affects the dipole moment of specific tubulin domains. C-terminal tails of tubulin are essential for interactions with microtubule-associated proteins, such as motor proteins 56 . Since the data in Fig. 3 suggested substantial effects on CTTs of tubulin, we next focused on a detailed analysis of CTTs. In panel B of Fig. 6, we found a strong polarization effect on the CTT of β-tubulin, with a mean dipolar shift for field strength ≥50 MV/m as high as 40 debye.  www.nature.com/scientificreports www.nature.com/scientificreports/ Apart from the EF effect on CTTs, we also asked a question whether a strong field can affect local electrostatic conditions of the tubulin dimer. To this end, we analyzed the shift of the dipole moment of every residue of the tubulin dimer. We plot the y′ component (internal tubulin coordinates) of the dipole moment in Fig. 7. It can be seen that some residues undergo a substantial change of the dipole moment when exposed to an EF of 100 MV/m. In particular, when considering both monomers, it is possible to appreciate dipolar shifts ranging from −10 up to almost +15 debye. Specifically, 22 out of the 451 residues of α-monomer residues exceed a ±5 debye shift, while for the β-monomer the concerned residues are 3 out of 445 ones. To understand if dipole moments in any functionally significant parts of the tubulin dimer are affected, we focused on selected important sites of the tubulin dimer: GTP binding/hydrolysis site, sites of tubulin longitudinal interactions, and the binding sites of the three most common tubulin drug families (paclitaxel, nocodazole/colchicine, and vinca alkaloids). The Methods section provides a rationale and identification of residues belonging to the tubulin sites based on energy considerations. In Fig. 8, we highlight the location of the important tubulin sites in a color-coded manner and also provide a list of the selected residues. We display the histograms of the dipole moments (y′ component) of six selected residues where we found the strongest effects. In all of these six cases, we see a field strength-dependent effect, shifting the y′ component of a residue dipole moment towards different values. Since the local electrostatic field is crucial in active sites of many enzymes enabling the process of catalysis, we speculate that influencing the local field could affect the GTP hydrolysis rate. We see that α-Asp 251, a crucial residue mediating GTP hydrolysis 57 , has its y′ component of the dipole moment (yDM in short) affected by almost 3 debye when comparing no  www.nature.com/scientificreports www.nature.com/scientificreports/ field and 100 MV/m conditions. Several residues mediate tubulin-tubulin longitudinal interactions. β-Arg 391 is one of them and has its yDM also influenced by the EF (Fig. 8). However, in this case, the field strength has opposite sign effect: 20 MV/m field tends to affect yDM in the opposite direction than the field with 50 and 100 MV/m field strengths, compared to no field condition. This is probably due to the position of the residue in the β-monomer (see Fig. 8, left-upper panel); in fact Arg-391 is located in the external part of the protein, presumably in contact with a layer of water; therefore its response to the external EF could be screened by water for lower field strength. The second most important residue for binding paclitaxel energy-wise, β-Val 23 (Fig. 8), also manifests the changed yDM up to few debyes when an EF of 100 MV/m is applied. Furthermore, residue β-Cys 239 belonging to the colchicine binding site has its yDM shifted towards zero with increasing field strength. We also found that an EF influences α-Pro 325 and β-Asp 177 belonging to the binding site of the vinca alkaloids family. In this case, however, the shift of yDM is only in the sub debye range.

Discussion
We have demonstrated that by using MD simulations combined with the Covariance Matrix method, it is possible to study with high precision the intrinsic structural and dipolar response of tubulin protein under the action of nanosecond scale EFs. The findings presented here cover several aspects of such tubulin/EFs interaction mechanisms. First, we quantified the unusually high magnitude of structural charge and dipole electric moment of tubulin family proteins in the absence of EFs and, specifically, the response of tubulin to nanosecond scale EFs. Second, we studied structural/conformation changes and possible unfolding effect on the whole protein. Finally, we provided evidence of dipolar coupling with the EF of residues forming the binding pockets for the most important tubulin-modulating pharmacological agents (taxane, vinblastine, and colchicine) as well as tubulin CTTs. In the following, we discuss the significance of such results, as well as the limitations and point out some possible extensions for future work. electric and dipolar properties of tubulin. Regarding the charge and dipolar properties of tubulin family, we showed that an average (across the list from various species) tubulin monomer has an approximately 4-5 fold higher structural electric charge and electric dipole moment than an average protein (Fig. 1). This result corroborates our further findings: tubulin seems to be more susceptible to EFs compared to other proteins analyzed earlier 11,12,19 , in a way that lower field strength, i.e., 20 MV/m, can attain a 50% increase of tubulin dipole. We should add that our bioinformatic analysis of the PISCES proteins (unique chains) is based only on the proteins which do www.nature.com/scientificreports www.nature.com/scientificreports/ have their structure determined. While there is no better and more straightforward approach, it strongly underrepresents membrane proteins and intrinsically disordered proteins.
In a recent paper 58 , a detailed energy balance calculation was provided for the stability of a microtubule with a seam (representing a so-called type B microtubule lattice). The calculations included the contributions from dipole-dipole interactions between tubulin dimers, solvent accessible surface area, van der Waals and electrostatics (generalized Born approximation) to demonstrate that the balance energy is such that a single GTP hydrolysis event can trigger a microtubule disassembly because when the seam is closed with GTP molecules attached to the β monomers, the net free energy is −9 kcal/mol. The dipole-dipole energy is positive (destabilizing) and amounts to 27 kcal/mol. When the seam becomes open due to GTP hydrolysis, the net free energy becomes vastly positive. These calculations used a dipole moment of approximately 4,500 debye for a tubulin dimer as an average value over various isotypes. In the present paper, we have shown that strong electric fields can substantially increase the dipole moment of tubulin. Here, we have shown the dipole moment to be close to 3,000 debye in zero field, which would translate into dipole-dipole energy for a microtubule lattice of 12 kcal/mol and net free energy of −24 kcal/mol, hence slightly more stable than the calculation provided in 58 . However, at strong electric field values (≥50 MV/m), the corresponding dipole moments were found to increase to as much as 6,000 debye, which translates into dipole-dipole energy of 48 kcal/mol and net free energy of +12 kcal/mol making the lattice unstable. This quantitatively supports our hypothesis that sufficiently strong electric fields disrupt a microtubule lattice by increasing the dipole moments of tubulin and contributing positive energy that cannot be balanced by the remaining contributions.
We need to add that the time scale of these simulations and the time scales of structural changes in microtubules are several orders of magnitude apart from each other and hence it is not possible to make a direct connection between the dynamics of microtubules occurring on a time scale of seconds to minutes and the nanosecond duration of pulses and the length of MD simulations. Our ongoing experimental work suggests that the train of a few hundred electric pulses of comparable field strength and duration as in our simulations modify tubulin up to the several minutes and they affect tubulin assembly to microtubules.
Structural/conformation changes on tubulin induced by electric fields. By approximating the tubulin with the ellipsoid as given by the covariance Matrix method, we were able to appreciate the actual shape changes under the action of the EFs. The protein slightly reduces its major axis, while increasing its medium one. This effect seems of particular interest since it is a consequence of protein polarization. The tubulin dimer undergoes a packing transition, rather than an elongation, since its dipole moment aligns along tubulin medium axis (see Fig. 3). This non-trivial effect is due to the high charge of tubulin CTTs. Tubulin seems to have a lower threshold for the unfolding transition 59 . For example, unfolding effects tend to appear at lower field strengths or within a shorter time frame (unfolding at 250 MV/m started at approximately 160 ns in insulin (see Fig. 1c of 11 ) compared to only a few ns at 200 MV/m for tubulin). Such faster unfolding is reasonable since the higher the charge and dipole moment of the structure, the higher the electric force and torque, respectively, acting on the protein.
Dipolar coupling of specific tubulin domains with the electric field. We found that a primary target of the electric field in the tubulin heterodimer are C-terminal tails since they carry a substantial amount of electric charge. Indeed, they manifest a strong polarization effect on the CTT of β-tubulin, with a mean dipolar shift for field strength ≥50 MV/m as high as 40 debye. We might ask a question: how biologically general is the effect on the CTT we observe, or put another way, is the CTT sequence we used common in other biological species?
An exact tubulin sequence and structure varies across biological species and tissues 35 . While the tubulin core is rather conserved across species, CTTs display higher variability across species and are also a target of post-translational modifications forming a so-called "tubulin code" 60,61 . The β-tubulin CTT sequence we used in our simulations has a 100% sequence identity to specific tubulin isotypes found for example in pig Sus scrofa (common experimental source of tubulin) tubulin gene TBB and to many tubulin isotypes of more than 20 species, see SI 3. It also has a 90% sequence similarity to CTT of TBB2B gene in the human, the only difference is that two immediate neighbor residues Gly 440 and Glu 441 are swapped, so the total charge of the CTT is the same as that of the CTT in our tubulin structure, see SI 3.
As a general result, we observed a widespread EF effect on protein residues. In particular, when analyzing specific tubulin domains, e.g. the GTP binding/hydrolysis site, sites of tubulin longitudinal interactions, and the binding sites of the three most common tubulin-binding drug classes (paclitaxel, nocodazole/colchicine, and vinca alkaloids), we obtained unexpected dipolar coupling even in the presence of a 20 MV/m acting for 30 ns. In this context it is interesting to note that at present, generators able to produce intense (>20 MV/m) ultra-short (30 ns) electric field pulses, are commercially available, hence the interesting perspective to control and modulate protein functions is becoming realistic. For instance, one can envisage that pulsed electric field could be used in synergy with taxane-based cancer treatment protocols to modulate the drug binding to tubulin and hence the dynamics of microtubules. It is worth noting that the most negatively charged part of the tubulin dimer, when embedded in a microtubule, is the cytosol-facing outer surface while the surface facing the lumen is less negatively charged. Also, the dipole moment of the dimer also includes a contribution from the CTT, which is highly variable due to the flexibility of this motif. Finally, the dielectric breakdown, which would pose a limitation on the strength of the applied field in practice, depends on the ion strength of the solution and the duration of the applied pulse 54,55 . Limitations. In discussing the implications of our modeling effort for experiments on microtubules, it could be argued that field strengths >100 MV/m are not readily experimentally attainable yet. However, one has to keep in mind that MD simulations are commonly best suitable for providing a relative assessment of various effects and not necessarily absolute values directly translatable the experimental situation. For example, the binding free energies of ligands interacting with proteins are typically off by a factor as large as two or more 62 . Hence, effects such as unfolding, which are computationally predicted to occur at large fields, may, in fact, require field strengths lower than predicted for tubulin 59 . Moreover, rapid technological advances in the field may well lead to the development of engineering solutions with sufficiently high electric field strengths attainable in the clinical setting. The ultimate physical limit, however, is the field strength at which the dielectric breakdown of the exposed biological material occurs.
Obviously, electronic degrees of freedom are not included in classical molecular dynamics simulations, and hence there is indeed an aspect that is not revealed by our present modeling effort in the context of electric field effects. This aspect requires the use of more sophisticated methodology such as density functional approaches or QM/MM simulations, both of which are not yet suitable for an extensive molecular system such as a tubulin dimer. However, as suggested in a paper dedicated to this problem 63 a simple use of a scaling factor may introduce a suitable correction to molecular dynamics results. Concerning the possible use of a polarizable force field, recently biological force fields have heavily relied on experimental NMR data to refine their parameters for example for nucleic acids; 64 despite this, deriving accurate force field parameters is challenging because of the large parameter space, non-linear interdependencies of parameters and limitation in the amount and quality of experimental and ab initio reference data 65 and there are cases where polarizable force fields are comparable to or even worse than non-polarizable force fields, due to the poor quality of parameters in the polarizable force fields 66 . The adopted force field with fixed atomic charges can provide the polarizability of the molecule due to semi-classical motions (i.e., structural deformations), missing only the possible and typically much weaker effects on electronic distribution induced by the electrostatic interactions. However, these subtle effects, that are somewhat limited for an electron in ground state conditions are difficult to be included in atomistic simulation procedures and often polarizable force fields aiming to model such effects might not always provide more accurate and/or reliable results as they require the use of further empirical/adjustable parameters.
Challenges for future work. In our future work we intend to better calibrate the simulation parameters in order to make the model quantitatively predictable so that such characteristics as field strength, frequency and duration may be tailor-designed to elicit a specific response of the target, namely the tubulin dimer or an entire www.nature.com/scientificreports www.nature.com/scientificreports/ microtubule. This modeling effort may need to be extended to include quantum effects through the use of quantum mechanics/molecular mechanics calculations especially to be able to elucidate covalent bond response to electric fields, for example at the GTP binding site.
Additional work is ongoing to address further questions of how the electric field affects the interaction of microtubule-associated proteins, such as motor proteins docked on the tubulin subunits of a microtubule. In general, it is very likely that the interaction of heads of the motor proteins will alter the currently observed effect on the C-terminus on the β-tubulin since the balance of the electrostatic forces will be altered. For instance, in the case of kinesin motor protein, the overall kinesin net charge depends on the kinesin type and can vary from −2 to +5 e 67 . Since the tubulin C-terminus is negatively charged, a kinesin head could tend to either attract or repulse the tubulin C-terminus. X-ray crystallographic studies directly confirmed the presence of direct interaction between tubulin's C-termini and kinesin involving electrostatic attraction between the oppositely charged domains on the two proteins 68 . A direct determination of the tubulin's C-termini involvement was provided by cleavage assays using subtilisin. These results show that the C-termini directly modulate the motor-tubulin interface and the binding properties of motors. Consequently, C-termini cleavage increases motor protein's binding stability so that kinesin adopts a binding conformation close to the nucleotide-free configuration and hence loses its processivity 56 . Various tubulin isoforms differ in the sequences of C-termini, which, together with post-translational modifications of C-termini is known to serve to modulate kinesin processivity 69 .
To conclude, our results allow us to draw biologically-relevant consequences in terms of microtubule stability and the changes in the strength of binding for tubulin-binding drugs under the influence of EF. Therefore, we believe the present paper provides new and vital quantitative insights into the effects of EF pulses of nanosecond duration on tubulin and microtubules. The knowledge of residue-specific EF realignments allows extrapolating the computer simulation results in terms of their consequences on the behavior of other tubulin isoforms and tubulin mutants under a both exogenous and intrinsic molecular electric field.

Methods
We follow a gapless numbering convention of tubulin residues in the current paper, in contrast to a gapped numbering in the Protein Data Bank entries 1JFF and 1TUB, which contain 2 gaps after β-Leu 44 and 8 gaps after β-Pro 360.
Comparison of tubulin family proteins to pIsCes set of proteins. Charges and dipole moments of proteins from the tubulin family were obtained from our earlier work 35 , which includes 246 isotypes of tubulin, mostly αand β-tubulin, from various species. PISCES set, a representative set of all proteins in the RSCB PDB database, was obtained using PISCES server 70 . We used following settings to run the server query (date 2018-05-18): sequence percentage identity threshold ≤90, Resolution: 0.0-4.0, R-factor: 0.5, sequence length: 40-10,000, Non X-ray entries: Included, CA-only entries: Excluded, Cull PDB by entry, Cull chains within entries: No. We obtained a list with 36,331 entries, see SI 4. Charges and dipole moment data for our PISCES set of proteins were obtained from Protein Dipole Moments Server of Weizmann Institute https://dipole.weizmann.ac.il/index. html 71 . Each protein chain was evaluated separately which leads to a total number of 63,110 records for both protein chain charge and dipole moment, see SI 5. The effective pH considered is 7 71 .
Note that the methods for charge and dipole moment calculation used by Tuszynski et al. 35 (GROMOS96 43a1 force field) and the dipole moment server 71  tubulin structure. The structure of the GMPCPP-bound tubulin was obtained from the Protein Data Bank under the PDB ID: 3J6E 72 . The cofactor GMPCPP, a non-hydrolyzable analogue of GTP, was modified into GTP. Furthermore, the CTT, which are usually missing in PDB structures of tubulin, were added to both αand β-tubulin subunits. This was achieved using the Molecular Operating Environment software (MOE, Chemical Computing Group Inc.) 73 . MOE was also used for the addition of hydrogen atoms and the prediction of ionization states. The CTTs were added in an extended conformation, and we depended on molecular dynamics simulations later to help to restructure the tails in the correct conformation in solution. The CTT is defined here as the last 16 residues of the α-tubulin (GVDSVEGEGEEEGEEY) and the last 20 residues of the β-tubulin (QDATADEQGEFEEEGEEDEA). The CTT sequence identity was tested using the BLAST tool on 74 , see the.fasta format result in SI 3.
Residues of drug binding sites. The paclitaxel, colchicine/nocodazole and vinca binding sites were identified based on the RSCB PDB structures 1JFF (paclitaxel is named as residue TA1, site identifier AC5) 75 , 1SA0 (colchicine is named as residue CN2, site identifier AC8) 76 and 5BMV (vinblastine is named as residue VLB, site identifier AE1) 77 , respectively. The residues belonging to the respective binding sites are identified there using an algorithm from 78 and are available in a respective.pdb file. For paclitaxel and colchicine, we additionally included the residues which contribute 75% binding energy between paclitaxel and tubulin based on the analysis in 79 and 80 , respectively. Those results provided additional β-tubulin residues LEU 215, GLN 280, and LEU 361 on top of the list from RSCB PDB 1JFF structure for paclitaxel and α-ALA 180 and β-LEU 246 on top of the list from RSCB PDB 1SA0 structure for colchicine. See the complete list of residues of individual interaction and binding sites in Fig. 8.

Molecular dynamics simulation.
We carried out MD simulations of a tubulin protein in water solution using the GROMACS package v. 4.6.5 81 . The simulated system consisted of a rectangular box (12 × 12 × 15 nm 3 ), where we placed a single tubulin heterodimer, 67,087 TIP3P (transferable intermolecular potential 3P) 82 www.nature.com/scientificreports www.nature.com/scientificreports/ water molecules and 47 sodium ions, resulting in a typical density of 1,000 kg/m 3 . Note that, to describe tubulin behavior properly, it was necessary to simulate a box of water molecules large enough to reproduce both the first hydration shells and a significant amount of bulk water. The construction of the simulation box, we followed the methodology we used earlier 59 , where the protein-water system was properly built to match the pressure of a pure water system with the same volume. Specifically, the pure TIP3 system at 55.32 mol/L and 300 K, corresponding to the typical experimental liquid water condition, served to set the reference pressure to be matched by the tubulin-water system. The latter system was built by solvating a single tubulin protein and then adjusting the number of water with this procedure; i) add few water molecules; ii) run a few ns MD simulation; iii) perform a pressure check (pressure is at equilibrium? Is it the same as that of pure water?); iv) repeat the previous steps, where appropriate. Following an energy minimization and subsequent solvent relaxation, the system was gradually heated from 50 K to 300 K using short (typically 50 ps) MD simulations. A first trajectory was propagated up to 150 ns in the NVT ensemble using an integration step of 2 fs and removing the tubulin center of mass translation but with no constraints on its related rotation. The temperature was kept constant at 300 K by the Berendsen thermostat 83 with the relaxation time equal to the simulation time step, hence virtually equivalent to the isothermal coupling 84 which provides consistent statistical mechanical behavior. All bond lengths were constrained using the LINCS algorithm 85 . Long-range electrostatics were computed by the Particle Mesh Ewald method 86 with 34 wave vectors in each dimension and a 4 th order cubic interpolation. The amber03 force field 87 parameters were adopted. Once obtained an exhaustive equilibrated-unexposed trajectory we evaluated possible effects induced by a high exogenous electric field ranging from 10 MV/m up to 300 MV/m (see Table 2), acting in the simulation box as explained in 88 for 30 ns in each exposure condition and applied along the z-axis of the Cartesian reference of frame. More precisely the application of the electric field takes place in continuity at the last frame of the unexposed simulation, thus allowing direct evaluation of the characteristic response over time of the system due to the switch on of the exogenous perturbation. Taking advantages of the previous findings in 19,59,89 on the predictable protein polarization process, for the highest field strengths (200 MV/m and 300 MV/m) the simulation box was enlarged along the z-direction up to 30 nm, resulting in a rectangular box of 12 × 12 × 30 nm 3 .

Covariance matrix method.
A simple method to characterize the geometry of a complex system like a protein is to approximate its overall geometrical structure, at each MD frame, to the ellipsoid given by the eigenvectors of the Covariance Matrix as obtained by the distribution of the x, y, and z coordinates of the protein atoms (the geometrical matrix  C) 19,90 . Therefore, for a system defined by N-atoms the elements of the 3 × 3 geometrical matrix at each time frame are given by

Data Availability
See SI 6 for initial and configuration files as well as molecular dynamics trajectories and scripts for post-processing.  Table 2. List of molecular dynamics simulations with the corresponding conditions. *The starting configuration comes from an equilibrium 150 ns trajectory (regular box, 12 × 12 × 15 nm 3 ) and **The starting configuration comes from an equilibrium 30 ns trajectory (big box, 12 × 12 × 30 nm 3 ).