Molecular Dynamics Simulations of Tofacitinib-JAK Complexes

Each tofacitinib/JAK complex was prepared by all-atom MD simulations
with three different initial velocities using the pmemd CUDA in AMBER16^{45 (link)} in periodic boundary conditions with the isobaric-isothermal
(NPT) ensemble. Details of simulations were as described for other
biological systems.^{46 (link)−48 (link)} The FF14SB^{49 (link)} and GAFF^{50 (link)} force fields were applied for protein and tofacitinib,
respectively. To generate the partial charges of inhibitor, the 3D
structure of tofacitinib was optimized by the HF/6-31(d) level of
theory as per previous studies^{51 (link)−53 (link)} using the Gaussian09 program.⁵⁴ The electrostatic potential (ESP) charges and
restrained ESP (RESP) charges of the ligand were generated using parmchk
of AMBER16. The systems were soaked in the boxes of explicit water
using the TIP3P model (∼10 750 molecules for JAK1, ∼10 696
molecules for JAK2, and ∼10 439 molecules for JAK3).
The time step was set as 2 fs at a constant pressure of 1 atm.^{55 (link)} The short-range cutoff of 12 Å was used
for nonbonded interactions, while long-range electrostatic interactions
were treated by Ewald’s method.^{56 (link)} Temperature and pressure were controlled by the Berendsen algorithm.^{57 (link)} The SHAKE algorithm was used to constrain all
covalent bonds involving hydrogen atoms.^{58 (link)} The simulated models were heated up to 310 K with the relaxation
time for 100 ps. The temperature was controlled by a Langevin thermostat
with a collision frequency of 2.0 ps. Finally, the unrestrained NPT
simulation was performed for 500 ns. The MD trajectories were recorded
every 500 steps for analysis. The RMSD analysis of each system was
performed using all atoms. The intermolecular HB occupation, SASA,
the number of contacts, and the motion of proteins were evaluated
using the CPPTRAJ module.^{59 (link)} Besides, the
MM-PB(GB)SA and QM/MM-GBSA ΔG_bind and ΔG_bind,residue were calculated
by the MM/PBSA.py module.^{60 (link)} Note that in
the QM/MM approach, tofacitinib was quantum-mechanically treated by
the semiempirical method PM3 and the self-consistent-charge density-functional
tight-binding method (SCC-DFTB),^{61 (link)} whereas
the remaining region was described by molecular mechanics using the
FF14SB force field. The same sets of MD snapshots were used to predict
the ΔG_bind based on the solvated
interaction energy (SIE) method.^{62 (link)} SIE
is an end-point physics-based scoring function for predicting binding
affinities in aqueous solution, which is calculated by an interaction
energy contribution, desolvation free energy contribution, electrostatic
component, and nonpolar component.^{62 (link)}