Noncovalent Interactions in Specific Recognition Motifs of Protein

Dec 19, 2016 - In view of the importance of protein–DNA interactions in biological processes, we extracted from the Protein Data Bank several one-to...
2 downloads 12 Views 1019KB Size
Subscriber access provided by UNIV OF WESTERN ONTARIO

Article

Noncovalent interactions in specific recognition motifs of protein–DNA complexes Olga A. Stasyuk, David Jakubec, Jiri Vondrasek, and Pavel Hobza J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00775 • Publication Date (Web): 19 Dec 2016 Downloaded from http://pubs.acs.org on December 19, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Noncovalent interactions in specific recognition motifs of protein–DNA complexes Olga A. Stasyuk,† David Jakubec,†,‡ Jiří Vondrášek,*,† and Pavel Hobza,*,†, †Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, Flemingovo nám. 2, 166 10 Prague, Czech Republic ‡Department of Physical and Macromolecular Chemistry, Faculty of Science, Charles University in Prague, Albertov 6, 128 43 Prague, Czech Republic Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Palacký University, 771 46 Olomouc, Czech Republic amino acid • nucleotide • hydrogen bond • interaction energy • implicit solvent

ABSTRACT: In view of the importance of protein–DNA interactions in biological processes, we extracted from the Protein Data Bank several one-to-one complexes of amino acids with nucleotides that matched certain geometric and energetic specificity criteria and investigated them using quantum chemistry methods. The CCSD(T)/CBS interaction energies were used as a benchmark to compare the performance of the MP2.5, MP2-F12, DFT-D3, and PM6-D3H4 methods. All methods yielded good agreement with the reference values, with declining accuracy

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 32

from MP2.5 to PM6-D3H4. Regardless of the site of interaction, the minima found after full optimization in implicit solvent with high dielectric constant were close to the structures experimentally detected in protein–DNA complexes. According to DFT-SAPT analysis, the nature of noncovalent interactions strongly depends on the type of amino acid. The negatively charged sugar-phosphate backbone of DNA heavily influences the strength of interactions and must be included in the computational model, especially in the case of interactions with charged amino acids.

Introduction Protein–DNA interactions play an essential role in biological processes including DNA replication, DNA repair, and regulation of gene expression. Over the years, two modes of sequence-specific recognition of DNA by proteins have been identified: direct recognition and indirect recognition.1 Direct recognition involves a series of hydrogen bonding interactions between amino acids and a complementary pattern of hydrogen bond donor/acceptor groups within the nucleobases of the recognised DNA sequence. For example, a statistical analysis of interactions within a set of 139 protein–DNA and 49 protein–RNA non-homologous complexes extracted from the Protein Data Bank (PDB) revealed that hydrogen bonds are the most common interactions (occurring in more than 50% of cases), followed by van der Waals, hydrophobic, and electrostatic interactions.2 Direct readout of DNA sequences is predominantly realised in the major groove of the DNA double-helix due to the greater variability of the hydrogen bond donor/acceptor groups offered. In 1976, Seeman et al.3 suggested that bidentate hydrogen bonds between certain amino acids and nucleobases are essential for specific base pair recognition in a

ACS Paragon Plus Environment

2

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

one-to-one fashion, whereas single hydrogen bonds are not sufficient for such unique identification. Since then, many other preferential interactions between amino acids and nucleobases have been identified.4 However, no generally applicable recognition code featuring one-to-one correspondences between individual amino acids and DNA residues has been described to date.5,6 In contrast, the indirect mode of DNA sequence recognition involves readout by the protein of certain structural and dynamic characteristics of the DNA molecule, including its flexibility, deformability, and intrinsic curvature, as well as the topologies of the major and minor grooves and sugar-phosphate backbone. Moreover, some deformation of the DNA helix is required to form hydrogen bonds or other noncovalent interactions between protein and DNA.7 Thus, most DNA-binding proteins utilize a combination of direct and indirect recognition mechanisms to achieve DNA binding specificity. A study of protein–DNA interactions at an atomic level8 revealed that two-thirds of hydrogen bonds with nucleobases involve bidentate and more complex interactions, which provide the greatest specificity, whereas van der Waals contacts are mainly used for stabilization of complexes without any notable sequence preferences. Additionally, some combinations of single hydrogen bonds, van der Waals contacts, and watermediated interactions can be essential for specificity in particular complexes.8 Analysis of the data in the Amino Acid-Nucleotide Interaction Database9 (AANT) helped reveal the significance of the sugar-phosphate part of the nucleotide in the specificity of protein–nucleic acid interactions. Various in vitro, in vivo, in silico, and biophysical techniques have been utilised for the study of protein–DNA interactions.10 For an accurate and quantitative description, application of computational approaches together with structural analysis seems to be reasonable. Because

ACS Paragon Plus Environment

3

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32

high-level quantum chemistry methods are limited to small- or medium-sized systems, a reductionist approach involving dissection of complex biological systems into their constituent parts is commonly utilised. This method relies on the assumption that the small isolated parts of a whole system and interactions between them reflect the characteristics of the real biological complex fairly well. Following this approach, interactions between selected pairs of amino acids and DNA bases have been examined extensively over the last two decades.4,8,11,12,13,14,15,16,17 Most of these studies were performed in the gas phase, although a few also took into account solvent effects.11,13,17 Recently, binding preferences for all 80 potential amino acid–DNA base combinations were studied by calculating the interaction energy profiles of the respective pairs.14 Using empirical potentials and DFT-D calculations, the researchers found that the geometries of several pairs featuring bidentate hydrogen bonds correspond to unique minima on the potential energy surface of the amino acid–DNA residue pairs. Subsequently, a benchmark study of the interaction energies was performed for 272 chosen amino acid‒base pairs.15 The results showed that all of the tested computational methods (MP2.5/CBS, MP2.X/CBS, MP2-F12, DFT-D3, PM6, and Amber force field) performed well for neutral complexes, while analysis of charged systems yielded poor agreement with the reference CCSD(T)/CBS data. To gain insight into the nature of interactions in hydrogen-bonded complexes of amino acid side chains and nucleobases, decomposition of the interaction energy was carried out.12 This revealed that for neutral complexes, stabilization is provided by delocalization energy associated with electron density deformation during the interaction, whereas for charged complexes, the electrostatic component is more important.

ACS Paragon Plus Environment

4

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Indeed, for a complete description of DNA sequence recognition, the role of the charged sugarphosphate DNA backbone should be analysed, because backbone contacts can also be important for specificity.18 The contribution of the backbone to the recognition process was recently studied in experimental structures of protein–DNA complexes using molecular mechanics (MM).17 This work showed that in the case of adenosine, the sugar-phosphate moiety can form a set of van der Waals contacts with aliphatic lysine or threonine side chains, increasing the stability of the protein–DNA complex. For pyrimidine nucleotides, single hydrogen bonds between hydroxyl or amide groups of the amino acid side chains and the phosphate group of the nucleotide were detected. Most of these contacts meet some specificity criteria17 in various environments and can be important for the direct recognition process. In our work, we focus on the amino acid–nucleotide pairs identified in the aforementioned study17 as contributors to the direct recognition of DNA sequences by proteins. This set of complexes contains different types of interactions, including single and bidentate hydrogen bonds and CH∙∙∙π contacts. In addition, the set of complexes includes both electrically neutral and charged systems (q = -1 or -2). To investigate the interaction specificity, we checked the geometry of the complexes upon optimization in different environments to see if “specific” complexes correspond to energy minima. The interaction energies were calculated using various computational methods, including CCSD(T)/CBS, MP2.5, MP2-F12, DFT-D3, and PM6-D3H4, to test the performance of these methods for the studied complexes. The greatest attention was paid to the applicability of the DFT-D3 and semiempirical PM6-D3H4 methods as potential methods for investigation of large protein–DNA binding motifs consisting of amino acids interacting with several DNA bases. Importantly, both of these techniques can be combined with the continuum solvent model, thus allowing calculation of interaction energy not only in the gas

ACS Paragon Plus Environment

5

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32

phase but also in solvent. To obtain new information about the nature of the observed noncovalent interactions and to study the influence of the sugar-phosphate backbone on the energetics of sequence-specific recognition, we analysed different components of the interaction energy using the DFT-SAPT approach.

Figure 1. Structures of the complexes analysed in this study (dAMP - deoxyadenosine monophosphate,

dGMP

-

deoxyguanosine

monophosphate,

dCMP

-

deoxycytidine

monophosphate, TMP – thymidine monophosphate). Computational Details Data set

ACS Paragon Plus Environment

6

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

A set of 1,584 X-ray structures of protein–DNA complexes with resolution better than 2.5 Å and R-factor no worse than 0.25 published in the Protein Data Bank (PDB) prior to April 2014 was dissected into the building blocks of each of the biopolymers to yield 3D distributions of amino acids around the four nucleobases (adenine, guanine, thymine, and cytosine). Only contacts found at the protein–DNA interface were considered for this process. Because of the directional nature of hydrogen bonds,19 amino acids tend to cluster in particular spatial regions around the nucleobases.8,14 We analysed each of these clusters based on geometric criteria and designated the contact with the greatest number of other contacts within the RMSD threshold of 1.5 Å as the cluster’s representative pair. For a more detailed derivation and description of the data set, please refer to previous work.14,17 Within each amino acid–nucleotide pair, the amino acid was reduced to its Cα representation to exclude nonspecific interactions between the DNA residue and the protein backbone.20 To this end, the carbonyl and amide groups of the peptide bond were replaced with hydrogen atoms. Interaction energies (including the deformation energy required to create the dimer from monomers) were calculated by a MM method for each amino acid–nucleotide pair to establish a link between the geometric conservation of the motifs and their energetic characteristics.17 MM interaction energy calculations utilising the generalised Born (GB) solvation model were also performed to study the effect of the surroundings. Three dielectric constants (ε = 4, 16, and 80) were chosen to simulate the properties of the protein core, protein– nucleic acid interface, and water, respectively. In the last step, specificity criteria were derived and applied to all amino acid–nucleotide pairs with the goal of identifying the most promising contacts for the direct recognition process. In summary, these criteria state that: (i) the orientation of the amino acid relative to the nucleotide must be found within one of the detected

ACS Paragon Plus Environment

7

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

geometric clusters; (ii) the interaction energy for a certain pair must be negative and unique (i.e., only a few contacts not belonging to the “specific” cluster can have similar interaction energy); and (iii) the stabilization energy of this complex must be greater than that of any other amino acid–nucleotide pair involving the same amino acid. Twelve complexes (Figure 1) that meet these conditions completely and follow some additional requirements were chosen for further investigation in this study. This set of complexes contains different types of amino acids (polar, non-polar, and charged) that are bound to DNA bases or phosphate groups via hydrogen bond(s) or CH∙∙∙π contacts. A detailed description of the selection of specific contacts is provided in the recent paper by Jakubec et. al.17

Interaction energy Calculations were performed for two types of geometry without any symmetry constraints: (i) after optimization of hydrogen atoms at the BLYP-D3/def2-QZVP level and (ii) after optimization of whole dimers at the BLYP-D3/def2-TZVPP level. The interaction energy was calculated according to the supermolecular method, in which energies of monomers are subtracted from the energy of the complex: ∆ =   −  −   

(1)

The CCSD(T) method extrapolated to the complete basis set (CBS) limit was used as a benchmark to calculate the reference values of the interaction energies. A composite CCSD(T)/CBS interaction energy can be written as: ()

∆

!" &'%"$% () = ∆   + ∆   − ∆ &'% )()$*∗∗(,.%.,,.$.) "$% + (∆

(2)

The first term, Hartree-Fock energy, is computed using the largest basis set. Two other terms account for correlation energy. For the second term, we used the MP2-F12 method21 instead of

ACS Paragon Plus Environment

8

Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

an extrapolation scheme to the CBS for MP2 method. This trick allowed us to reduce the number of calculations without much loss in accuracy, as shown in previous work.15 The last term, CCSD(T) correction, was calculated using the 6-31G**(0.25,0.15) basis set, because the size of our complexes varies from 44 to 58 atoms and the use of a larger basis set is not feasible. Moreover, the difference in correlation energies for the CCSD(T) and MP2 methods depends only slightly on the size of the basis set.22 Together with accurate CCSD(T) calculations, the interaction energies of the dimers were computed using several less expensive methods. Specifically, we tested the performance of the MP2-F12, MP2.5, DFT-D3, and PM6-D3H4 methods. The explicitly correlated MP2-F12 method, with a greatly improved rate of convergence compared to conventional MP2 theory, yields results close to the CBS limit even for relatively small basis sets.23 In our calculations, we used the cc-pVDZ-F12 basis set,24 which is sufficient for molecules involving elements of the first and second rows. The next step in perturbation theory, MP3 does not bring improvement over MP2. In contrast to MP2, MP3 underestimates the stabilization energy. Interestingly, the magnitude of these errors correlates well across complexes of different nature. This observation has led us to the introduction of MP2.5 method25,26 where the correlation energy is taken as average of MP2 and MP3 energies. As was previously shown for S22 and S66 data sets, this method can be considered as an accurate and computationally feasible alternative to the CCSD(T) method in the case of noncovalently bound systems and can be applied to large systems.27,28 The DFT-D3 method is a reasonable choice for fairly large weakly bound complexes. The DFT energy was calculated using the BLYP functional with the def2-QZVP basis set. London dispersion interactions were treated using empirical Grimme’s D3 correction29 with Becke-

ACS Paragon Plus Environment

9

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

Johnson (BJ) damping.30 According to the results for the S22 and S66 data sets,26,31 the BLYP functional shows a large improvement in combination with D3(BJ) correction and is one of the best GGA functionals with the RMSE (in kcal/mol) of 0.33 and 0.25 for the S22 and S66 data sets. The B3LYP/6-31G* and TPSS/TZVP methods provide for these data sets much larger values of RMSE (2.68 and 3.63; 0.69 and 0.58 respectively).31 Three-body dispersion contribution, which is important for biomolecules,32 was also taken into account. To facilitate investigation of more extended fragments of DNA–protein complexes, we tested the performance of the semiempirical PM6 method augmented with empirical corrections for dispersion and hydrogen bond (D3H4).33 The D3H4 correction is parametrized on a large set of benchmark data including different types of noncovalent interactions, and considerably improves results obtained by semiempirical quantum mechanical methods. The PM6-D3H4 calculations were carried out in MOPAC2012.34 All other calculations were performed in Turbomole V6.635 using the resolution of identity (RI) approximation36,37,38 to speed up the calculations. Except for semiempirical and DFT results, all interaction energies were corrected for BSSE using the counterpoise procedure proposed by Boys and Bernardi.39 The calculations were automated using the computational chemistry framework Cuby.40 Solvent effects were treated in an implicit fashion using the conductor-like screening model (COSMO)41,42 implemented in Turbomole V6.6 and MOPAC2012. The accuracy and performance of various implicit models based on quantum mechanical (QM) electronic densities (SMD, MST, COSMO-RS) as well as on molecular mechanics partial charges (GB, PB) were investigated in our recent paper,43 and it was shown that QM-based models provided best agreement with experimental data. We applied dielectric constants of 4, 16, and 78.5 to simulate the properties of the protein core, protein–DNA interface, and water, respectively. The

ACS Paragon Plus Environment

10

Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

interaction energy calculations in solvents were performed by DFT-D3 and PM6-D3H4 methods using the DFT-D3 geometries obtained after optimization of the complexes in implicit solvent. To estimate the amount of energy required to deform the equilibrium geometries of amino acids and nucleotides (E') into their geometries in the one-to-one complex (E), the deformation energy was calculated at the BLYP-D3/def2-QZVP level of theory: 1 1 ) ∆0 = ( −  + (   −    )

(3)

The total energy of interactions was obtained from the two aforementioned components: the interaction (∆Eint) and deformation (∆Edef) energies. (4)

∆  = ∆ + ∆0

Energy decomposition analysis The interaction energy in the complexes of amino acids and nucleotides was examined by the density fitting DFT-SAPT method44 performed with the aug-cc-pVDZ basis set, which allows quantitative decomposition of energy into several physically meaningful components: ($)

($)

(%)

(%)

(%)

(%)

(5)

 =  2 +  3 +  +  3 + 2 +  32 + 4(56) ($)

($)

The first term,  2 , corresponds to the electrostatic interaction between monomers;  3 (%)

(%)

contains the exchange-repulsion contribution; the next two terms  and  3 can be combined into one term called an induction part, which also covers the charge-transfer energy; (%)

(%)

the sum (2 +  32 ) accounts for the dispersion contribution; and the last term, δ(HF), contains the interaction energy contributions of third and higher order, which are computed using Hartree-Fock approximation and can be grouped with the induction part.45 For the monomer DFT calculations, the recommended treatment44,46 based on the asymptotically corrected LPBE0 potential with aug-cc-pVDZ basis set was used as implemented

ACS Paragon Plus Environment

11

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

in the Molpro 2010.1 program package.47 As it is known that DFT canonical virtual orbital energies are less than optimal for use in a perturbative scheme, their values were corrected before SAPT treatment using a gradient-controlled shift procedure. This step uses the difference between the exact vertical ionisation potential (IP) and the energy of the highest occupied molecular orbital (HOMO) of each monomer obtained by additional calculations of a neutral system and its cation.46 Different energy terms with the aug-cc-pVDZ basis set practically converged, with the exception of dispersion energy, which is underestimated by approximately 10% with this basis set. The dispersion term was therefore scaled using results from the aug-ccpVDZ and aug-cc-pVTZ basis sets.48 Scaling factors for each studied complex were calculated using model systems with a reduced number of atoms. To construct the model systems, we removed part of the complexes that did not participate directly in the interactions (the sugarphosphate backbone for complexes 1-7 and the nucleobase for complexes 8-12). Results and Discussion The complexes studied here (Figure 1) can be divided into groups in three different ways: (i) based on the site of interaction—complexes in which the amino acid interacts with the nucleobase (1-7) or phosphate group (8-12); (ii) based on the type of interaction—complexes with one hydrogen bond (3, 4, 8-12), with two hydrogen bonds (1, 2, 5, 6), and with CH∙∙∙π contacts (7); and (iii) based on the type of amino acid—complexes with polar (1-3, 8-10, 12), non-polar (7, 11), and charged amino acids (4-6).

Testing of methods Recently, the performance of several computational methods for calculation of the interaction energies between amino acid side chains and nucleobases was investigated in a variety of

ACS Paragon Plus Environment

12

Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

representative complexes.15 According to the results and following our expectations, the largest errors were observed for charged systems. In this study, the performance of various semiempirical, DFT and ab initio methods was compared for the complexes between amino acids and nucleotides, i.e., containing charged sugar-phosphate backbone; therefore, the applicability of different lower-level computational methods to complexes of such type needs to be carefully investigated. We paid greatest attention to the performance of the PM6-D3H4 and DFT-D3 methods as potential methods for energy calculations in more extended protein–DNA motifs. Geometries taken from crystal structures with optimized hydrogen atoms at the BLYPD3/def2-QZVP level of theory were used in this step. Interaction energies obtained at various levels of theory for all studied complexes are listed in Table S1. The reference CCSD(T)/CBS values are found in a large range from -5.7 kcal/mol (for the dCMP-Ile complex) to -85.8 kcal/mol (for the dGMP-Arg complex). As expected, the strengths of interaction between oppositely charged systems (negatively charged nucleotide and positively charged amino acid) are the largest, as in complexes 4 and 5. These are followed by complexes of nucleotides with neutral amino acids, and then complex 6 (nucleotide with a negatively charged amino acid). The weakest interaction is in stacked complex 7. All applied methods yielded well-correlated results, with correlation coefficients (cc) greater than 0.997. As our attention focused on the performance of the DFT-D3 and PM6-D3H4 methods as potential methods for investigation of large protein–DNA motifs, Figure 2 shows the correlations between the results obtained by these two methods and the reference values. Correlations for other methods are shown in Figure S1. A comparison of the results is presented in Figure 3 as MAD and root-mean-square deviation (RMSD) from the reference CCSD(T)/CBS values, which characterize the overall accuracy of methods; more details are given in the Supporting

ACS Paragon Plus Environment

13

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 32

Information (Table S2, Figure S2). The performance of MP2.5/CBS was the best among the methods tested, with MAD of 0.27 kcal/mol and RMSD of 0.31 kcal/mol. The error cancellation between the overestimated MP2 and underestimated MP3 results is mostly responsible for the high accuracy of the MP2.5 method.25 For all studied complexes, the MP2.5 results were less negative than the CCSD(T) values, which indicates underestimation of the interaction energy. These results are in good agreement with the results obtained for large noncovalent complexes.49 The largest deviations from the benchmark data (ca. 0.5 kcal/mol) were observed for the following three complexes: dAMP with lysine (4); dGMP with arginine (5), and dCMP with isoleucine (7). This is most likely a consequence of different magnitudes of the errors produced by MP2 and MP3, which cannot compensate for each other. This imbalance can be larger for interactions between charged monomers and in the case of dispersion interaction.25

Figure 2. Correlations between interaction energy values (in kcal/mol) obtained by the DFT-D3 and PM6-D3H4 methods and the reference CCSD(T) values.

ACS Paragon Plus Environment

14

Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 3. MAD and RMSD values of tested methods for interaction energy over the studied set of complexes, compared with the CCSD(T)/CBS method. Deviations for the results obtained at the MP2-F12/cc-pVDZ-F12 level of theory are almost identical in magnitude to those obtained with the MP2.5/CBS method, although systematic underestimation was not observed. MP2-F12 is second best of all tested methods. This method performed poorly only for the complex of TMP with Tyr (i.e., with an aromatic amino acid). For the TMP-Tyr complex (11), the interaction energy was overestimated by 1.0 kcal/mol compared to the reference value. This could be due to partial stacking between TMP and Tyr; it is known that the MP2 method overestimates stabilization energy for stacking of aromatic systems.50 For other complexes, the results were reproduced very well. Taking into consideration its high accuracy and relatively small computational cost, this method can be considered in the future as a benchmark method in situations in which standard high-level quantum chemical methods such as CCSD(T) and MP2.5 are unreachable and no stacking between subsystems occurs. As mentioned, the DFT and PM6 methods with corrections for noncovalent interactions are the most promising methods for investigation of large complexes between DNA and proteins. The DFT-D3 method yielded systematically overestimated interaction energy, which is in line with previous results.15,51 The greatest deviations (ca. 2 kcal/mol) were found for complexes of

ACS Paragon Plus Environment

15

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 32

positively charged amino acids and nucleotides (i.e., for complexes between oppositely charged monomers). However, the relative error in these cases is not significant, because the interactions are very strong (-∆Eint > 75 kcal/mol). A similar tendency was previously reported for oppositely charged pairs of amino acids.52 Interestingly, the strength of the interaction in complex 6, which includes a negatively charged amino acid, was underestimated by 0.9 kcal/mol. We do not yet have an explanation for this observation. The semiempirical PM6 method with corrections for dispersion (D3) and hydrogen bonding (H4) can be applied to very large systems (up to 10,000 atoms), providing acceptable accuracy.33 For our complexes, the PM6-D3H4 method reproduces the reference CCSD(T)/CBS interaction energies (Figure 2b) quite well but with the largest errors among all the methods studied. The largest deviation (5.6 kcal/mol) was found in the dGMP-Arg complex, and substantial deviations were also found in complexes in which amino acids interact with the sugar-phosphate part of the nucleotide. The reason for this is likely the correction for hydrogen bonding, which was parametrized for charge-assisted hydrogen bonds but not including phosphate groups.

Effect of solvent A common technique for modelling different environments in biomolecules is use of an implicit solvent model.53 We first tested the dependency of calculated solvation energies on DFT functionals and basis sets. Forty solvation free energies of various small systems were taken from Ref. [54] and they were recalculated with the following functionals – BLYP, TPSS and BP-86, and using SVP and TZVP basis sets. The RMSEs of all these calculated energies differ by less than 10% showing their small dependence on DFT functional and/or basis set. Therefore, to be consistent we have performed calculations by DFT-D3 method at the BLYP-D3/def2-

ACS Paragon Plus Environment

16

Page 17 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

QZVP level of theory, as well as by PM6-D3H4 method. To calculate the interaction energies between amino acids and nucleotides in the protein interior, we chose a dielectric constant of 4, which is commonly used in continuum solvent methods.55 In reality, the protein surface consists of parts with very different dielectric properties; therefore, the dielectric constant of proteins is not a “constant” and depends on structural characteristics such as packing and the presence of internal cavities and their shapes.56 The average dielectric constant for each type of amino acid varies from 11.0 to 25.6. Additionally, the dielectric constant of DNA spans a wide range from 8 to 16 or even more.57,58 Thus, we applied a dielectric constant of 16 to simulate the properties of the DNA–protein interface. In the last step, we investigated the influence of the water environment on the interaction energy using a dielectric constant value of 78.5. We used the geometries from crystal structures in which only hydrogen atoms were optimized at the BLYPD3/def2-QZVP level of theory with an implicit solvent model. As expected, the absolute values of the interaction energy decreased linearly when passing from the gas phase to higher dielectric constants due to solvent attenuation of electrostatics between interacting partners (Table 1). The only exception was the dGMP-Asp complex, in which interaction between two negatively charged monomers occurs. In this case, introduction of solvent reduces unfavorable electrostatic repulsion between partners; thus, the interaction becomes stronger. For complexes between oppositely charged partners (4 and 5 in Figure 1), we observed a dramatic reduction of the interaction energy in all solvents, but these interactions were still stabilizing, even in the water environment where the energy loss was nearly 90%. The strength of the interaction in the stacked dCMP-Ile complex decreased to a much lesser extent. Correlation between the results obtained by the PM6-D3H4 and DFT-D3 methods worsened with increasing dielectric constant

ACS Paragon Plus Environment

17

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32

(Figure S3). The greatest outliers are complexes with charged amino acids, suggesting that the problem may stem from a poor description of charge distribution by the PM6 method. Table 1. Interaction energy in solution (in kcal/mol) obtained at the BLYP-D3/def2-QZVP level of theory using COSMO calculations. Complex

gas phase ε=4

ε=16

ε=78.5

dAMP-Gln dAMP-Asn dAMP-Thr dAMP-Lys

-17.5 -17.9 -10.2 -75.1

-8.0 -7.2 -5.5 -12.3

-7.1 -6.1 -5.0 -7.8

-11.0 -10.5 -6.8 -28.4

dGMP-Arg -86.9 dGMP-Asp -5.9

-35.6 -16.4 -10.6 -10.8 -11.3 -11.3

dCMP-Ile dCMP-Gln

-5.8 -23.7

-3.0 -2.0 -12.6 -8.1

-1.7 -6.7

TMP-Thr TMP-Gln TMP-Tyr TMP-Ser

-22.1 -23.0 -28.5 -18.0

-11.2 -12.4 -15.8 -9.3

-5.6 -6.7 -9.9 -4.8

-6.9 -8.1 -11.3 -5.8

Nature of noncovalent interactions To understand why the environment affects the stabilization of the studied complexes to varying degrees, we decomposed the interaction energies into physically meaningful components using the DFT-SAPT approach.44 The interaction energies obtained by DFT-SAPT with the aug-ccpVDZ basis set are systematically smaller that the reference CCSD(T)/CBS values. It is known that the dispersion term is mostly underestimated;48 therefore, we scaled this term by a factor of 1.10 for complexes 1-7 (amino acid interacts with nucleobase) and by a factor of 1.13 for complexes 8-12 (amino acid interacts with phosphate group). Both scaling factors are close to widely used value of 1.1 which corresponds to 10% underestimation of dispersion energy when

ACS Paragon Plus Environment

18

Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

aug-cc-pVDZ basis set is used. Evidently, the individual scaling for different cluster type is not required and recommended value of 1.1 can be safely used. Recalculated interaction energies accounting for the new dispersion term correlated well with the reference values (Figure S4). The accuracy of the results was close to that of the DFT-D3 method (Table S2). The percentage contribution of all attractive terms, namely electrostatic, induction, and dispersion, to their total value is shown in Figure 4. We found that for complexes with positively charged amino acids (Lys, Arg), the attractive interaction is provided mostly by the electrostatic component (more than 60%), followed by induction and dispersion. The same was observed for complexes with polar amino acids. However, for complexes with a negatively charged amino acid (Asp), the electrostatic term is the least important, and the induction term prevails. For the stacked dCMP-Ile complex, the dispersion term was more important, as expected. Thus, we can conclude that complexes with a high prevalence of electrostatic interactions are more susceptible to solvent effects, regardless of the site of interaction (either the nucleobase or the phosphate group). The strength of the interaction between charged monomers decreases with increasing value of the dielectric constant. Interestingly, the behavior of the dCMP-Ile complex (prevailing role of dispersion) in solvents is similar to the trend found for hydrogen-bonded complexes of amino acids with the phosphate group of nucleotides.

ACS Paragon Plus Environment

19

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32

Figure 4. Relative contributions of electrostatic attraction (∆Eelst), induction (∆Eind), and dispersion interaction (∆Edisp) to all bonding interactions (∆Eelst + ∆Eind + ∆Edisp) for the studied complexes in the gas phase. Another benefit of DFT-SAPT decomposition is the opportunity to study the effect of more distant parts of complexes, those that do not participate directly in the interactions, on the nature of bonding. For complexes 1-7 (interaction site is the nucleobase), we studied the effect of the charged sugar-phosphate backbone on the nature of the interaction. For complexes 8-12 (interaction site is centred on the phosphate group), we estimated the contribution of the nucleobase to the energy components. We found that the role of the sugar-phosphate part is important only for the complexes with charged amino acids (Lys, Arg, and Asp), whereas nucleobases do not affect the character of interactions with the phosphate group (Figure S5). This means that for realistic modelling of interactions in protein–DNA complexes, the negatively charged sugar-phosphate backbone has to be considered, especially in the case of interactions with charged amino acids. However, the neutral distant parts of complexes could be eliminated to simplify the analysis.

ACS Paragon Plus Environment

20

Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Geometry optimization To check the specificity and importance of the studied motifs, we performed a full optimization of their geometries by the DFT-D3 method. As the strength of noncovalent interactions is dependent on the environment, we also tested the geometry preservation of the studied complexes upon optimization in different solvents. The most suitable indicator of geometrical changes is the RMSD value of atom positions in a fully optimized structure and their positions in a crystal. To compare the structures after constrained (only the H atoms) and full optimization, we calculated RMSD values for complexes using the VMD program.59 Figure 5 shows the RMSD values for all complexes studied in different environments. These deviations are largest for the gas phase, as we expected, but taking solvent into account led to the final structures, which are very similar geometrically to those found in crystals. In most cases, a higher solvent dielectric constant resulted in high similarity between the experimental and optimized structures. This indicates that the closest local minima found in solvents with high dielectric constants are more likely to represent a biologically relevant geometric orientation of the nucleotide with respect to the amino acid. It is important to note that during the optimization process all existing hydrogen bonds are not broken, even in the structures with RMSD greater than 2.5 Å (Figure 6). The only exception is the dGMP-Arg complex. Full optimization in the gas phase resulted in a bowl-shaped structure that cannot normally be adopted in protein–DNA complexes. Comparison of crystal and fully optimized structures for other complexes is provided in Table S3.

ACS Paragon Plus Environment

21

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32

Figure 5. RMSD values (in Å) of atom positions of nucleotides (fully optimized complex vs. crystal structure) in different environments.

Figure 6. Structures of complexes with the largest RMSD values in the gas phase. Crystal structures are shown in red and fully optimized structures in green.

ACS Paragon Plus Environment

22

Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

We are aware that the potential energy surface of a biomolecular cluster is characterized by a large amount of energy minima, while optimization mostly terminates in the nearest local minimum. To avoid this geometry bias, we started full optimization from structures generated by elongation of the distance between interacting atoms by 1.0 Å along the corresponding axis. It is worth noting that all disrupted interactions were completely recovered after the optimization. However, for several complexes, new minima found in the gas phase and three types of solvent are lower-lying than those previously discussed and determined experimentally. In the gas phase, the maximum energy difference is only -0.2 kcal/mol for the dCMP-Gln complex, whereas in high-dielectric solvents, the new found minima are even lower, up to -4.9 kcal/mol (for the TMPThr complex in water). Generally, the largest energy differences between the new and previously found minima are observed for complexes with the interaction located at the phosphate group (complexes 8-12 in Figure 1). Structures with two hydrogen bonds formed between the nucleobase and amino acid are more rigid and can change their geometry only marginally, whereas the sugar-phosphate part of nucleotides is more flexible and can form more stable complexes with amino acids than those found the strongest in crystals.

Deformation energy In real protein–DNA complexes, the energy required to deform DNA from its native conformation to the conformation in a protein-bound complex was found to be responsible for indirect recognition for a series of DNA-binding proteins.60 Therefore, the amount of energy needed to deform the monomers (amino acid and nucleotide) from their equilibrium conformations to the geometries in the complex can serve as an additional factor affecting the

ACS Paragon Plus Environment

23

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 32

recognition process. Calculation of the deformation energy revealed that for structures with hydrogen atoms optimized, its value does not exceed 2.0 kcal/mol in the gas phase and 1.2 kcal/mol in solvent. In all environments, the largest nucleotide deformation was observed in the dGMP-Asp complex. For fully optimized complexes, the deformation energy was even higher, with a maximum value of 8.8 kcal/mol in the gas phase and 3.9 kcal/mol in water. Gas-phase calculations showed the largest deformation of nucleotide in complexes with the hydrogen bond(s) between nucleobases and charged amino acids, whereas in solvent, the nucleotide deformed more in complexes with the hydrogen bond between the phosphate group and amino acid. Importantly, the total energy calculated as a sum of interaction and deformation energies was negative, indicating stabilizing interactions in all complexes (Table S4). The only exception was that the total energy for the fully optimized structure of the dGMP-Asp complex which was positive (+1.7 kcal/mol) in the gas phase. This indicates that deformation of the monomers exceeds the interaction energy between them, which is not surprising in this case because two negatively charged monomers interact with each other.

Conclusions Noncovalent interactions were analysed in twelve one-to-one amino acid–nucleotide complexes that matched geometric and energy criteria of specificity.17 Testing of several computational methods, including MP2.5, MP2-F12, DFT-D3, and PM6D3H4, showed that all these methods can be applied to systems containing a charged sugarphosphate backbone and different types of amino acids. They reproduce the reference CCSD(T)/CBS interaction energies quite well, with the largest errors found for the PM6-D3H4

ACS Paragon Plus Environment

24

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

method. DFT-D3 provides sufficient accuracy, but systematically overestimates the interaction energy. MP2-F12/cc-pVDZ-F12 is as accurate as MP2.5/CBS and could be used as a benchmark method in cases in which standard high-level quantum chemical methods such as CCSD(T) and MP2.5 are unreachable and complexes do not exhibit stacking features. Full optimization of the structures in different environments revealed that inclusion of solvents with high dielectric constants leads to minima that are close to structures experimentally detected in the protein–DNA complexes. Due to attenuation of electrostatic attraction between interacting partners by solvents, the absolute values of the interaction energy decrease monotonically when passing from the gas phase to water. The strength of interaction between charged monomers decreases with increasing dielectric constant due to stronger stabilization of monomers in the solvent. Decomposition of the interaction energies using the DFT-SAPT approach helps explain the nature of stabilizing interactions and assess the role of the nucleic acid backbone. For positively charged and polar amino acids, the attractive interaction is provided mostly by the electrostatic component (more than 60%) followed by induction and dispersion. For negatively charged amino acids, the induction term is prevalent, whereas in a stacked complex with a non-polar amino acid, dispersion contributes significantly. For realistic modelling of interactions in protein–DNA complexes, the negatively charged sugar-phosphate backbone is essential for complexes with charged amino acids (Lys, Arg, and Asp), whereas inclusion of nucleobases in the model does not affect the nature of interactions with the phosphate group and could be eliminated to simplify the analysis.

ASSOCIATED CONTENT

ACS Paragon Plus Environment

25

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 32

Supporting Information Gas phase interaction energies obtained by different computational methods and their errors (Table S1 and S2); comparison of crystal and fully optimized structures in different environments (Table S3); total energies for the studied complexes (Table S4). Correlations for interaction energies obtained by different computational methods and signed errors (Figure S1, S2 and S4); correlations for DFT-D3 and PM6-D3H4 interaction energies in different environments (Figure S3); DFT-SAPT decomposition for the reduced model systems (Figure S5). This material is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION Corresponding Authors * E-mail: [email protected] (J.V.). * E-mail: [email protected] (P.H.). ACKNOWLEDGMENT This work was part of the Research Project RVO: 61388963 of the Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic. This work was also supported by the Czech Science Foundation [P208/12/G016].

The authors gratefully

acknowledge support from project LO1305 of the Ministry of Education, Youth and Sports of the Czech Republic.

REFERENCES

ACS Paragon Plus Environment

26

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1

) Steffen, N.R.; Murphy, S.D.; Tolleri, L.; Hatfield, G.W.; Lathrop, R.H. Bioinformatics 2002, 18, S22-S30. 2

) Lejeune, D.; Delsaux, N.; Charloteaux, B.; Thomas, A.; Brasseur, R. Proteins 2005, 61, 258–271. 3

) 808.

Seeman, N.C.; Rosenberg, J.M.; RICH, A. Proc. Natl. Acad. Sci. U.S.A. 1976, 73, 804-

4

Cheng, A.C.; Chen, W.W.; Fuhrmann, C.N.; Frankel, A.D. J. Mol. Biol. 2003, 327, 781–

) 796. 5

)

Sarai, A.; Kono, H. Annu. Rev. Biophys. Biomol. Struct. 2005, 34, 379–398.

6

)

Siggers, T.; Gordan, R. Nucleic Acids Res. 2013, 1–13.

7

) Rohs, R.; Jin, X.; West, S.M.; Joshi, R.; Honig, B.; Mann, R.S. Annu Rev Biochem. 2010, 79, 233–269. 8

) Luscombe, N.M., Laskowski, R.A., Thornton, J.M. Nucleic Acids Res. 2001, 29, 2860– 2874. 9

) Hoffman, M.M.; Khrapov, M.A.; Cox, J.C.; Yao, J.; Tong, L.; Ellington, A.D. Nucleic Acids Res. 2004, 32, D174-D181. 10

) Dey, B.; Thukral, S.; Krishnan, S.; Chakrobarty, M.; Gupta, S.; Manghani, C.; Rani, V. Mol. Cell. Biochem. 2012, 365, 279–299. 11

) Rutledge, L.R;. Durst, H.F.; Wetmore, S.D. Phys. Chem. Chem. Phys. 2008, 10, 28012812. 12

) Czyżnikowska, Ż.; Lipkowski, P.; Góra, R.W.; Zaleśny, R.; Cheng, A. C. J. Phys. Chem. B 2009, 113, 11511–11520 13

)

de Ruiter, A.; Zagrovic, B. Nucleic Acids Res. 2015, 43, 708-718.

14

) Jakubec, D.; Hostaš, J.; Laskowski, R.A.; Hobza, P.; Vondrášek, J. J. Chem. Theory Comput. 2015, 11, 1939–1948. 15

) Hostaš, J.; Jakubec, D.; Laskowski, R.A.; Gnanasekaran, R.; Řezáč, J.; Vondrášek, J.; Hobza, P. J. Chem. Theory Comput. 2015, 11, 4086–4092. 16

) González, J.; Baños, I.; León, I.; Contreras-García, J.; Cocinero, E.J.; Lesarri, A.; Fernández, J.A.; Millán, J. J. Chem. Theory Comput. 2016, 12, 523–534. 17

)

Jakubec, D.; Laskowski, R.A.; Vondrasek, J. PLoS One 2016, 11 (7), e0158704

ACS Paragon Plus Environment

27

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

18

Page 28 of 32

) Luscombe, N.M.; Thornton, J.M. J Mol Biol. 2002, 320, 991–1009.

19

) Arunan, E.; Desiraju, G.R.; Klein, R.A.; Sadlej, J.; Scheiner, S.; Alkorta, I.; Clary, D.C.; Crabtree, R.H.; Dannenberg, J.J.; Hobza, P.; Kjaergaard, H. G.; Legon, A.C.; Mennucci, B.; Nesbitt, D.J. Pure Appl. Chem. 2011, 83, 1637–1641. 20

) Berka, K.; Laskowski, R.A.; Riley, K.E.; Hobza, P.; Vondrasek, J. J. Chem. Theory Comput. 2009, 5, 982−992 21

) Bachorz, R.A.; Bischoff, F.A.; Glöß, A.; Hättig, C.; Höfener, S.; Klopper, W.; Tew, D.P. J. Comput. Chem. 2011, 32, 2492–2513. 22

)

Jurečka, P.; Hobza, P. Chem. Phys. Lett. 2002, 365, 89–94.

23

) Riley, K.E.; Platts, J.A.; Řezac, J.; Hobza, Hill, J.G. J. Phys. Chem. A 2012, 116, 4159– 4169. 24

)

Peterson, K.A.; Adler, T.B.; Werner, H.-J. J. Chem. Phys. 2008, 128, 084102.

25

) Pitoňák, M., Neogrády, P., Černý, J., Grimme, S.; Hobza, P. ChemPhysChem 2009, 10, 282–289. 26

) Řezáč, J.; Hobza, P. Chem. Rev. 2016, 116, 5038–5071.

27

) Sedlak, R.; Riley, K.E.; Řezáč, J.; Pitoňák, M.; Hobza, P. ChemPhysChem 2013, 14, 698–707. 28

) Brauer, B.; Kesharwani, M.K.; Kozuch, S.; Martin, J.M.L. Phys. Chem. Chem. Phys. 2016, DOI: 10.1039/C6CP00688D. 29

)

Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104.

30

)

Becke, A.D.; Johnson, E.R. J. Chem. Phys. 2005, 122, 154104.

31

) Goerigk, L.; Kruse, H.; Grimme, S. ChemPhysChem 2011, 12, 3421-3433.

32

)

von Lilienfeld, O.A.; Tkatchenko, A. J. Chem. Phys. 2010, 132, 234109.

33

)

Řezáč, J.; Hobza, P. J. Chem. Theory Comput. 2012, 8, 141–151.

34

) Stewart, J.J.P. MOPAC2012; Stewart Computational Chemistry: Colorado Springs, CO, USA, 2012; http://OpenMOPAC.net.

ACS Paragon Plus Environment

28

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

35

) TURBOMOLE V6.6, http://www.turbomole.com.

TURBOMOLE

GmbH:

Karlsruhe,

Germany,

2014;

36

) Eichkorn, K.; Treutler, O.; Öhm, H.; Häser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 242, 652-660. 37

)

Weigend, F.; Häser, M. Theor. Chem. Acc. 1997, 97, 331-340.

38

)

Weigend, F.; Haser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143−152.

39

)

Boys, S.F.; Bernardi, F. Mol. Phys. 1970, 19, 553−566.

40

)

Řezáč, J. J. Comput. Chem. 2016, 37, 1230-1237.

41

)

Klamt, A.; Schüürmann, G. J. Chem. Soc. Perkin Trans. 2 1993, 799-805.

42

)

Klamt, A. J. Phys. Chem. 1995, 99, 2224-2235.

) Kolář, M.; Fanfrlík, J.; Lepšík, M.; Forti, F.; Luque, F.J.; Hobza, P. J. Phys. Chem. B 2013, 117, 5950–5962. 43

44

)

Heßelmann, A.; Jansen, G.; Schütz, M. J. Chem. Phys. 2005, 122, 014103.

45

)

Hohenstein, E.G.; Sherrill, C.D. J. Chem. Phys. 2010, 133, 014101.

) Grüning, M.; Gritsenko, O.V.; van Gisbergen, S.J.A.; Baerends, E.J. J. Chem. Phys. 2001, 114, 652-660. 46

47

) MOLPRO v. 2010.1, a package of ab initio programs, Werner, H.-J.; Knowles, P.J;. Knizia, G.; Manby, F.R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K.R.; Adler, T.B.; Amos, R.D.; Bernhardsson, A.; Berning, A.; Cooper, D.L.; Deegan, M.J.O.; Dobbyn, A.J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; Köppl, C.; Liu, Y.; Lloyd, A.W.; Mata, R.A.; May, A.J.; McNicholas, S.J.; Meyer, W.; Mura, M.E.; Nicklass, A.; O'Neill, D.P.; Palmieri, P.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A.J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; Wolf, A.; http://www.molpro.net. 48

)

Řezáč, J.; Hobza, P. J. Chem. Theory Comput. 2011, 7, 685–689.

ACS Paragon Plus Environment

29

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 32

49

) Sedlak, R.; Janowski, T.; Pitoňák, M.; Řezáč, J.; Pulay, P.; Hobza, P. J. Chem. Theory Comput. 2013, 9, 3364−3374. 50

)

Riley, K.E.; Pitoňák, M.; Černý, J.; Hobza, P. J. Chem. Theory Comput. 2010, 6, 66–80.

51

)

Antony, J.; Grimme, S.; Liakos, D.G.; Neese, F. J. Phys. Chem. A 2011, 115, 11210–

11220. 52

)

Antony, J.; Grimme, S. Phys. Chem. Chem. Phys. 2006, 8, 5287–5293.

53

) Roux, B. In Computational Biochemistry and Biophysics; Becker, O.M; Mackerell, Jr., A.D.; Roux, B.; Watanabe, M., Ed.; Marcel Dekker: New York, 2001; pp. 133–152. 54

) Kongsted, J.; Soderhjelm, P.; Ryde, U. J. Comp. Aided Mol. Des. 2009, 23, 395-409.

55

) Dwyer, J.J.; Gittis, A.G.; Karp, D.A.; Lattman, E.E.; Spencer, D.S.; Stites, W.E.; GarcıaMoreno, B.E. Biophys. J. 2000, 79, 1610–1620. 56

)

Li, L.; Li, C.; Zhang, Z.; Alexov, E. J. Chem. Theory Comput. 2013, 9, 2126−2136.

57

) Cuervo, A.; Dans, P.D.; Carrascosa, J.L.; Orozco, M.; Gomila, G.; Fumagalli, L. Proc. Natl. Acad. Sci. U.S.A. 2014, 111, E3624-30. 58

)

Yang, L.; Weerasinghe, S.; Smith, P.E.; Pettitt, B.M. Biophys. J. 1995, 69, 1519–1527.

59

)

Humphrey, W., Dalke, A.; Schulten, K. J. Molec. Graphics 1996, 14, 33-38.

60

) Aeling, K.A.; Steffen, N.R.; Johnson, M.; Hatfield, G.W.; Lathrop, R.H.; Senear, D.F. IEEE/ACM Trans Comput Biol Bioinform. 2007, 4, 117-25.

ACS Paragon Plus Environment

30

Page 31 of 32

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

31 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 32

Table of Contents/ Abstract Graphic

32 ACS Paragon Plus Environment