To address the aforementioned issue, this study presents a new application of molecular simulation for pKa prediction using the polarizable AMOEBA force field to improve over traditional empirical models like PROPKA. In particular, all proteins containing residues with experimentally determined pKa values were downloaded from the Protein Data Bank (PDB). These files were then cleaned via a multi-step preparation protocol that include: first, PDBFixer (via OpenMM7) was used to add missing atoms and residues and convert non-standard residues to canonical amino acids. Next, ChimeraX was used to strip extraneous water molecules while preserving essential metal ions (e.g., Cu, Fe). For problematic structures (e.g., 1I0E), SWISS–PDBViewer was employed to correct anomalies. The cleaned structures were then converted into .xyz files compatible with the Tinker molecular modeling package using the pdbxyz tool that utilized the AMOEBABio18 polarizable force field. Proteins were solvated in pre-equilibrated or custom-generated cubic water boxes with explicit water molecules. Systems were neutralized to achieve net charge balance, and energy minimization was performed using Tinker. A sphere centered at the residue's Cα was cropped from the whole protein to isolate the residue from the whole protein, 5 different radii from 7A to 11A were tested to optimize the representation of the residue's local environment and to better capture the unbounded Repulsion-dispersion (vdW) effects based on the Halgren’s buffered 14−7 function employed by AMOEBA force field.
The protein features extracted from the aforementioned molecular modeling workflow were used as inputs for three advanced graph-based deep learning models: Graph Convolutional Network (GCN), Graph Isomorphism Network (GIN), and Graph Attention Network (GAT). These models were selected, constructed, and tested for pKa prediction due to the inherently graph-structured nature of molecular data. Each protein residue was represented as a graph, with atoms as nodes and chemical bonds as edges. The same set of 15 node feature descriptors derived from molecular simulations was used across all models, including induced dipole moment vectors, normalized atomic coordinates, hydrogen bonding counts, solvent accessibility, and local heavy atom density. Adjacency matrices encoded bond connectivity and enabled the aggregation of information from neighboring atoms. Compared to the traditional PROPKA method, the graph neural network models—particularly GAT—demonstrated significantly improved accuracy in pKa prediction by capturing detailed spatial and electrostatic features that empirical approaches often overlook.
Reference