The starting models for simulations were built from the respective XFEL-based time-resolved structures. We removed the DARPin protein and reconstructed the βM-loop (res 272-288) using the one in chain D of PDB 4I4T as a template by employing the Prime65 (link) module of the Schrödinger 2020-4 suite. Analogously, Mg2+ ion in chain B (β-tubulin) was added by superimposition with the one of chain D of 4I4T structure. Calcium ions were removed from the model.
The final model consisted of an αβ-tubulin heterodimer in complex with azo-CA4, GTP, and GDP molecules bound to the α- and β-tubulin subunits, respectively, as well as their associated Mg2+ ions and their coordinating water molecules. The resulting protein structure possessed 437 out of 451 residues of α-tubulin (UniProtKB ID P81947) and 431 out of 445 residues of β-tubulin (UniProtKB ID Q6B856). Missing residues belonging to the intrinsically disordered C-terminal tails of α- and β-tubulin were not modeled and C-termini were capped with N-methyl amide (NME) groups. Residue protonation states were evaluated at pH 7.0 using the Protein Preparation Wizard tool66 (link) implemented in the Schrödinger 2020-4 suite67 . The αβ-tubulin heterodimer structure was solvated with the TIP3P-model68 (link) for water molecules in a truncated octahedron box using 12 Å as the minimum distance between the protein and the box edges. The system was neutralized by adding Na+ ions resulting in a total of about 122k atoms. The atomistic force field Amber-ff14SB69 (link) was used for all simulations. Parameters for Mg2+ ions and the GTP and GDP molecules were developed by Allner et al.70 (link) and Meagher et al.71 (link), respectively.
Concerning the ligand, partial atomic charges were obtained with the Restrained Electrostatic Potential (RESP) method on a DFT-optimized trans conformation at the B3LYP/6-31 g* level. The geometric optimization was carried out with Terachem v1.94. The General Amber Force Field (GAFF) was used, with custom parameters for the three main dihedral angles: 1-1a-1b-1’, 6-1-1a-1b, 1a-1b-1’-6’.
To obtain such optimized parameters, we fitted the GAFF energy curves associated with a rigid torsion of the dihedral angles to the corresponding DFT energy curves. For each dihedral, we proceeded as follows: we started out with the geometrically optimized trans structure. We then performed a clockwise rigid rotation with 15° steps followed by an analogous counter-clockwise rigid rotation. For each rotation step we performed a DFT geometric optimization restraining the rotated dihedral, and we evaluated the energy. We retained the minimum energy among the clockwise and anticlockwise rotation cycles to overcome hysteresis. The three resulting DFT energy profiles consisting of 24 evaluation points were used as a target for the fitting procedure. DFT calculations were performed with Terachem v1.94, at the B3LYP/6-31 g* level. The fitting procedure was carried out separately for the three energy profiles with mdgx, part of AmberTools21. The αβ-tubulin heterodimer system was assembled with the LEaP tool implemented in the AmberTools21 software package67 .
The designed models for QM calculations consisted of protein residues α(99-101, 178-181), β(237-242, 248-259, 314-321, 349-354, 378-380), with the respective N- and C- termini capped with acetyl (ACE) and NME groups, the azo-CA4 molecule, and its two closest water molecules. They resulted in 743 atoms.
Free full text: Click here