Extending the perturbed matrix method beyond the dipolar approximation: comparison of different levels of theory

Conformational heterogeneity and protonation equilibria shape the photocycle branching in channelrhodopsin-2

Channelrhodopsin-2 is a photoactive membrane protein serving as an ion channel, gathering significant interest for its applications in optogenetics. Despite extensive investigation, several aspects of its photocycle remain elusive and continue to be subjects of ongoing debate. Of particular interest are the localization of the P480 intermediate within the photocycle and the timing of the deprotonation of glutamic acid E90, a critical residue for ChR2 functioning. In this study, we explore the possibility of an early-P480 state, formed directly upon photoillumination of the dark-adapted state, where E90 is deprotonated, as hypothesized in a previous work [Kuhne et al. Proc. Natl. Acad. Sci. 116.19 (2019): 9380]. Employing extended molecular dynamics simulations, deprotonation free energy calculations, and the computation of the infrared band associated with E90, we provide support to the photocycle model proposed by Kuhne et al. Furthermore, our findings show that E90 protonation state is influenced by diverse interconnected variables and provide molecular detail insights that connect E90 interaction pattern with its deprotonation propensity. Our data demonstrate in fact that both protonated and deprotonated E90 are possible in P480 depending on E90 hydrogen bonding pattern and explaining the molecular mechanism at the basis of P480 accumulation under continuous illumination.

Modulating Charge Transfer Kinetics along Poly Adenine: Chemical Modifications, Temperature, and Conformational Effects

The charge transfer (CT) reactions in nucleic acids are crucial for genome damage and repair and nanoelectronics using DNA as a molecular conductor. Previous experimental and theoretical works underlined the significance of nucleic acid structural dynamics on CT kinetics, requiring models that incorporate the dynamics of the nucleic acid, solvents, and counterions. Here, we investigated hole transfer kinetics in poly adenine single and double strands at various temperatures and the rate enhancement due to adenine-to-7-deazaadenine mutation by means of a QM/MM approach. We found that the hole transfer rate in poly adenine double strands increases with temperature while the helix conformation is retained, whereas single strands exhibit the opposite thermal response. Additionally, the positive charge migrates more efficiently in poly-7-deazaadenine double strands. Our results, consistent with experimental data, suggest that a thermally induced hopping model can accurately describe CT kinetics in these sequences. The approach is transferable for studying CT reactions in other nucleic acid strands.

Unveiling the Cleavage Mechanism of an RNA Model Compound on the whole pH Scale: Computations Meet Experiments in the Determination of Reaction Rates

The knowledge of the mechanism of reactions occurring in solution is a primary research line both in the context of theoretical-computational chemistry and in the field of organic and bio-organic chemistry. Given the importance of the hydrolysis of nucleic acids in life-related phenomena, here we present a combined experimental and computational study on the cleavage of an RNA model compound. This phosphodiester features a cleavage rate strictly dependent on the pH with three different dependence domains. Such experimental evidence, highlighted by an in-depth kinetic investigation, unequivocally suggests a change in the reaction mechanism along the pH scale. In order to interpret the data and to explain the experimental behavior, we have applied a theoretical-computational procedure, involving a hybrid quantum/classical approach, able to model chemical reactions in complex environments, i. e. in solution. This study turns out to quantitatively reproduce the experimental data with accuracy and, in addition, provides useful mechanistic insight into the transesterification process of the investigated compound. The study indicates that the cleavage can occur through anA N D N ${A_N D_N }$ , anA N + D N ${A_N + D_N }$ , and aD N A N ${D_N A_N }$ mechanism depending on the pH values.

Unraveling redox pathways of the disulfide bond in dimethyl disulfide: Ab initio modeling

CONTEXT

In cellular environments, the reduction of disulfide bonds is pivotal for protein folding and synthesis. However, the intricate enzymatic mechanisms governing this process remain poorly understood. This study addresses this gap by investigating a disulfide bridge reduction reaction, serving as a model for comprehending electron and proton transfer in biological systems. Six potential mechanisms for reducing the dimethyl disulfide (DMDS) bridge through electron and proton capture were explored. Thermodynamic and kinetic analyses elucidated the sequence of proton and electron addition. MD-PMM, a method that combines molecular dynamics simulations and quantum-chemical calculations, was employed to compute the redox potential of the mechanism. This research provides valuable insights into the mechanisms and redox potentials involved in disulfide bridge reduction within proteins, offering an understanding of phenomena that are challenging to explore experimentally.

METHODS

All calculations used the Gaussian 09 software package at the MP2/6-311 + g(d,p) theory level. Visualization of the molecular orbitals and electron densities was conducted using Gaussview6. Molecular dynamics simulations were performed using GROMACS with the CHARMM36 force field. The PyMM program (Python Program for QM/MM Simulations Based on the Perturbed Matrix Method) is used to apply the Perturbed Matrix Method to MD simulations.

Modeling the temperature dependence of the fluorescence properties of Indole in aqueous solution

In a recent paper, we proposed a scheme to describe the relaxation mechanism of the excited Indole in aqueous solution, involving the fluctuations among the diabatic electronic states 1Lb, 1La and 1πσ∗. Such a theoretical and computational model reproduced accurately the available experimental data at room temperature. Following these results, in the present work, we model the complex temperature dependence of the fluorescence properties of Indole in aqueous solution, with results further validating the proposed relaxation scheme. This scheme is able to explain the temperature effects on the fluorescence behavior indicating the water fluctuations as the main cause of (i) the stabilization of the dark state (1πσ∗) and (ii) the increase in temperature of the kinetics of the irreversible transition towards such a state.

A Theoretical-Computational Study of Phosphodiester Bond Cleavage Kinetics as a Function of the Temperature

The hydrolysis of the phosphodiester bond is an important chemical reaction involved in several biological processes. Here, we study the cleavage of this bond by means of a theoretical-computational method in a model system, the dineopentyl phosphate. By such an approach, we reconstructed the kinetics and related thermodynamics of this chemical reaction along an isochore. In particular, we evaluated the kinetic constants of all the reaction steps within a wide range of temperatures, mostly corresponding to conditions where no experimental measures are available due to the extremely slow kinetics. Our results, in good agreement with the experimental data, show the robustness of our theoretical-computational methodology which can be easily extended to more complex systems.

Molecular Modeling of the Deamidation Reaction in Solution: A Theoretical-Computational Study

In this work, a theoretical-computational method is applied to study the deamidation reaction, a critical post-translational modification in proteins, using a simple model molecule in solution. The method allows one to comprehensively address the environmental effect, thereby enabling one to accurately derive the kinetic rate constants for the three main steps of the deamidation process. The results presented, in rather good agreement with the available experimental data, underline the necessity for a rigorous treatment of environmental factors and a precise kinetic model to correctly assess the overall kinetics of the deamidation reaction.

Benchmarking Aspects of Ab Initio Fragment Models for Accurate Excimer Potential Energy Surfaces

While Coupled-Cluster methods have been proven to provide an accurate description of excited electronic states, the scaling of the computational costs with the system size limits the degree for which these methods can be applied. In this work different aspects of fragment-based approaches are studied on noncovalently bound molecular complexes with interacting chromophores of the fragments, such as π-stacked nucleobases. The interaction of the fragments is considered at two distinct steps. First, the states localized on the fragments are described in the presence of the other fragment(s); for this we test two approaches. One method is founded on QM/MM principles, only including the electrostatic interaction between the fragments in the electronic structure calculation with Pauli repulsion and dispersion effects added separately. The other model, a Projection-based Embedding (PbE) using the Huzinaga equation, includes both electrostatic and Pauli repulsion and only needs to be augmented by dispersion interactions. In both schemes the extended Effective Fragment Potential (EFP2) method of Gordon et al. was found to provide an adequate correction for the missing terms. In the second step, the interaction of the localized chromophores is modeled for a proper description of the excitonic coupling. Here the inclusion of purely electrostatic contributions appears to be sufficient: it is found that the Coulomb part of the coupling provides accurate splitting of the energies of interacting chromophores that are separated by more than 4 Å.

Unveiling the Excited State Dynamics of Indole in Solution

In this paper, we reconstruct in detail the dynamics of the emitting electronic excited state of aqueous indole, investigating its relaxation mechanism and kinetics to be related to the time-dependent fluorescence signal. Taking advantage of the results shown in a very recent paper, we were able to model the relaxation process in solution in terms of the transitions between two gas-phase singlet electronic states (¹L_a and ¹L_b), subsequently irreversibly relaxing to the gas-phase singlet dark state (¹πσ*). A comparison of the results with the available experimental data shows that the relaxation mechanism we obtain by our theoretical-computational model is reliable, reproducing rather accurately all the experimental observables.

Theoretical-Computational Modeling of CD Spectra of Aqueous Monosaccharides by Means of Molecular Dynamics Simulations and Perturbed Matrix Method

The electronic circular dichroism (ECD) spectra of aqueous d-glucose and d-galactose were modeled using a theoretical-computational approach combining molecular dynamics (MD) simulations and perturbed matrix method (PMM) calculations, hereafter termed MD-PMM. The experimental spectra were reproduced with a satisfactory accuracy, confirming the good performances of MD-PMM in modeling different spectral features in complex atomic-molecular systems, as already reported in previous studies. The underlying strategy of the method was to perform a preliminary long timescale MD simulation of the chromophore followed by the extraction of the relevant conformations through essential dynamics analysis. On this (limited) number of relevant conformations, the ECD spectrum was calculated via the PMM approach. This study showed that MD-PMM was able to reproduce the essential features of the ECD spectrum (i.e., the position, the intensity, and the shape of the bands) of d-glucose and d-galactose while avoiding the rather computationally expensive aspects, which were demonstrated to be important for the final outcome, such as (i) the use of a large number of chromophore conformations; (ii) the inclusion of quantum vibronic coupling; and (iii) the inclusion of explicit solvent molecules interacting with the chromophore atoms within the chromophore itself (e.g., via hydrogen bonds).

Modelling Complex Bimolecular Reactions in a Condensed Phase: The Case of Phosphodiester Hydrolysis

(1) Background: the theoretical modelling of reactions occurring in liquid phase is a research line of primary importance both in theoretical-computational chemistry and in the context of organic and biological chemistry. Here we present the modelling of the kinetics of the hydroxide-promoted hydrolysis of phosphoric diesters. (2) Method: the theoretical-computational procedure involves a hybrid quantum/classical approach based on the perturbed matrix method (PMM) in conjunction with molecular mechanics. (3) Results: the presented study reproduces the experimental data both in the rate constants and in the mechanistic aspects (C-O bond vs. O-P bond reactivity). The study suggests that the basic hydrolysis of phosphodiesters occurs through a concerted ANDN mechanism, with no formation of penta-coordinated species as reaction intermediates. (4) Conclusions: the presented approach, despite the approximations, is potentially applicable to a large number of bimolecular transformations in solution and therefore leads the way to a fast and general method to predict the rate constants and reactivities/selectivities in complex environments.

Origins of the pH-Responsive Photoluminescence of Peptide-Functionalized Au Nanoclusters

Ultrasmall peptide-protected gold nanoclusters are a promising class of bioresponsive material exhibiting pH-sensitive photoluminescence. We present a theoretical insight into the effect peptide-ligand environment has on pH-responsive fluorescence, with the aim of enhancing the rational design of gold nanoclusters for bioapplications. Employing a hybrid quantum/classical computational methodology, we systematically calculate deprotonation free energies of N-terminal cysteine amine groups in proximity to the inherently fluorescent core of Au₂₅(Peptide)₁₈ nanoclusters. We find that subtle changes in hexapeptide sequence alter the electrostatic environment and significantly shift the conventional N-terminal amine pK_a expected for amino acids free-in-solution. Our findings provide an insight into how the deprotonation equilibrium of N-terminal amine and side chain carboxyl groups cooperatively respond to solution pH changes, explaining the experimentally observed, yet elusive, pH-responsive fluorescence of peptide-functionalized Au₂₅ clusters.

A Simplified Treatment for Efficiently Modeling the Spectral Signal of Vibronic Transitions: Application to Aqueous Indole

In this paper, we introduce specific approximations to simplify the vibronic treatment in modeling absorption and emission spectra, allowing us to include a huge number of vibronic transitions in the calculations. Implementation of such a simplified vibronic treatment within our general approach for modelling vibronic spectra, based on molecular dynamics simulations and the perturbed matrix method, provided a quantitative reproduction of the absorption and emission spectra of aqueous indole with higher accuracy than the one obtained when using the existing vibronic treatment. Such results, showing the reliability of the approximations employed, indicate that the proposed method can be a very efficient and accurate tool for computational spectroscopy.

Theoretical Evaluation of Sulfur-Based Reactions as a Model for Biological Antioxidant Defense

Sulfur-containing amino acids, Methionine (Met) and Cysteine (Cys), are very susceptible to Reactive Oxygen Species (ROS). Therefore, sulfur-based reactions regulate many biological processes, playing a key role in maintaining cellular redox homeostasis and modulating intracellular signaling cascades. In oxidative conditions, Met acts as a ROS scavenger, through Met sulfoxide formation, while thiol/disulfide interchange reactions take place between Cys residues as a response to many environmental stimuli. In this work, we apply a QM/MM theoretical-computational approach, which combines quantum-mechanical calculations with classical molecular dynamics simulations to estimate the free energy profile for the above-mentioned reactions in solution. The results obtained, in good agreement with experimental data, show the validity of our approach in modeling sulfur-based reactions, enabling us to study these mechanisms in more complex biological systems.

PyMM: An Open-Source Python Program for QM/MM Simulations Based on the Perturbed Matrix Method

Quantum mechanical/molecular mechanics (QM/MM) methods are important tools in molecular modeling as they are able to couple an extended phase space sampling with an accurate description of the electronic properties of the system. Here, we describe a Python software package, called PyMM, which has been developed to apply a QM/MM approach, the perturbed matrix method, in a simple and efficient way. PyMM requires a classical atomic trajectory of the whole system and a set of unperturbed electronic properties of the ground and electronic excited states. The software output includes a set of the most common perturbed properties, such as the electronic excitation energies and the transitions dipole moments, as well as the eigenvectors describing the perturbed electronic states, which can be then used to estimate whatever electronic property. The software is composed of a simple and complete command-line interface, a set of internal input validation, and three main analyses focusing on (i) the perturbed eigenvector behavior, (ii) the calculation of the electronic absorption spectrum, and (iii) the estimation of the free energy differences along a reaction coordinate.

Modeling Charge Transfer Reactions by Hopping between Electronic Ground State Minima: Application to Hole Transfer between DNA Bases

In this paper, we extend the previously described general model for charge transfer reactions, introducing specific changes to treat the hopping between energy minima of the electronic ground state (i.e., transitions between the corresponding vibrational ground states). We applied the theoretical-computational model to the charge transfer reactions in DNA molecules which still represent a challenge for a rational full understanding of their mechanism. Results show that the presented model can provide a valid, relatively simple, approach to quantitatively study such reactions shedding light on several important aspects of the reaction mechanism.

Integration of Quantum Chemistry, Statistical Mechanics, and Artificial Intelligence for Computational Spectroscopy: The UV-Vis Spectrum of TEMPO Radical in Different Solvents

The ongoing integration of quantum chemistry, statistical mechanics, and artificial intelligence is paving the route toward more effective and accurate strategies for the investigation of the spectroscopic properties of medium-to-large size chromophores in condensed phases. In this context we are developing a novel workflow aimed at improving the generality, reliability, and ease of use of the available computational tools. In this paper we report our latest developments with specific reference to unsupervised atomistic simulations employing non periodic boundary conditions (NPBC) followed by clustering of the trajectories employing optimized feature spaces. Next accurate variational computations are performed for a representative point of each cluster, whereas intracluster fluctuations are taken into account by a cheap yet reliable perturbative approach. A number of methodological improvements have been introduced including, e.g., more realistic reaction field effects at the outer boundary of the simulation sphere, automatic definition of the feature space by continuous perception of solute-solvent interactions, full account of polarization and charge transfer in the first solvation shell, and inclusion of vibronic contributions. After its validation, this new approach has been applied to the challenging case of solvatochromic effects on the UV-vis spectra of a prototypical nitroxide radical (TEMPO) in different solvents. The reliability, effectiveness, and robustness of the new platform is demonstrated by the remarkable agreement with experiment of the results obtained through an unsupervised approach characterized by a strongly reduced computational cost as compared to that of conventional quantum mechanics and molecular mechanics models without any accuracy reduction.

Effects of Environmental and Electric Perturbations on the pKa of Thioredoxin Cysteine 35: A Computational Study

Here we present a theoretical-computational study dealing with the evaluation of the pKa of the Cysteine residues in Thioredoxin (TRX) and in its complex with the Thioredoxin-interacting protein (TXNIP). The free energy differences between the anionic and neutral form of the Cysteine 32 and 35 have been evaluated by means of the Perturbed Matrix Method with classical perturbations due to both the environment and an exogenous electric field as provided by Molecular Dynamics (MD) simulations. The evaluation of the free energies allowed us to show that the effect of the perturbing terms is to lower the pKa of Cysteine 32 and Cysteine 35 with respect to the free amino-acid. On the other hand, in the complex TRX-TXNIP, our data show an enhanced stabilization of the neutral reduced form of Cys 35. These results suggest that external electric stimuli higher than 0.02 V/nm can modulate the Cysteine pKa, which can be connected to the tight regulation of the TRX acting as an antioxidant agent.

Computation of Oxidation Potentials of Solvated Nucleobases by Static and Dynamic Multilayer Approaches

The determination of the redox properties of nucleobases is of paramount importance to get insight into the charge-transfer processes in which they are involved, such as those occurring in DNA-inspired biosensors. Although many theoretical and experimental studies have been conducted, the value of the one-electron oxidation potentials of nucleobases is not well-defined. Moreover, the most appropriate theoretical protocol to model the redox properties has not been established yet. In this work, we have implemented and evaluated different static and dynamic approaches to compute the one-electron oxidation potentials of solvated nucleobases. In the static framework, two thermodynamic cycles have been tested to assess their accuracy against the direct determination of oxidation potentials from the adiabatic ionization energies. Then, the introduction of vibrational sampling, the effect of implicit and explicit solvation models, and the application of the Marcus theory have been analyzed through dynamic methods. The results revealed that the static direct determination provides more accurate results than thermodynamic cycles. Moreover, the effect of sampling has not shown to be relevant, and the results are improved within the dynamic framework when the Marcus theory is applied, especially in explicit solvent, with respect to the direct approach. Finally, the presence of different tautomers in water does not affect significantly the one-electron oxidation potentials.

Rationalizing Sequence and Conformational Effects on the Guanine Oxidation in Different DNA Conformations

The effect of the environment on the guanine redox potential is studied by means of a theoretical–computational approach. Our data, in agreement with previous experimental findings, clearly show that the presence of consecutive guanine bases in both single- and double-stranded DNA oligomers lowers their reduction potential. Such an effect is even more marked when a G-rich quadruplex is considered, where the oxidized form of guanine is particularly stabilized. To the best of our knowledge, this is the first computational study reporting on a quantitative estimate of the dependence of the guanine redox potential on sequence and conformational effects in complex DNA molecules, ranging from single-stranded DNA to G-quadruplex.

Evidence for a high pKa of an aspartic acid residue in the active site of CALB by a fully atomistic multiscale approach

Candida antarctica Lipase B (CALB) is a paradigm for the family of lipases. At pH 7, the optimal pH for catalysis, the protonation state of an aspartic acid of the active site (Asp134) could not be conclusively assigned. In fact, the pK_a estimate provided by a widely used computational tool, namely PropKa, that predicts pKa values of ionizable groups in proteins based on the crystallographic structure, is only slightly above 7 (pKa = 7.25). This, along with the lack of an experimental evaluation, makes the assignment of its protonation state at neutral pH challenging. Here, we calculate the pK_a of Asp134 by means of a fully atomistic multiscale computational approach based on classical molecular dynamics (MD) simulation and the perturbed matrix method (PMM), namely the MD-PMM approach. MD-PMM is able to take into account the dynamics of the system and, at the same time, to treat the deprotonation step at the quantum level. The calculations provide a pK_a value of 8.9 ± 1.1, hence suggesting that Asp134 in CALB should be protonated at neutral, and even at slightly basic, pH.

Absorption behavior of doxorubicin hydrochloride in visible region in different environments: a combined experimental and computational study

The experimental absorption measurements in the interval 350-600 nm (Vis), molecular dynamics simulations, quantum-mechanics calculations and an advanced molecular treatment of simulation data are here combined to provide a complete picture of the absorption behavior in the visible portion of the electromagnetic spectrum of the doxorubicin hydrochloride (DX) molecule in different environments. By such an approach, we have shown that it is possible to characterize the effect of the environment on the DX absorption behavior - including the vibronic contributions - as well as to interpret such differences in terms of molecular electronic excited states, which are found to be strongly influenced by the environment.

Theoretical Modeling of Redox Potentials of Biomolecules

The estimation of the redox potentials of biologically relevant systems by means of theoretical-computational approaches still represents a challenge. In fact, the size of these systems typically does not allow a full quantum-mechanical treatment needed to describe electron loss/gain in such a complex environment, where the redox process takes place. Therefore, a number of different theoretical strategies have been developed so far to make the calculation of the redox free energy feasible with current computational resources. In this review, we provide a survey of such theoretical-computational approaches used in this context, highlighting their physical principles and discussing their advantages and limitations. Several examples of these approaches applied to the estimation of the redox potentials of both proteins and nucleic acids are described and critically discussed. Finally, general considerations on the most promising strategies are reported.

Stationary and Time-Dependent Carbon Monoxide Stretching Mode Features in Carboxy Myoglobin: A Theoretical-Computational Reappraisal

The stationary and time-dependent infrared spectrum (IR) of the CO stretching mode (ν_CO) in carboxymyoglobin (MbCO), a longstanding problem of biophysical chemistry, has been modeled through a theoretical-computational method specifically designed for simulating quantum observables in complex atomic-molecular systems and based on a combined application of long time scale molecular dynamics simulations and quantum-chemical calculations. This study is basically focused on two aspects: (i) the origin of the stationary IR substates (termed as A₀, A₁, and A₃) and (ii) the modeling and the interpretation of the ν_CO energy relaxation. The results, strengthened by a more than satisfactory agreement with the experimental data, concisely indicate that (i) the conformational His64-FeCO relevant substates, i.e., characterized by the formation-disruption of the H-bond between the above moieties, are the main responsible of the presence of two distinct and well separated (A₀ and A₁/A₃) spectroscopic regions; (ii) the characteristic bimodal shape of the A₁/A₃ spectral region, according to our model, is the result of the fluctuation of the electric field pattern as provided by the protein-solvent framework perturbing the bound His64-CO-Heme complex; and (iii) the electric field pattern, in conjunction with the relatively high density of MbCO vibrational states, is also the main determinant of the ν_CO energy relaxation, characterizing its kinetic efficiency.

Segregation on the nanoscale coupled to liquid water polyamorphism in supercooled aqueous ionic-liquid solution

The most intriguing hypothesis explaining many water anomalies is a metastable liquid-liquid phase transition (LLPT) at high pressure and low temperatures, experimentally hidden by homogeneous nucleation. Recent infrared spectroscopic experiments showed that upon addition of hydrazinium trifluoroacetate to water, the supercooled ionic solution undergoes a sharp, reversible LLPT at ambient pressure, possible offspring of that in pure water. Here, we calculate the temperature-dependent signature of the OH-stretching band, reporting on the low/high density phase of water, in neat water and in the same experimentally investigated ionic solution. The comparison between the infrared signature of the pure liquid and that of the ionic solution can be achieved only computationally, providing insight into the nature of the experimentally observed phase transition and allowing us to investigate the effects of ionic compounds on the high to low density supercooled liquid water transition. We show that the experimentally observed crossover behavior in the ionic solution can be reproduced only if the phase transition between the low- and high-density liquid states of water is coupled to a mixing-unmixing transition between the water component and the ions: at low temperatures, water and ions are separated and the water component is a low density liquid. At high temperatures, water and ions get mixed and the water component is a high-density liquid. The separation at low temperatures into ion-rich and ion-poor regions allows unveiling the polyamorphic nature of liquid water, leading to a crossover behavior resembling that observed in supercooled neat water under high pressure.

Integration of theory, simulation, artificial intelligence and virtual reality: a four-pillar approach for reconciling accuracy and interpretability in computational spectroscopy

The established pillars of computational spectroscopy are theory and computer based simulations. Recently, artificial intelligence and virtual reality are becoming the third and fourth pillars of an integrated strategy for the investigation of complex phenomena. The main goal of the present contribution is the description of some new perspectives for computational spectroscopy, in the framework of a strategy in which computational methodologies at the state of the art, high-performance computing, artificial intelligence and virtual reality tools are integrated with the aim of improving research throughput and achieving goals otherwise not possible. Some of the key tools (e.g., continuous molecular perception model and virtual multifrequency spectrometer) and theoretical developments (e.g., non-periodic boundaries, joint variational-perturbative models) are shortly sketched and their application illustrated by means of representative case studies taken from recent work by the authors. Some of the results presented are already well beyond the state of the art in the field of computational spectroscopy, thereby also providing a proof of concept for other research fields.

A computational insight into the relationship between side chain IR line shapes and local environment in fibril-like structures

Infrared spectroscopy is a widely used technique to characterize protein structures and protein mediated processes. While the amide I band provides information on proteins' secondary structure, amino acid side chains are used as infrared probes for the investigation of protein reactions and local properties. In this paper, we use a hybrid quantum mechanical/classical molecular dynamical approach based on the perturbed matrix method to compute the infrared band due to the C=O stretching mode of amide-containing side chains. We calculate, at first, the infrared band of zwitterionic glutamine in water and obtain results in very good agreement with the experimental data. Then, we compute the signal arising from glutamine side chains in a microcrystal of the yeast prion Sup35-derived peptide, GNNQQNY, with a fibrillar structure. The infrared bands obtained by selective isotopic labeling of the two glutamine residues, Q4 and Q5, of each peptide were experimentally used to investigate the local hydration in the fibrillar microcrystal. The experimental spectra of the two glutamine residues, which experience different hydration environments, feature different spectral signals that are well reproduced by the corresponding calculated spectra. In addition, the analysis of the simulated spectra clarifies the molecular origin of the experimentally observed spectroscopic differences that arise from the different local electric field experienced by the two glutamine residues, which is, in turn, determined by a different hydrogen bonding pattern.

Tuning Proton Transfer Thermodynamics in SARS-CoV-2 Main Protease: Implications for Catalysis and Inhibitor Design

The catalytic reaction in SARS-CoV-2 main protease is activated by a proton transfer (PT) from Cys145 to His41. The same PT is likely also required for the covalent binding of some inhibitors. Here we use a multiscale computational approach to investigate the PT thermodynamics in the apo enzyme and in complex with two potent inhibitors, N3 and the α-ketoamide 13b. We show that with the inhibitors the free energy cost to reach the charge-separated state of the active-site dyad is lower, with N3 inducing the most significant reduction. We also show that a few key sites (including specific water molecules) significantly enhance or reduce the thermodynamic feasibility of the PT reaction, with selective desolvation of the active site playing a crucial role. The approach presented is a cost-effective procedure to identify the enzyme regions that control the activation of the catalytic reaction and is thus also useful to guide the design of inhibitors.

Theoretical Characterization of the Reduction Potentials of Nucleic Acids in Solution

Here, we present the theoretical-computational modeling of the oxidation properties of four DNA nucleosides and nucleotides and a set of dinucleotides in solutions. Our estimates of the vertical ionization energies and reduction potentials, close to the corresponding experimental data, show that an accurate calculation of the molecular electronic properties in solutions requires a proper treatment of the effect of the environment. In particular, we found that the effect of the environment is to stabilize the oxidized state of the nucleobases resulting in a remarkable reduction-up to 6.6 eV-of the energy with respect to the gas phase. Our estimates of the aqueous and gas-phase vertical ionization energies, in good agreement with photoelectron spectroscopy experiments, also show that the effect on the reduction potential of the phosphate group and of the additional nucleotide in dinucleotides is rather limited.

Cooperative protein-solvent tuning of proton transfer energetics: carbonic anhydrase as a case study

We investigate the coupling between the proton transfer (PT) energetics and the protein-solvent dynamics using the intra-molecular PT in wild type (wt) human carbonic anhydrase II and its ten-fold faster mutant Y7F/N67Q as a test case. We calculate the energy variation upon PT, and from that we also calculate the PT reaction free energy, making use of a hybrid quantum mechanics/molecular dynamics approach. In agreement with the experimental data, we obtain that the reaction free energy is basically the same in the two systems. Yet, we show that the instantaneous PT energy is on average lower in the mutant possibly contributing to the faster PT rate. Analysis of the contribution to the PT energetics of the solvent and of each protein residue, also not in the vicinity of the active site, provides evidence for electrostatic tuning of the PT energy arising from the combined effect of the solvent and the protein environment. These findings open up a way to the more general task of the rational design of mutants with either enhanced or reduced PT rate.

Tuning Proton Transfer Thermodynamics in SARS-Cov-2 Main Protease: Implications for Catalysis and Inhibitor Design

In this comutational work a hybrid quantum mechanics/molecular mechanics approach, the MD-PMM approach, is used to investigate the proton transfer reaction the activates the catalytic activity of SARS-CoV-2 main protease. The proton transfer thermodynamics is investigated for the apo ensyme (i.e., without any bound substrate or inhibitor) and in the presence of a inhibitor, N3, which was previously shown to covalently bind SARS-CoV-2 main protease.

Exploiting Machine Learning to Efficiently Predict Multidimensional Optical Spectra in Complex Environments

The excited-state dynamics of chromophores in complex environments determine a range of vital biological and energy capture processes. Time-resolved, multidimensional optical spectroscopies provide a key tool to investigate these processes. Although theory has the potential to decode these spectra in terms of the electronic and atomistic dynamics, the need for large numbers of excited-state electronic structure calculations severely limits first-principles predictions of multidimensional optical spectra for chromophores in the condensed phase. Here, we leverage the locality of chromophore excitations to develop machine learning models to predict the excited-state energy gap of chromophores in complex environments for efficiently constructing linear and multidimensional optical spectra. By analyzing the performance of these models, which span a hierarchy of physical approximations, across a range of chromophore-environment interaction strengths, we provide strategies for the construction of machine learning models that greatly accelerate the calculation of multidimensional optical spectra from first principles.

Computational Spectroscopy in Solution by Integration of Variational and Perturbative Approaches on Top of Clusterized Molecular Dynamics

Multiscale QM/MM approaches have become the most suitable and effective methods for the investigation of spectroscopic properties of medium- or large-size chromophores in condensed phases. On these grounds, we are developing a novel workflow aimed at improving the generality, reliability, and ease of use of the available tools. In the present paper, we report the latest developments of such an approach with specific reference to a general workplan starting with the addition of acetonitrile to the panel of solvents already available in the General Liquid Optimized Boundary (GLOB) model enforcing nonperiodic boundary conditions (NPBC). Next, the solvatochromic shifts induced by acetonitrile on both rigid (uracil and thymine) and flexible (thyrosine) chromophores have been studied introducing in our software a number of new features ranging from rigid-geometry NPBC molecular dynamics based on the quaternion formalism to a full integration of variational (ONIOM) and perturbative (perturbed matrix method (PMM)) approaches for describing different solute–solvent topologies and local fluctuations, respectively. Finally, thymine and uracil have been studied also in methanol to point out the generality of the computational strategy. While further developments are surely needed, the strengths of our integrated approach even in its present version are demonstrated by the accuracy of the results obtained by an unsupervised approach and coupled to a computational cost strongly reduced with respect to that of conventional QM/MM models without any appreciable accuracy deterioration.

CPL Spectra of Camphor Derivatives in Solution by an Integrated QM/MD Approach

We extend a recently proposed computational strategy for the simulation of absorption spectra of semi-rigid molecular systems in condensed phases to the emission spectra of flexible chromophores. As a case study, we have chosen the CPL spectrum of camphor in methanol solution, which shows a well-defined bisignate shape. The first step of our approach is the quantum mechanical computation of reference spectra including vibrational averaging effects and taking bulk solvent effects into account by means of the polarizable continuum model. In the present case, the large amplitude inversion mode is explicitly treated by a numerical approach, whereas the other small-amplitude vibrational modes are taken into account within the harmonic approximation. Next, the snapshots of classical molecular dynamics computations are clusterized and one representative configuration from each cluster is used to compute a reference spectrum. In the present case, different clusters correspond to the two stable conformers of camphor in the S1 excited electronic state and, for each of them, to different numbers of strong solute-solvent hydrogen bonds. Finally, local fluctuation effects within each cluster are taken into account by means of the perturbed matrix model. The overall procedure leads to good agreement with experiment for absorption and emission spectra together with their chiral counterparts, thus allowing to analyze the role of different effects (stereo-electronic, vibrational, environmental) in tuning the overall experimental spectra.

Fully Atomistic Multiscale Approach for pKa Prediction

The ionization state of titratable amino acids strongly affects proteins structure and functioning in a large number of biological processes. It is therefore essential to be able to characterize the pK_a of ionizable groups inside proteins and to understand its microscopic determinants in order to gain insights into many functional properties of proteins. A big effort has been devoted to the development of theoretical approaches for the prediction of deprotonation free energies, yet the accurate theoretical/computational calculation of pK_a values is recognized as a current challenge. A methodology based on a hybrid quantum/classical approach is here proposed for the computation of deprotonation free energies. The method is applied to calculate the pK_a of formic acid, methylammonium, and methanethiol, providing results in good agreement with the corresponding experimental estimates. The pK_a is also calculated for aspartic acid and lysine as single residues in solution and for three aspartic/glutamic acids inside a well-characterized protein: hen egg white lysozyme. While for small molecules the method is able to deal with multiple protonation states of all titratable groups, this becomes computationally very expensive for proteins. The calculated pK_a values for the single amino acids (except for the zwitterionic aspartic acid) and inside the protein display a systematic shift with respect to the experimental values that suggests that the fine balance between hydrophobic and polar interactions might be not accurately reproduced by the usual classical force-fields, thus affecting the computation of deprotonation free energies. The calculated pK_a shifts inside the protein are in good agreement with the corresponding experimental ones (within 1 pK_a unit), well reproducing the pK_a changes due to the protein environment even in the case of large pK_a shifts.

The challenge of non-covalent interactions: theory meets experiment for reconciling accuracy and interpretation

In the past decade, many gas-phase spectroscopic investigations have focused on the understanding of the nature of weak interactions in model systems. Despite the fact that non-covalent interactions play a key role in several biological and technological processes, their characterization and interpretation are still far from being satisfactory. In this connection, integrated experimental and computational investigations can play an invaluable role. Indeed, a number of different issues relevant to unraveling the properties of bulk or solvated systems can be addressed from experimental investigations on molecular complexes. Focusing on the interaction of biological model systems with solvent molecules (e.g., water), since the hydration of the biomolecules controls their structure and mechanism of action, the study of the molecular properties of hydrated systems containing a limited number of water molecules (microsolvation) is the basis for understanding the solvation process and how structure and reactivity vary from gas phase to solution. Although hydrogen bonding is probably the most widespread interaction in nature, other emerging classes, such as halogen, chalcogen and pnicogen interactions, have attracted much attention because of the role they play in different fields. Their understanding requires, first of all, the characterization of the directionality, strength, and nature of such interactions as well as a comprehensive analysis of their competition with other non-covalent bonds. In this review, it is shown how state-of-the-art quantum-chemical computations combined with rotational spectroscopy allow for fully characterizing intermolecular interactions taking place in molecular complexes from both structural and energetic points of view. The transition from bi-molecular complex to microsolvation and then to condensed phase is shortly addressed.

Ergodicity Breaking in Thermal Biological Electron Transfer? Cytochrome C Revisited II

It was recently suggested that cytochrome c operates in an ergodicity-breaking regime characterized by unusually large energy gap thermal fluctuations and associated reorganization free energies for heme oxidation of up to 3.0 eV. The large fluctuations were reported to lower activation free energy for oxidation of the heme cofactor by almost a factor of 2 compared to the case where ergodicity is maintained. Our group has recently investigated this claim computationally at several levels of theory and found no evidence for such large energy gap fluctuations. Here we address the points of our earlier work that have raised criticism and we also extend our previous investigation by considering a simple linear polarizability model for cytochrome c oxidation. We find very consistent results among all our computational approaches, ranging from classical molecular dynamics, to the linear polarizability model to QM(PMM)/MM to full QM(DFT)/MM electrostatic emdedding. None of them support the notion of very large energy gap fluctuations or ergodicity breaking. The deviation between our simulations and the ones reported in [ J. Phys. Chem. B 2017, 121, 4958] is traced back to rather large electric fields at the Fe site of the heme c cofactor in that study, not seen in our simulations, neither with the AMBER nor with the CHARMM force field. While ergodicity breaking effects may well occur in other biological ET, our numerical evidence suggests that this is not the case for cytochrome c.

The ONIOM/PMM Model for Effective Yet Accurate Simulation of Optical and Chiroptical Spectra in Solution: Camphorquinone in Methanol as a Case Study

This paper deals with the development and first validation of a composite approach for the simulation of chiroptical spectra in solution aimed to strongly reduce the number of full QM computations without any significant accuracy loss. The approach starts from the quantum mechanical computation of reference spectra including vibrational averaging effects and taking average solvent effects into account by means of the polarizable continuum model. Next, the snapshots of classical molecular dynamics computations are clusterized and one reference configuration from each cluster is used to compute a reference spectrum. Local fluctuation effects within each cluster are then taken into account by means of the perturbed matrix model. The performance of the proposed approach is tested on the challenging case of the optical and chiroptical spectra of camphorquinone in methanol solution. Although further validations are surely needed, the results of this first study are quite promising also taking into account that agreement with experimental data is reached by just a couple of full quantum mechanical geometry optimizations and frequency computations.

Modeling amino-acid side chain infrared spectra: the case of carboxylic residues

Infrared (IR) spectroscopy is commonly utilized for the investigation of protein structures and protein-mediated processes.

Modeling multidimensional spectral lineshapes from first principles: application to water-solvated adenine

We theoretically describe spectral lineshape from first principles, providing insight into solvent–solute interactions in terms of static and dynamic disorder and how these shape experimental signals in linear and non-linear optical spectroscopies.

Modelling vibrational relaxation in complex molecular systems

In this paper we show how it is possible to treat the quantum vibrational relaxation of a chromophore, embedded in a complex atomic-molecular environment, via the explicit solution of the time-dependent Schroedinger equation once using a proper separation between quantum and semiclassical degrees of freedom. The rigorous theoretical framework derived, based on first principles and making use of well defined approximations/assumptions, is utilized to construct a general model for the kinetics of the vibrational relaxation as obtained by the direct evaluation of the density matrix for all the relevant quantum state transitions. Application to (deuterated) N-methylacetamide (the typical benchmark used as a model for the amino acids) shows that the obtained theoretical-computational approach captures the essential features of the experimental process, unveiling the basic relaxation mechanism involving several vibrational state transitions.

Ergodicity-Breaking in Thermal Biological Electron Transfer? Cytochrome C Revisited

It was recently suggested that certain redox proteins operate in an ergodicity-breaking regime to facilitate biological electron transfer (ET). A signature for this is a large variance reorganization free energy (several electronvolts) but a significantly smaller Stokes reorganization free energy due to incomplete protein relaxation on the time scale of the ET event. Here we investigate whether this picture holds for oxidation of cytochrome c in aqueous solution, at various levels of theory including classical molecular dynamics with two additive and one electronically polarizable force field, and QM/MM calculations with the QM region treated by full electrostatic DFT embedding and by the perturbed matrix method. Sampling the protein and energy gap dynamics over more than 250 ns, we find no evidence for ergodicity-breaking effects. In particular, the inclusion of electronic polarizability of the heme group at QM/MM levels did not induce nonergodic effects, contrary to previous reports by Matyushov et al. The well-known problem of overestimation of reorganization free energies with additive force fields is cured when the protein and solvent are treated as electronically polarizable. Ergodicity-breaking effects may occur in other redox proteins, and our results suggest that long simulations, ideally on the ET time scale, with electronically polarizable force fields are required to obtain strong numerical evidence for them.