With just a simple alphabet of natural amino acids and common metals, enzymes perform a spectacular number of reactions necessary for life on earth with enviable ease – with neither the extreme conditions nor chemical waste seen from many industrial reactors in human-developed catalysis. Principally, protein structure and related function is not well understood beyond being decisive for its unique selectivity and efficiency. A reliable treatment of larger protein movements/sampling coupled with precise quantum mechanical treatment of metals and bonds breaking/forming (multiscale) is still an open problem in Computational Chemistry. Metalloenzymes can be particularly challenging to model so we present some of the latest in multiscale modeling techniques. We have developed methods sensitive enough to study “selection” of similar metals in enzymes such as HDAC8, where previous literature failed to conclude which metal is active during in vivo catalysis, while fast enough to study larger protein movements including fold stability.

Eminently, multiscale modeling opens the discussion of engineering enzymes to cater to modern day needs in catalysis. Society has not only developed needs in catalysis outside the scope of natural evolution (e.g., drug synthesis, energy conversion) but also has access to more of the periodic table than nature has had a chance to explore. Thus, combining the efficacy of natural enzymes with modern catalytic processes has enormous potential. We have studied iridium, a rare metal with thermodynamic advantages over other metals in promoting catalysis of hydroamination and transfer hydrogenation. In the context of an enzyme we predict catalytic rates near and even exceeding existing organometallic catalysts with further design for specificity available. We are in collaboration with an experimental group in designing a protein fold which both directly accommodates and uses iridium in catalysis. We hope these foundations will support metalloenzyme design efforts towards novel chemical transformations performed as efficiently and environmentally soundly as nature has shown us is possible.

Molecular dynamics (MD) is a powerful tool to study atomic scale changes in proteins underpinning biological and catalytic pathways. It is routinely used to identify small molecule binding sites and study protein folding. MD simulations struggle to achieve quantitatively accurate energies and traditionally sample a fixed chemical state – they do not allow covalent bonds to form or break. However, chemical change underlies many important biochemical processes, including shifts in metal ligation that occur during metalloenzyme activity or when metal ions bind to metalloproteins, or the protonation and deprotonation of amino acids that dictate pH-dependent activity and stability. Understanding these processes at the atomic scale demands computational tools that operate in conjunction with MD to calculate more accurate energies.

Hybrid quantum mechanical/molecular mechanics approaches (QM/MM) can capture the dynamics of chemical change as they provide accurate energies for small regions of interest, such as the active site. We have developed one such hybrid method, QM/DMD, which combines density functional theory (DFT) calculations with discrete molecular dynamics simulations for rapid metalloprotein conformational sampling that can even capture shifts in metal ligation. With QM/DMD we identified the orientational preferences of a cofactor in phenylalanine hydroxylase that distinguishes the disease state of phenylketonuria. To study some forms of chemical change further theory even beyond QM/DMD is necessary. We created a competitive metal affinity (CMA) method, a semi-empirical thermodynamic cycle, to calculate relative metal binding affinities to metalloproteins. With our CMA method and QM/DMD we assessed the ability of human serum transferrin to transport non-native metals and therefore its role in metal cytotoxicity. We also identified the different mechanisms by which Li⁺ and Be�⁺ inhibit glycogen synthase kinase 3β, which could inform future drug design targeting that protein.

In contrast to the detailed QM treatment of a small region needed to study metal behavior, pH-dependent behavior requires rapid sampling of chemical changes across the whole protein in the form of protonation and deprotonation reactions. We therefore developed a constant-pH molecular dynamics method, called titr-DMD, that stochastically updates protonation states based on the efficient DMD method and the semi-empirical electrostatic method Propka. We successfully benchmarked titr-DMD on experimentally verified pH dependent conformational changes in a staphylococcal nuclease mutant. Our work demonstrates the utility of properly modified molecular dynamics, and QM/DMD in particular, to study many forms of chemical change in proteins with good accuracy and speed.

The relative genotoxicity of different polycyclic aromatic hydrocarbons (PAHs) is widely thoughtto be a function of the structural and thermodynamic features of their corresponding PAH-DNA adducts. As a result, accurate parameters for molecular mechanics force fields are crucial to the study of such systems via molecular dynamics (MD). While transferability of parameters among structurally similar molecular systems is frequently a goal when parameterizing novel residues for the CHARMM force field, we will show that planar bay region and non-planar fjord region PAH- DNA adduct systems require distinct dihedral terms to accurately model the torsional potential energy surface of the adduct covalent bond that links a PAH-diol-epoxide and adenine, despite identical atomic connectivity. We then examine the use of the Truncated Singular Value Decom- position and Tikhonov Regularization in standard form to address ill-posed least squares problems Ax = b that frequently arise in molecular mechanics force field parameter optimization. Utilizing the Discrete Picard Condition and/or a well-defined gap in the singular value spectrum when A has a well-determined numerical rank, we are able to systematically determine truncation and in turn regularization parameters that are correspondingly used to produce truncated and regularized solutions to the ill-posed least squares problem at hand. These solutions in turn result in optimized force field dihedral terms that accurately parameterize the torsional energy landscape. As the solu- tions produced by this approach are unique, it has the advantage of avoiding the multiple iterations and guess and check work often required to optimize molecular mechanics force field parameters. With optimized parameters for bay and fjord region PAH-DNA adduct systems developed, we con- duct alchemical free energy perturbation calculations over closed thermodynamic cycles in order to gauge the relative genotoxicities of several IARC Group 2A/B and 3 PAHs in the NRAS(Q61) DNA sequence context. These calculations reveal that the fjord region PAHs examined in this work as well as other IARC Group 2A/B and 3 PAHs exhibit greater relative binding affinity as compared to the IARC Group 1 known human carcinogen B[a]P. These PAHs are also less likely to form productive PAH-DNA protein binding complexes required in the recognition step of global genomic - nucleotide excision repair, indicating that they are more likely to persist and induce mutations in subsequent DNA replication cycles. Further examination reveals that the intercalated conformation and structural differences among PAH-DNA adducts have an impact on stabilizing van der Waals interactions and hydrogen bonding between nucleobase pairs in NRAS(Q61) that are generally associated with trends in relative binding free energies.

Realistic modelling of catalysis requires the incorporation of complexity in the form of fluxionality. This means the inclusion of multiple thermally accessible isomers in the starting ensemble, and along potential reaction pathways. Furthermore, the role of fluxionality on fundamental physical behavior of the system must also be accounted for. This work focuses on extending this fluxionality paradigm both along reaction mechanisms to capture the full complexity of experiment, as well as part of the fundamental physical processes of cluster catalyst deactivation such as sintering or poisoning. Fluxionality is incorporated into DFT calculations via global optimization of cluster catalyst structures and relevant reaction intermediates. It essential to account for isomeric diversity for improved interpretation of experimental results, which enables the identification of novel size-dependent sintering behavior that has subsequently been confirmed experimentally. Furthermore, this allows us to go beyond simple interpretation of experimental results to more complex predictions of changes in the structure and composition of the system. For example, this enabled the pre- diction of self-limiting coke formation in Pt4Ge systems, and aided in understanding why Pt4 in contrast rapidly continues to coke.

WB4.2is one of the hardest metals known. Though not harder than diamond and cubic boron nitride, it surpasses these established hard materials in being cheaper, easier to produce and process, and also more functional. Metal impurities have been shown to affect and in some cases further improve the intrinsic hardness of WB4.2, but the mechanism of hardening remained elusive. In this work we first theoretically elucidate the preferred placements of Ti, V, Cr, Mn, Zr, Nb, Mo, Hf, Ta in the WB4.2 structure, and show these metals to preferentially replace W in two competing positions with respect to the partially occupied B3-cluster site. The impurities avoid the void position in the structure. Next, we analyze the chemical bonding within these identified doped structures, and propose two different mechanisms of strengthening the material, afforded by these impurities, and dependent on their nature. Smaller impurity atoms (Ti, V, Cr, Mn) with deeply lying valence atomic orbitals cause the inter-layer compression of WB4.2, which strengthens the Bhex-Bclusterbonding slightly. Larger impurities (Zr, Nb, Mo, Hf, Ta) with higher-energy valence orbitals, while expanding the structure and negatively impacting the Bhex-Bclusterbonding, also form strong BclusterMbonds. The latter effect is an order of magnitude more substantial than the effect on the Bhex-Bcluster bonding. We conclude that the effect of the impurities on the boride hardness does not simply reduce to structure interlocking due to the size difference between M and W, but instead, has a significant electronic origin.