Tubulin as a sensitive target of nanosecond-scale intense electric field: quantitative insights from molecular dynamics simulations

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 intra-cellular components such as cytoskeletal proteins could be directly influenced by electric pulses within biomedically-attainable parameters. 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 β-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 novel therapeutic protocols. Author summary α/β-tubulin heterodimers are the basic building blocks of microtubules, that form diverse cellular structures responsible for essential cell functions such as cell division and intracellular transport. The ability of tubulin protein to adopt distinct conformations contributes to control the architecture of microtubule networks, microtubule-associated proteins, and motor proteins; moreover, it regulates microtubule growth, shrinkage, and the transitions between these states. Previous recent molecular dynamics simulations demonstrated that the interaction of the tubulin protein macrodipole with external electric field modifies orientation and conformations of key loops involved in lateral contacts: as a result, the stability of microtubules can be modulated by such fields. In this study, we seek to exploit these findings by investigating the possibility of fine-tuning the dipolar properties of binding sites of major drugs, by means of the action of electric fields. This may open the way to control tubulin-drug interactions using electric fields, thus modulating and altering the biological functions relative to the molecular vectors of microtubule assembly or disassembly. The major finding of our study reveals that intense (> 20 MV/m) ultra-short (30 ns) electric fields induce changes in the major residues of selected binding sites in a field strength-dependent manner.

properties, were proposed to be involved in endogenous electrodynamic processes in 78 cells [43][44][45]. However, all-atom molecular simulations of external EF effect on tubulin 79 have been carried out only recently [46,47]. They specifically investigated the EF effects 80 on protein mechanics but did not include the C-terminal tail. The C-terminal tail is a 81 highly flexible unstructured domain of tubulin, which (i) is essential for tubulin-protein 82 interactions [48,49], (ii) is the main site of protein mutation and post-translational 83 modifications [49], and (iii) greatly contributes to the overall electric properties of 84 tubulin accounting for approximately 30-40% of the total charge [33]. A very recent 85 study [50] presented results of the molecular dynamics simulation analysis of the electric 86 field effect on tubulin. This latter study investigated electric field strengths in the range 87 between 50 and 750 MV/m, which overlaps with the values used in our simulations, and 88 carried out simulations over 10 ns while our study reports on simulations that ran up to 89 30 ns. The previous study only examined displacement effects on key secondary 90 structure motifs such as tubulin's C-termini and alpha helices. In contrast, in the 91 present study, we aim to unravel the mechanisms of interaction of intense EFs on the 92 tubulin protein family, given their important and attractive role as drug targets for 93 cancer therapy, due to their involvement as key-players in cellular self-organization. In 94 fact, it is already known that various stabilizing/destabilizing tubulin-binding drugs 95 such as taxanes, colchicines, and vinca alkaloids, bind to different sites on the tubulin 96 dimer, modulating microtubule-based processes [51]. This type of binding is assumed to 97 be based on purely electrostatic interactions like those employed in the recognition 98 between proteins [2]. Conversely, what is presently unexplored, is the possibility to 99 modulate such electrostatic environment by means of EF. Therefore, we first employ 100 bioinformatics tools to systematically compare electric charges and dipole moments of 101 the various isoforms of tubulin to the so-called PISCES set, which represents all unique 102 chains in the whole protein database. Then, we quantitatively evaluate the effects of 103 intense nanosecond EFs on the overall shape of the tubulin dimer as well as the 104 time-dependent evolution of its dipole moment. Finally, we closely analyze the dipole 105 moments of individual residues with special attention to those forming the binding 106 pockets for the most important tubulin-modulating pharmacological agents (e.g. 107 paclitaxel or vinca alkaloids) and for GTP/GDP molecules. Thus, in this paper, we cast 108 light on the mechanism of direct EF action on tubulin at the molecular level. In this study we first compared the structural charges and dipole moments of proteins 113 from the tubulin family to those of the PISCES set of proteins (all unique protein chains 114 available in the RSCB protein database). It can be seen from Fig. 1A that tubulin 115 proteins possess a generally much higher structural electric charge (with a mean value of 116 -22 e per monomer) than most PISCES proteins (mean value of -4 e). This is likely due 117 to an excess of acidic over basic amino acid residues present among the tubulin proteins 118 compared to the PISCES set of proteins, see S1-1TUB-aa content.xlsx. Taking as an 119 example the 1TUB tubulin structure, we found that while a relative content of acidic 120 residues is similar in tubulin versus the average from the PISCES set, tubulin contains 121 fewer basic residues than an average protein. It is worth noting that a large fraction (on 122 average around 34% across the tubulin dataset) of electric charge of tubulin is due to the 123 unstructured and highly flexible C-terminal tail (CTT). This is remarkable, since the 124 CTT of tubulins is rather short -having on average around 20 residues. A rather high 125 value of the electric charge of tubulin CTT is due to a large content of acidic residues. 126 Fig. 1B shows that tubulins possess a much higher dipole moment than PISCES 127 proteins: with a mean of (2,166 debye) versus a mean of (555 debye), respectively. The 128 high dipole moment arises from an asymmetric electric charge distribution, i.e. an 129 asymmetric distribution of acidic versus basic residues. This asymmetric charge 130 distribution is partially due to the fact that a large fraction of charge is located on the 131 C-terminal tail. However, even the tubulin structures where the CTT is not resolved 132 possess a rather high dipole moment, e.g. both 1TUB and 1JFF crystal structures 133 result in a dipole moment above 2,000 debye. The results in Fig. 1 clearly demonstrate 134 that the tubulin proteins have exceptional electric properties, hence it is reasonable to 135 analyze in detail if these properties can be exploited to manipulate tubulin's structure 136 and hence its function using an external EF with specifically designed characteristics. Rounding effect on tubulin shape induced by electric field 138 We have employed molecular dynamics (MD) simulations, see Methods for details, to 139 analyze effects of an EF acting on tubulin. The tubulin structure used in these 140 simulations consists of α-and β-monomers, which form a stable noncovalently bonded 141 heterodimer, see Methods for details.

142
In the MD simulation performed, we have added an EF as an external Coulomb 143 force acting on every atom in the system. We analyzed the effect of the EF strength 144 starting from 20 MV/m up to 300 MV/m. The rationale for the selection of this range 145 of field strength is both theoretical and experimental. Theoretically, the interaction 146 energy U = p · E, where p is the tubulin dipole moment vector and E the electric field 147 vector, exceeds thermal energy for field strengths in the range > MV/m. The range 148 < 70 MV/m is also experimentally attainable since it is below the field strength of 149 dielectric breakdown of water-like media with an exact value depending on the ionic 150 strength and the electric pulse duration [52]. The length of the simulation was selected 151 to be up to 30 ns. This time scale is short enough to be computationally tractable and 152 long enough to be comparable to experimental data. Furthermore, electric pulses of this 153 duration typically do not cause any appreciable heating in experiments conducted even 154 for the experimentally attainable field strengths mentioned above.

155
The most basic integral parameter of protein is its shape. Hence, we first analyzed 156 how the EF affects the shape of the tubulin dimer, see Fig. 2 (note that no water New coordinate system and dimensions of the tubulin approximated by an ellipsoid. The electric field affects not only the orientation of the tubulin dimer but also its overall shape. The distributions on the right are from 2,500 frames from the last 5 ns of the MD simulation (sampling rate 2 ps). process, bypassing the trivial roto-translational effects taking place due to the external 161 EF. At first, we observed that the effective shape of the equivalent tubulin ellipsoid is 162 affected by a 100 MV/m EF in a way that the medium axis is elongated by 11% and 163 the long axis is shortened by 2%. We also observed that tubulin undergoes rotation -164 this will be discussed in more detail in the next section as the rotation is related to the 165 dipole properties of the tubulin analyzed there. Orientation effect of electric field on tubulin dipole moment 167 As the next step in this investigation, we focused on the effect on the electric dipole 168 moment of the tubulin dimer, see Fig. 3. It is readily seen that the tubulin dipole 169 moment under zero-field conditions has an average value of around 2,500 debye. 170 However, in the presence of a 100 MV/m EF, the dipole moment is increased to more 171 than 6,500 debye, i.e. it more than doubles in the process. In the lower part of the 172 figure it is also visible the polarized C-terminal tails extending from the structure. To 173 get deeper insight into time evolution and field strength dependence, we further 174 analyzed the tubulin y dipole component, the most affected one (i.e. with the highest 175 polarization change) among the three dipole components shown in Fig. 3. The external 176 electric field vector was oriented in the z direction (Cartesian reference system, see  Fig. 4). This stable dipole value evolution is due to the coordinate system 182 transformation adopted (i.e. from the external Cartesian system to the internal protein 183 system), which removes the dipole roto-translational effects, highlighting internal 184 polarization effects. In fact, the EF tends to act by exerting a torque on the tubulin's The fundamental and functionally-crucial characteristics of any protein structure are the 197 type and location of its secondary structure motifs. Therefore, we have focused next on 198 the analysis of how an EF affects the count of tubulin residues being part of a certain 199 secondary structure motif. The results of this analysis can be seen in Fig. 5, which 200 displays the count of residues in alpha helices, beta-sheets, coils and turns as averaged 201 over the last 5 ns of the corresponding trajectories. The field strength of 100 MV/m 202 seems not to cause any effect on the structure -just slightly increasing the content of 203 the tubulin residues present in coils at the expense of beta-sheets, alpha helices and 204 turns. Only 300 MV/m induces a strong effect on tubulin secondary structures, where 205 more than 70 residues are converted from alpha-helices to coils. This is a sign of 206 unfolding, when charged CTTs are experiencing an electrostatic force, which is strong 207 enough to undergo pulling out from the tubulin body. Since there are two alpha helices 208 (H12 and H11) close to CTTs, those are the first ones which are unfolded. Furthermore, 209 we selected the CTT of β-tubulin since it is longer and carries more charge than the 210 α-tubulin CTT for a more detailed view. In Fig. 6A, we quantified the count of CTT 211 residues being in coil structures within the 30 ns time scale for 0, 20, 50, and 212 100 MV/m field strength conditions. We can see that the 20 MV/m condition does not 213 tend to change the secondary structure of CTT compared to the no field condition. being in a coil structure at the expense of α-helices and turns. In these cases, the higher 216 the field strength, the shorter the time needed to achieve the maximum amount of 217 residues in a coil structure.  Apart from the EF effect on CTTs, we also asked a question whether an intense field 226 can affect local electrostatic conditions of the tubulin dimer. To this end, we analyzed 227 the shift of the dipole moment of every residue of the tubulin dimer. We plot the y 228 component (internal tubulin coordinates) of the dipole moment in Fig. 7. It can be seen 229 that some residues undergo substantial change of the dipole moment when exposed to 230 an EF of 100 MV/m. In particular, when considering both monomers it is possible to 231 appreciate dipolar shifts ranging from -10 up to almost +15 debye. Specifically, 22 out 232 of the 451 residues of α-monomer residues exceed a ±5 debye shift, while for the The Methods section provides a rationale and identification of residues belonging to the 239 tubulin sites based on energy considerations. In Fig. 8, we highlight the location of the 240 important tubulin sites in a color-coded manner and also provide a list of the selected 241 residues. We display the histograms of the dipole moments (y component) of six 242 selected residues where we found the strongest effects. In all of these six cases, we see a 243 field strength-dependent effect, shifting the y component of a residue dipole moment 244 towards different values. Since the local electrostatic field is crucial in active sites of 245 many enzymes enabling the process of catalysis, we speculate that influencing the local 246 field could affect the GTP hydrolysis rate. We see that α-Asp 251, a crucial residue 247 mediating GTP hydrolysis [54], has its y component of the dipole moment (yDM in 248 short) affected by almost 3 debye when comparing no field and 100 MV/m conditions. 249 Tubulin-tubulin longitudinal interactions are mediated by several residues. β-Arg 391 is 250 one of them and has its yDM also influenced by the EF (Fig. 8). However, in this case 251 the field strength has opposite sign effect. 20 MV/m field tends to affect yDM in an 252 opposite direction than fields with 50 and 100 MV/m field strengths, compared to no 253 field condition. This is probably due to the position of the residue in the β-monomer 254 (see Fig. 8, left-upper panel); in fact Arg-391 is located in the external part of the 255 protein, presumably in contact with a layer of water, therefore its response to the 256 external EF could be screened by water for lower field strength. The second most 257 important residue for binding paclitaxel energy-wise, β-Val 23 (Fig. 8), also manifests 258 the changed yDM up to few debye when an EF of 100 MV/m is applied. Furthermore, 259 residue β-Cys 239 belonging to the colchicine binding site has its yDM shifted towards 260 zero with an increasing field strength. We also found that α-Pro 325 and β-Asp 177  (Fig. 1). This result corroborates our further findings: tubulin seems to be more 282 susceptible to EFs compared to other proteins analyzed earlier [11,12,19], in a way that 283 lower field strength, i.e. 20 MV/m, is able to attain a 50% increase of tubulin dipole.

284
In a recent paper [55] a detailed energy balance calculation was provided for the 285 stability of a microtubule with a seam (representing a so-called type B microtubule 286 lattice). The calculations included the contributions from dipole-dipole interactions 287 between tubulin dimers, solvent accessible surface area, van der Waals and electrostatics 288 (generalized Born approximation) to demonstrate that the balance energy is such that a 289 single GTP hydrolysis event can trigger a microtubule disassembly because when the 290 Fig 8. Electric field effect on selected important tubulin residues; in the upper panel a graphical representation of known αβ-tubulin sites is given together with a look-up table -see text criteria and for references where the residue numbers were obtained. In the lower panel, the histogram representation of electric field effects (ranging from zero up to 100 MV/m) on a dipole moment of selected relevant residues is presented calculated from the frames of the last 5 ns of simulation. seam is closed with GTP molecules attached to the β monomers, the net free energy is 291 −9 kcal/mol. The dipole-dipole energy is positive (destabilizing) and amounts to 292 27 kcal/mol. When the seam becomes open due to GTP hydrolysis, the net free energy 293 becomes vastly positive. These calculations used a dipole moment of approximately 294 4,500 debye for a tubulin dimer as an average value over various isotypes. In the present 295 paper we have shown that strong electric fields can substantially increase the dipole 296 moment of tubulin. Here, we have shown the dipole moment to be close to 3,000 debye 297 in zero field, which would translate into a dipole-dipole energy for an microtubule 298 lattice of 12 kcal/mol and a net free energy of -24 kcal/mol, hence slightly more stable 299 than the calculation provided in [55]. However, at strong electric field values 300 (≥ 50 MV/m), the corresponding dipole moments were found to increase to as much as 301 6,000 debye, which translates into a dipole-dipole energy of 48 kcal/mol and a net free 302 energy of +12 kcal/mol making the lattice unstable. This quantitatively supports our 303 hypothesis that sufficiently strong electric fields disrupt a microtubule lattice by  Structural/conformation changes on tubulin induced by electric 307 fields 308 By approximating the tubulin with the ellipsoid as given by the covariance Matrix 309 method, we were able to appreciate the actual shape changes under the action of the 310 EFs. The protein slightly reduces its major axis, while clearly increasing its medium one. 311 This effect seems of particular interest since it is a consequence of protein polarization. 312 In fact, the tubulin dimer undergoes a packing transition, rather than an elongation, 313 since its dipole moment aligns along tubulin medium axis (see Fig. 3). This non-trivial 314 effect is due to the high charge of tubulin CTTs. Tubulin seems to have a lower 315 threshold for the unfolding transition [56]. For example, unfolding effects tend to appear 316 at lower field strengths or within a shorter time frame (unfolding at 250 MV/m started 317 at approximately 160 ns in insulin [11, Fig.1c] compared to only a few ns at 200 MV/m 318 for tubulin). This is reasonable since the higher the charge and dipole moment of the 319 structure, the higher the electric force and torque, respectively, acting on the protein. general is the effect on the CTT we observe, or put another way, is the CTT sequence 327 we used common in other biological species?

328
An exact tubulin sequence and structure varies across biological species and 329 tissues [33]. While the tubulin core is rather conserved across species, CTTs display 330 higher variability across species and are also a target of post-translational modifications 331 forming a so-called "tubulin code" [57,58]. The β-tubulin CTT sequence we used in our 332 simulations has a 100% sequence identity to certain tubulin isotypes found for example 333 in pig Sus scrofa (common experimental source of tubulin) tubulin gene TBB and to 334 many tubulin isotypes of more than 20 species, see S3-CTT alignment.fasta. It has also 335 a 90% sequence similarity to CTT of TBB2B gene in the human, the only difference is 336 that two immediate neighbor residues Gly 440 and Glu 441 are swapped, so the total 337 charge of the CTT is the same as that of the CTT in our tubulin structure, see S3-CTT 338 alignment.fasta.

339
As a general result, we observed a widespread EF effect on protein residues. In embedded in a microtubule, is the cytosol-facing outer surface while the surface facing 352 the lumen is less negatively charged. Also, the dipole moment of the dimer also includes 353 a contribution from the CTT which is highly variable due to the flexibility of this motif. 354 Finally, the dielectric breakdown, which would pose a limitation on the strength of the 355 applied field in practice, depends on the ion strength of the solution and the duration of 356 the applied pulse [52].

358
In discussing the implications of our modeling effort for experiments on microtubules, it 359 could be argued that field strengths > 100 MV/m are not readily experimentally 360 attainable yet. However, one has to keep in mind that MD simulations are commonly 361 best suitable for providing relative assessment of various effects and not necessarily 362 absolute values directly translatable the experimental situation. For example, the 363 binding free energies of ligands interacting with proteins are typically off by a factor as 364 large as two or more [59]. Hence, effects such as unfolding, which are computationally 365 predicted to occur at large fields, may in fact require field strengths lower than predicted 366 for tubulin [56]. Moreover, rapid technological advances in the field may well lead to the 367 development of engineering solutions with sufficiently high electric field strengths 368 attainable in the clinical setting. The ultimate physical limit, however, is the field 369 strength at which the dielectric breakdown of the exposed biological material occurs.

370
Challenges for future work 371 In our future work we intend to better calibrate the simulation parameters in order to 372 make the model quantitatively predictable so that such characteristics as field strength, 373 frequency and duration may be tailor-designed to elicit specific response of the target, 374 namely the tubulin dimer or an entire microtubule. This modeling effort may need to 375 be extended to include quantum effects through the use of quantum 376 mechanics/molecular mechanics calculations especially to be able to elucidate covalent 377 bond response to electric fields, for example at the GTP binding site.

378
To conclude, our results allow to draw biologically-relevant consequences in terms of 379 microtubule stability and the changes in the strength of binding for tubulin-binding 380 drugs under the influence of EF. Therefore, we believe the present paper provides new 381 and important quantitative insights into the effects of EF pulses of nanosecond duration 382 on tubulin and microtubules. The knowledge of residue-specific EF realignments allows 383 to extrapolate the computer simulation results in terms of their consequences on the 384 behavior of other tubulin isoforms and tubulin mutants under both exogenous and 385 intrinsic molecular electric field.

387
We follow a gapless numbering convention of tubulin residues in the current paper, in 388 contrast to the Protein Data Bank entries 1JFF and 1TUB, which contain 2 gaps after 389 β-Leu 44 and 8 gaps after β-Pro 360. Charge and dipole moment of proteins from tubulin family was obtained from our 393 earlier work [33], which includes 246 isotypes of tubulin, mostly α-and β-tubulin, from 394 various species. PISCES set, a representative set of all proteins in the RSCB PDB 395 database, was obtained using PISCES server [60]. We used following settings to run the 396 server query (date 2018-05-18): sequence percentage identity threshold ≤ 90, Resolution: 397 0.0 -4.0, R-factor: 0.5, sequence length: 40 -10,000, Non X-ray entries: Included, 398 CA-only entries: Excluded, Cull PDB by entry, Cull chains within entries: No. We 399 obtained a list with 36,331 entries, see S4-PDB entries.zip. Charges and dipole moment 400 data for our PISCES set of proteins were obtained from Protein Dipole Moments Server 401 of Weizmann Institute https://dipole.weizmann.ac.il/index.html [61]. Each 402 protein chain was evaluated separately which leads to a total number of 63,110 records 403 for both protein chain charge and dipole moment, see S5-PISCES-Charge,Dipole.zip.

405
Note that the methods for charge and dipole moment calculation used by Tuszynski 406 et al. [33] (GROMOS96 43a1 force field) and the dipole moment server [61]  The structure of the GMPCPP-bound tubulin was obtained from the Protein Data Bank 413 under the PDB ID: 3J6E [62]. The cofactor GMPCPP, a non-hydrolyzable analogue of 414 GTP, was modified into GTP. Furthermore, the CTT, which are usually missing in PDB 415 structures of tubulin, were added to both α-and β-tubulin subunits. This was achieved 416 using the Molecular Operating Environment software (MOE, Chemical Computing 417 Group Inc.) [63]. MOE was also used for the addition of hydrogen atoms and the 418 prediction of ionization states. The CTTs were added in an extended conformation and 419 we depended on molecular dynamics simulations later to help to restructure the tails in 420 the correct conformation in solution. The CTT is defined here as the last 16 residues of 421 the α-tubulin (GVDSVEGEGEEEGEEY) and the last 20 residues of the β-tubulin 422 (QDATADEQGEFEEEGEEDEA). The CTT sequence identity was tested using BLAST 423 tool on [64], see the .fasta format result in S3-CTT alignment.fasta.

424
Residues of drug binding sites 425 The paclitaxel, colchicine/nocodazole and vinca binding sites were identified based on 426 the RSCB PDB structures 1JFF (paclitaxel is named as residue TA1, site identifier 427 AC5) [65], 1SA0 (colchicine is named as residue CN2, site identifier AC8) [66] and 428 5BMV (vinblastine is named as residue VLB, site identifier AE1) [67], respectively. The 429 residues belonging to the respective binding sites are identified there using an algorithm 430 from [68] and are available in a respective .pdb file. For paclitaxel and colchicine, we 431 additionally included the residues which contribute 75% binding energy between 432 paclitaxel and tubulin based on the analysis in [69] and [70], respectively. Those results 433 provided β-tubulin additional residues LEU 215, GLN 280, and LEU 361 on top of the 434 list from RSCB PDB 1JFF structure for paclitaxel and α-ALA 180 and β-LEU 246 on 435 top of the list from RSCB PDB 1SA0 structure for colchicine. See the complete list of 436 residues of individual interaction and binding sites in Fig. 8.

437
Molecular dynamics simulation 438 We carried out MD simulations of a tubulin protein in water solution using the 439 GROMACS package v. 4.6.5 [71] Table 1. List of molecular dynamics simulations with the 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 ).
hydration shells and a significant amount of bulk water. Following an energy 445 minimization and subsequent solvent relaxation, the system was gradually heated from 446 50 K to 300 K using short (typically 50 ps) MD simulations. A first trajectory was 447 propagated up to 150 ns in the NVT ensemble using an integration step of 2 fs and 448 removing the tubulin center of mass translation but with no constraints on its related 449 rotation. The temperature was kept constant at 300 K by the Berendsen thermostat [73] 450 with the relaxation time equal to the simulation time step, hence virtually equivalent to 451 the isothermal coupling [74] which provides consistent statistical mechanical behavior. 452 All bond lengths were constrained using LINCS algorithm [75]. Long range  Table 1 ), acting in the simulation box as 458 explained in [78] for 30 ns in each exposure condition and applied along the z-axis of the 459 Cartesian reference of frame. More precisely the application of the electric field takes 460 place in continuity at the last frame of the unexposed simulation, thus allowing a direct 461 evaluation of the characteristic response over time of the system due to the switch on of 462 the exogenous perturbation. Taking advantages of the previous findings in [19,56,79] on 463 the predictable protein polarization process, for the highest field strengths (200 MV/m 464 and 300 MV/m) the simulation box was enlarged along the z-direction up to 30 nm, 465 resulting in a rectangular box of 12×12×30 nm 3 .

466
Covariance matrix method 467 A simple method to characterize the geometry of a complex system like a protein is to 468 approximate its overall geometrical structure, at each MD frame, to the ellipsoid given 469 by the eigenvectors of the Covariance Matrix as obtained by the distribution of the x, y, 470 and z coordinates of the protein atoms (the geometrical matrixC) [19,80]. Therefore, 471 for a system defined by N-atoms the elements of the 3×3 geometrical matrix at each