Characterization of the Phenoale-Keto Oxyluciferin/Lucifrase Interactions in the S1 State by QM/MM Energy Decomposition Analysis

Supporting Information for Characterization of the Phenoale-Keto Oxyluciferin/Lucifrase Interactions in the S1State by QM/MM Energy Decomposition Analysis Henar Mateo-delaFuente,†Nuria Anguita-Ortiz,†Marcos Mandado,‡and Juan J. Nogueira∗,†,¶ †Chemistry Department,Universidad Autónoma de Madrid,Calle Francisco Tomás y Valiente, 7, 28049 Madrid, Spain ‡Department of Physical Chemistry, University of Vigo, Lagoas-Marcosende S/N, S36310 Vigo, Galicia, Spain ¶IADCHEM, Institute for Advanced Research in Chemistry, Universidad Autónoma de Madrid, Calle Francisco Tomás y Valiente, 7, 28049 Madrid, Spain E-mail: [email protected] 1 Computational Details of the Interaction Energy Decompositions To study the nature of the total interaction energy in the S1state, the last 600 ns of the classical molecular dynamics (CMD) simulation were used for the conformational sampling of the system. In a first stage, these 600 ns were used to calculate, every 3 ns, the binding free energy by means of the 1A-MM-GBSA approach1,2 using the MMPBSA.py tool.3This binding free energy was decomposed into (i) interaction (electrostatic and van der Waals) and solvation (polar and non-polar) terms and (ii) amino acid contributions to characterize the nature of the interactions and to determine which amino acids present a higher contribution to the interaction energy, as shown in Figure 1(B). In order to reduce the computational effort, the decomposition was performed for the 40 residues that are located at shortest distances from OLU along the CMD. Consequently, in order to obtain a more accurate description of the interaction energy, the classical sampling was employed to perform an Energy Decomposition Analysis (EDA) at quantum mechanics (QM)/molecular mechanics (QM/MM) level along the simulation, using a scheme based on electron deformation densities4–6 and treating the interaction between the QM and MM regions by an electrostatic embedding scheme. To do so, (i) the DFT functional and basis set are required to be suitable for the calculation and (ii) the number of amino acids included in the QM has to be converged. For that, a representative structure of the 600 ns of the trajectory was obtained by means of a hierarchical agglomerative (bottomup) clustering analysis using the cpptraj 7module of the AmberTools22.8Starting from an initial random frame, this algorithm analyzed the trajectory every 20 snapshots, finishing the clustering when at least the average distance between members of two clusters is greater than 50.0 Å. This convergence criteria ensures that only one structure was obtained. In the case of the selection of a suitable functional and basis set, monomer 1, treated at QM level, corresponded to the OLU chromophore while monomer 2, treated at QM/MM 2 level, represents the protein. Moreover, only PHE249—identified as the largest contributor to binding energy as shown in Figure 1(B) —was included in the QM region of monomer 2, and the rest of the protein was described by a force field within a MM electrostatic embedding scheme. The QM region was described with different combinations of the CAMB3LYP9–12 and the M062X13 DFT functionals and 6-31G(2d,p),14–18 6-31+G(2d,p),14–20 6-311+G(2d,p)15,19–22 and 6-31G(2df,p)14–16,18,21,23 basis sets with and without empirical dispersion treated by the D3BJ version of Grimme’s dispersion correction24 (for CAM-B3LYP) or by the D3 version of Grimme’s dispersion correction with the original D3 damping function25 (for M062X) due to implementation constrains. As discussed in the following sections (Sections 2 and 3), CAM-B3LYP is not a suitable DFT functional for the calculation of interaction energies and the most suitable functional/basis set combination for calculating binding energies is the M062X functional with the 6-31G(2d,p), as the inclusion of fpolarization functions only increases the computational time without improving the results, and the use of diffuse functions leads not only to a completely wrong decomposition of polarization into dispersive and inductive terms, but also to convergence problems. Once a suitable DFT functional was selected, the convergence with the number of amino acids in the QM region of the protein was studied by increasing the number of amino acids in the QM region of the luciferase from one to twenty and performing the EDA. The amino acids were added following the energetic ordering obtained from the 1A-MM-GBSA approach. The convergence for all the energy terms was achieved with fifteen amino acids. Subsequently, the convergence with respect to the number of frames was studied. For this, fifteen amino acids were included in the QM region of the protein, and the number of equidistant frames to perform the EDA was increased up to 250, achieving convergence after 125 frames. All the QM calculations were performed with the Gaussian09 software package26 using the counterpoise correction to overcome the basis set superposition error.27,28 Moreover, in order to calculate the interaction energies in the S1state, the electron densities of monomer 1 and complex in the first electronic excited state were described with TD-DFT. On the flip 3 side, the electron density of monomer 2 (the protein) was directly obtained at DFT level, using the same functional as for the chromophore, since the protein does not participate in the excited state of the system. This is clearly reflected on the analysis of the one-electron transition density matrix between the ground state and the S1states using the TheoDORE software.29 The input preparation for the QM calculations was automatically performed by means of the MoBioTools toolkit (https://github.com/mobiochem/MoBioTools),30 which employs the link atom approach to describe the QM/MM interface,31,32 and the EDA was carried out using the QM/MM version of the EDA-NCI program (https://github.com/marcosmandado/EDA-NCI).33 Benchmark of Functionals and Basis Sets Employing the representative structure of the trajectory, not only the DFT functional and the basis set were selected but also the convergence with the number of amino acids in the QM region of the protein was studied, as explained in the previous section. In order to determine the most suitable DFT functional/basis set pair for the QM/MMEDA calculation, different combinations of the CAM-B3LYP and the M062X DFT functionals and 6-31G(2d,p), 6-31+G(2d,p), 6-311+G(2d,p) and 6-31G(2df,p) basis sets with and without empirical dispersion described by the D3BJ version of Grimme’s dispersion with Becke-Johnson damping in the case of the CAM-B3LYP functional or by the D3 version of Grimme’s dispersion with the original D3 damping function in the case of the M062X functional due to implementation constrains, have been studied. As it can be observed in Figure S1(A-D), CAM-B3LYP, which has not been parameterized to properly describe dispersion, presents very consistent results. However, when empirical dispersion is included by means of Grimme’s dispersion with Becke-Johnson damping (Figure S1(D)), the polarization term gains participation to the total interaction energy due to a better description of dispersion (Figure S2(D)). In the case of the M062X functional (see Figure S1(E-G)), all the studied 4 basis set present very similar contributions to the total interaction energy of the EDA components. However, when diffuse functions are included in the calculation, an underestimate of the contribution of dispersion energy to the attractive energy components of the interaction energy is observed (Figure S2/(E-F)). (A) CAM-B3LYP (B) CAM-B3LYP 6-311+G(2d,p) 6-311G(2d,p) (C) CAM-B3LYP (D) CAM-B3LYP 6-31+G(2d,p) 6-31+G(2d,p) D3BJ (E) M062X (F) M062X (G) M062X 6-31+G(2d,p) D3 6-31G(2d,p) D3 6-31G(2df,p) D3 Figure S1: Percentages of contribution to the total interaction energy of the EDA components relative to the sum of the absolute values of the three components (electrostatic, pauli and polarization) calculated on the representative geometry. 5 (A) CAM-B3LYP (B) CAM-B3LYP 6-311+G(2d,p) 6-311G(2d,p) (C) CAM-B3LYP (D) CAM-B3LYP 6-31+G(2d,p) 6-31+G(2d,p) D3BJ (E) M062X (F) M062X (G) M062X 6-31+G(2d,p) D3 6-31G(2d,p) D3 6-31G(2df,p) D3 Figure S2: Percentages of contribution to the attractive terms of the interaction energy of the EDA components relative to the three attractive components (electrostatic, dispersion and induction) calculated on the representative geometry. 6 Convergence with the Number of Amino Acids For the convergence of the EDA with the number of amino acids, the twenty amino acids with the largest binding Gibbs free energy have been included in the QM region of the TD-DFT calculation one by one in decreasing order of contribution according to the 1AMM-GBSA analysis. These calculations were performed using the M062X functional and the 6-31+G(2d,p), 6-31G(2d,p) and 6-31G(2df,p) basis sets, together with the D3 version of Grimme’s dispersion with the original D3 damping function. Upon the inclusion of diffuse functions in the basis set, any of the EDA components seem to be close to convergence when the first fourteen amino acids are included in the QM region, as shown in Figure S3(A). Moreover, the dispersion term becomes divergent reaching an unrealistic positive value even though it is an attractive term. Thus, including diffuse functions in the basis set leads to a completely wrong description of the interaction energy. When increasing the number of polarization functions in the heavy atoms (see Figure S3(B,C,E,F)), no significant differences are observed in the energy components but the computational cost increases by a factor of 2.5. Thus, it is more suitable to use the 6-31G(2d,p) basis set in order to make the QM/MMEDA calculations more efficient. 7 (A) (B) (C) (D) (E) (F) Figure S3: (A-C) Convergence of the interaction energy and the EDA components and (D-F) error of the total interaction energy with respect to the calculation with twenty amino acids in the QM region as the number of amino acids in the QM region increases with the M062X functional and (A,D) 6-31+G(2d,p), (B,E) 6-31G(2d,p) and (C,F) 6-31G(2df,p) basis sets. Localization of the Excitation along the MD Trajectory Figure S4: Number of transferred electrons between the two monomers (green) and averaged position of the orbital before and after the electronic transition (orange) along the MD trajectory. 8 References (1) Srinivasan, J.; Miller, J.; Kollman, P. A.; Case, D. A. Continuum solvent studies of the stability of RNA hairpin loops and helices. J. Biomol. Struct. Dyn. 1998,16, 671–682. (2) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; others Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 2000,33, 889–897. (3) Miller III, B. R.; McGee Jr, T. D.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E. MMPBSA.py: an efficient program for end-state free energy calculations. J. Chem. Theory Comput. 2012,8, 3314–3321. (4) Mandado, M.; Hermida-Ramón, J. M. Electron density based partitioning scheme of interaction energies. J. Chem. Theory Comput. 2011,7, 633–641. (5) Ramos-Berdullas, N.; Pérez-Juste, I.; Van Alsenoy, C.; Mandado, M. Theoretical study of the adsorption of aromatic units on carbon allotropes including explicit (empirical) DFT dispersion corrections and implicitly dispersion-corrected functionals: the pyridine case. Phys. Chem. Chem. Phys. 2015,17, 575–587. (6) Mandado, M. Electron Density-Based Energy Decomposition Analysis: From QM to QM/MM Calculations; Springer, 2025. (7) Roe, D. R.; Cheatham III, T. E. PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data. J. Chem. Theory Comput. 2013,9, 3084–3095. (8) Case, D. A.; Aktulga, H. M.; Belfon, K.; Ben-Shalom, I. Y.; Berryman, J. T.; Brozell, S. R.; Cerutti, D. S.; Cheatham III, T. E.; Cisneros, G. A.; Cruzeiro, V. W. D.; others Amber 2022; University of California, San Francisco, 2022. 9