4.3. In Silico Analyses

The models of full-length of toxin ζ models bound to UNAG and ATP were built using I-TASSER with ab initio structure prediction [28]. Pockets in the protein surface were identified using 3 V [29]. Initial ligand–protein conformations were assigned by transference from related structures using Dali [33] and UCSF Chimera [34]. The complexes were initially minimized using OpenBabel [35] with the UFF, MMFF94S and GAFF force fields (4000 minimization steps), and further optimized using GROMACS with the AMBER and OPLS-AA force fields [64,65,66]. In silico alanine scanning of toxin ζ and its variants were generated with TRITON [43], and Modeller v9 [44] using a very large MD refinement.

The structures of toxin ζ substrates (ATP, ADP, UNAG and UNAG-P) were optimized with a SQM conformational search using Gabedit [36] and OpenMopac [37] with the PM7 parameter set and an implicit solvent model.

Additional protein–ligand conformations were generated using molecular docking with DOCK6 and AutoDock VINA [38], and the complexes were optimized as described. Ligand-protein affinities of selected conformers were estimated using Xscore [40] and DrugScoreDSX [41].

For MD simulations, the target temperature and pressure were 310 °K and 1 bar. Structures were subjected to an initial optimization step (up to 10,000 steps), followed by NVT ensemble (200 ps, in 10⁵ 2fs steps) using the V-rescale thermostat and an in NPT equilibrium for 200 ps (in 10⁵ 2fs steps) with the V-rescale thermostat and the isotropic Parrinello-Raman barostat. The MD simulations were also run (4 ns, in 25·10⁵ 2fs steps) using the Nosé-Hoover thermostat and the isotropic Parrinello-Raman barostat using GROMACS [39] with the AMBER and OPLS-AA force fields, and the ligand parameters were generated with ACPYPE [66] using AM1-BCC charges [64].

Detailed reaction simulations were performed using GTKDynamo [50] and pDynamo [51]. The most promising models were first refined using QM with OpenMopac following a protocol consisting of H assignment, H optimization and self-referenced optimization of the protein structure, and then subjected to a QM/MM simulation. The complete systems consisted of the ligands alone and in solution and of the protein and ligands surrounded by 10 Å of water or NaCl solution, partitioned into a QM region (complete ligands, nearby residues and relevant solvent molecules), a MM region (defined as the whole protein and all solvent molecules within 4 Å of the protein, and a fixed region (defined as solvent molecules beyond the MM limit). The reaction was modeled in the QM region approaching the reactants in 0.05 Å steps, using SQM with the PM6 parameter set and considering the appropriate charges for the system. The MM region used the AMBER force field, with parameters for the ligands being generated with GTKDynamo [50] and Ambertools [65]. Results were visualized using TRITON, UCSF Chimera, PyMol (www.pymol.org) [67], Jmol [68] and POVray (www.povray.org).