Deborah A.
Drost
and
Christian
Merten
*
Ruhr University Bochum, Faculty of Chemistry and Biochemistry, Organic Chemistry II, Universitätsstraße 150, 44801 Bochum, Germany. E-mail: www.mertenlab.de; christian.merten@ruhr-uni-bochum.de
First published on 4th June 2024
Recording VCD spectra of aqueous solution poses a particular challenge as water is a strong infrared absorber. Likewise, the computational analysis of VCD spectra by means of DFT-based spectral calculations requires the consideration of explicit solvent molecules, thus posing an even greater challenge. Several studies suggested that by modeling the solvent environment with a few water molecules in a micro-solvation approach would be sufficient to describe experimental spectra. For example, using proline at different pH values, we herein show that a change in the relative spatial orientation of a single water molecule in five-fold solvated structures strongly affects the computed VCD spectral signatures and that Boltzmann-weighted spectra do not correctly reproduce the experiment. We thus explored an approach based on molecular dynamics and subsequent DFT-calculations, in which we considered 30 water molecules (about 1.5 solvation shells). Once again, it was found that the Boltzmann-weighted spectra obtained on the basis of several hundred structures did not correctly reproduce experimental signatures, and a simple averaging scheme resulted in well-matching spectra with comparable bandwidths. The rationale behind the procedure was that sampling the configurational space of the solvent molecules is as equally important as the conformational sampling of the solute. For conformationally more flexible molecules, it is assumed that a much larger set of structures will have to be computed in order to properly sample the conformational space of both solute and solvent.
Recording VCD spectra for samples in aqueous solution, however, poses a particular challenge. Vibrational modes of water strongly overlap with the fingerprint region,1 thus requiring very short pathlengths of about 6–8 μm and high solute concentrations. Heavy water (D2O) allows for slightly longer pathlengths and thus lower concentrations, as the bending mode occurs below 1250 cm−1 instead of in the range 1700–1500 cm−1. However, using D2O instead of H2O also results in H/D exchange at acidic positions and consequently a shift of X–H related vibrational modes out of the spectral window of the measurement. Probing hydrogen bonding networks involving solvent molecules thus becomes challenging and must be done indirectly through the vibrational modes remaining accessible/visible.
Analysing VCD spectra recorded for aqueous solutions also requires huge computational efforts. While solute–solvent interactions with organic solvents occasionally need to be considered explicitly in the computational analysis, their modelling is somewhat straightforward as ACN-d3 and DMSO-d6 only act as hydrogen bond acceptors.2 Even methanol, which can act as both a donor and an acceptor, can easily be modelled as its interactions are typically weak. Water, however, is a strong donor and acceptor, forming tightly bound solvation shells around solutes. As shown for small organic molecules, it can also insert into intramolecular hydrogen bonds and give rise to induced VCD phenomena.3–7a
First, VCD studies in aqueous solution were mostly focussed on biopolymers and targeted qualitative structure–spectra relationships.8–12 Later, combined experimental and computational studies on amino acids in aqueous solution (D2O) and under varying pH conditions were reported. Among the investigated systems were serine,13 leucine,14 cysteine,15 and N-acetyl-cysteine16 for which the manual or molecular dynamics-guided placement of 4–6 D2O molecules was found to sufficiently resemble the key spectral features.
The apparent need to consider several water molecules and to sample their configurational space around the solute stimulated the development of methods to compute VCD spectra based on molecular dynamics (MD) approaches. A numerical approach to obtain VCD spectra directly from semi-empirical QM/MM MD simulations was tested for (Ala)2,17 Ac-Pro-NH218,19 and Ac-Ala-NHCH320 in D2O. With the implementation of nuclear velocity perturbation theory into a plane-wave electronic structure code,21,22 VCD spectral calculations directly from the trajectories of ab initio molecular dynamics (AIMD) simulations became possible.23,24 As an alternative approach for computing VCD spectra directly from MD trajectories, it is also possible to post-process snapshots on a static DFT level.25,26 This has been demonstrated for amphetamine and methamphetamine in D2O, for which clusters with an average of 18 water molecules were extracted from the trajectory, categorized by conformers and subsequently fitted against the experimental spectra.27 For this fitting procedure, the water-based molecular property tensors were set to zero, so that all vibrational modes of the solvent were suppressed. Using the quantum-cluster growth method, solvation shells could also be built up systematically.7,28,29
In this contribution, we investigate the IR and VCD spectra of the amino acid proline (Pro) in H2O at different pH values. Despite its unique structural properties compared to other amino acids, the experimental VCD spectra of Pro were only recorded in small regions of the fingerprint range.30 Theoretical studies on Pro-H2O clusters have been reported, but no comparison with the experiment has been established.31 Interestingly, also Raman optical activity (ROA) data on proline in water were not complemented with computational studies going beyond continuum solvation.32,33 Hence, we selected it as an example to address the fundamental issues of modelling VCD spectra of solutes in aqueous solution, especially regarding the necessity for a comprehensive configurational sampling of the water hydrogen bonding network around the solute. We evaluated how the orientation of a single water molecule in a Pro-(H2O)n cluster affects VCD spectral signatures and we have shown how a QM/MM-MD can serve as basis for VCD spectra calculations without additional clustering or processing of the computed spectra.
Baseline corrections of the VCD spectra were carried out by subtracting those of the racemic mixture recorded under identical conditions. This procedure gave good mirror-image quality spectra, yet they exhibited some noise and a certain degree of baseline drift. Hence, for better comparison with the computed spectra, the spectra presented in the main text are baseline corrected by subtracting the half-sum spectrum of the enantiomers, that is, the average of the VCD spectra recorded for (R)- and (S)-proline. Comparison with the racemic-corrected spectra ensured that no artificial signals were generated through this procedure (cf. Fig. S1, ESI†). The IR spectra are presented as raw spectra.
Fig. 2 Structures of two zwitterionic conformers and the lowest energy anionic and protonated conformers of Pro. |
Geometry optimizations of these clusters were tedious, and due to convergence issues, we failed to optimize all 64 structures. Nonetheless, a clear picture can be drawn from the successfully converged 46 ZwPro-(H2O)5, of which 26 are based on ZwPro-c1 and 20 on ZwPro-c2. From an energy perspective, it is interesting to note that the structures based on ZwPro-c2 are generally lower in energy than those based on ZwPro-c1, so the explicit solvation with water flipped the conformational preferences compared to the only implicitly solvated structures. Within the set of ZwPro-(H2O)5 structures of each conformer, the zero-point corrected energies (ΔEZPC) are very similar and almost all are within a 0.2 kcal mol−1 range. Flipping the orientation of a single water molecule apparently does not cost much energy, resulting in evenly distributed Boltzmann weights.
The simulated VCD and IR spectra of ZwPro-(H2O)5 match the experimental signatures notably better than the spectra of non-solvated ZwPro (Fig. 4A). It is possible to make band assignments, and the characteristic signature found around 1400 cm−1 could be considered to be present in the simulated VCD spectrum. Yet, in terms of relative intensities, the spectral signatures are not well reproduced.
Fig. 4 VCD and IR spectra of zwitterionic ZwPro. (A) Comparison of the experimental spectra of (S)-proline in water (2 M, 6 μm pathlength) with the Boltzmann-averaged computed spectra of isolated ZwPro and solvated clusters ZwPro-(H2O)5. (B) Overlay of all spectra of ZwPro-c2-(H2O)5; the offset spectra belong to the structures shown in Fig. 3. |
As the conformer weights of the different solvated structures were quite similar, we investigated whether the water configuration around the solute affects the spectral signatures or if it may be sufficient to consider only one solvated structure for each conformer. Not unexpectedly, the spectral region above 1500 cm−1, where the CO stretching modes of ZwPro overlap with the water bending modes, exhibits strong differences among the clusters (Fig. 4B). However, the characteristic region of 1500–1300 cm−1 also changes notably, and the two offset examples demonstrate that the changes are significant. The two bands at 1410 and 1375 cm−1, which change in sign and magnitude, can both be ascribed to a combined motion of CH rocking, NH2 twisting, and CH2 wagging (displacement vectors are shown in Fig. S2, ESI†). The change in the orientation of a single water molecule causes a re-distribution of intensities between these two modes. Interestingly, none of the five water molecules contribute to the vibrational motion, which suggests that the influence of the water shell traces back solely to polarizing effects. The observation that the orientation of a single water molecule has such a tremendous effect on the computed spectral signatures led to a key conclusion: all or at least many water configurations need to be considered and it cannot be ruled out that higher clusters are required to better resemble the spectra.
Due to this conclusion, we did not attempt the micro-solvation approach for other cluster sizes of ZwPro-(H2O)n or any clusters with ProH(+) and Pro(−). The efforts to generate various cluster sizes, and even more so to ensure that a large number of structures also converge, did not seem feasible. Instead, we evaluated another route as described in the next section.
In the next step, we extracted 100 uncorrelated snapshots from each of the MD runs and removed all but the closest 30 water molecules from each of them. The decision to keep about 1.5 solvation shells was a compromise between the expected computer time required for the intended subsequent DFT-based geometry optimizations and having enough solvent molecules to stabilize the first shell during the optimization. The full structures were then optimized without any geometric constraints. All but a few structures converged, so that VCD and IR spectra calculations were subsequently carried out for 190 structures of ZwPro-(H2O)30. Expectedly, many structures showed small negative eigenvalues, but when they accounted for less than 150 cm−1, we still considered them for the following analysis. Removing all structures with negative eigenvalues gave very similar spectral features, but the spectra were less smooth and appeared “noisy”.
The individual VCD spectra of the ZwPro-(H2O)30 structures continue the trend seen already in ZwPro-(H2O)5 structures with VCD bands varying in intensity and sign and spanning the entire investigated fingerprint range (Fig. 5a). Interestingly, when considering the zero-point corrected energies of all clusters, the Boltzmann distribution is dominated by only two structures with weights of 68 and 21%, respectively. The rather even distribution observed for ZwPro-(H2O)5 thus breaks down when the cluster grows. This could be related to differences in the hydrogen bonding networks of the clusters, e.g., the number and types of hydrogen bonds, as both were kept identical in the evaluated ZwPro-(H2O)5 structures. The conformational space of ZwPro in the 30-fold solvated structures is also not represented by two sharply defined minima c1 and c2, but rather exhibits a conformational distribution (Fig. S4, ESI†). Most importantly, however, the resulting Boltzmann weighted spectra in fact exhibit poorer matching than those of ZwPro-(H2O)5 (Fig. 5b). Instead of relying on the DFT-computed energies, which may be inaccurate for various reasons, we simply computed the arithmetic average spectra over all computed structures of ZwPro-(H2O)30. This approach assumes that the MD simulation would deliver a reliable conformational distribution for ZwPro itself and the large number of samples should properly account for the configurational ensemble of water molecules. The VCD and IR spectra obtained from such a simple averaging procedure are in surprisingly good agreement with the experimental data (Fig. 5c). It becomes particularly evident in the IR spectrum, that the averaging reproduced band positions, shapes and line broadening well, and it even captured the onset of the libration band at the low-energy side of the spectrum. Likewise, the predicted VCD spectral pattern resembles the experimental signatures and relative intensities are reproduced quiet well.
The main observations described for ZwPro were made again for the other two protonation states. Examining the sets of hundreds of single spectra revealed bands over the entire investigated spectral range with changing signs and magnitudes as well as frequencies (cf. Fig. S5a andb, ESI†). Again, only few structures were favoured based on zero-point corrected energies and the Boltzmann-weighted spectra of ProH(+)-(H2O)30 do not match well with the experimentally observed spectra (cf. Fig. S5c, ESI†). For the anionic form, the Boltzmann weighted spectra of Pro(−)-(H2O)30 in fact resemble the spectral pattern, but significant line broadening has to be applied to match the overall spectra shapes (cf. Fig. S5d, ESI†). Once again, the simulated spectra obtained through simple averaging over the computed hundreds of single cluster spectra exhibit high agreement with the experimental VCD and IR spectral features of both ProH(+) and Pro(−) (Fig. 6). They closely resemble both the spectral pattern and band shapes, as well as the broadening.
Fig. 6 Comparison of the experimental VCD and IR spectra of (A) ProH(+) and (B) Pro(−) with the computed average spectra. |
To better account for the large configurational space of water molecules around ZwPro, we applied a QM/MM-MD based procedure to compute the VCD spectrum. Based on hundreds of snapshots, we obtained optimized geometries of ZwPro-(H2O)30 clusters. The vibrational spectra of ZwPro-(H2O)30 predicted based on zero-point corrected energies and the corresponding Boltzmann weights did not correctly resemble even the most prominent spectral features. In contrast, a simple averaging over all structures yielded VCD and IR spectra that were in excellent agreement with the experimental data. This observation led to the conclusion that the energies of such large clusters are not the appropriate basis for spectral simulations as these require a broad basis of structures with a statistical distribution of solvent orientations. This conclusion was further strengthened by the fact that the spectra of the protonated form of proline, ProH(+), and the anionic form, Pro(−), were equally well described by averaging of several hundreds of spectra computed for clusters with 30 water molecules.
The results of this study have shown that VCD spectra may be even more sensitive to the actual shape of the solvation shell than expected. However, while the presented approach to compute the VCD spectra of proline at different pH values appears simple to apply, we must stress that the conformational space of the solute is extremely small. Hence, rather few optimized snapshots may already suffice to cover its conformational distribution, so that most snapshots sample solely the possible solvent configurations. For molecules with greater degrees of conformational freedom, the number of snapshots and most likely also the number of solvent molecules in the clusters need to be adjusted, which results in a potentially drastic increase in the required number of snapshots. Future studies will thus evaluate the transferability of the approach to more conformationally complex molecules.
Footnote |
† Electronic supplementary information (ESI) available: Additional spectra plots, conformational distributions, conformer energies and Cartesian coordinates. See DOI: https://doi.org/10.1039/d4cp01768d |
This journal is © the Owner Societies 2024 |