Electro-opening of a microtubule lattice in silico.

Modulation of the structure and function of biomaterials is essential for advancing bio-nanotechnology and biomedicine. Microtubules (MTs) are self-assembled protein polymers that are essential for fundamental cellular processes and key model compounds for the design of active bio-nanomaterials. In this in silico study, a 0.5 μs-long all-atom molecular dynamics simulation of a complete MT with approximately 1.2 million atoms in the system indicated that a nanosecond-scale intense electric field can induce the longitudinal opening of the cylindrical shell of the MT lattice, modifying the structure of the MT. This effect is field-strength- and temperature-dependent and occurs on the cathode side. A model was formulated to explain the opening on the cathode side, which resulted from an electric-field-induced imbalance between electric torque on tubulin dipoles and cohesive forces between tubulin heterodimers. Our results open new avenues for electromagnetic modulation of biological and artificial materials through action on noncovalent molecular interactions.

Introduction

The ability to control the structure and function of biomolecular matter [1], [2] is crucial for the development of novel biomedical therapeutic methods and bioinspired nanotechnology [3], as well as new methods to understand biological processes [4]. Cytoskeletal proteins and their assemblies, including microtubules (MTs), are an excellent example of such biomolecular matter because they are a multipurpose platform that can be used to investigate the effects of external physicochemical modulators at the molecular level. Here, we demonstrate in silico that a nanosecond-scale intense electric field, as a new class of an oriented external electric field [5], affects the stability of noncovalent bonds in MTs, leading to deformation of the supramolecular lattice and ultimately to its opening.

MTs are attractive components for the development of biologically inspired active and self-repairing matter [6], [7]. MTs are self-assembled protein tubes, 25 nm in diameter and lengths of up to several tens of micrometers, that are present in virtually all cells and enable cell division and intracellular transport [8]. The ability of MTs to generate force and micro- and nanoscopic transport has attracted significant interest in molecular biophysics [9] and bioinspired engineering of self-assembling [10], self-repairing [11], and active matter [1], [2] and nanostructures [12].

MT-based systems are often controlled by optical forces [13]; control of the motion of individual MT-related proteins [14] and MT-based active matter has been demonstrated [1]. Electric forces have also been employed to act on the charges of a MT as it glides on a nanomotor-covered substrate and for steering and sorting of molecules [15]. However, changes in MTs induced by external electric forces have not been explored at the molecular level, either experimentally or in silico.

The effects of short (nanosecond range) intense electric pulses on MTs have been demonstrated in cells: A MT network was disrupted, either immediately [16] or with some delay [17], after electric pulses were applied. A recent study demonstrated the possibility to remodel the MT cytoskeleton in cells without a complete de-polymerization phase [18]. However, it is not clear whether the electric pulses were acting directly on the MTs or the effect was just a consequence of downstream signaling from other processes such as membrane permeabilization [19], action on membrane channels [20], or changes in intracellular ion concentration [21].

The direct effects of an intense external electric field on single proteins have been only recently explored. Most of the investigations have been performed in silico [22], where they explored either membrane [23] or protein systems. The effect of the electric field was mostly on the secondary structure, conformation, and orientation of various proteins [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], including cytoskeletal system proteins such as tubulin [43], [44], [45], [46], [47] and kinesin [48], even leading to the unfolding of some proteins [49], [50], [51]. Most of this work explored electric field effects either in a single protein or membrane-bound proteins systems. Our current study goes much beyond these works, as it explores in MD the electric field effects on a large multimeric non-membrane protein system for the first time.

Although the effect of a microsecond-scale intense electric field on the enzymatic activity and structure of proteins has been demonstrated experimentally [52], [53], the effects of the nanosecond regime have only been recently explored. Hekstra et al. [54] used time-resolved X-ray crystallography to show that the conformation of the PDZ domain of a protein can be modified by applying a 300 ns, 100 MV/m electric pulse. Although providing exquisite atomistic details, that experiment was performed on a protein crystal, so its relevance to proteins under physiological condition (in solution) is limited. In solution, Urabe et al. have shown that a train (500) of 5 ns, 3 MV/m electric pulses can physically denature urease [55] and a train (1000) of 1 ns, 200 MV/m pulses promotes the disintegration of transthyretin amyloid [56]. In addition, in our recent pioneering studies, we showed that an intense electric field has multiple effects on a single tubulin in vitro and in silico. In silico [46], we showed that the major effect of the electric field was on the highly electrically charged C-terminus of tubulin, and the shape and dipole moment were also affected. Experimentally, in addition to the structural modifications and electrical properties changes of tubulin, we demonstrated the ability of tubulin to reversibly or irreversibly self-assemble into MTs with different shapes depending on the number of applied electric pulses [3]. However, the direct effect of nanosecond-scale intense electric pulses on MTs has not been investigated in silico or in vitro. To explore this case and predict potential experimental scenarios, we performed an all-atom molecular dynamics simulation of the effect of an electric field on an MT lattice in solution.

Materials and methods

Molecular structure

Our system consisted of a ring of 13 tubulin heterodimers that were solvated and equilibrated. The ring and the solvent were constructed to exploit symmetry and utilize the effect of periodic boundary conditions to multiply the ring in all directions, resulting in an arbitrary number of arbitrarily long MTs in the B-lattice form (see Fig. 1 and Section SI for details). This system was in accordance with our previous work [57], in which we determined the binding energy between adjacent tubulin heterodimers utilizing a similar setting. Furthermore, the C-terminal tail, which is not resolved in the crystal structure of tubulin, was added to each tubulin heterodimer because it makes an important contribution to the overall tubulin electric charge [58] and affects the response of tubulin to external electric fields [3], [46]. The specific amino acid sequences of the C-terminal tails were GVDSVEGEGEEEGEEY for -tubulin and QDATADEQGEFEEEGEEDEA for -tubulin. The resulting protein sequence totaled 451 and 445 residues for and tubulin, respectively. We chose to simulate GTP-bound tubulin because it represents the protective cap of stable MTs.

Fig. 1

The constitution of the simulated microtubule (MT) system. a) The size and content of the simulated unit cell, which consists of one ring of 13 tubulin heterodimers bound to guanosine triphosphate (GTP) and a water box containing ions to neutralize the charge of the system and simulate a specific (close to physiological) ionic strength. b) The system is repeated through the application of periodic boundary conditions to form arbitrarily long MTs. c) Numbering of the tubulin heterodimers in a MT ring. Black radial lines represent relative binding energies along the lateral interfaces between neighboring tubulin heterodimers in the absence of any electric field. The longer the line, the stronger the binding between the neighboring dimers. -tubulin is colored in blue and -tubulin is colored in red.

Molecular dynamics details

We attached tubulin C-terminal tails as unstructured chains to the previously well-equilibrated MT ring and solvated the system in a water box. Ions were added to neutralize the system and simulate normal (1296 Na+ and 633 Cl−), high (1929/1266), and low (663/0) salt concentration environments (see Fig. S1); the total number of atoms in the three systems was 1235213, 1231535, and 1233947, respectively, which included 1929, 3195, and 663 ions, respectively, 1053312, 1048368, and 1053312 water atoms, respectively, and the same number of protein atoms (179972, including GTP cofactor) in each system. After energy minimization with the steepest descent algorithm followed by the conjugate gradient algorithm, we slowly heated the systems to a final temperature of 310 K with temperature and pressure coupling. This was followed by another 5 ns run under constant temperature to let the system relax into a suitable energy state. Before each production run there was a 500 ps pre-run without an applied E-field. This pre-run was performed under the same simulation conditions as the production run, in which the initial velocities were pseudorandomly assigned to all atoms from the Maxwell–Boltzmann distribution corresponding to a temperature of 310 K. We ran the main production run at constant volume with temperature scaled by a stochastic Langevin thermostat with a collision frequency of 2 . All bonds involving hydrogen were constrained using the SHAKE algorithm and a 2-fs time step was used for numerical integration of the Newton equation of motion. Periodic boundary conditions were imposed on all edges of our cubic system and the particle-mesh Ewald method was employed for treating long-range electrostatics. A non-bonded interaction cut-off of 10  was employed. The rest of the parameters were set to Amber defaults. All the simulations were run on the CUDA-accelerated version of AMBER software [59], [60] using the Amber ff14SB force field in combination with the TIP3P water model for solvation [61], and Joung and Cheatham parameters for ions [62]. The GTP cofactor was parameterized utilizing the parameter set developed by Meagher et al. [63]. The effect of an external electric field was implemented through force acting on each charged atom i in simulation box: , where is charge of that atom and is the vector of the electric field.

Distance analysis

For the analysis of distance we took the Euclidean distance of the centers of mass of individual tubulin heterodimers.

Angle analysis

We used a simple geometric analysis to characterize the rotation of the tubulin heterodimer. The resulting angle of rotation is the angle between two vectors: The xy projection of the dipole moment vector of the tubulin heterodimer and the xy projection of the vector connecting the center of mass of the tubulin heterodimer and the center of mass of the whole MT ring.

Interaction energy analysis

The energies of interaction between adjacent pairs of tubulin heterodimers along the lateral interfaces were estimated using an approach developed in a previous study [57]. The trajectory for each pair of adjacent tubulin heterodimers was separated from the original molecular dynamics trajectory, resulting in 13 different trajectories for the 13 different tubulin heterodimer pairs (tubulin 1–2, tubulin 2–3, …, tubulin 12–13, tubulin 13–1). The interaction energy between the pair was estimated using the molecular mechanics/generalized Born surface area (MM/GBSA) method, as implemented in AmberTools via the MMPBSA.py.MPI 14.0 script [64], employing default parameters and salt concentrations matching those concentrations in the original molecular dynamics simulation trajectories. In addition, another MMGBSA run was done on the trajectory of -X 10 MV/m, under the same conditions as the previous runs except that the decomposition of the contribution of each residue to the binding energy was estimated. Multiple scripts were written to automate the process and select only interfacial residues with significant contribution to the binding energy between tubulin dimer number 2 and number 3, the dimers that get detached due to the application of an electric field during the course of the simulations.

Results

Dependence on the electric field orientation

First, we investigated whether the orientation of the electric field had any effect on the MT structure. Initially, the tubulin dipole moments parallel and antiparallel to the electric field abruptly increased and decreased, respectively, in the sub-nanosecond time scale; this trend continued, although at a slower rate (Fig. 2 and SI (Section 7), resulting in a change to the polarization of the whole MT ring (see SI Section 59). The snapshots in Fig. 2 at 5–10 ns demonstrate our major finding: The electric field deformed the MT lattice and opened a gap between neighboring tubulin heterodimers, which consistently happened at the cathode side.

Fig. 2

Snapshots of the MT ring a) when no electric field is applied, the 100 MV/m electric field is applied in the b) X, c) -X, d) Y, and e) -Y direction. Black arrows represent the dipole moment vector of each tubulin heterodimer. The four arrows on the left side depict the electric field vector orientation. The cathode and anode side of the MT ring correspond to the end and beginning of the arrow, respectively. The snapshots b)–e) are from trajectory 1 for each condition.

The opening at the cathode can be intuitively understood, considering that the dipole moment of individual tubulin dimers point to the center of the MT cross-section (Fig. 2). Such an orientation renders the dipole of the tubulin closest to the cathode side (MT ring “pole”) to be antiparallel to the electric field vector, which is the least energetically favourable orientation. However, the torque on a dipole antiparallel to the electric field vector is zero, so why does the tubulin with an antiparallel dipole constitute an opening in the MT lattice? To obtain an integrative qualitative understanding of the MT lattice opening, we formulated a coarse-grained model, see SI (Section 2) and Fig. S2. Therein, each tubulin heterodimer is represented by a single grain and only Lennard–Jones-type forces and tubulin electrostatic dipole-moment-based forces are present in the model. The balance of these types of forces underlies the MT lattice stability [65]. An external electric field acts by exerting a torque on the tubulin dipoles, most strongly on the MT ring “equator”. This torque causes a strain of the whole MT ring, which accumulates at the MT ring poles (see the Fig. 2, 5—10 ns), ultimately leading to the opening of the MT lattice on the cathode side. See SI, Section 2, for more details.

However, which tubulin–tubulin interface opens up also depends on other factors. To determine how the heterogeneity of the energetics of the MT lattice might modulate the effect of the electric field, we calculated the binding energy at the lateral interface between each pair of neighboring tubulin heterodimers utilizing the MM/GBSA method on the molecular dynamics trajectory. Although the absolute binding energies obtained from the MM/GBSA approach are not very reliable, those binding energies can reliably be used for relative comparison within the same model [66]. We found that the tubulin–tubulin interface that preferentially opens up for each electric field orientation has a lower binding energy than the neighboring interfaces (see Fig. 1c and Fig. 2) in the unperturbed state. Over time, the gap in the MT lattice widened and the tubulin dimers at the gap tended to turn to align their dipole moments with the electric field (Fig. 2). This is understandable because the electric torque acts on the tubulin dipole. As soon as the tubulin breaks free from one of its lateral contacts, it starts to turn. Close inspection of Fig. 2 shows that the widest opening in the ring occurs in the case of the -X electric field orientation. This difference is more observable upon inspecting the distances between neighboring tubulin dimers, as can be seen when comparing Fig.  S5 to Fig.  S4, S6 and S7. This behavior is expected because the -X electric field orientation coincides with the interface between dimer 2 and 3, which shows the weakest binding of all 13 interfaces in Fig. 1C. This could also be the same reason why in Fig. S44, S45, S46, and S47, the -X electric field orientation shows the fastest opening of the ring at the interface between dimer 2 and 3, at which the opening is achieved in less than 2 ns. Based on that, we can infer that the opening of the ring is faster and more efficient when the field direction coincides with a weak linkage such as that between dimer 2 and 3.

Beyond the single gap formation between the two neighbor dimers, other consequent effects can be observed in some trajectories and field directions. For example, the evolution of the distance between dimers in individual trajectories for 100 MV/m, X field direction, trajectories 1 and 3 (Fig. S4a,c), indicates that there is not only a single opening occurring at the 8–9 interface, but simultaneously at the 9–10 interface, making the whole dimer (protofilament) leave the MT ring. Similar effect can also be seen for the 100 MV/m, Y field direction, trajectory 3 (Fig. S6c).

Effect of temperature and field strength

To quantitatively assess the kinetics of the opening process, we selected the -X electric field orientation and analyzed the time evolution of (i) the distance between the tubulins that are delimiting the gap in the MT lattice, (ii) the angle between the MT ring radius vector and the dipole moment vector of the tubulin at the gap, and (iii) the energy of the interface forming the gap. To understand the physical effects that might modulate the MT electro-opening, we performed an analysis of these three observables for two values of the electric field strength, temperature (Fig. 3), and ionic strength (SI Figs. S13, S14, S25, S26, S53, S54). See figures in SI (Sections 4,5,8) for all data on these three observables. The kinetics of the opening was field-strength-dependent: the width of the gap increased faster, the turning of the tubulin dipole was more substantial, and zero binding energy was reached much more rapidly with the 100 MV/m than with the 50 MV/m field strength. We did not analyze lower field strengths, but we expect that the electric field effect would become negligible when the interaction energy ( is the tubulin dipole moment and the electric field) is much smaller than the tubulin–tubulin binding energy . For instance, will be equal to  kcal/mol (see specific values in Fig. 3e,f and Fig. S43) for  MV/m and tubulin heterodimer of approximately 4000 debye (Fig. S31). However, beyond a quantitative difference, a qualitative difference between the two analyzed field strengths was evident. The slope of the gap widening and dipole turning kinetics was relatively linear and reproducible for 100 MV/m, whereas the kinetics of these observables for 50 MV/m displayed a biphasic behavior with a crossover point between 15 and 25 ns (Fig. 3a,c). The behavior of the observables at 100 MV/m was also more deterministic, whereas at 50 MV/m there was a much larger spread of values. This can be understood considering that the lower the interaction energy , the more significant the contribution of the stochastic thermal motion.

Fig. 3

Time evolution of the tubulin–tubulin distance a)–b) and tubulin heterodimer axial angle c)–d) and binding energy between tubulin heterodimer 2 and 3 e)–f). All data for all of the -X field direction. a), c), e) show the effect of the electric field strength for 100 MV/m and 50 MV/m, at the standard (310 K) temperature and standard ion concentration (1296 Na+ atoms (0.17 M) and 633 Cl− atoms (0.08 M). b), d), f) show the effect of temperature at 100 MV/m and standard ion concentration. Mean values (thick line) and the standard deviation (shaded envelope) are from N = 3 trajectories. See Figs. S3, S15,S43 for data at no-field condition.

This result led us to analyze the effect of temperature on the MT electro-opening process. We compared the kinetics of our observables at two different temperature values: room (297 K) and body (310 K) temperature (Fig. 3b,d,f). The temperature effect on the kinetics of the gap widening was negligible, whereas the turning of the tubulin dipole was significantly promoted at the higher temperature. This might be because polymers, including biopolymers such as MTs, become effectively softer at higher temperatures, enabling a higher flexibility for the rotation of weakly bound tubulin at the gap. The increased temperature also seems to affect the energetics by decreasing the binding energy between tubulin dimers.

Molecular details of the electro-opening

Analysis of the per-residue binding energy contribution at the tubulin–tubulin lateral interface of interest, which is the one between tubulin number 2 and 3, was also performed to give an insight into the atomistic events that take place during the electric-field-driven separation (Fig. 4). Snapshots from the first 2 ns of the simulation show that hydrogen bonds and other interactions quickly fade away under the effect of the electric field (Fig.4a–h). Beyond the first 2 ns, complete separation has already happened and no significant interactions between the tubulin dimers are apparent. The decline of the per-residue binding energy contribution over the simulation time for the residues with the most significant interfacial interactions is plotted in Fig.4i–j. It is clear that within 1–3 ns, all significant interactions have completely vanished. Interestingly, the residues that contributed the most to the binding energy were in perfect agreement with the residues that we identified in our previous studies [67], [57]. The integral scheme describing the mechanism of the electro-opening of the MT wall is in Fig.5.

Fig. 4

Snapshots of the opening detail (all from the first trajectory of -X 100 MV/m data): a)–d) between -tubulin 2 and 3, e)–h) between -tubulin 2 and 3. i)–j) Time evolution of the binding energy contribution of six highly contributing selected amino acids. Mean values (thick line) and the standard deviation (shaded envelope) are from N = 3 trajectories.

Fig. 5

The mechanism of electro-opening of the MT lattice. a) The electric field (grey arrows in the background) acts by torque (blue rounded arrows show direction of turning) on tubulin dipoles (red arrows). The magnitude of the force depends on the angle of the dipole vector with respect to the electric field vector. b) The torque on individual dipoles causes a deformation of the ring. c) The deformation accumulates at the poles of the ring. The orientation of the tubulin dipole parallel (on the right) to the electric field vector is stabilized (being at the minimum of the energy potential), whereas the orientation of the dipole anti-parallel (on the left) to the electric field is unstable (the dipole is at the peak of the energy profile). d) The opening occurs near the unstable tubulin at the interface with the lower binding energy, see Fig.1c. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

Discussion and conclusion

The deformation and opening of the MT lattice in this study are supported by the experimental observations of Chafai et al. [3], in which electrically-engineered tubulins self-assemble to a deformed and open MT structure. Based on that study, the nanosecond pulses transiently modified the effective tubulin charge, thus different tubulin subunits might have heterogeneous dipoles. The results of our current work confirm that modification of the electric dipole of tubulin is the key feature that leads to different deformation of the MT structure. We hypothesize that these findings can have future applications in controlling MT-based systems in artificial (nanotechnological) settings as well as in natural settings in patient’s tissue for potential novel treatment protocols. First, the opening might initiate MT “breathing” (fast opening and closing of the MT lattice), which is believed to facilitate transport of drugs [68] and enzymes [69] that bind to and modify the intralumenal side of the MT. Second, extensive opening of the MT lattice might promote the destabilization of the MT and consequent de-polymerization. Third, opening might induce the formation of defects in the MT lattice, hence affecting MT functionality [7].

In the current work, we aimed to investigate the effect of an electric field on the GTP cap of MTs. We selected the GTP-model over the GDP-model for two main reasons: First, it is believed that the GTP cap protects MTs from depolymerization [70]. In addition, it is more strongly bound together than regular GDP tubulin [57]. Therefore, we were interested in investigating the effect of the electric field on a GTP MT because if the electric field is capable of breaking apart the protective GTP cap of a MT, it should be more capable of breaking apart the rest of the GDP-bound subunits along the length of a MT. Second, the size of the GTP-cap in vitro ranges from 1 layer to approximately 15 layers, whereas in vivo it has been estimated to be over 50 layers (ca. 650 tubulin dimers) and extending over 1 μm [71], [18]. Therefore, modeling a complete GTP MT is a good approximation of the GTP-cap, especially when non-bonding cutoff in MD simulations would prevent any subunit from sensing subunits that are a few layers away. Similarly, we followed the same approach and modeled a GTP MT in previous research [57]. Moreover, some experimental studies and X-ray crystallography have been routinely performed on MTs bound to GMPCPP, which is a non–hydrolyzable GTP analogue [72], [73]. It would be intriguing to hypothesize what effect the electric field would have on a GDP-lattice compared with the GTP-lattice analyzed in this paper. We expect that the electric-field-mediated opening of the GDP lattice would occur faster or at lower field strength values than those for the GTP-lattice because some lateral contacts in the GDP lattice (e.g. 8–9, 1–13) have an even lower binding energy than the weakest lateral interface in the GTP-lattice [57]. It is well-known that a pulsed electric field can be used to open a cell membrane, in a process termed electroporation or electropermeabilization [74]. In this paper, for the first time, we show that a hollow protein structure with a nanoscale cross-section can be opened by applying an electric field. Although there seems to be a similarity in the phenomena, the underlying mechanisms are distinct: Cell electroporation is believed to occur owing to the formation of aqueous pores in the lipid membrane because of a large induced transmembrane voltage [75]. However, in our case, the major mechanism is the torque action on the tubulin dipole, which is a few-fold higher than in most other proteins [46]. Furthermore, in contrast to the intense microsecond and millisecond electric pulses that have been explored in bio-electrochemistry for several decades [76], the nanosecond electric pulses represent a much more recent area of research [77] and offer several distinct new features [17]. The feature with the strongest impact on the chemistry is that in the nanosecond pulse regime an electric field strength of a few tens of MV/m (approaching the strength of local molecular fields [78], [79]) can be attained throughout the volume exposed to the electric field without appreciable heating [17], [3]. Therefore, the reactivity and molecular structure can be modulated. We foresee new avenues of future research including investigation of the reversibility of the MT electro-opening after the field is switched off, the effect of the Z-electric field direction on the MT structure, the effect of C-terminus charge (which is modulated by various posttranslational modifications) on the MT response to electric field, and the effect of the nanosecond electric field on the MT tip dynamics (larger systems and time scales — coarse grained and larger computational power will be required). Our results open new potential applications in bio-nanotechnology and biomedical approaches that rely on employing or regulating noncovalent interactions in tubulin or other protein and supramolecular assemblies.

CRediT authorship contribution statement

Jiří Průša: Data curation, Formal analysis, Investigation (lead), Methodology, Software, Visualization, Writing - original draft (supporting), Writing - review & editing (supporting). Ahmed Taha Ayoub: Formal analysis, Investigation (supporting), Methodology, Visualization (supporting), Writing - original draft (supporting), Writing - review & editing (supporting). Djamel Eddine Chafai: Conceptualization, Validation, Writing - review & editing (supporting). Daniel Havelka: Conceptualization (supporting), Visualization (supporting), Resources, Writing - review & editing (supporting). Michal Cifra: Conceptualization, Formal analysis (supporting), Funding acquisition, Investigation (supporting), Project administration, Resources, Supervision, Validation, Visualization (supporting), Writing - original draft (lead), Writing - review & editing (lead).

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.