Force Field X: A computational microscope to study genetic variation and organic crystals using theory and experiment.

Force Field X (FFX) is an open-source software package for atomic resolution modeling of genetic variants and organic crystals that leverages advanced potential energy functions and experimental data. FFX currently consists of nine modular packages with novel algorithms that include global optimization via a many-body expansion, acid-base chemistry using polarizable constant-pH molecular dynamics, estimation of free energy differences, generalized Kirkwood implicit solvent models, and many more. Applications of FFX focus on the use and development of a crystal structure prediction pipeline, biomolecular structure refinement against experimental datasets, and estimation of the thermodynamic effects of genetic variants on both proteins and nucleic acids. The use of Parallel Java and OpenMM combines to offer shared memory, message passing, and graphics processing unit parallelization for high performance simulations. Overall, the FFX platform serves as a computational microscope to study systems ranging from organic crystals to solvated biomolecular systems.

Assessing variants of uncertain significance implicated in hearing loss using a comprehensive deafness proteome.

Hearing loss is the leading sensory deficit, affecting ~ 5% of the population. It exhibits remarkable heterogeneity across 223 genes with 6328 pathogenic missense variants, making deafness-specific expertise a prerequisite for ascribing phenotypic consequences to genetic variants. Deafness-implicated variants are curated in the Deafness Variation Database (DVD) after classification by a genetic hearing loss expert panel and thorough informatics pipeline. However, seventy percent of the 128,167 missense variants in the DVD are "variants of uncertain significance" (VUS) due to insufficient evidence for classification. Here, we use the deep learning protein prediction algorithm, AlphaFold2, to curate structures for all DVD genes. We refine these structures with global optimization and the AMOEBA force field and use DDGun3D to predict folding free energy differences (∆∆GFold) for all DVD missense variants. We find that 5772 VUSs have a large, destabilizing ∆∆GFold that is consistent with pathogenic variants. When also filtered for CADD scores (> 25.7), we determine 3456 VUSs are likely pathogenic at a probability of 99.0%. Of the 224 genes in the DVD, 166 genes (74%) exhibit one or more missense variants predicted to cause a pathogenic change in protein folding stability. The VUSs prioritized here affect 119 patients (~ 3% of cases) sequenced by the OtoSCOPE targeted panel. Approximately half of these patients previously received an inconclusive report, and reclassification of these VUSs as pathogenic provides a new genetic diagnosis for six patients.

A generalized Kirkwood implicit solvent for the polarizable AMOEBA protein model.

Computational simulation of biomolecules can provide important insights into protein design, protein-ligand binding interactions, and ab initio biomolecular folding, among other applications. Accurate treatment of the solvent environment is essential in such applications, but the use of explicit solvents can add considerable cost. Implicit treatment of solvent effects using a dielectric continuum model is an attractive alternative to explicit solvation since it is able to describe solvation effects without the inclusion of solvent degrees of freedom. Previously, we described the development and parameterization of implicit solvent models for small molecules. Here, we extend the parameterization of the generalized Kirkwood (GK) implicit solvent model for use with biomolecules described by the AMOEBA force field via the addition of corrections to the calculation of effective radii that account for interstitial spaces that arise within biomolecules. These include element-specific pairwise descreening scale factors, a short-range neck contribution to describe the solvent-excluded space between pairs of nearby atoms, and finally tanh-based rescaling of the overall descreening integral. We then apply the AMOEBA/GK implicit solvent to a set of ten proteins and achieve an average coordinate root mean square deviation for the experimental structures of 2.0 Å across 500 ns simulations. Overall, the continued development of implicit solvent models will help facilitate the simulation of biomolecules on mechanistically relevant timescales.

Genomic Landscape and Mutational Signatures of Deafness-Associated Genes.

The classification of genetic variants represents a major challenge in the post-genome era by virtue of their extraordinary number and the complexities associated with ascribing a clinical impact, especially for disorders exhibiting exceptional phenotypic, genetic, and allelic heterogeneity. To address this challenge for hearing loss, we have developed the Deafness Variation Database (DVD), a comprehensive, open-access resource that integrates all available genetic, genomic, and clinical data together with expert curation to generate a single classification for each variant in 152 genes implicated in syndromic and non-syndromic deafness. We evaluate 876,139 variants and classify them as pathogenic or likely pathogenic (more than 8,100 variants), benign or likely benign (more than 172,000 variants), or of uncertain significance (more than 695,000 variants); 1,270 variants are re-categorized based on expert curation and in 300 instances, the change is of medical significance and impacts clinical care. We show that more than 96% of coding variants are rare and novel and that pathogenicity is driven by minor allele frequency thresholds, variant effect, and protein domain. The mutational landscape we define shows complex gene-specific variability, making an understanding of these nuances foundational for improved accuracy in variant interpretation in order to enhance clinical decision making and improve our understanding of deafness biology.

In Silico and In Vivo Analysis of Amino Acid Substitutions That Cause Laminopathies.

Mutations in the LMNA gene cause diseases called laminopathies. LMNA encodes lamins A and C, intermediate filaments with multiple roles at the nuclear envelope. LMNA mutations are frequently single base changes that cause diverse disease phenotypes affecting muscles, nerves, and fat. Disease-associated amino acid substitutions were mapped in silico onto three-dimensional structures of lamin A/C, revealing no apparent genotype-phenotype connections. In silico analyses revealed that seven of nine predicted partner protein binding pockets in the Ig-like fold domain correspond to sites of disease-associated amino acid substitutions. Different amino acid substitutions at the same position within lamin A/C cause distinct diseases, raising the question of whether the nature of the amino acid replacement or genetic background differences contribute to disease phenotypes. Substitutions at R249 in the rod domain cause muscular dystrophies with varying severity. To address this variability, we modeled R249Q and R249W in Drosophila Lamin C, an orthologue of LMNA. Larval body wall muscles expressing mutant Lamin C caused abnormal nuclear morphology and premature death. When expressed in indirect flight muscles, R249W caused a greater number of adults with wing posturing defects than R249Q, consistent with observations that R249W and R249Q cause distinct muscular dystrophies, with R249W more severe. In this case, the nature of the amino acid replacement appears to dictate muscle disease severity. Together, our findings illustrate the utility of Drosophila for predicting muscle disease severity and pathogenicity of variants of unknown significance.

Nanophthalmos patient with a THR518MET mutation in MYRF, a case report.

BACKGROUND: Nanophthalmos has a significant genetic background and disease-causing mutations have been recently been reported in the myelin regulatory factor (MYRF) gene. We report clinical features in a patient with nanophthalmos and a Thr518Met MYRF mutation. CASE PRESENTATION: A three-year-old male was discovered to have nanophthalmos after first presenting to the emergency department for a frontal headache, eye pain, emesis, and lethargy. Imaging studies (CT and MRI) were negative except for increased posterior fossa cerebrospinal fluid. Subsequent examinations revealed nanophthalmos (short axial eye lengths 18.1 mm OD and 18.3 mm OS), microcornea, and a large crystalline lens. Peripheral chorioretinal pigment abnormalities were also observed. He experienced episodes of marked ocular hypertension (53 mmHg OD and 60 mmHg) likely due to intermittent angle closure precipitated by nanophthalmos. The ocular hypertension was responsive to topical medicines. Genetic analysis of known nanophthalmos genes MFRP and TMEM98 were negative, while a novel mutation, Thr518Met was detected in MYRF. The Thr518Met mutation was absent from 362 matched normal controls and was extremely rare in a large population database, allele frequency of 0.000024. The Thr518Met mutation altered a highly conserved amino acid in the MYRF protein and three of four algorithms suggested that this mutation is likely pathogenic. Finally, molecular modeling showed that the Thr518Met mutation is damaging to MYRF structure. Together these data suggest that the Thr518Met mutation causes nanophthalmos. CONCLUSIONS: Nanophthalmos may present at an early age with features of angle closure glaucoma and a Thr518Met mutation in MYRF was detected in a patient with nanophthalmos. Prevalence data, homology data, mutation analysis data, and protein modeling data suggest that this variant is pathogenic and may expand the phenotypic range of syndromic nanophthalmos caused by MYRF mutations to include central nervous system abnormalities (increased posterior fossa cerebrospinal fluid).

Characterization of a TP53 Somatic Variant of Unknown Function From an Ovarian Cancer Patient Using Organoid Culture and Computational Modeling.

In our proof-of-concept study of 1 patient with stage IIIC carcinosarcoma of the ovary, we discovered a rare mutation in the tumor suppressor, TP53, that results in the deletion of N131. Immunofluorescence imaging of the organoid culture revealed hyperstaining of p53 protein. Computational modeling suggests this residue is important for maintaining protein conformation. Drug screening identified the combination of a proteasome inhibitor with a histone deacetylase inhibitor as the most effective treatment. These data provide evidence for the successful culture of a patient tumor and analysis of drug response ex vivo.

Structural Insights into Hearing Loss Genetics from Polarizable Protein Repacking.

Hearing loss is associated with ∼8100 mutations in 152 genes, and within the coding regions of these genes are over 60,000 missense variants. The majority of these variants are classified as "variants of uncertain significance" to reflect our inability to ascribe a phenotypic effect to the observed amino acid change. A promising source of pathogenicity information is biophysical simulation, although input protein structures often contain defects because of limitations in experimental data and/or only distant homology to a template. Here, we combine the polarizable atomic multipole optimized energetics for biomolecular applications force field, many-body optimization theory, and graphical processing unit acceleration to repack all deafness-associated proteins and thereby improve average structure MolProbity score from 2.2 to 1.0. We then used these optimized wild-type models to create over 60,000 structures for missense variants in the Deafness Variation Database, which are being incorporated into the Deafness Variation Database to inform deafness pathogenicity prediction. Finally, this work demonstrates that advanced polarizable atomic multipole force fields are efficient enough to repack the entire human proteome.

Hypomorphic mutations in TRNT1 cause retinitis pigmentosa with erythrocytic microcytosis.

Retinitis pigmentosa (RP) is a highly heterogeneous group of disorders characterized by degeneration of the retinal photoreceptor cells and progressive loss of vision. While hundreds of mutations in more than 100 genes have been reported to cause RP, discovering the causative mutations in many patients remains a significant challenge. Exome sequencing in an individual affected with non-syndromic RP revealed two plausibly disease-causing variants in TRNT1, a gene encoding a nucleotidyltransferase critical for tRNA processing. A total of 727 additional unrelated individuals with molecularly uncharacterized RP were completely screened for TRNT1 coding sequence variants, and a second family was identified with two members who exhibited a phenotype that was remarkably similar to the index patient. Inactivating mutations in TRNT1 have been previously shown to cause a severe congenital syndrome of sideroblastic anemia, B-cell immunodeficiency, recurrent fevers and developmental delay (SIFD). Complete blood counts of all three of our patients revealed red blood cell microcytosis and anisocytosis with only mild anemia. Characterization of TRNT1 in patient-derived cell lines revealed reduced but detectable TRNT1 protein, consistent with partial function. Suppression of trnt1 expression in zebrafish recapitulated several features of the human SIFD syndrome, including anemia and sensory organ defects. When levels of trnt1 were titrated, visual dysfunction was found in the absence of other phenotypes. The visual defects in the trnt1-knockdown zebrafish were ameliorated by the addition of exogenous human TRNT1 RNA. Our findings indicate that hypomorphic TRNT1 mutations can cause a recessive disease that is almost entirely limited to the retina.

LADD syndrome with glaucoma is caused by a novel gene.

PURPOSE: Lacrimo-auriculo-dento-digital (LADD) syndrome is an autosomal dominant disorder displaying variable expression of multiple congenital anomalies including hypoplasia or aplasia of the lacrimal and salivary systems causing abnormal tearing and dry mouth. Mutations in the FGF10, FGFR2, and FGFR3 genes were found to cause some cases of LADD syndrome in prior genetic studies. The goal of this study is to identify the genetic basis of a case of LADD syndrome with glaucoma and thin central corneal thickness (CCT). METHODS: Whole exome sequencing was performed, and previously described disease-causing genes (FGF10, FGFR2, and FGFR3) were first evaluated for mutations. Fifty-eight additional prioritized candidate genes were identified by searching gene annotations for features of LADD syndrome. The potential pathogenicity of the identified mutations was assessed by determining their frequency in large public exome databases; through sequence analysis using the Blosum62 matrix, PolyPhen2, and SIFT algorithms; and through homology analyses. A structural analysis of the effects of the top candidate mutation in tumor protein 63 (TP63) was also conducted by superimposing the mutation over the solved crystal structure. RESULTS: No mutations were detected in FGF10, FGFR2, or FGFR3. The LADD syndrome patient's exome data was searched for mutations in the 58 candidate genes and only one mutation was detected, an Arg343Trp mutation in the tumor protein 63 (TP63) gene. This TP63 mutation is absent from the gnomAD sequence database. Analysis of the Arg343Trp mutation with Blosum62, PolyPhen2, and SIFT all suggest it is pathogenic. This arginine residue is highly conserved in orthologous genes. Finally, crystal structure analysis showed that the Arg343Trp mutation causes a significant alteration in the structure of TP63's DNA binding domain. CONCLUSIONS: We report a patient with no mutations in known LADD syndrome genes (FGF10, FGFR2, and FGFR3). Our analysis provides strong evidence that the Arg343Trp mutation in TP63 caused LADD syndrome in our patient and that TP63 is a fourth gene contributing to this condition. TP63 encodes a transcription factor involved in the development and differentiation of tissues affected by LADD syndrome. These data suggest that TP63 is a novel LADD syndrome gene and may also influence corneal thickness and risk for open-angle glaucoma.

High-Throughput Genetic Testing for Thrombotic Microangiopathies and C3 Glomerulopathies.

The thrombotic microangiopathies (TMAs) and C3 glomerulopathies (C3Gs) include a spectrum of rare diseases such as atypical hemolytic uremic syndrome, thrombotic thrombocytopenic purpura, C3GN, and dense deposit disease, which share phenotypic similarities and underlying genetic commonalities. Variants in several genes contribute to the pathogenesis of these diseases, and identification of these variants may inform the diagnosis and treatment of affected patients. We have developed and validated a comprehensive genetic panel that screens all exons of all genes implicated in TMA and C3G. The closely integrated pipeline implemented includes targeted genomic enrichment, massively parallel sequencing, bioinformatic analysis, and a multidisciplinary conference to analyze identified variants in the context of each patient's specific phenotype. Herein, we present our 1-year experience with this panel, during which time we studied 193 patients. We identified 17 novel and 74 rare variants, which we classified as pathogenic (11), likely pathogenic (12), and of uncertain significance (68). Compared with controls, patients with C3G had a higher frequency of rare and novel variants in C3 convertase (C3 and CFB) and complement regulator (CFH, CFI, CFHR5, and CD46) genes (P<0.05). In contrast, patients with TMA had an increase in rare and novel variants only in complement regulator genes (P<0.01), a distinction consistent with differing sites of complement dysregulation in these two diseases. In summary, we were able to provide a positive genetic diagnosis in 43% and 41% of patients carrying the clinical diagnosis of C3G and TMA, respectively.

Toward polarizable AMOEBA thermodynamics at fixed charge efficiency using a dual force field approach: application to organic crystals.

First principles prediction of the structure, thermodynamics and solubility of organic molecular crystals, which play a central role in chemical, material, pharmaceutical and engineering sciences, challenges both potential energy functions and sampling methodologies. Here we calculate absolute crystal deposition thermodynamics using a novel dual force field approach whose goal is to maintain the accuracy of advanced multipole force fields (e.g. the polarizable AMOEBA model) while performing more than 95% of the sampling in an inexpensive fixed charge (FC) force field (e.g. OPLS-AA). Absolute crystal sublimation/deposition phase transition free energies were determined using an alchemical path that grows the crystalline state from a vapor reference state based on sampling with the OPLS-AA force field, followed by dual force field thermodynamic corrections to change between FC and AMOEBA resolutions at both end states (we denote the three step path as AMOEBA/FC). Importantly, whereas the phase transition requires on the order of 200 ns of sampling per compound, only 5 ns of sampling was needed for the dual force field thermodynamic corrections to reach a mean statistical uncertainty of 0.05 kcal mol-1. For five organic compounds, the mean unsigned error between direct use of AMOEBA and the AMOEBA/FC dual force field path was only 0.2 kcal mol-1 and not statistically significant. Compared to experimental deposition thermodynamics, the mean unsigned error for AMOEBA/FC (1.4 kcal mol-1) was more than a factor of two smaller than uncorrected OPLS-AA (3.2 kcal mol-1). Overall, the dual force field thermodynamic corrections reduced condensed phase sampling in the expensive force field by a factor of 40, and may prove useful for protein stability or binding thermodynamics in the future.

Calculating binding free energies of host-guest systems using the AMOEBA polarizable force field.

Molecular recognition is of paramount interest in many applications. Here we investigate a series of host-guest systems previously used in the SAMPL4 blind challenge by using molecular simulations and the AMOEBA polarizable force field. The free energy results computed by Bennett's acceptance ratio (BAR) method using the AMOEBA polarizable force field ranked favorably among the entries submitted to the SAMPL4 host-guest competition [Muddana, et al., J. Comput.-Aided Mol. Des., 2014, 28, 305-317]. In this work we conduct an in-depth analysis of the AMOEBA force field host-guest binding thermodynamics by using both BAR and the orthogonal space random walk (OSRW) methods. The binding entropy-enthalpy contributions are analyzed for each host-guest system. For systems of inordinate binding entropy-enthalpy values, we further examine the hydrogen bonding patterns and configurational entropy contribution. The binding mechanism of this series of host-guest systems varies from ligand to ligand, driven by enthalpy and/or entropy changes. Convergence of BAR and OSRW binding free energy methods is discussed. Ultimately, this work illustrates the value of molecular modelling and advanced force fields for the exploration and interpretation of binding thermodynamics.

Dead-End Elimination with a Polarizable Force Field Repacks PCNA Structures.

A balance of van der Waals, electrostatic, and hydrophobic forces drive the folding and packing of protein side chains. Although such interactions between residues are often approximated as being pairwise additive, in reality, higher-order many-body contributions that depend on environment drive hydrophobic collapse and cooperative electrostatics. Beginning from dead-end elimination, we derive the first algorithm, to our knowledge, capable of deterministic global repacking of side chains compatible with many-body energy functions. The approach is applied to seven PCNA x-ray crystallographic data sets with resolutions 2.5-3.8 Å (mean 3.0 Å) using an open-source software. While PDB_REDO models average an Rfree value of 29.5% and MOLPROBITY score of 2.71 Å (77th percentile), dead-end elimination with the polarizable AMOEBA force field lowered Rfree by 2.8-26.7% and improved mean MOLPROBITY score to atomic resolution at 1.25 Å (100th percentile). For structural biology applications that depend on side-chain repacking, including x-ray refinement, homology modeling, and protein design, the accuracy limitations of pairwise additivity can now be eliminated via polarizable or quantum mechanical potentials.

PDZD7 and hearing loss: More than just a modifier.

Deafness is the most frequent sensory disorder. With over 90 genes and 110 loci causally implicated in non-syndromic hearing loss, it is phenotypically and genetically heterogeneous. Here, we investigate the genetic etiology of deafness in four families of Iranian origin segregating autosomal recessive non-syndromic hearing loss (ARNSHL). We used a combination of linkage analysis, homozygosity mapping, and a targeted genomic enrichment platform to simultaneously screen 90 known deafness-causing genes for pathogenic variants. Variant segregation was confirmed by Sanger sequencing. Linkage analysis and homozygosity mapping showed segregation with the DFNB57 locus on chromosome 10 in two families. Targeted genomic enrichment with massively parallel sequencing identified causal variants in PDZD7: a homozygous missense variant (p.Gly103Arg) in one family and compound heterozygosity for missense (p.Met285Arg) and nonsense (p.Tyr500Ter) variants in the second family. Screening of two additional families identified two more variants: (p.Gly228Arg) and (p.Gln526Ter). Variant segregation with the hearing loss phenotype was confirmed in all families by Sanger sequencing. The missense variants are predicted to be deleterious, and the two nonsense mutations produce null alleles. This report is the first to show that mutations in PDZD7 cause ARNSHL, a finding that offers addition insight into the USH2 interactome. We also describe a novel likely disease-causing mutation in CIB2 and illustrate the complexity associated with gene identification in diseases that exhibit large genetic and phenotypic heterogeneity.

Biomolecular electrostatics and solvation: a computational perspective.

An understanding of molecular interactions is essential for insight into biological systems at the molecular scale. Among the various components of molecular interactions, electrostatics are of special importance because of their long-range nature and their influence on polar or charged molecules, including water, aqueous ions, proteins, nucleic acids, carbohydrates, and membrane lipids. In particular, robust models of electrostatic interactions are essential for understanding the solvation properties of biomolecules and the effects of solvation upon biomolecular folding, binding, enzyme catalysis, and dynamics. Electrostatics, therefore, are of central importance to understanding biomolecular structure and modeling interactions within and among biological molecules. This review discusses the solvation of biomolecules with a computational biophysics view toward describing the phenomenon. While our main focus lies on the computational aspect of the models, we provide an overview of the basic elements of biomolecular solvation (e.g. solvent structure, polarization, ion binding, and non-polar behavior) in order to provide a background to understand the different types of solvation models.

Endothelial lipase variant T111I does not alter inhibition by angiopoietin-like proteins.

High levels of HDL-C are correlated with a decreased risk of cardiovascular disease. HDL-C levels are modulated in part by the secreted phospholipase, endothelial lipase (EL), which hydrolyzes the phospholipids of HDL and decreases circulating HDL-C concentrations. A 584C/T polymorphism in LIPG, the gene which encodes EL, was first identified in individuals with increased HDL levels. This polymorphism results in a T111I point mutation the EL protein. The association between this variant, HDL levels, and the risk of coronary artery disease (CAD) in humans has been extensively studied, but the findings have been inconsistent. In this study, we took a biochemical approach, investigating how the T111I variant affected EL activity, structure, and stability. Moreover, we tested whether the T111I variant altered the inhibition of phospholipase activity by angiopoietin-like 3 (ANGPTL3) and angiopoietin-like 4 (ANGPTL4), two known EL inhibitors. We found that neither the stability nor enzymatic activity of EL was altered by the T111I variant. Moreover, we found no difference between wild-type and T111I EL in their ability to be inhibited by ANGPTL proteins. These data suggest that any effect this variant may have on HDL-C levels or cardiovascular disease are not mediated through alterations in these functions.

METTL23 Variants and Patients With Normal-Tension Glaucoma.

Importance: This research confirms and further establishes that pathogenic variants in a fourth gene, METTL23, are associated with autosomal dominant normal-tension glaucoma (NTG). Objective: To determine the frequency of glaucoma-causing pathogenic variants in the METTL23 gene in a cohort of patients with NTG from Iowa. Design, Setting, and Participants: This case-control study took place at a single tertiary care center in Iowa from January 1997 to January 2024, with analysis occurring between January 2023 and January 2024. Two groups of participants were enrolled from the University of Iowa clinics: 331 patients with NTG and 362 control individuals without glaucoma. Patients with a history of trauma; steroid use; stigmata of pigment dispersion syndrome; exfoliation syndrome; or pathogenic variants in MYOC, TBK1, or OPTN were also excluded. Main Outcomes and Measures: Detection of an enrichment of METTL23 pathogenic variants in individuals with NTG compared with control individuals without glaucoma. Results: The study included 331 patients with NTG (mean [SD] age, 68.0 [11.7] years; 228 [68.9%] female and 103 [31.1%] male) and 362 control individuals without glaucoma (mean [SD] age, 64.5 [12.6] years; 207 [57.2%] female and 155 [42.8%] male). There were 5 detected instances of 4 unique METTL23 pathogenic variants in patients with NTG. Three METTL23 variants-p.Ala7Val, p.Pro22Arg, and p.Arg63Trp-were judged to be likely pathogenic and were detected in 3 patients (0.91%) with NTG. However, when all detected variants were evaluated with either mutation burden analysis or logistic regression, their frequency was not statistically higher in individuals with NTG than in control individuals without glaucoma (1.5% vs 2.5%; P = .27). Conclusion and Relevance: This investigation provides evidence that pathogenic variants in METTL23 are associated with NTG. Within an NTG cohort at a tertiary care center, pathogenic variants were associated with approximately 1% of NTG cases, a frequency similar to that of other known normal-tension glaucoma genes, including optineurin (OPTN), TANK-binding kinase 1 (TBK1), and myocilin (MYOC). The findings suggest that METTL23 pathogenic variants are likely involved in a biologic pathway that is associated with glaucoma that occurs at lower intraocular pressures.

Crystal Polymorph Search in the NPT Ensemble via a Deposition/Sublimation Alchemical Path.

The formulation of active pharmaceutical ingredients involves discovering stable crystal packing arrangements or polymorphs, each of which has distinct pharmaceutically relevant properties. Traditional experimental screening techniques utilizing various conditions are commonly supplemented with in silico crystal structure prediction (CSP) to inform the crystallization process and mitigate risk. Predictions are often based on advanced classical force fields or quantum mechanical calculations that model the crystal potential energy landscape but do not fully incorporate temperature, pressure, or solution conditions during the search procedure. This study proposes an innovative alchemical path that utilizes an advanced polarizable atomic multipole force field to predict crystal structures based on direct sampling of the NPT ensemble. The use of alchemical (i.e., nonphysical) intermediates, a novel Monte Carlo barostat, and an orthogonal space tempering bias combine to enhance the sampling efficiency of the deposition/sublimation phase transition. The proposed algorithm was applied to 2-((4-(2-(3,4-dichlorophenyl)ethyl)phenyl)amino)benzoic acid (Cambridge Crystallography Database Centre ID: XAFPAY) as a case study to showcase the algorithm. Each experimentally determined polymorph with one molecule in the asymmetric unit was successfully reproduced via approximately 1000 short 1 ns simulations per space group where each simulation was initiated from random rigid body coordinates and unit cell parameters. Utilizing two threads of a recent Intel CPU (a Xeon Gold 6330 CPU at 2.00 GHz), 1 ns of sampling using the polarizable AMOEBA force field can be acquired in 4 h (equating to more than 300 ns/day using all 112 threads/56 cores of a dual CPU node) within the Force Field X software (https://ffx.biochem.uiowa.edu). These results demonstrate a step forward in the rigorous use of the NPT ensemble during the CSP search process and open the door to future algorithms that incorporate solution conditions using continuum solvation methods.

Constant-pH Simulations with the Polarizable Atomic Multipole AMOEBA Force Field.

Accurately predicting protein behavior across diverse pH environments remains a significant challenge in biomolecular simulations. Existing constant-pH molecular dynamics (CpHMD) algorithms are limited to fixed-charge force fields, hindering their application to biomolecular systems described by permanent atomic multipoles or induced dipoles. This work overcomes these limitations by introducing the first polarizable CpHMD algorithm in the context of the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field. Additionally, our implementation in the open-source Force Field X (FFX) software has the unique ability to handle titration state changes for crystalline systems including flexible support for all 230 space groups. The evaluation of constant-pH molecular dynamics (CpHMD) with the AMOEBA force field was performed on 11 crystalline peptide systems that span the titrating amino acids (Asp, Glu, His, Lys, and Cys). Titration states were correctly predicted for 15 out of the 16 amino acids present in the 11 systems, including for the coordination of Zn2+ by cysteines. The lone exception was for a HIS-ALA peptide where CpHMD predicted both neutral histidine tautomers to be equally populated, whereas the experimental model did not consider multiple conformers and diffraction data are unavailable for rerefinement. This work demonstrates the promise polarizable CpHMD simulations for pKa predictions, the study of biochemical mechanisms such as the catalytic triad of proteases, and for improved protein-ligand binding affinity accuracy in the context of pharmaceutical lead optimization.