Tale of Two Polymorphs: Investigating the Structural Differences and Dynamic Relationship between Nirmatrelvir Solid Forms (Paxlovid).

Two anhydrous polymorphs of the novel antiviral medicine nirmatrelvir were discovered during the development of Paxlovid, Pfizer's oral Covid-19 treatment. A comprehensive experimental and computational approach was necessary to distinguish the two closely related polymorphs, herein identified as Forms 1 and 4. This approach paired experimental methods, including powder X-ray diffraction and single-crystal X-ray diffraction, solid-state experimental methods, thermal analysis, solid-state nuclear magnetic resonance and Raman spectroscopy with computational investigations comprising crystal structure prediction, Gibbs free energy calculations, and molecular dynamics simulations of the polymorphic transition. Forms 1 and 4 were ultimately determined to be enantiotropically related polymorphs with Form 1 being the stable form above the transition temperature of ∼17 °C and designated as the nominated form for drug development. The work described in this paper shows the importance of using highly specialized orthogonal approaches to elucidate the subtle differences in structure and properties of similar solid-state forms. This synergistic approach allowed for unprecedented speed in bringing Paxlovid to patients in record time amidst the pandemic.

Leveraging Language Model Multitasking To Predict C-H Borylation Selectivity.

C-H borylation is a high-value transformation in the synthesis of lead candidates for the pharmaceutical industry because a wide array of downstream coupling reactions is available. However, predicting its regioselectivity, especially in drug-like molecules that may contain multiple heterocycles, is not a trivial task. Using a data set of borylation reactions from Reaxys, we explored how a language model originally trained on USPTO_500_MT, a broad-scope set of patent data, can be used to predict the C-H borylation reaction product in different modes: product generation and site reactivity classification. Our fine-tuned T5Chem multitask language model can generate the correct product in 79% of cases. It can also classify the reactive aromatic C-H bonds with 95% accuracy and 88% positive predictive value, exceeding purpose-developed graph-based neural networks.

Kinetic Modeling of API Oxidation: (2) Imipramine Stress Testing.

Gauging the chemical stability of active pharmaceutical ingredients (APIs) is critical at various stages of pharmaceutical development to identify potential risks from drug degradation and ensure the quality and safety of the drug product. Stress testing has been the major experimental method to study API stability, but this analytical approach is time-consuming, resource-intensive, and limited by API availability, especially during the early stages of drug development. Novel computational chemistry methods may assist in screening for API chemical stability prior to synthesis and augment contemporary API stress testing studies, with the potential to significantly accelerate drug development and reduce costs. In this work, we leverage quantum chemical calculations and automated reaction mechanism generation to provide new insights into API degradation studies. In the continuation of part one in this series of studies [Grinberg Dana et al., Mol. Pharm. 2021 18 (8), 3037-3049], we have generated the first ab initio predictive chemical kinetic model of free-radical oxidative degradation for API stress testing. We focused on imipramine oxidation in an azobis(isobutyronitrile) (AIBN)/H2O/CH3OH solution and compared the model's predictions with concurrent experimental observations. We analytically determined iminodibenzyl and desimipramine as imipramine's two major degradation products under industry-standard AIBN stress testing conditions, and our ab initio kinetic model successfully identified both of them in its prediction for the top three degradation products. This work shows the potential and utility of predictive chemical kinetic modeling and quantum chemical computations to elucidate API chemical stability issues. Further, we envision an automated digital workflow that integrates first-principle models with data-driven methods that, when actively and iteratively combined with high-throughput experiments, can substantially accelerate and transform future API chemical stability studies.

Kinetic Modeling of API Oxidation: (1) The AIBN/H2O/CH3OH Radical "Soup".

Stress testing of active pharmaceutical ingredients (API) is an important tool used to gauge chemical stability and identify potential degradation products. While different flavors of API stress testing systems have been used in experimental investigations for decades, the detailed kinetics of such systems as well as the chemical composition of prominent reactive species, specifically reactive oxygen species, are unknown. As a first step toward understanding and modeling API oxidation in stress testing, we investigated a typical radical "soup" solution an API is subject to during stress testing. Here we applied ab initio electronic structure calculations to automatically generate and refine a detailed chemical kinetics model, taking a fresh look at API oxidation. We generated a detailed kinetic model for a representative azobis(isobutyronitrile) (AIBN)/H2O/CH3OH stress-testing system with a varied cosolvent ratio (50%/50%-99.5%/0.5% vol water/methanol) for 5.0 mM AIBN and representative pH values of 4-10 at 40 °C that was stirred and open to the atmosphere. At acidic conditions, hydroxymethyl alkoxyl is the dominant alkoxyl radical, and at basic conditions, for most studied initial methanol concentrations, cyanoisopropyl alkoxyl becomes the dominant alkoxyl radical, albeit at an overall lower concentration. At acidic conditions, the levels of cyanoisopropyl peroxyl, hydroxymethyl peroxyl, and hydroperoxyl radicals are relatively high and comparable, while, at both neutral and basic pH conditions, superoxide becomes the prominent radical in the system. The present work reveals the prominent species in a common model API stress testing system at various cosolvent and pH conditions, sets the stage for an in-depth quantitative API kinetic study, and demonstrates the usage of novel software tools for automated chemical kinetic model generation and ab initio refinement.

Identifying Conformational Isomers of Organic Molecules in Solution via Unsupervised Clustering.

We present a systematic approach for the identification of statistically relevant conformational macrostates of organic molecules from molecular dynamics trajectories. The approach applies to molecules characterized by an arbitrary number of torsional degrees of freedom and enables the transferability of the macrostates definition across different environments. We formulate a dissimilarity measure between molecular configurations that incorporates information on the characteristic energetic cost associated with transitions along all relevant torsional degrees of freedom. Such metric is employed to perform unsupervised clustering of molecular configurations based on the Fast Search and Find of Density Peaks algorithm. We apply this method to investigate the equilibrium conformational ensemble of Sildenafil, a conformationally complex pharmaceutical compound, in different environments including the crystal bulk, the gas phase, and three different solvents (acetonitrile, 1-butanol, and toluene). We demonstrate that while Sildenafil can adopt more than 100 metastable conformational configurations, only 12 are significantly populated across all of the environments investigated. Despite the complexity of the conformational space, we find that the most abundant conformers in solution are the closest to the conformers found in the most common Sildenafil crystal phase.

Solvation of the Glycyl Radical.

The effect of adding explicit water molecules to the neutral (N) and zwitterionic (Z) forms of the glycyl radical has been examined. The results show that a minimum of three water molecules is required to stabilize the Z radical as a local minimum, with an energy gap of 123 kJ mol-1 between the N and Z forms at this point, in favor of the N form. Increasing the number of water molecules to ∼20 leads to a converged Z-N energy difference of ∼50 kJ mol-1 still in favor of the N form, even though the radical is not considered fully solvated from a structural point of view. Thus, energetic convergence is determined mainly by solvation of the polar functional groups, and a complete coverage of the entire molecule is not necessary. Because aqueous closed-shell glycine exists as a zwitterion while aqueous glycyl radical prefers the neutral form, the conversion between the two necessitates a change along the hydrogen-abstraction reaction pathway. In this regard, the transition structure for α-hydrogen abstraction by the ·OH radical has greater resemblance to glycine than to the glycyl radical. Overall, the barrier for hydrogen abstraction from Z glycine is larger than that from the N isomer, and this might act to provide some protection against radical damage to the free amino acid in the (aqueous) biological environment.

Drug-Excipient Interactions in the Solid State: The Role of Different Stress Factors.

Understanding properties and mechanisms that govern drug degradation in the solid state is of high importance to ensure drug stability and safety of solid dosage forms. In this study, we attempt to understand drug-excipient interactions in the solid state using both theoretical and experimental approaches. The model active pharmaceutical ingredients (APIs) under study are carvedilol (CAR) and codeine phosphate (COP), which are known to undergo esterification with citric acid (CA) in the solid state. Starting from the crystal structures of two different polymorphs of each compound, we calculated the exposure and accessibility of reactive hydroxyl groups for a number of relevant crystal surfaces, as well as descriptors that could be associated with surface stabilities using molecular simulations. Accelerated degradation experiments at elevated temperature and controlled humidity were conducted to assess the propensity of different solid forms of the model APIs to undergo chemical reactions with anhydrous CA or CA monohydrate. In addition, for CAR, we studied the solid state degradation at varying humidity levels and also under mechano-activation. Regarding the relative degradation propensities, we found that variations in the exposure and accessibility of molecules on the crystal surface play a minor role compared to the impact of molecular mobility due to different levels of moisture. We further studied drug-excipient interactions under mechano-activation (comilling of API and CA) and found that the reaction proceeded even faster than in physical powder mixtures kept at accelerated storage conditions.

Molecular Investigation of the Mechanism of Non-Enzymatic Hydrolysis of Proteins and the Predictive Algorithm for Susceptibility.

A number of potential degradation routes can limit the shelf life of a biotherapeutic. While these are experimentally measurable, the tests to do so require a substantial investment in both time and material, resources rarely available early in the drug development process. To address the potential degradation route of non-enzymatic hydrolysis, we performed a molecular modeling analysis, together with an experimental study, to gain detailed insight into the reaction. On the basis of the mechanism, an algorithm for predicting the likely cleavage sites of a protein has been created. This algorithm measures four key properties during a molecular dynamics simulation, which relate to the key steps of the hydrolysis mechanism, in particular the rate-determining step (which can vary depending on the local environment). The first two properties include the secondary structure and the surface exposure of the amide bond, both of which help detect if the addition of the proton to the amide bond is possible. The second two properties relate to whether the side chain can cyclize and form a furane ring. These two properties are the orientation of the side chain relative to the amide bond and the number of hydrogen bonds between the side chain and the surrounding protein. Overall, the algorithm performs well at identifying reactive versus nonreactive bonds. The algorithm correctly classifies nearly 90% of all amide bonds following an aspartic or glutamic acid residue as reactive or nonreactive.

Mechanistic Insights into Radical-Mediated Oxidation of Tryptophan from ab Initio Quantum Chemistry Calculations and QM/MM Molecular Dynamics Simulations.

An assessment of the mechanisms of (•)OH and (•)OOH radical-mediated oxidation of tryptophan was performed using density functional theory calculations and ab initio plane-wave Quantum Mechanics/Molecular Mechanics (QM/MM) molecular dynamics simulations. For the (•)OH reactions, addition to the pyrrole ring at position 2 is the most favored site with a barrierless reaction in the gas phase. The subsequent degradation of this adduct through a H atom transfer to water was intermittently observed in aqueous-phase molecular dynamics simulations. For the (•)OOH reactions, addition to the pyrrole ring at position 2 is the most favored pathway, in contrast to the situation in the model system ethylene, where concerted addition to the double bond is preferred. From the (•)OOH position 2 adduct QM/MM simulations show that formation of oxy-3-indolanaline occurs readily in an aqueous environment. The observed transformation starts from an initial rupture of the O-O bond followed by a H atom transfer with the accompanying loss of an (•)OH radical to solution. Finally, classical molecular dynamics simulations were performed to equate observed differential oxidation rates of various tryptophan residues in monoclonal antibody fragments. It was found that simple parameters derived from simulation correlate well with the experimental data.

Reactions of benzene and 3-methylpyrrole with the •OH and •OOH radicals: an assessment of contemporary density functional theory methods.

A high-level quantum chemistry investigation has been carried out for the addition and abstraction reactions by the radicals (•)OH and (•)OOH to and from the model alkenes 3-methylpyrrole and benzene. These models were chosen to reflect the functionalities contained in the side chain of the amino acid tryptophan. The W1BD procedure was used to calculate benchmark barriers and reaction energies for the smaller model system of (•)OOH addition to ethylene. It was found that the CBS-QB3 methodology compares best with the W1BD benchmark, demonstrating a mean absolute deviation (MAD) from W1BD of 3.9 kJ mol(-1). For the reactions involving the (•)OH radical and benzene or 3-methylpyrrole, addition is favored over abstraction in all cases. In particular the CBS-QB3 calculations suggest a barrierless addition reaction of the (•)OH radical to position two of 3-methylpyrrole. For the analogous addition and abstraction reactions involving the (•)OOH radical, the same order of reactivity was found, albeit with higher barriers. A number of other processes involving the addition of the (•)OOH radical were also investigated. The main findings of these studies determined that the initial (•)OOH barrier of stepwise addition to 3-methylpyrrole (+18.8 kJ mol(-1)) is significantly smaller than the concerted addition barrier (+71.5 kJ mol(-1)). This conclusion contrasts starkly with the situation for ethylene in which it is well established that the concerted process has the smaller barrier. A considerable variety of contemporary density functional theory procedures have been tested to examine their accuracy in predicting the CBS-QB3 results. It was found that the best overall performing method was UBMK with an MAD of 7.3 kJ mol(-1). A number of other functionals additionally performed well. They included UM06, RM06, UXYG3 and RXYG3, all of which have MADs of less than 8 kJ mol(-1).

Insights and Challenges in Correcting Force Field Based Solvation Free Energies Using a Neural Network Potential.

We present a comprehensive study investigating the potential gain in accuracy for calculating absolute solvation free energies (ASFE) using a neural network potential to describe the intramolecular energy of the solute. We calculated the ASFE for most compounds from the FreeSolv database using the Open Force Field (OpenFF) and compared them to earlier results obtained with the CHARMM General Force Field (CGenFF). By applying a nonequilibrium (NEQ) switching approach between the molecular mechanics (MM) description (either OpenFF or CGenFF) and the neural net potential (NNP)/MM level of theory (using ANI-2x as the NNP potential), we attempted to improve the accuracy of the calculated ASFEs. The predictive performance of the results did not change when this approach was applied to all 589 small molecules in the FreeSolv database that ANI-2x can describe. When selecting a subset of 156 molecules, focusing on compounds where the force fields performed poorly, we saw a slight improvement in the root-mean-square error (RMSE) and mean absolute error (MAE). The majority of our calculations utilized unidirectional NEQ protocols based on Jarzynski's equation. Additionally, we conducted bidirectional NEQ switching for a subset of 156 solutes. Notably, only a small fraction (10 out of 156) exhibited statistically significant discrepancies between unidirectional and bidirectional NEQ switching free energy estimates.

Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation.

Biomolecular simulations have become an essential tool in contemporary drug discovery, and molecular mechanics force fields (FFs) constitute its cornerstone. Developing a high quality and broad coverage general FF is a significant undertaking that requires substantial expert knowledge and computing resources, which is beyond the scope of general practitioners. Existing FFs originate from only a limited number of groups and organizations, and they either suffer from limited numbers of training sets, lower than desired quality because of oversimplified representations, or are costly for the molecular modeling community to access. To address these issues, in this work, we developed an AMBER-consistent small molecule FF with extensive chemical space coverage, and we provide Open Access parameters for the entire modeling community. To validate our FF, we carried out benchmarks of quantum mechanics (QM)/molecular mechanics conformer comparison and free energy perturbation calculations on several benchmark data sets. Our FF achieves a higher level of performance at reproducing QM energies and geometries than two popular open-source FFs, OpenFF2 and GAFF2. In relative binding free energy calculations for 31 protein-ligand data sets, comprising 1079 pairs of ligands, the new FF achieves an overall root-mean-square error of 1.19 kcal/mol for ΔΔG and 0.92 kcal/mol for ΔG on a subset of 463 ligands without bespoke fitting to the data sets. The results are on par with those of the leading commercial series of OPLS FFs.

Templated nucleation of acetaminophen on spherical excipient agglomerates.

We investigated the effect of spherical agglomeration of heterogeneous crystalline substrates on the nucleation of acetaminophen (AAP). Optical and electron microscopy showed that the surface morphologies of single crystal triclinic lactose and D-mannitol differed significantly from their counterparts formed via spherical agglomeration. Spherical agglomerates of lactose were shown to enhance the nucleation rate of acetaminophen (AAP) by a factor of 11 compared to single crystal lactose; however, no such enhancement was observed for D-mannitol. X-ray powder diffraction identified the presence of new crystal faces of lactose present only in the spherical agglomerates However, D-mannitol did not show any significant change in crystal morphology. The new crystal faces of triclinic lactose were analyzed using geometric lattice matching software and molecular dynamics simulations to establish any new and significant epitaxial matches between lactose and AAP. A coincident lattice match and a large favorable energy interaction from hydrogen bonding were observed between the (141¯) and (001) crystal faces of lactose and AAP, respectively. The enhanced nucleation kinetics, X-ray data, and computational studies indicated that the spherical crystallization of lactose exposed the (141¯) face on the surface of the agglomerates, which subsequently enhanced the nucleation rate of AAP through geometric lattice matching and molecular functionality. This study highlights the importance of exploring different heterogeneous substrate morphologies for enhancing nucleation kinetics.

Investigating the Role of Solvent in the Formation of Vacancies on Ibuprofen Crystal Facets.

Surface defects play a crucial role in the process of crystal growth, as incorporation of growth units generally takes place on undercoordinated sites on the growing crystal facet. In this work, we use molecular simulations to obtain information on the role of the solvent in the roughening of three morphologically relevant crystal faces of form I of racemic ibuprofen. To this aim, we devise a computational strategy to evaluate the energetic cost associated with the formation of a surface vacancy for a set of ten solvents, covering a range of polarities and hydrogen bonding propensities. We find that the mechanism as well as the work of defect formation are markedly solvent and facet dependent. Based on Mean Force Integration and Well Tempered Metadynamics, the methodology developed in this work has been designed with the aim of capturing solvent effects at the atomistic scale while maintaining the computational efficiency necessary for implementation in high-throughput in-silico screenings of crystallization solvents.

Crystal structures of anhydrous and hydrated ceftibuten.

Ceftibuten, C15H14N4O6S2, with the systematic name (6R,7R)-7-{[(Z)-2-(2-amino-1,3-thia-zol-4-yl)-4-carb-oxy-but-2-eno-yl]amino}-8-oxo-5-thia-1-aza-bicyclo-[4.2.0]oct-2-ene-2-carb-oxy-lic acid, is a third generation, orally administered cephalosporin anti-biotic with broad anti-microbial activity and stability against extended spectrum ß-lactamases. Ceftibuten can exist in various hydration states and to better understand the location of the water mol-ecules of crystallization and their effect on the structure, the crystal structures of anhydrous (I) and hydrated (II) ceftibuten were determined and both occur as zwitterions with proton transfer from the carboxyl-ate group adjacent to the ß-lactam ring to the N atom of the thia-zole ring. The ß-lactam ring in (I) is almost planar but the equivalent grouping in (II) is slightly buckled. In the extended structure of (I), O-H⋯O and N-H⋯O hydrogen bonds link the mol-ecules into a three-dimensional network. In (II), O-H⋯Oc, N-H⋯Oc, O-H⋯Ow, N-H⋯Ow and Ow-H⋯Ow (c = ceftibuten, w = water) hydrogen bonds link the components into a three-dimensional network. A large void space is present within the anhydrous crystal structure that can accommodate between two and three mol-ecules of water.

Binary isobaric phase diagrams of stearyl alcohol-poloxamer 407 formulations in the molten and solid state.

Multiparticulate formulations allow for the design of specialized pharmaceutical dosage forms that cater to the needs of a wide range of patient demographics, such as pediatric and geriatric populations, by affording control over the release rate and facilitating the formulation of fixed-dose combination drugs. Melt spray-congealing (MSC) is a method for preparing multiparticulate dosage forms from a suspension or solid solution of active pharamaceutical ingredients (API) and a molten carrier matrix. Stearyl alcohol and poloxamer 407 mixtures are widely used as carrier matrices in MSC microsphere formulations. In this report, the phase equilibria of stearyl alcohol-poloxamer 407 mixtures were investigated by generating binary phase diagrams of composition, i.e. weight/weight percent of poloxamer 407 in stearyl alcohol, and temperature in the molten form and the solid state. The phase equilibria of the molten state were characterized by 1H NMR measurements. The miscibility curves of stearyl alcohol-poloxamer 407 molten mixtures revealed that stearyl alcohol and poloxamer 407 are not miscible in all proportions and that miscibility substantially increases with temperature. The phase equilibria of the solid state were characterized by DSC and PXRD experiments. The phase diagrams of the solid state indicate that stearyl alcohol and poloxamer 407 crystallize and melt separately and, thus, do not form a eutectic or a single phase. The phases equilibria of the bulk mixtures were compared to the phases observed in placebo MSC microspheres and it was determined that the microspheres consist of a mixture of thermodynamically stable and metastable stearyl alcohol crystals immediately after manufacture.

Hydrogen abstraction by chlorine atom from small organic molecules containing amino acid functionalities: an assessment of theoretical procedures.

A high-level quantum chemistry investigation has been carried out for the abstraction by chlorine atom of hydrogen from methane and five monosubstituted methanes, chosen to reflect the chemical functionalities contained in amino acids and peptides. A modified W1' procedure is used to calculate benchmark barriers and reaction energies for the six reactions. The reactions demonstrate a broad range of barrier heights and reaction energies, which can be rationalized using curve-crossing and molecular orbital theory models. In addition, the performance of a range of computationally less demanding electronic structure methods is assessed for calculating the energy profiles for the six reactions. It is found that the G3X(MP2)-RAD procedure compares best with the W1' benchmark, demonstrating a mean absolute deviation (MAD) from W1' of 2.1 kJ mol(-1). The more economical RMP2/G3XLarge and UB2-PLYP/G3XLarge methods are also shown to perform well, with MADs from W1' of 2.9 and 3.0 kJ mol(-1), respectively.

Dynamics and Thermodynamics of Ibuprofen Conformational Isomerism at the Crystal/Solution Interface.

Conformational flexibility of molecules involved in crystal growth and dissolution is rarely investigated in detail and usually considered to be negligible in the formulation of mesoscopic models of crystal growth. In this work, we set out to investigate the conformational isomerism of ibuprofen as it approaches and is incorporated in the morphologically dominant {100} crystal face, in a range of different solvents: water, 1-butanol, toluene, cyclohexanone, cyclohexane, acetonitrile, and trichloromethane. To this end, we combine extensive molecular dynamics and well-tempered metadynamics simulations to estimate the equilibrium distribution of conformers, compute conformer-conformer transition rates, and extract the characteristic relaxation time of the conformer population in solution, adsorbed at the solid/liquid interface, incorporated in the crystal in contact with the mother solution, and in the crystal bulk. We find that, while the conformational equilibrium distribution is weakly dependent on the solvent, relaxation times are instead significantly affected by it. Furthermore, differences in the relaxation dynamics are enhanced on the crystal surface, where conformational transitions become slower and specific conformational transition pathways are hindered. This leads us to observe that the dominant mechanisms of conformational transition can also change significantly moving from the bulk solution to the crystal interface, even for a small molecule with limited conformational flexibility such as ibuprofen. Our findings suggest that understanding conformational flexibility is key to provide an accurate description of the solid/liquid interface during crystal dissolution and growth, and therefore, its relevance should be systematically assessed in the formulation of mesoscopic growth models.

Secondary Structure Assignment of Amyloid-ß Peptide Using Chemical Shifts.

The distinct conformational dependence of chemical shifts caused by α-helices and ß-sheets renders NMR chemical shift analysis a powerful tool for the structural determination of proteins. However, the time scale of NMR experiments can make a secondary structure assignment of highly flexible peptides or proteins, which may be converting between conformational substates, problematic. For instance the amyloid-ß monomer, according to NMR chemical shifts, adopts a predominately random coil structure in aqueous solution (with <3% α-helical content). Molecular dynamics simulations, on the other hand, suggest that α-helical content can be significant (10-25%). In this paper, we explore the possible reasons for this discrepancy and show that the different results from experiments and theory are not necessarily mutually exclusive but may reflect a general problem of secondary structure assignment of conformationally flexible biomolecules.

Nature of Glycine and Its α-Carbon Radical in Aqueous Solution: A Theoretical Investigation.

Quantum chemistry calculations and classical molecular dynamics simulations have been used to examine the equilibria in solution between the neutral and zwitterionic forms of glycine and also of the glycyl radical. The established preference (by 30 kJ mol(-1)) for the zwitterion of glycine was confirmed by both the quantum chemical calculations and the classical molecular dynamics simulations. The best agreement with experiment was derived from thermodynamic integration calculations of explicitly solvated systems, which gives a free energy difference of 36.6 ± 0.6 kJ mol(-1). In contrast, for the glycyl radical in solution, the neutral form is preferred, with a calculated free energy difference of 54.8 ± 0.6 kJ mol(-1). A detailed analysis of the microsolvation environments of each species was carried out by evaluating radial distribution functions and hydrogen bonding patterns. This analysis provides evidence that the change in preference between glycine and glycyl radical is due to the inherent gas-phase stability of the neutral α-carbon radical rather than to any significant difference in the solvation behavior of the constituent species.