Clare McCabe

The accurate modeling of metalloproteins remains a central challenge in computational chemistry, particularly when classical force fields fail to describe metal–ligand interactions with significant covalent character. Herein, we present a diagnostic methodology to delimit quantum mechanical regions in hybrid quantum mechanics/molecular mechanics (QM/MM) simulations, considering (i) the hydration of F^– and (ii) the zinc finger (ZF) domain of the Nuclear Protein Localization 4 (NPL4) as case studies. Our approach is grounded in methods of wave function analysis within quantum chemical topology. Specifically, we rely on the quantum theory of atoms in molecules and the interacting quantum atoms approaches to partition interaction energies between (i) the F^– anions and solvating water molecules and (ii) the Zn²⁺ and its surrounding amino acids within the ZF of NPL4 into classical (ionic) and exchange–correlation (covalent) components. The presented methodology provides also a rational criterion to determine either the inclusion or exclusion of water molecules within the QM region as Zn²⁺···OH₂ interactions might be relevant as indicated by electronic structure geometry optimizations. Our diagnostic offers a systematic way to define suitable QM/MM boundaries, and it is useful to decide on the inclusion of solvent or other species in the QM region of hybrid simulations when covalent contributions are uncertain or only partially present. This analysis reveals that (i) the first solvation layer of F^– and (ii) the first coordination shell of Zn²⁺ exhibit significant covalency in their interaction with the F^– and Zn²⁺ ions, respectively, justifying a minimal QM region comprising three cysteines, one histidine, and the metal center for the ZF in NPL4. The QM/MM hybrid simulations for the hydration of the F^– anion showed important differences with respect to classical molecular mechanics concerning the coordination number and maxima of F^–···O peaks in the radial distribution functions. Hybrid QM/MM simulations using this covalency-based region in the ZF of NPL4 recover accurate coordination geometries and bond lengths, correcting important distortions introduced by classical force fields, including spurious hydration and hypercoordination of the Zn²⁺ center. Both the enlargement and reduction of the QM region result in incorrect coordination geometries and less accurate Zn²⁺–ligand distance distributions. The diagnostic method put forward herein shows that metrics based on the chemical bonding scenario can inform the construction of transferable QM regions for the solvation of ions and metalloproteins. Overall, our results establish a general protocol to integrate wave function analysis with hybrid modeling, and they highlight the potential of such analysis to improve the accuracy of QM/MM simulations.

The stratum corneum (SC), the outermost layer of human skin, owes its barrier function to highly ordered lipid lamellae composed primarily of ceramides (CERs), cholesterol (CHOL), and free fatty acids (FFAs). Molecular dynamics simulations offer the opportunity to gain valuable insights into the structural organization of lipids, complementing experimental approaches. However, simulations using atomistic models are computationally expensive when studying the large, multilamellar structures characteristic of the SC. Coarse-grained (CG) models of SC lipids provide an efficient alternative but have largely been limited to CER NS. In this work, a previously developed Multistate Iterative Boltzmann Inversion (MS-IBI) CG model for CER NS has been extended to three additional CER subclassesCERs NP, AP, and ASwhich differ from CER NS in headgroup hydroxylation. By leveraging structural similarity and transferring nonbonded interaction parameters for CER NS hydroxyl groups, we have developed new models with minimal reparameterization. The models have been validated against atomistic simulations of both pure and mixed bilayers containing CHOL and FFA. To capture the multilamellar organization, six-leaflet multilayers were self-assembled. The resulting CG systems exhibited lamellar organization, chain order, and repeat distances consistent with the available experimental data. Comparisons across CER subclasses revealed that headgroup hydroxylation influences lipid packing, chain tilt, and whether the CER tails exhibit a hairpin or extended conformation. This work demonstrates the flexibility and transferability of the MS-IBI approach and provides CG models for key CER subclasses in human SC, enabling large-scale simulations of realistic SC lipid compositions beyond the reach of atomistic models.

Amphiphilic monolayers composed of end groups with distinct polar and nonpolar functional groups offer rapid and reversible interfacial adaptation in response to environmental stimuli such as a change in interfacial medium polarity. We have synthesized and characterized a suite of monolayers with functional groups of competing polarity designed to reconfigure their interfacial chemical composition in response to solvent polarity. In these films, the end group is designed to be able to reorient and expose the functional groups that minimize the interfacial free energy between the film and the environment. Using a combination of spectroscopic, computational, and wettability studies, we have investigated the responsive interfacial behavior of different end groups upon exposure to environments with varying polarities. Contact angle measurements across a series of polar and dispersive probe liquids reveal trends that reflect the underlying molecular flexibility and composition. Vibrational sum frequency generation (SFG) spectroscopy and atomistic molecular dynamics (MD) simulations confirm solvent-driven reorientation of the end groups, with restructuring observed at the interface. To quantify these effects, we have developed a surface energy calculation model that incorporates solvent-induced surface rearrangements into the estimations. Our findings reveal a strong dependence of surface energy and switching behavior on the length and flexibility of the functionalities in the end group, which affects the exposure of certain interfacial compositions under different solvents. These results offer new insights into the design of adaptive monolayers and provide a framework for evaluating solvent-responsive surfaces.

Shear viscosity is a crucial property for optimising various industrial processes. While experimental data is widely available for common fluids under ambient conditions, viscosity often needs to be estimated at state conditions far from ambient temperature and pressure or for solvent systems lacking experimental measurements. Although empirical correlations can be used to estimate viscosity, molecular simulation offers a molecular-level approach that not only provides accurate viscosity estimates but also yields insight into the underlying molecular mechanisms and can be readily extended to fluid mixtures. Of the available simulation approaches to calculating the viscosity, the non-equilibrium molecular dynamics (NEMD) periodic perturbation method (PPM) has demonstrated both computational efficiency and accuracy in predicting shear viscosity. However, implementing PPM requires significant trial and error in selecting a suitable acceleration amplitude or sampling at multiple acceleration amplitudes to determine an extrapolation scheme, both of which increase computational demands and reduce the appeal of the method. In this work, we demonstrate that Quentrec’s local order theory provides a superior fit for the dependence of viscosity on acceleration amplitude, enabling accurate extrapolation to zero-perturbation viscosity. The method is applied to a diverse set of 144 organic solvents and yields results that show good agreement with both experimental data and values obtained from equilibrium molecular dynamics simulations. We further show that data-driven models can accurately estimate the acceleration amplitude corresponding to a given relative deviation from zero-perturbation viscosity. By specifying a reasonable relative deviation, the estimated acceleration amplitude exhibits lower statistical noise while simultaneously enabling precise reproduction of the computed zero-perturbation viscosity after compensating for the deviation. This approach circumvents the need for sampling multiple acceleration amplitudes and thereby facilitates the implementation of the PPM method.

Functionalized alkylsilane films exhibit switching interfacial behavior in response to different solvents. Image shows exposed oxygen (red) and carbon atoms (green). Hydroxyl films switch more favorably than carboxyl counterparts.

Molecular simulations are increasingly used to predict thermophysical properties and explore molecular-level phenomena beyond modern imaging techniques. To make these tools accessible to nonexperts, several open-source molecular dynamics (MD) and Monte Carlo (MC) codes have been developed. However, using these tools is challenging, and concerns about the validity and reproducibility of the simulation data persist. In 2017, Schappals et al. reported a benchmarking study involving several research groups independently performing MD and MC simulations using different software to predict densities of alkanes using common molecular mechanics force fields [ J. Chem. Theory Comput. 2017, 4270−4280 28738147 10.1021/acs.jctc.7b00489 ]. Although the predicted densities were reasonably close (mostly within 1%), the data often fell outside of the combined statistical uncertainties of the different simulations. Schappals et al. concluded that there are unavoidable errors inherent to molecular simulations once a certain degree of complexity of the system is reached. The Molecular Simulation Design Framework (MoSDeF) is a workflow package designed to achieve TRUE (Transparent, Reproducible, Usable-by-others, and Extensible) simulation studies by standardizing the implementation of molecular models for various simulation engines. This work demonstrates that using MoSDeF to initialize a simulation workflow results in consistent predictions of system density, even while increasing model complexity.

The lipids located in the outermost layer of the skin, the stratum corneum (SC), play a crucial role in maintaining the skin barrier function. The primary components of the SC lipid matrix are ceramides (CERs), cholesterol (CHOL), and free fatty acids (FFAs). They form two crystalline lamellar phases: the long periodicity phase (LPP) and the short periodicity phase (SPP). In inflammatory skin conditions like atopic dermatitis and psoriasis, there are changes in the SC CER composition, such as an increased concentration of a sphingosine-based CER (CER NS) and a reduced concentration of a phytosphingosine-based CER (CER NP). In the present study, a lipid model was created exclusively forming the SPP, to examine whether alterations in the CER NS:CER NP molar ratio would affect the lipid organization. Experimental data were combined with molecular dynamics simulations of lipid models containing CER NS:CER NP at ratios of 1:2 (mimicking a healthy SC ratio) and 2:1 (observed in inflammatory skin diseases), mixed with CHOL and lignoceric acid as the FFA. The experimental findings show that the acyl chains of CER NS and CER NP and the FFA are in close proximity within the SPP unit cell, indicating that CER NS and CER NP adopt a linear conformation, similarly as observed for the LPP. Both the experiments and simulations indicate that the lamellar organization is the same for the two CER NS:CER NP ratios while the SPP NS:NP 1:2 model had a slightly denser hydrogen bonding network than the SPP NS:NP 2:1 model. The simulations show that this might be attributed to intermolecular hydrogen bonding with the additional hydroxide group on the headgroup of CER NP compared with CER NS.

We have developed a multi-input E(n) equivariant graph convolution-based model designed for the prediction of chemical properties that result from the interaction of heterogeneous molecular structures. By incorporating spatial features and constraining the functions learned from these features to be equivariant to E(n) symmetries, the interactional-equivariant graph neural network (IEGNN) can efficiently learn from the 3D structure of multiple molecules. To verify the IEGNN's capability to learn interactional properties, we tested the model's performance on three molecular data sets, two of which are curated in this study and made publicly available for future interactional model benchmarking. To enable the loading of these data sets, an open-source data structure based on the PyTorch Geometric library for batch loading multigraph data points is also created. Finally, the IEGNN's performance on a data set consisting of an unknown interactional relationship (the frictional properties resulting between monolayers with variable composition) is examined. The IEGNN model developed was found to have the lowest mean absolute percent error for the predicted tribological properties of four of the six data sets when compared to previous methods.