Fragment Molecular Orbital Calculations with Implicit Solvent Based

Subscriber access provided by Warwick University Library

B: Biophysical Chemistry and Biomolecules

Fragment Molecular Orbital Calculations with Implicit Solvent Based on the Poisson–Boltzmann Equation: Implementation and DNA Study Yoshio Okiyama, Tatsuya Nakano, Chiduru Watanabe, Kaori Fukuzawa, and Shigenori Tanaka J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b01172 • Publication Date (Web): 20 Mar 2018 Downloaded from http://pubs.acs.org on April 5, 2018

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 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 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.

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 51 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

The Journal of Physical Chemistry

Fragment Molecular Orbital Calculations with Implicit Solvent Based on the Poisson–Boltzmann Equation: Implementation and DNA Study Yoshio Okiyama,∗,†,∥ Tatsuya Nakano,†,‡ Chiduru Watanabe,† Kaori Fukuzawa,†,¶ and Shigenori Tanaka§ †Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan ‡National Institute of Health Sciences, 3-25-26 Tonomachi, Kawasaki-ku, Kawasaki, Kanagawa 210-9501, Japan ¶Faculty of Pharmaceutical Sciences, Hoshi University, 2-4-41 Ebara, Shinagawa-ku, Tokyo 142-8501, Japan §Graduate School of System Informatics, Kobe University, 1-1 Rokkodai, Nada-ku, Kobe 657-8501, Japan ∥Current address: National Institute of Health Sciences, 3-25-26 Tonomachi, Kawasaki-ku, Kawasaki, Kanagawa 210-9501, Japan E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

Abstract In this study, an ab initio fragment molecular orbital (FMO) methodology was developed to evaluate the solvent effects on electrostatic interactions, which make a significant contribution to the physical and chemical processes occurring in biological systems. Here, a fully polarizable solute consisting of the FMO electron density was electrostatically coupled with an implicit solvent based on the Poisson–Boltzmann (PB) equation; in addition, the nonpolar contributions empirically obtained from the molecular surface area (SA) were added. Interaction analysis considering solvent-screening and dispersion effects is now available as a powerful tool to determine the local stabilities inside solvated biomolecules. This methodology is applied to a deoxyribonucleic acid (DNA) duplex known as the Dickerson dodecamer. We found that excessively large electrostatic interactions inside the duplex are effectively damped by the screening, and the frontier molecular orbital energies also successfully lowered. These observations indicate the stability of highly charged DNA duplexes in solution. Moreover, the solvation free energies in the implicit model show fairly good agreement with those in the explicit model while avoiding the costly statistical sampling of the electrolyte distribution. Consequently, our FMO-PBSA approach could yield new insights into biological phenomena and pharmacological problems via this ab initio methodology.

2


Page 2 of 51

Page 3 of 51 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


1

Introduction

The fragment molecular orbital (FMO) method 1–5 is an efficient and practical quantumchemical approach for computing the electronic states of large molecular systems and accurately predicts whole molecular properties; this technique has enabled us to analyze proteins, nucleic acids, and their complexes in an ab initio manner. The fragmentation into small substructures not only drastically reduces the computational cost but also provides the interaction energy between the fragments as a byproduct, the so-called inter-fragment interaction energy (IFIE) or pairwise interaction energy (PIE). IFIE analysis is a powerful tool for evaluating the local stabilities inside large molecules (see refs 6,7 for a review), and its estimation in vacuo often leads to a substantial overestimation of the electrostatic interactions. 8 After the successful development and applications of the FMO approach, obtaining an understanding of the molecular properties in aqueous solution is the next essential issue for FMO-based analyses because chemical and biological processes related to life phenomena such as molecular recognition, reactions, transportation, regulation, and other activities are universally mediated by solvents around solutes: 9 Solvents can dramatically change the structures and functions of biomolecular systems both directly and indirectly. 10 Including solvent effects in high-level ab initio theory, such as electron-correlated FMO, is, thus, necessary, especially for the accurate prediction of ligand binding affinities 11 in drug discovery, where even small differences in the ligand structure can cause unexpected changes in the pharmacological activity landscape. 12 Because solutes and solvents interact with each other not only via electrostatic forces but also via several quantum effects such as dispersion forces or charge transfer over the solute–solvent interface, 9,10 solvents must be explicitly treated to obtain a correct description of their effects. For this reason, FMO calculations are also performed for each biomolecule as a solute surrounded by explicit water molecules and ions, where the solvent molecules are treated as fully quantum-mechanical (QM) fragments 13–29 or partially classical ones using the effective fragment potential (EFP) model in combination 3



with a molecular mechanical (MM) approach. 30–34 Explicit solvent effects, however, must be statistically investigated from an ensemble obtained from the multiple sampling of energetically feasible states corresponding to each solvent distribution. This results in several difficulties concerning how to generate and sample the distribution, and how to deal with the growth of the computational costs arising from the multiple FMO calculations for the ensemble. The implicit treatment of solvents is an alternative approach to cope with these difficulties by including the electrostatic effects as a statistical average. Implicit solvents can be incorporated into FMO calculations in two ways: Continuum dielectric models, such as the polarized continuum model (PCM) 35 or Poisson–Boltzmann (PB)-based models, 36,37 or the three-dimensional reference interaction site model (3D-RISM) based on the integral equation theory for molecular liquids; 38,39 the former models are the focus of this investigation and are used to estimate the electrostatic contribution of the solvents. One of the most commonly used continuum dielectric models employed in QM approaches is PCM, which was developed for small molecules. There are several variations of the original PCM model. 35 Fedorov et al. interfaced the conductor-like screening model (COSMO) 40 with the FMO method, 41–46 and this approach has been applied to several biological and biochemical targets. 26,47–51 Furthermore, a practical definition of the fragment interaction energy in the PCM solvent has been suggested, 52 and the electrostatic interactions are effectively reduced by solvent screening in their estimation. On the other hand, PB-based models, even including ionic effects in accordance with the Boltzmann distribution, are popularly used for macromolecules modeled using classical treatment 53–55 to study, for example, the molecular surface potential, 53,56 Brownian dynamics, 36,57 the electrostatic solvation free energy, 58,59 or pH dependency. 59,60 In particular, the solvation free energy is often estimated in combination with a nonpolar contribution empirically derived from molecular surface area (SA), and this methodology is known as the MM-PBSA approach, 61,62 as well as the MM-GBSA based on the generalized Born (GB) model, 63,64 which mimics the numerical solution of the PB equation, and is used as the de facto stan-

4


Page 4 of 51

Page 5 of 51 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


dard tool for evaluating protein–ligand and protein–protein binding affinities and nucleic acid stabilities. 55,65–68 Several programs for solving the PB equation have been developed so far with sophisticated numerical solutions and discretized formulations of the PB equation, 69,70 for example, DelPhi, 71–74 Uhbd, 75–77 Apbs, 78–80 Pbsa, 81,82 Pbeq, 83 Zap, 84,85 Mead 86 and Mibpb. 87 The PB solvent model has also been employed in quantum chemistry, coupled with the electron density obtained from molecular orbital (MO) theory for small molecules, and applied to the electronic energy 88–91 and analytical gradient 92 calculations. As for applications to large biomolecules, there have been some attempts, although with limitations, such as the fully polarizable but semiempirical approach 93,94 or ab initio but partially polarizable approach with QM/MM hybridization. 95,96 Watanabe et al. first incorporated PB solvent effects into ab initio FMO calculations for small proteins using Hartree–Fock (HF) theory 97 by employing the FMO-processing program Abinit-mp 3,5,98,99 and DelPhi, which were linked via scripts with input/output (I/O) file handling. 100 In this framework, a solute charge distribution given by the FMO electron density of the whole molecule is self-consistently polarized through the electrostatic interaction with the induced surface charges updated by solving the PB equation to obtain the FMO-PB electronic energy. To integrate the PB solvent model with the FMO method more comprehensively, we have developed a PB equation solver and implemented it into the Abinit-mp program, as well as extended the approximation level to second-order Møller–Plesset perturbation (MP2) theory. 101 These improvements enable us an IFIE analysis considering solvent-screening 52 and dispersion effects in solution. Moreover, by adding the empirical SA energy as a nonpolar term to FMO-PB electronic energy, the FMO-PBSA energy is finally obtained as a total form of the implicit solvent model used here. Note that the FMO-PBSA approach using a fully polarizable medium 28,29,102 is clearly distinct from a practical combination of the FMO electronic energy in vacuo and the MM-PBSA solvation free energy with a fixed solute charge distribution. 103–105

5



In the following sections, we first present the fundamental theories, implementation, and background of the FMO-PBSA methodology. We then apply the methodology to a deoxyribonucleic acid (DNA) duplex of 12 base pairs, the so-called Dickerson dodecamer, 106 to investigate its electrostatic character, fragment interaction energy, molecular orbital energy, and solvation free energy in the implicit model. Moreover, we validate the availability and efficiency of the FMO-PBSA calculation by comparing them with those obtained from the explicit model.

6


Page 6 of 51

Page 7 of 51 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


2

Methodology

Here, we summarize the theoretical background of the FMO method with a PB solvent model calculation. First, the methodology of the ab initio FMO calculation is introduced to obtain electronic properties of a whole molecular system with its associated interaction analysis, and the integration with the PB solvent model is shown. Next, details of the computational settings and datasets used in this study are described.

2.1 2.1.1

FMO Method in Cooperation with the PB Solvent Model Theoretical Basics

To solve the electronic state of the solute in the FMO method, 1–5 a large molecular system is divided into small fragments such as amino acid residues for proteins and nucleic acids for DNA. Alternatively, the units can comprise their subdivided functional groups. 102,107 Fragments comprise not only monomers but also dimers or multiple combinations, and these polymers are treated as the targets of the MO calculation in the FMO method. First, all monomer densities are iteratively calculated through self-consistent charge (SCC) procedure, where the electronic energy of each monomer is solved by self-consistent field (SCF) calculations including the environmental electrostatic potential (ESP) from the surrounding fragments at every step. Each dimer or higher-order polymer is additionally solved under the embedded potential generated by the other fragment densities converged in the monomer SCC according to the order of the many-body expansion. In the pairwise approximation of the FMO (abbreviated to FMO2), the total electronic energy of whole a system in vacuo is derived from the electronic energies (EX ) of the fragment

7



Page 8 of 51

monomers (X = I, J, K, ...) and dimers (X = IJ, JK, IK, ...) 3–5

E FMO =

∑ I

=

∑

EI +

=

I

(EIJ − EI − EJ )

I>J

∑[ ( )] ′ EI′ + EIJ − EI′ − EJ′ + Tr ∆P IJ V IJ

I

∑

∑

I>J

EI′ +

∑

e′ , ∆E IJ

(1)

I>J

where P X and V X are the electron density and the environmental ESP imposed on each frag′ ment, respectively. EX shows the internal energy with the embedded potential contribution

subtracted in the form of ( ) ′ EX = EX − Tr P X V X ,

(2)

e ′ yields the IFIE defined by and ∆E IJ ( ) e ′ = E ′ − E ′ − E ′ + Tr ∆P IJ V IJ ∆E IJ IJ I J

(3)

under ∆P IJ = P IJ − P I ⊕ P J . These formulations are easily extended to those in the framework of MP2 perturbation theory, for example, by adding the electronic correlation ′ 6,7 . energy to EX

2.1.2

Formulation in the Reaction Field

In the PB implicit solvent model, the external field potential from the induced charges on the molecular surface derived from solving the PB equation is additively imposed on the FMO calculation above for the solute. The computational formalism for solving the PB equation is represented in Section S1 and results in the generation of M points of induced charge, s qm , at position r sm (m = 1, ..., M ) on the solvent-excluded molecular surface (SES) 108,109

(see Figure 1). The induced charges increase the coulombic electrostatic solute–solvent interactions, which are given by a contribution of one-electron integrals, W , to the Fock

8


Page 9 of 51 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


matrix 41 X Wµν

M ⟨ ∑ qs =− µ ′ m s r −r

m

m=1

⟩ ν

(4)

with indices of atomic orbital µ, ν and electron coordinate r ′ . Considering the interaction with the nuclear charge ZA on atom A, the electrostatic solute–solvent interaction energy on a fragment is then given by [ ] N ∑ M s ∑ ( ) 1 Z q A m es EX = Tr P X W X + , s | 2 |r − r A m A∈X m=1

(5)

where N is the number of atoms belonging to the fragment. Accordingly, when the internal ′ ′′ es energy of each fragment can be divided as EX = EX + EX , the total FMO-PB solute energy

with the solute–solvent interaction is defined as

E FMO-PB =

∑

EI′ +

I

=

∑ I

+

∑[ )] ( ′ EIJ − EI′ − EJ′ + Tr ∆P IJ V IJ I>J

EI′′ +

∑

∑[ ( )] ′′ EIJ − EI′′ − EJ′′ + Tr ∆P IJ V IJ I>J

EIes +

I

=

∑(

∑ I>J

EI′′

+

es EI(I)

I

es (EIJ − EIes − EJes )

)

] ∑[ ′′ es es es e ∆EIJ + EI(J) + EJ(I) + ∆EIJ , +

(6)

I>J

where the electrostatic solute–solvent interaction energy on fragment I, EIes , is decomposed into each partial contribution with induced charges on surface belonging to fragment J es 6 denoted by EI(J) (see Figure 2) as

es EIes = EI(I) +

∑

es EI(J) .

(7)

I̸=J

e ′′ , is represented in the same form as eq 3: The internal IFIE, ∆E IJ ( ) e ′′ = E ′′ − E ′′ − E ′′ + Tr ∆P IJ V IJ , ∆E IJ IJ I J 9


(8)


Page 10 of 51

es and ∆EIJ shows the coupling term with induced charges arising from charge transfer in

dimer IJ ) 1 ( es es ∆EIJ = EIJ − EIes − EJes = Tr ∆P IJ W IJ . 2

(9)

On the other hand, the corresponding SA energy is independently written in the same way,

E SA =

∑

EInp +

I

=

∑

∑

np (EIJ − EInp − EJnp )

I>J

EInp +

∑

np , ∆EIJ

(10)

I>J

I

where a nonpolar contribution of each fragment is given by the relationship

np EX = γSX

(11)

for solvent surface tension, γ, and solvent accessible surface (SAS) area 110,111 belonging to a fragment, SX (see Figure 1). Finally, a complete form of the total FMO-PBSA energy in the pairwise approximation is given by combining the FMO-PB and SA energies

E FMO-PBSA =E FMO-PB + E SA ] ∑( ) ∑[ np np ′′ es ′′ es es es e = EI + EI(I) + EI + ∆EIJ + EI(J) + EJ(I) + ∆EIJ + ∆EIJ I

=

∑(

) solv

EI′′ + EI

I

=

∑ I

EIPBSA +

∑

+

∑(

I>J

e ′′ + ∆E solv ∆E IJ IJ

)

I>J PBSA ∆EIJ ,

(12)

EIPBSA =EI′′ + EIsolv ,

(13)

I>J

where

PBSA ′′ solv eIJ ∆EIJ =∆E + ∆EIJ

10


(14)

Page 11 of 51 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


represent the solvation monomer energy and IFIE with solvent screening by including their solvent components for the monomer energy and IFIE

es EIsolv =EI(I) + EInp , np solv es es es ∆EIJ =EI(J) + EJ(I) + ∆EIJ + ∆EIJ ,

(15) (16)

respectively. The difference in the total solute energies in vacuo and in the solvent gives, consequently, the solvation free energy. 88 Finally, it is worth discussing the origin of the solvent screening for IFIE in the implicit solvent model. When both fragments I and J have positive net charges, for example, their e ′′ , is electrostatically repulsive. At the same time, the solvent surface bare interaction, ∆E IJ charges are correspondingly induced around fragment J with negative charge, so their interes actions with fragment I, EI(J) , represent a totally attractive interaction and vice versa for es (Figure 2). These partial solute–solvent interactions screen out the bare interactions EJ(I)

and, thus, the IFIE in solvent is effectively relaxed. 6

2.1.3

Grand Iteration of Solute–Solvent Coupling

The FMO-PB solute energy is derived with self-consistency between the solute and inducedsurface charge distributions through grand iterations under their direct coupling. The PB equation is solved under the solute charge and dielectric constant distributions to determine the ESP field, whereas the solvent charges induced on the molecular surface are determined to reproduce the ESP field and are embedded in the next FMO calculation step to redistribute the solute charges under the electric induction as a solvent effect. Through the grand iteration, the HF density is used to obtain a solute charge distribution and HF energy to evaluate the convergence. The solute charge distribution is determined by an ESP-fitting scheme 112–115 within an FMO calculation and passed to the PB equation. In this paper, the Merz–Kollman scheme, 112,113 one of the most popular methods, is employed to

11



assign atomic charges to each atom in a solute molecule. Here, the atom-centered charges are fitted to each fragment-derived ESP and aggregated as a whole molecular charge distribution with the FMO convention; 116 that is, they are not fitted to whole molecular electrostatic potential (MEP) directly. 117,118 By using the converged MOs of the monomers obtained in the previous step as the initial guess in the next step, the convergence of the SCC procedure is improved substantially. At the same time, the MP2 electron correlation energies are obtained only in the initial gaseous and converged solvent phases to estimate the dispersion effects in the fragment-based interaction analysis by using the corresponding HF references; the Cholesky decomposition with adaptive metric (CDAM) approximation 119 were used for the time-consuming twoelectron repulsion integrals to efficiently accelerate the electron correlation calculations. A nonpolar contribution is then estimated from the SAS area of each fragment and its pair according to eq 11 after completion of the entire electronic calculation; this FMO-PBSA computational scheme (see Figure 3 in summary) has been implemented into Abinit-mp 120 in Fortran.

2.2 2.2.1

Computational Details Parameters for Solving the PB Equation

The PB equation in the finite difference approach is solved with the following configuration. The size of a solvated box is defined so that a solute occupies 60% length at the center of each side of the box with a grid spacing of 0.2 ˚ A. An exterior dielectric constant, εout , of 80 was used for the water solvent, which is depicted as a spherical probe with a radius of 1.4 ˚ A and no ionic strength. A set of atomic radii 88 shown in Table 1 was employed with the solvent probe radius for defining the SES as the dielectric boundary between the solute and solvent. These atomic radii are also used to define the SAS area for the estimation of the ˚2 ). The algebraic equation nonpolar energy with a surface tension, γ, of 0.0072 kcal/(mol A shown in eq S6 uses the successive over relaxation (SOR) method 72 and is iteratively solved 12


Page 12 of 51

Page 13 of 51 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


until the root mean square error of the ESP of the entire lattice space is converged to within 10−5 a.u. between two successive steps.

2.2.2

FMO Configuration and Related Analysis

The FMO calculation on the solute is performed with the following configuration. We impose some ESP approximations for efficient computation between two sufficiently-separated fragments. 3,5 The electron densities of fragments in near and far regions, which provide the embedded potential for an SCF calculation of each fragment, are approximately replaced by the Mulliken atomic orbital charge (esp-aoc) populations and further-reduced point charge (esp-ptc) populations, respectively. 5 At the same time, the dimer electronic energy for separated fragments is also approximately estimated by summing the corresponding monomer electronic energies and their electrostatic coulombic interaction energies (dimer-es). 3 Each approximation is applied for each fragment pair whose distance between the nearest atoms is larger than the product of a multiplier constant and the contact distance, which is the summation of van der Waals radii of the atoms. The multipliers are set to zero for the ‘esp-aoc’ approximation and two for the ‘esp-ptc’ and ‘dimer-es’ approximations, as usual. In this study, the grand iteration where solute and surface charge distributions were alternately updated in each step and fully repeated until the total energy of the solute was converged to be within 10−5 a.u., although a practical truncation of the iteration has also been proposed to reduce the computational cost while including the many-body effects efficiently. 41 All the electronic energies and densities were obtained with the 6-31G* basis set. The atomic charge distribution using natural population analysis (NPA) 121,122 was also computed based on the HF density to estimate the amount of charge transfer.

13



2.3

Datasets

We applied the FMO-PBSA methodology to a realistic biomolecule, a DNA duplex, to estimate its unbound properties in solution. DNA duplexes are carriers of hereditary information and the genetic codes comprise four nucleobases, adenine, guanine, cytosine, and thymine, which form a sequence of Watson–Crick (WC) base pairs attached to antiparallel sugar–phosphate (SP) backbones with highly negative charges. 123 In this work, we employed the same DNA duplex structure as that used in the previous studies considering explicit solvent effects; 28,107 the duplex is in the B-type conformation, known as the Dickerson dodecamer, 106 d(CGCGAATTCGCG)2 , as shown in Figure 4. The model is constructed from a crystal structure (PDB ID: 355D 124 ) and solvated with explicit water molecules and a sufficient number of counterions (sodium cations (Na+ )) for charge neutralization. The whole system was annealed by the AMBER99 force field. 125 See ref 28 for more details of system preparation. Finally, the solvent molecules were removed and the remaining DNA duplex was used as the solute in the FMO-PBSA calculations. For comparison with the explicit solvent model, DNA configurations solvated by explicit solvent shells with different solvent layer thicknesses were prepared from ref 28, as discussed later. The DNA duplex is divided into fragments per nucleobase, which includes nitrogenous bases with purine or pyrimidine scaffolds, and the SP backbone, which consists of the fivecarbon 2’-deoxyribose sugar and a phosphate group. They are separated at the N -glycosidic bond in each nucleotide unit and the C4’–C5’ bond on the backbone (Figure 5).

14


Page 14 of 51

Page 15 of 51 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


3

Results and Discussion

To evaluate our FMO-PBSA methodology, we applied it to estimate the influence of the solvent on the electronic properties of the DNA duplex. These properties include the internal polarization, local stabilization, frontier molecular orbitals, and solvation free energy in comparison to those obtained used an explicit solvent model. 28,107

3.1

Electrostatics

Electrostatic interactions make a significant contribution to the specific recognition or binding between molecules and this depends on the molecular electrostatic properties including the internal charge distributions and externally generated ESPs. Solvent greatly affects these properties through the electrostatic induction of solute molecules. Thus, we first investigated the changes of molecular electrostatic properties resulting from solvent effects. As shown in Figure 6a, the ESP generated by a solute charge distribution in vacuo is mapped on the molecular surface defined by the isoelectron density equal to 0.001 a.u. Note that an ESP-fitted charge set that reduces a spatially distributed electron density to points on nuclei is here used as the solute charge distribution because the potentials generated by different electron densities in vacuo and solvent are compared on the same molecular surface. Because each phosphate group constituting the DNA backbone, PO4 – , has a negative charge of −1e, the negative potential is broadly distributed on the molecular surface of the duplex with a total charge of −22e (cf. the ESP study for an estrogen response element of DNA with −32e in ref 126). The difference in the potential under the reaction field generated by PB solvent model is shown in Figure 6b. Because of the solvation, the intermediate and both end regions on the sequence of the duplex are, in total, polarized to be negative and positive as a difference, respectively. That is, globally, electrons move from the edges to the middle of the sequence. In particular, the negative polarization is emphasized inside the major groove of the duplex.

15



The net charge variation in each fragment is estimated from NPA measurements. The distribution has mirror symmetry on the forward and reverse strands per fragment type (Figure S1a). Focusing on the base fragments of each strand, electrons are likely to gather on the guanine 28 or adenine bases; focusing on the SP fragments, electrons are assigned to 3’end side of the backbone with some oscillation. The summation of the charges on each base pair also supports the electron transfer to the intermediate region mentioned above (Figure S1d). The fluctuations at both ends of the duplex are consequently increased to modulate the fragment interactions, as discussed later. Because the fluctuation does not arise in the explicit treatment, 28 we recognize it as an edge effect originating from a polarized strain peculiar to the implicit solvent model, where electrons are not relaxed by outflow into the solvent. Finally, the induced solvent charge distribution as a source of the reaction field is shown in Figure 6c; the electrostatic induction opposite in sign to the molecular charge distribution stabilizes the molecular polarization of the solute. Overall, the distribution reflects the local electrostatic properties, such as the positive charge around phosphate groups or negative charge around amine groups, for example, showing a clear patch-like pattern. Note that the solvent charge induction does not necessarily correspond to the molecular potential variation, which represents the global polarization of the solute.

3.2

Fragment Interaction Energy

The highly structured DNA duplex is realized by WC hydrogen bondings between complementary nucleobases and π-stackings between their piles, and the charged SP groups support these interactions in the backbone while repelling each other. One of our concerns is how and to what degree the solvation affects the local stability between such partial moieties inside the DNA duplex. The overall distribution of IFIEs for all the fragment pairs inside the 12base-paired DNA duplex is shown in Figure 7 as a function of the interfragment distance by comparison of the results in vacuo with those in the implicit solvent. The distance is defined 16


Page 16 of 51

Page 17 of 51 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


as the separation of the nearest atoms belonging to different fragments. These pairs are classified into six categories: WC hydrogen-bonding base pairs, intrastrand-stacking base pairs, interstrand-stacking base pairs, other base pairs, base and SP backbone pairs, and other SP backbone pairs. For each type of proximate base pair included in the three former categories, the IFIEs are mapped as a peak spectrum of distances on the plots in Figures 7a–c with moderately stabilized on the whole by solvation. For these pairs, a detailed breakdown of the energetic components of each IFIE is shown in Figure 8. As an overall trend of the IFIEs, the solvation affects the HF energy component the most (by several kilocalories per mole), but the solvation affects the MP2 correlation energy little. For the WC hydrogen-bonding pairs, all the IFIEs of A–T and G–C pairs in the middle region of the duplex show few changes in the magnitude; however, those of G–C pairs closer to the termini show more ambiguous and larger fluctuations (Figure 8a). The intrastrand stacking interactions mainly consist of the MP2 correlation energy, and this is gradually strengthened towards the intermediate region; in contrast, this is weakened towards the termini (Figure 8b). This trend corresponds to the ESP variation in Figure 6b, where electron transfer to the middle region is induced by the solvation. On the other hand, the polarization has influence only on the HF energy component of IFIE and hardly on the MP2 correlation component. The HF interactions are originally repulsive or attractive in vacuo are electrostatically relaxed by the solvent-screening effect. Therefore, the electron enrichment in the intermediate region directly arises not in the π-orbitals but in the solvent-exposed orbitals of the lone pairs (see Figure S2), and no enhancement in the dispersion force between intrastrand-stacking base pairs occurs (cf. ref 127). The same can be said for the interstrand-stacking base pairs, where almost all the pairs have positive IFIE values for each HF component in vacuo and the solvation mitigates this repulsion (Figure 8c). This IFIE relaxation in solvent can be explained by a screening effect, including the attractive interactions of lone pairs in a base with a positive surface charges induced around the lone pairs in another base and vice versa.

17



The screening increases especially for base pairs that have a large number of lone pairs on the sides exposed to the grooves of a duplex (e.g., guanine). Both the penultimate G2–G3’ and G10–G11’ base pairs from the ends are stabilized by more than 5 kcal/mol; however, this may also include the edge effects mentioned above. It would be worth comparing solvent effects on these interactions with an explicit solvent model of the same DNA duplex covered by a solvent shell. Here, we used an explicit model with a 12-˚ A solvent shell, which is the maximum in ref 28. For WC hydrogen bonding, the IFIE values of A–T base pairs are steady and those of G–C all drop by about 3 kcal/mol when the duplex is immersed in the explicit solvent. This stabilization is explained as the result of a modulation in electrostatic interactions between the base pairs. 28 However, the modulation is due to the loss of electrons that flow out into solvent to reduce the charge–charge repulsion. In principle, such a picture can never appear in the implicit solvent model. In this manner, the difference in the restriction of electron relaxation in each solvent model results in some inconsistency in the response to the solvent. For example, the polarization of each WC hydrogen-bonding base pair arises not in the WC hydrogen-bonding atom pairs in the explicit model but in the solvent-exposed polar atoms including lone pairs in the implicit model (Figure S2). Meanwhile, the total intrastacking and interstacking arrangements between each two base pairs are in turn destabilized, unlike those in the implicit solvent, which are overall stabilized (Figure S3). In contrast with these inconsistencies, the stabilization and destabilization mainly arise from the modulation of electrostatic interactions, 28 which is a factor common to both solvent models. Thus, the difference in the IFIE fluctuation between the explicit and implicit models can be summarized as follows. In the explicit model, the electron outflow to the solvent, especially through guanine, brings about a relaxation in the repulsion in the WC hydrogen bonding and a destabilization of π-stackings. 28 In the implicit model, on the other hand, a screening effect basically introduces a relaxation of the electrostatic interactions. The IFIEs of other base pairs illustrate a discrete pattern of distances on the plot in

18


Page 18 of 51

Page 19 of 51 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


Figure 7d. Because each of these contains mutually uncharged fragments, the screening effect yields a relaxation of only 2 kcal/mol at most, even for the close pairs, and has little effect on the separated pairs whose IFIE is thoroughly damped in vacuo. For the base and SP backbone pairs as interactions between neutral and charged fragments, the overall IFIEs are broadly distributed, as in Figure 7e. The energies of the proximate pairs, which range from −20 to 10 kcal/mol in vacuo, are moderately relaxed to half of this range after solvation. On the other hand, those of the distant pairs are partly scattered toward the attractive directions by approximately by 5 kcal/mol. This is because of the attractive interaction between the SP backbone and positive surface charges induced around the base lone pairs, which contributes to the stabilization. Finally, the distribution of IFIEs for the SP backbone pairs in vacuo can be clustered into two groups in Figure 7f. One is a group consisting of the pairs of mutually charged fragments, and the repulsive IFIEs are mapped on a curve that decays with increasing distance. In this group, the IFIE values ranging 10–35 kcal/mol in vacuo are successfully screened out by the solvation and reduced to 5 kcal/mol or less, being nearly zero in the far region. The picture is consistent with IFIE between solvated two ions in ref 6. The other is a group consisting of the pairs including an uncharged fragment at the ends of a duplex. These interactions result in IFIE values of around zero in vacuo, but all become slightly repulsive after solvation. Because a conspicuous coulomb interaction between the charged fragments remains in explicit model no matter how many solvent molecules is included, an intuitive analysis of the IFIEs would be advantageous in the implicit model.

3.3

Frontier Molecular Orbital Energy

The frontier MO energies consisting of the highest occupied and lowest unoccupied MOs, HOMO and LUMO, and the difference between them, the so-called band gap, are deeply related to the reactivity, conductivity, and stability of the DNA duplex because it is a highly-charged molecular system. 28 As is apparent from previous studies for the Dickerson 19



dodecamer, 28,94 quantum chemical calculations in vacuo yield the same result; that is, the DNA duplex is unstable with a positive HOMO energy, whether treated as a whole molecular system or divided into fragments. Thus, a proper description of the frontier MOs is crucial for understanding the biological mechanism for the damage, protection, and engineering of functional materials. 94 Figure 9 summarizes the mean changes in the frontier MO energies on each fragment type of a 12 base-paired DNA duplex between the in vacuo and in implicit solvent cases. The details are shown in Table S1, where each entry is averaged over the same type of fragments excluding those in terminal base pairs and shown with its standard deviation. As shown, the HOMO energies for all the fragment types decrease to negative values when the reaction field is imposed on the duplex; this, in short, guarantees that the molecule is stable. For example, the HOMO energy, which is 18.21 eV for the adenine base fragment in vacuo, decreases by 22.85 eV and becomes −4.64 eV in solvent. Similarly, in the case of the LUMO with the same decrement, the MO energy of 27.86 eV in vacuo falls to 6.40 eV after solvation, and the band gap between the HOMO and LUMO energies is kept maintained, irrespective of the influence of the external field. These features are common to other fragment types, whereas, with respect to the SP backbone, the band gap widens by 0.74 eV because the decrease in the MO energy for the HOMO is larger than that for LUMO. In addition, for the fragment types of guanine, cytosine, and the SP backbone, which receive edge effects because of their closeness to the termini of the duplex, the frontier MO energies in vacuo have large deviations reduced by approximately 25% in solvent. Here, we make a comparison with the explicit solvent model in ref 28. By including a sufficient number of water molecules and counterions, the frontier MO energies on each nucleobase and SP backbone fragment are successfully reduced to be comparable to those of the implicit model. From the similarity of the influence on the MO energies, we found that the induced surface charges of the implicit model have a similar effect to the counterions of the explicit model. These results are also consistent with the previous study involving the

20


Page 20 of 51

Page 21 of 51 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


semi-empirical QM Hamiltonian and PB implicit solvent model, where the HOMO energy of the whole DNA duplex changed from 15 to −8 eV on solvation. 94 The frontier MOs on a cytosine base and its associated SP backbone fragments, constituting deoxycytosine at the ninth position of the forward strand (dC9), are illustrated in Figure 10 with a comparison between the in vacuo and implicit solvent cases. In the cytosine base fragment, the population of the HOMO on the pyrimidine ring in vacuo is somewhat transferred onto the amine group in solvent, and a p-orbital vanishes from the nitrogen atom at the 3-position of the pyrimidine ring, while another appears the nitrogen atom of the amine group. The corresponding LUMO, in contrast, shows little change. Similarly, in the SP backbone fragment, the electrons move out from the phosphate group to the sugar ring, and the p-orbitals disappear from the two oxygen atoms of the phosphate group. The LUMO still shows little change, although some of the electron density in the bond between the two SP backbone fragments moves from the bond to the cytosine base fragment. These responses of the frontier MOs to the reaction field are also consistent in their shape with those in explicit model with counterions (see Figure S4), whose effect on the MOs of the modeled base pairs is discussed in detail in ref 107.

3.4

Solvation Free Energy

To investigate the solvent effects of a whole DNA duplex energetically, the enthalpic contribution to the solvation free energy in the implicit solvent model was compared with those in the explicit solvent model where the solute is covered by explicit solvent shells with different thicknesses and numbers of water molecules and counterions. The definition of the solvation free energy is given in Appendix A.1, where it consists of the internal solute energy, E in , as the difference between the energy in vacuo and in solvent, ∆E in , and the solute–solvent interaction energy, E es . Note that E es in the explicit solvent model is the sum of solute–solvent IFIEs by half in the definition, while the IFIEs are sometimes used as is when not using the concept of solute free energy. 13,22,28,29 The values 21



of these two types of energies, ∆E in and E es , were obtained from the previous study 28 using an explicit solvent model with a solvent thickness of 0–12 ˚ A. Here, E es does not seem to fully converge with a shell thickness of 12 ˚ A yet, whereas ∆E in is already converged (see Table S2). To achieve convergence, we prepared 13, 14, and 20-˚ A shell models of the same snapshot in ref 28 and carried out regular FMO calculations for them; the solvated DNA system is completely neutralized with 22 Na+ ions and a thickness of 20 ˚ A. The solvation free energies in the explicit solvent models above are plotted with the number of Na+ ions as a function of the shell thickness in Figure 11. Through neutralization, the DNA duplex is internally destabilized by 555.5 kcal/mol because of its polarization and the outflow of electrons to solvent molecules of around −7e. 28 Meanwhile, the solute–solvent interaction becomes increasingly attractive, corresponding to the number of counterions in the solvent; thus, the DNA duplex is totally stabilized. On the other hand, the solvation free energy in the implicit model, whose definition is given in eq 5, is also drawn in Figure 11. The induced surface charges in the implicit model impose a reaction field on the solute and also polarize it. There are no explicit water molecules or counterions to which electrons of solute can outflow in the implicit model; however, the solvation free energies of both models finally approach each other. These values are also comparable with that of the semi-empirical MO calculation for the whole duplex, −6235.8 kcal/mol. 93 To obtain a more reliable estimation of a solvation free energy, multiple solvated configurations of a solute via multiple sampling are required for an explicit treatment because there is a fluctuation of about ±150 kcal/mol in the solvation free energy, even for a neutral solute. 13 In addition, because the energetic contribution of the solute–solvent interaction is significant, neutralization with a solvent including counterions, which requires an increase in the thickness of a solvent shell, is essential for a charged solute. Considering these costs, the PB solvent model is a good solution for conveniently estimating the solvation free energy, even targeting a highly charged molecule such as a DNA duplex.

22


Page 22 of 51

Page 23 of 51 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


4

Conclusions

An FMO methodology coupled with PB implicit solvent model was developed to analyze the electronic properties of large biomolecules in solution. We applied the FMO-PBSA methodology to a DNA duplex and estimated the influence of the solvent reaction field potential imposed on the solute by comparison with those of the previous explicit solvent model. In the fragment-based interaction study, solvent screening induced by the interaction with the induced surface charges effectively reduced the excessive electrostatic interactions between charged fragments. Such an electrostatic relaxation also often occurs between proximate nucleobase pairs depending on the number of lone pairs exposed to the solvent while hardly affecting the electron-correlated interactions. In the frontier MO study, the energy levels of both of the HOMO and LUMO on each fragment are perturbed by the reaction field and successfully shifted down to those guaranteeing stable electronic states. The pictures along with the shape of the MOs show a good agreement with those in the explicit treatment. At the same time, concerning the responsiveness of the entire solute, the solvation free energy obtained in explicit solvent asymptotically approaches that of the implicit solvent through its neutralization with increasing explicit solvent shell thickness. These findings indicate that the induced surface charges work as an explicit electrolyte including counterions without any costly statistical sampling of the solvent configurations. On the other hand, some inconsistencies with the explicit model certainly exist because electron flow cannot occur over the solute–solvent interface after the polarization of solute in the implicit model. This may cause a slightly different interpretation of the response to the solvent on the internal charge distribution or interaction mechanisms, for example, with ambiguous edge effects. Despite the limitations, the FMO-PBSA calculation provides a reliable description of not only local electronic states and interactions but also the whole molecular properties, even for highly charged systems. These merits should be emphasized in practical use. Because, in explicit treatment, addition of solvent molecules until convergence of the electronic properties of a solute results in increasing the number of fragments enormously; 23



this often causes serious problems in writing out a large amount of results, depended to network and disk I/O performance, or formidable obstacles in processing and visualizing such large size data in analysis. The FMO-PBSA approach would also have great potential in relative free energy and IFIE-based analyses for proteins and their binding systems such as protein–ligand, protein–protein, and protein–DNA complexes. These verifications will be carried out in the future.

Supporting Information Available In Section S1, the essential formulation and the numerical solution of the PB equation using the finite difference method is briefly reviewed as the implicit solvent model implemented in this work. Section S2 contains the additional figures and tables listed below: Figure S1 shows the amount of charge transferred via solvation, focusing on the fragment unit of each base pair. Figure S2 shows atomic charge distribution in vacuo and the difference with that in solvent for C3–G3’ and A6–T6’ base pairs. Figure S3 gives the summation of the IFIEs on stacking interactions between each base pair individually, as shown in Figure 8. In Figure S4, exactly the same comparison of the frontier MOs on fragments of dC9 in vacuo and in solvent is illustrated by using the explicit solvent. Table S1 shows a detailed breakdown of the frontier MO energies shown in Figure 9. Table S2 also shows a detailed breakdown of the solvation free energies shown in Figure 11.

Acknowledgement The authors thank Prof. Yuji Mochizuki at Rikkyo University for supplying computational resources for development and validation. We thank Dr. Yuto Komeiji at the National Institute of Advanced Industrial Science and Technology (AIST) for providing the structure of the DNA duplex. We also thank Dr. Hirofumi Watanabe at Kobe University for fruitful discussion on the implementation of FMO-PBSA. This research was performed in “Research 24


Page 24 of 51

Page 25 of 51 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


and Development of Innovative Simulation Software (RISS)” project supported by Research and Development for Next-Generation Information Technology of Ministry of Education, Culture, Sports, Science, and Technology (MEXT).

25



Page 26 of 51

Appendix A.1

Solvation Free Energy Definitions in Explicit and Implicit Solvent Models

Here, we define the solute free energy in explicit/implicit solvent, E solute , as the summation of the internal solute energy component, E in , and the electrostatic solute-solvent interaction in component, E es , as the enthalpic contribution. In the explicit model, Eexplicit is obtained

from the energetic summation in FMO convention with limited fragment space to the solute es alone, 13,22,28 whereas Eexplicit is obtained by summing up the IFIEs by half for all solute–

solvent fragment pairs because their interaction energies must be equally shared with the solvent fragments. 128

solute in es Eexplicit = Eexplicit + Eexplicit ∑ ∑ in EI′ + Eexplicit =

e′ ∆E IJ

(18)

I>J I,J∈solute

I∈solute es Eexplicit =

(17)

∑ I∈solute J∈solvent

1 e′ ∆EIJ 2

(19)

in Note that any kind of energy type (either HF or MP2) can be adapted to estimate Eexplicit , es whereas only the HF energy is conveniently used here for Eexplicit because of the assumption

of an electrostatic interaction between the solute and solvent; several quantum effects other than electrostatic effects are also included in the proximity interaction, even if HF is used. On the other hand, the two energy components above are explicitly separated in the es implicit model (see eq 6) and Eimplicit already comprises the bisection concept 129,130 in the es definition of EX , as in eq 5. Thus, these are simply written down in FMO formulae as

26


Page 27 of 51 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


follows.

solute in es Eimplicit = Eimplicit + Eimplicit ≡ E FMO-PB ∑ ∑ in e ′′ Eimplicit = EI′′ + ∆E IJ I

es = Eimplicit

∑ I

(20) (21)

I>J

EIes +

∑

es ∆EIJ

(22)

I>J

Finally, the solvation free energy, G, can be defined as the difference in the solute free energy in vacuo and in solvent, 93 E0solute and E solute , respectively, and written as G = E solute − E0solute = E in − E0in + E es = ∆E in + E es

(23)

where ∆E in = E in − E0in and E0solute = E0in because the solute in vacuo has no interaction with the solvent.

27



References (1) Kitaura, K.; Sawai, T.; Asada, T.; Nakano, T.; Uebayasi, M. Pair interaction molecular orbital method: An approximate computational method for molecular interactions. Chem. Phys. Lett. 1999, 312, 319–324. (2) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment molecular orbital method: An approximate computational method for large molecules. Chem. Phys. Lett. 1999, 313, 701–706. (3) Nakano, T.; Kaminuma, T.; Sato, T.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment molecular orbital method: Application to polypeptides. Chem Phys Lett 2000, 318, 614–618. (4) Kitaura, K.; Sugiki, S. I.; Nakano, T.; Komeiji, Y.; Uebayasi, M. Fragment molecular orbital method: Analytical energy gradients. Chem. Phys. Lett. 2001, 336, 163–170. (5) Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment molecular orbital method: Use of approximate electrostatic potential. Chem. Phys. Lett. 2002, 351, 475–480. (6) Fedorov, D. G.; Nagata, T.; Kitaura, K. Exploring chemistry with the fragment molecular orbital method. Phys. Chem. Chem. Phys. 2012, 14, 7562. (7) Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. Electroncorrelated fragment-molecular-orbital calculations for biomolecular and nano systems. Phys. Chem. Chem. Phys. 2014, 16, 10310–10344. (8) Tanaka, S.; Watanabe, C.; Okiyama, Y. Statistical correction to effective interactions in the fragment molecular orbital method. Chem. Phys. Lett. 2013, 556, 272–277. (9) Chaplin, M. F. Do we underestimate the importance of water in cell biology? Nat. Rev. Mol. Cell Biol. 2006, 7, 861–866. 28


Page 28 of 51

Page 29 of 51 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


(10) Wiggins, P. M. Role of water in some biological processes. Microbiol. Rev. 1990, 54, 432–49. (11) Ryde, U.; Söderhjelm, P. Ligand-binding affinity estimates supported by quantummechanical methods. Chem. Rev. 2016, 116, 5520–5566. (12) Stumpfe, D.; Bajorath, J. Exploring activity cliffs in medicinal chemistry. J. Med. Chem. 2012, 55, 2932–2942. (13) Komeiji, Y.; Ishida, T.; Fedorov, D. G.; Kitaura, K. Change in a protein’s electronic structure induced by an explicit solvent: An ab initio fragment molecular orbital study of ubiquitin. J. Comput. Chem. 2007, 28, 1750–1762. (14) Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y.; Ishikawa, T.; Nakano, T. How does an SN2 reaction take place in solution? Full ab initio MD simulations for the hydrolysis of the methyl diazonium ion. J. Am. Chem. Soc. 2008, 130, 2396–2397. (15) Komeiji, Y.; Ishikawa, T.; Mochizuki, Y.; Yamataka, H.; Nakano, T. Fragment molecular orbital method-based molecular dynamics (FMO-MD) as a simulator for chemical reactions in explicit solvation. J. Comput. Chem. 2009, 30, 40–50. (16) Komeiji, Y.; Mochizuki, Y.; Nakano, T.; Fedorov, D. G. Fragment molecular orbitalbased molecular dynamics (FMO-MD), a quantum simulation tool for large molecular systems. J. Mol. Struct. THEOCHEM 2009, 898, 2–7. (17) Kistler, K. A.; Matsika, S. Solvatochromic shifts of uracil and cytosine using a combined multireference configuration interaction/molecular dynamics approach and the fragment molecular orbital method. J. Phys. Chem. A 2009, 113, 12396–12403. (18) Ishikawa, T.; Kuwata, K. Interaction analysis of the native structure of prion protein with quantum chemical calculations. J. Chem. Theory Comput. 2010, 6, 538–547.

29



(19) Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y.; Nakano, T. Does amination of formaldehyde proceed through a zwitterionic intermediate in water? Fragment molecular orbital molecular dynamics simulations by using constraint dynamics. Chem. - A Eur. J. 2010, 16, 6430–6433. (20) Mori, H.; Hirayama, N.; Komeiji, Y.; Mochizuki, Y. Differences in hydration between cis- and trans-platin: Quantum insights by ab initio fragment molecular orbital-based molecular dynamics (FMO-MD). Comput. Theor. Chem. 2012, 986, 30–34. (21) Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y. FMO-MD simulations on the hydration of formaldehyde in water solution with constraint dynamics. Chem. - A Eur. J. 2012, 18, 9714–9721. (22) Ishikawa, T.; Burri, R. R.; Kamatari, Y. O.; Sakuraba, S.; Matubayasi, N.; Kitao, A.; Kuwata, K. A theoretical study of the two binding modes between lysozyme and triNAG with an explicit solvent model based on the fragment molecular orbital method. Phys. Chem. Chem. Phys. 2013, 15, 3646–54. (23) Okamoto, A.; Yano, A.; Nomura, K.; Higai, S.; Kurita, N. Stable conformation of fulllength amyloid-β (1–42) monomer in water: Replica exchange molecular dynamics and ab initio molecular orbital simulations. Chem. Phys. Lett. 2013, 577, 131–137. (24) Watanabe, C.; Fukuzawa, K.; Tanaka, S.; Aida-Hyugaji, S. Charge clamps of lysines and hydrogen bonds play key roles in the mechanism to fix helix 12 in the agonist and antagonist positions of estrogen receptor α: Intramolecular interactions studied by the ab initio fragment molecular orbital method. J. Phys. Chem. B 2014, 118, 4993–5008. (25) Yano, A.; Okamoto, A.; Nomura, K.; Higai, S.; Kurita, N. Difference in dimer conformation between amyloid-β (1–42) and (1–43) proteins: Replica exchange molecular dynamics simulations in water. Chem. Phys. Lett. 2014, 595, 242–249. 30


Page 30 of 51

Page 31 of 51 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


(26) Koyama, Y.; Ueno-Noto, K.; Takano, K. Affinity of HIV-1 antibody 2G12 with monosaccharides: A theoretical study based on explicit and implicit water models. Comput. Biol. Chem. 2014, 49, 36–44. (27) Okamoto, A.; Yano, A.; Nomura, K.; Higai, S.; Kurita, N. Effect of D23N mutation on the dimer conformation of amyloid β-proteins: Ab initio molecular simulations in water. J. Mol. Graph. Model. 2014, 50, 113–124. (28) Fukuzawa, K.; Kurisaki, I.; Watanabe, C.; Okiyama, Y.; Mochizuki, Y.; Tanaka, S.; Komeiji, Y. Explicit solvation modulates intra- and inter-molecular interactions within DNA: Electronic aspects revealed by the ab initio fragment molecular orbital (FMO) method. Comput. Theor. Chem. 2015, 1054, 29–37. (29) Tokuda, K.; Watanabe, C.; Okiyama, Y.; Mochizuki, Y.; Fukuzawa, K.; Komeiji, Y. Hydration of ligands of influenza virus neuraminidase studied by the fragment molecular orbital method. J. Mol. Graph. Model. 2016, 69, 144–153. (30) Nagata, T.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. A combined effective fragment potential–fragment molecular orbital method. I. The energy expression and initial applications. J. Chem. Phys. 2009, 131, 024101. (31) Nagata, T.; Fedorov, D. G.; Sawada, T.; Kitaura, K.; Gordon, M. S. A combined effective fragment potential–fragment molecular orbital method. II. Analytic gradient and application to the geometry optimization of solvated tetraglycine and chignolin. J. Chem. Phys. 2011, 134, 034110. (32) Nagata, T.; Fedorov, D. G.; Kitaura, K. Analytic gradient and molecular dynamics simulations using the fragment molecular orbital method combined with effective potentials. Theor. Chem. Acc. 2012, 131, 1136. (33) Nagata, T.; Fedorov, D. G.; Sawada, T.; Kitaura, K. Analysis of solute–solvent interactions in the fragment molecular orbital method interfaced with effective fragment 31



potentials: Theory and application to a solvated griffithsin–carbohydrate complex. J. Phys. Chem. A 2012, 116, 9088–9099. (34) Steinmann, C.; Fedorov, D. G.; Jensen, J. H. Mapping enzymatic catalysis using the effective fragment molecular orbital method: Towards all ab initio biochemistry. PLoS One 2013, 8, e60602. (35) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999–3093. (36) Davis, M. E.; McCammon, J. A. Electrostatics in biomolecular structure and dynamics. Chem. Rev. 1990, 90, 509–521. (37) Sharp, K. A.; Honig, B. H. Calculating total electrostatic energies with the nonlinear Poisson–Boltzmann equation. J. Phys. Chem. 1990, 94, 7684–7692. (38) Ratkova, E. L.; Palmer, D. S.; Fedorov, M. V. Solvation thermodynamics of organic molecules by the molecular integral equation theory: Approaching chemical accuracy. Chem. Rev. 2015, 115, 6312–6356. (39) Yoshida, N. Efficient implementation of the three-dimensional reference interaction site model method in the fragment molecular orbital method. J. Chem. Phys. 2014, 140, 214118. (40) Klamt, A.; Schuurmann, G. COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc. Perkin Trans. 2 1993, 799–805. (41) Fedorov, D. G.; Kitaura, K.; Li, H.; Jensen, J. H.; Gordon, M. S. The polarizable continuum model (PCM) interfaced with the fragment molecular orbital method (FMO). J. Comput. Chem. 2006, 27, 976–985.

32


Page 32 of 51

Page 33 of 51 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


(42) Chiba, M.; Fedorov, D. G.; Kitaura, K. Polarizable continuum model with the fragment molecular orbital-based time-dependent density functional theory. J. Comput. Chem. 2008, 29, 2667–2676. (43) Li, H.; Fedorov, D. G.; Nagata, T.; Kitaura, K.; Jensen, J. H.; Gordon, M. S. Energy gradients in combined fragment molecular orbital and polarizable continuum model (FMO/PCM) calculation. J. Comput. Chem. 2010, 31, 778–790. (44) Nagata, T.; Fedorov, D. G.; Li, H.; Kitaura, K. Analytic gradient for second order Møller–Plesset perturbation theory with the polarizable continuum model based on the fragment molecular orbital method. J. Chem. Phys. 2012, 136, 204112. (45) Nakata, H.; Fedorov, D. G.; Kitaura, K.; Nakamura, S. Extension of the fragment molecular orbital method to treat large open-shell systems in solution. Chem. Phys. Lett. 2015, 635, 86–92. (46) Nishimoto, Y.; Fedorov, D. G. The fragment molecular orbital method combined with density-functional tight-binding and the polarizable continuum model. Phys. Chem. Chem. Phys. 2016, 18, 22047–22061. (47) Sawada, T.; Fedorov, D. G.; Kitaura, K. Role of the key mutation in the selective binding of avian and human influenza hemagglutinin to sialosides revealed by quantummechanical calculations. J. Am. Chem. Soc. 2010, 132, 16862–16872. (48) He, X.; Fusti-Molnar, L.; Cui, G.; Merz, K. M., Jr. Importance of dispersion and electron correlation in ab initio protein folding. J. Phys. Chem. B 2009, 113, 5290– 300. (49) Jensen, J. H.; Willemoës, M.; Winther, J. R.; De Vico, L. In silico prediction of mutant HIV-1 proteases cleaving a target sequence. PLoS One 2014, 9, e95833.

33



(50) Simoncini, D.; Nakata, H.; Ogata, K.; Nakamura, S.; Zhang, K. Y. Quality assessment of predicted protein models using energies calculated by the fragment molecular orbital method. Mol. Inform. 2015, 34, 97–104. (51) Prato, G.; Silvent, S.; Saka, S.; Lamberto, M.; Kosenkov, D. Thermodynamics of binding of di- and tetrasubstituted naphthalene diimide ligands to DNA G-quadruplex. J. Phys. Chem. B 2015, 119, 3335–47. (52) Fedorov, D. G.; Kitaura, K. Energy decomposition analysis in solution based on the fragment molecular orbital method. J. Phys. Chem. A 2012, 116, 704–719. (53) Honig, B. H.; Nicholls, A. Classical electrostatics in biology and chemistry. Science 1995, 268, 1144–1149. (54) Fogolari, F.; Brigo, A.; Molinari, H. The Poisson–Boltzmann equation for biomolecular electrostatics: A tool for structural biology. J. Mol. Recognit. 2002, 15, 377–392. (55) Antosiewicz, J. M.; Shugar, D. Poisson–Boltzmann continuum-solvation models: Applications to pH-dependent properties of biomolecules. Mol. Biosyst. 2011, 7, 2923– 2949. (56) Fogolari, F.; Zuccato, P.; Esposito, G.; Viglino, P. Biomolecular electrostatics with the linearized Poisson–Boltzmann equation. Biophys. J. 1999, 76, 1–16. (57) Sharp, K. A. Incorporating solvent and ion screening into molecular dynamics using the finite-difference Poisson–Boltzmann method. J. Comput. Chem. 1991, 12, 454– 468. (58) Honig, B.; Sharp, K.; Yang, A.-S. S. Macroscopic models of aqueous solutions: Biological and chemical applications. J. Phys. Chem. 1993, 97, 1101–1109. (59) Schaefer, M.; Sommer, M.; Karplus, M.; Biophysique, L. D. C.; Bel, I.; Uni, V.;

34


Page 34 of 51

Page 35 of 51 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


Pasteur, L.; Pascal, B. pH-dependence of protein stability: Absolute electrostatic free energy differences. J. Phys. Chem. B 1997, 5647, 1663–1683. (60) Antosiewicz, J. M.; McCammon, J. A.; Gilson, M. K. Prediction of pH-dependent properties of proteins. J. Mol. Biol. 1994, 238, 415–436. (61) Srinivasan, J.; Cheatham, T. E., III; Cieplak, P.; Kollman, P. A.; Case, D. A. Continuum solvent studies of the stability of DNA, RNA, and phosphoramidate–DNA helices. J. Am. Chem. Soc. 1998, 120, 9401–9409. (62) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W. et al. Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33, 889–897. (63) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990, 112, 6127–6129. (64) Onufriev, A.; Case, D. A.; Bashford, D. Effective Born radii in the generalized Born approximation: The importance of being perfect. J. Comput. Chem. 2002, 23, 1297– 1304. (65) Gilson, M. K.; Zhou, H.-X. Calculation of protein–ligand binding affinities. Annu. Rev. Biophys. Biomol. Struct. 2007, 36, 21–42. (66) Bertonati, C.; Honig, B. H.; Alexov, E. Poisson–Boltzmann calculations of nonspecific salt effects on protein–protein binding free energies. Biophys. J. 2007, 92, 1891–1899. (67) Hou, T.; Wang, J.; Li, Y.; Wang, W. Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations. J. Chem. Inf. Model. 2010, 51, 69–82. 35



(68) Genheden, S.; Ryde, U. The MM/PBSA and MM/GBSA methods to estimate ligandbinding affinities. Expert Opin. Drug Discov. 2015, 10, 449–461. (69) Lu, B. Z.; Zhou, Y. C.; Holst, M. J.; McCammon, J. A. Recent progress in numerical methods for the Poisson–Boltzmann equation in biophysical applications. Commun. Comput. Phys. 2008, 3, 973–1009. (70) Li, L.; Li, C.; Zhang, Z.; Alexov, E. On the dielectric “constant” of proteins: Smooth dielectric function for macromolecular modeling and its implementation in DelPhi. J. Chem. Theory Comput. 2013, 9, 2126–2136. (71) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K. A.; Honig, B. H. Focusing of electric fields in the active site of Cu–Zn superoxide dismutase: Effects of ionic strength and amino-acid modification. Proteins Struct. Funct. Genet. 1986, 1, 47–59. (72) Nicholls, A.; Honig, B. H. A rapid finite-difference algorithm, utilizing successive overrelaxation to solve the Poisson–Boltzmann equation. J. Comput. Chem. 1991, 12, 435–445. (73) Rocchia, W.; Alexov, E.; Honig, B. H. Extending the applicability of the nonlinear Poisson–Boltzmann equation: Multiple dielectric constants and multivalent ions. J. Phys. Chem. B 2001, 105, 6507–6514. (74) Rocchia, W.; Sridharan, S.; Nicholls, A.; Alexov, E.; Chiabrera, A.; Honig, B. H. Rapid grid-based construction of the molecular surface and the use of induced surface charge to calculate reaction field energies: Applications to the molecular systems and geometric objects. J. Comput. Chem. 2002, 23, 128–137. (75) Davis, M. E.; Madura, J. D.; Luty, B. A.; McCammon, J. A.; Lutty, B. A.; Mac Cammon, A. J. Electrostatics and diffusion of molecules in solution: Simulations with the University of Houston Brownian Dynamics program. Comput. Phys. Commun. 1991, 62, 187–197. 36


Page 36 of 51

Page 37 of 51 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


(76) Luty, B. A.; Davis, M. E.; McCammon, J. A. Solving the finite-difference non-linear Poisson–Boltzmann equation. J. Comput. Chem. 1992, 13, 1114–1118. (77) Madura, J. D.; Davis, M. E.; Gilson, M. K.; Wade, R. C.; Luty, B. a.; McCammon, J. A.; Wades, R. C.; Luty, B. a.; McCammon, J. A.; Wade, R. C. et al. Biological applications of electrostatic calculations and Brownian dynamics simulations. Rev. Comput. Chem. 1994, 5, 229–267. (78) Holst, M. J.; Baker, N. A.; Wang, F. Adaptive multilevel finite element solution of the Poisson–Boltzmann equation I. Algorithms and examples. J. Comput. Chem. 2000, 21, 1319–1342. (79) Baker, N. A.; Holst, M. J.; Wang, F. Adaptive multilevel finite element solution of the Poisson–Boltzmann equation II. Refinement at solvent-accessible surfaces in biomolecular systems. J. Comput. Chem. 2000, 21, 1343–1352. (80) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of nanosystems: Application to microtubules and the ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037–10041. (81) Wang, J.; Luo, R. Assessment of linear finite-difference Poisson–Boltzmann solvers. J. Comput. Chem. 2010, 31, 1689–1698. (82) Cai, Q.; Hsieh, M. J.; Wang, J.; Luo, R. Performance of nonlinear finite-difference Poisson–Boltzmann solvers. J. Chem. Theory Comput. 2010, 6, 203–211. (83) Jo, S.; Vargyas, M.; Vasko-Szedlar, J.; Roux, B.; Im, W. PBEQ-Solver for online visualization of electrostatic potential of biomolecules. Nucleic Acids Res. 2008, 36, W270–W275. (84) Prabhu, N. V.; Zhu, P.; Sharp, K. A. Implementation and testing of stable, fast

37



implicit solvation in molecular dynamics using the smooth-permittivity finite difference Poisson–Boltzmann method. J. Comput. Chem. 2004, 25, 2049–2064. (85) Prabhu, N. V.; Panda, M.; Yang, Q.; Sharp, K. A. Explicit ion, implicit water solvation for molecular dynamics of nucleic acids and highly charged molecules. J. Comput. Chem. 2008, 29, 1113–1130. (86) Bashford, D. An object-oriented programming suite for electrostatic effects in biological molecules An experience report on the MEAD project. In Sci. Comput. ObjectOriented Parallel Environ. First Int. Conf. ISCOPE 97 Mar. del Rey, California, USA December 8–11, 1997 Proc.; Ishikawa, Y., Oldehoeft, R. R., Reynders, J. V. W., Tholburn, M., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 1997; pp 233–240. (87) Chen, D.; Chen, Z.; Chen, C.; Geng, W.; Wei, G.-W. MIBPB: A software package for electrostatic analysis. J. Comput. Chem. 2011, 32, 756–770. (88) Tannor, D. J.; Marten, B.; Murphy, R. B.; Friesner, R. A.; Sitkoff, D.; Nicholls, A.; Ringnalda, M. N.; Goddard, W. A., III; Honig, B. H. Accurate first principles calculation of molecular charge distributions and solvation energies from ab initio quantum mechanics and continuum dielectric theory. J. Am. Chem. Soc. 1994, 116, 11875– 11882. (89) Marten, B.; Kim, K.; Cortis, C. M.; Friesner, R. A.; Murphy, R. B.; Ringnalda, M. N.; Sitkoff, D.; Honig, B. H. New model for calculation of solvation free energies: Correction of self-consistent reaction field continuum dielectric theory for short-range hydrogen-bonding effects. J. Phys. Chem. 1996, 100, 11775–11788. (90) Wang, M.; Wong, C. F. Calculation of solvation free energy from quantum mechanical charge density and continuum dielectric theory. J. Phys. Chem. A 2006, 110, 4873– 4879.

38


Page 38 of 51

Page 39 of 51 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


(91) Roos, K.; Viklund, J.; Meuller, J.; Kaspersson, K.; Svensson, M. Potency prediction of β-secretase (BACE-1) inhibitors using density functional methods. J. Chem. Inf. Model. 2014, 54, 818–825. (92) Cortis, C. M.; Langlois, J.-M.; Beachy, M. D.; Friesner, R. A. Quantum mechanical geometry optimization in solution using a finite element continuum electrostatics method. J. Chem. Phys. 1996, 105, 5472. (93) Gogonea, V.; Merz, K. M., Jr. Fully quantum mechanical description of proteins in solution. Combining linear scaling quantum mechanical methodologies with the Poisson– Boltzmann equation. J. Phys. Chem. A 1999, 103, 5171–5188. (94) Westerhoff, L. M.; Merz, K. M., Jr. Quantum mechanical description of the interactions between DNA and water. J. Mol. Graph. Model. 2006, 24, 440–455. (95) Kaukonen, M.; Soderhjelm, P.; Heimdal, J.; Ryde, U. QM/MM–PBSA method to estimate free energies for reactions in proteins. J. Phys. Chem. B 2008, 112, 12537– 12548. (96) Hayik, S. A.; Liao, N.; Merz, K. M., Jr. A combined QM/MM Poisson–Boltzmann approach. J. Chem. Theory Comput. 2008, 4, 1200–1207. (97) Szabo, A.; Ostlund, N. S. Modern quantum chemistry: Introduction to advanced electronic structure theory. In Introd. to Adv. Electron. Struct. Theory; Dover Publications: Mineola, N. Y., 1996; p 480. (98) Nakano, T.; Mochizuki, Y.; Fukuzawa, K.; Amari, S.; Tanaka, S. Developments and applications of ABINIT-MP software based on the fragment molecular orbital method. In Mod. Methods Theor. Phys. Chem. Biopolym.; Starikov, E. B., Lewis, J. P., Tanaka, S., Eds.; Elsevier: Amsterdam, 2006; Chapter 2, pp 39–52.

39



(99) Mochizuki, Y.; Yamashita, K.; Murase, T.; Nakano, T.; Fukuzawa, K.; Takematsu, K.; Watanabe, H.; Tanaka, S. Large scale FMO-MP2 calculations on a massively parallelvector computer. Chem. Phys. Lett. 2008, 457, 396–403. (100) Watanabe, H.; Okiyama, Y.; Nakano, T.; Tanaka, S. Incorporation of solvation effects into the fragment molecular orbital calculations with the Poisson–Boltzmann equation. Chem. Phys. Lett. 2010, 500, 116–119. (101) Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. Phys. Rev. 1934, 46, 618–622. (102) Watanabe, C.; Fukuzawa, K.; Okiyama, Y.; Tsukamoto, T.; Kato, A.; Tanaka, S.; Mochizuki, Y.; Nakano, T. Three- and four-body corrected fragment molecular orbital calculations with a novel subdividing fragmentation method applicable to structurebased drug design. J. Mol. Graph. Model. 2013, 41, 31–42. (103) Mazanetz, M. P.; Ichihara, O.; Law, R. J.; Whittaker, M. Prediction of cyclindependent kinase 2 inhibitor potency using the fragment molecular orbital method. J. Cheminform. 2011, 3, 1–15. (104) Shigemitsu, Y. Quantum chemical study on molecular-level affinity of DJ-1 binding compounds. Int. J. Quantum Chem. 2013, 113, 574–579. (105) Watanabe, C.; Watanabe, H.; Fukuzawa, K.; Parker, L. J.; Okiyama, Y.; Yuki, H.; Yokoyama, S.; Nakano, H.; Tanaka, S.; Honma, T. Theoretical analysis of activity cliffs among benzofuranone-class Pim1 inhibitors using the fragment molecular orbital method with molecular mechanics Poisson–Boltzmann surface area (FMO+MMPBSA) approach. J. Chem. Inf. Model. 2017, 57, 2996–3010. (106) Wing, R.; Drew, H.; Takano, T.; Broka, C.; Tanaka, S.; Itakura, K.; Dickerson, R. E. Crystal structure analysis of a complete turn of B-DNA. Nature 1980, 287, 755–758.

40


Page 40 of 51

Page 41 of 51 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


(107) Fukuzawa, K.; Watanabe, C.; Kurisaki, I.; Taguchi, N.; Mochizuki, Y.; Nakano, T.; Tanaka, S.; Komeiji, Y. Accuracy of the fragment molecular orbital (FMO) calculations for DNA: Total energy, molecular orbital, and inter-fragment interaction energy. Comput. Theor. Chem. 2014, 1034, 7–16. (108) Richards, F. M. Areas, volumes, packing, and protein structure. Annu. Rev. Biophys. Bioeng. 1977, 6, 151–176. (109) Connolly, M. L. Analytical molecular surface calculation. J. Appl. Crystallogr. 1983, 16, 548–558. (110) Lee, B.; Richards, F. M. The interpretation of protein structures: Estimation of static accessibility. J. Mol. Biol. 1971, 55, 379–400. (111) Shrake, A.; Rupley, J. A. Environment and exposure to solvent of protein atoms. Lysozyme and insulin. J. Mol. Biol. 1973, 79, 351–371. (112) Singh, U. C.; Kollman, P. A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 1984, 5, 129–145. (113) Besler, B. H.; Merz, K. M., Jr.; Kollman, P. A. Atomic charges derived from semiempirical methods. J. Comput. Chem. 1990, 11, 431–439. (114) Breneman, C. M.; Wiberg, K. B. Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis. J. Comput. Chem. 1990, 11, 361–373. (115) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model. J. Phys. Chem. 1993, 97, 10269–10280. (116) Fedorov, D. G.; Slipchenko, L. V.; Kitaura, K. Systematic study of the embedding

41



potential description in the fragment molecular orbital method. J. Phys. Chem. A 2010, 114, 8742–8753. (117) Okiyama, Y.; Watanabe, H.; Fukuzawa, K.; Nakano, T.; Mochizuki, Y.; Ishikawa, T.; Tanaka, S.; Ebina, K. Application of the fragment molecular orbital method for determination of atomic charges on polypeptides. Chem. Phys. Lett. 2007, 449, 329–335. (118) Okiyama, Y.; Watanabe, H.; Fukuzawa, K.; Nakano, T.; Mochizuki, Y.; Ishikawa, T.; Ebina, K.; Tanaka, S. Application of the fragment molecular orbital method for determination of atomic charges on polypeptides. II. Towards an improvement of force fields used for classical molecular dynamics simulations. Chem. Phys. Lett. 2009, 467, 417–423. (119) Okiyama, Y.; Nakano, T.; Yamashita, K.; Mochizuki, Y.; Taguchi, N.; Tanaka, S. Acceleration of fragment molecular orbital calculations with Cholesky decomposition approach. Chem. Phys. Lett. 2010, 490, 84–89. (120) Mochizuki, Y.; Nakano, T.; Okiyama, Y.; Sakakura, K.; Yamamoto, J.; Komeiji, Y.; Yamashita, K.; Murase, T.; Ishikawa, T.; Tsukamoto, T. et al. ABINIT-MP 7.0 & BioStation Viewer 16.0 ; Center for Research on Innovative Simulation Software (CISS), Insititute of Industrial Science, the University of Tokyo: Tokyo, 2013. (121) Reed, A. E.; Weinhold, F. Natural bond orbital analysis of near-Hartree–Fock water dimer. J. Chem. Phys. 1983, 78, 4066–4073. (122) Reed, A. E.; Weinstock, R. B.; Weinhold, F. Natural population analysis. J. Chem. Phys. 1985, 83, 735. (123) Travers, A.; Muskhelishvili, G. DNA structure and function. FEBS J. 2015, 282, 2279–2295.

42


Page 42 of 51

Page 43 of 51 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


(124) Shui, X.; McFail-Isom, L.; Hu, G. G.; Williams, L. D. The B-DNA dodecamer at high resolution reveals a spine of water on sodium. Biochemistry 1998, 37, 8341–8355. (125) Wang, J.; Cieplak, P.; Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 2000, 21, 1049–1074. (126) Watanabe, T.; Inadomi, Y.; Fukuzawa, K.; Nakano, T.; Tanaka, S.; Nilsson, L.; Nagashima, U. DNA and estrogen receptor interaction revealed by fragment molecular orbital calculations. J. Phys. Chem. B 2007, 111, 9621–9627. (127) Martinez, C. R.; Iverson, B. L. Rethinking the term “pi-stacking”. Chem. Sci. 2012, 3, 2191–2201. (128) Fedorov, D. G.; Kitaura, K. Subsystem analysis for the fragment molecular orbital method and its application to protein–ligand binding in solution. J. Phys. Chem. A 2016, 120, 2218–2231. (129) Karelson, M. M.; Zerner, M. C. Theoretical treatment of solvent effects on electronic spectroscopy. J. Phys. Chem. 1992, 96, 6949–6957. (130) Szafran, M.; Karelson, M. M.; Katritzky, A. R.; Koput, J.; Zerner, M. C. Reconsideration of solvent effects calculated by semiempirical quantum chemical methods. J. Comput. Chem. 1993, 14, 371–377.

43



Table 1: Definition of the Atomic Radius for each Elementa element radius/˚ A H 1.15 C 1.90 N 1.60 O 1.60 P 2.00 S 1.90 a

See ref 88.

Figure 1: Definitions of molecular surface, the solvent accessible surface (SAS) and solvent exclusive surface (SES). These are drawn by a spherical probe, which approximates the effective size of a solvent molecule, rolling over a solute molecule with atomic radius volume. The former is defined as a path of the center of the probe. The latter is, on the other hand, defined as a boundary between the probe and atomic spheres. The boundary consists of convex- and concave-shaped regions; the former is the contact region exposed to solvent, where the probe contacts only one atom; the latter is the reentrant region curved inward, where the probe contacts two atoms or more simultaneously.

44


Page 44 of 51

Page 45 of 51 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


Figure 2: Fragment-based solute–solvent interactions.

Figure 3: (F)MO-PB iterations with the MP2 correlation energy calculations.

45



Figure 4: Structure and sequence of the 12-base-paired DNA duplex.

Figure 5: Definition of fragmentation for DNA. Fragments are dividing at the positions of each bond detached atom (BDA) connected to bond attached atom (BAA), where an occupied single bond is disconnected from the BDA on a fragment and assigned to the BAA on the other.

46


Page 46 of 51

Page 47 of 51 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


Figure 6: Electrostatic potentials in a.u. are mapped on the molecular surface of the 12base-paired DNA duplex generated by its atomic charges: (a) that in vacuo and (b) the difference with that in implicit solvent. The isoelectron density equals 0.001 a.u. in vacuo and defines the surface. Moreover, (c) the induced surface charges on the dielectric boundary defined in the PB equation is also shown in the unit of elementary charge.

47



Figure 7: IFIEs of all the fragment pairs inside 12 base-paired DNA duplex compared in vacuo to those in implicit solvent as a function of interfragment distance. The distance is defined as the one between the nearest atoms in different fragments. These pairs are categorized into (a) Watson–Crick hydrogen-bonding base pairs, (b) intrastrand-stacking base pairs, (c) interstrand-stacking base pairs, (d) other base pairs, (e) base and sugar– phosphate backbone pairs, and (f) other sugar–phosphate backbone pairs, respectively.

48


Page 48 of 51

Page 49 of 51 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


Figure 8: IFIEs of the proximate base-fragment pairs inside the 12-base-paired DNA duplex compared in vacuo to those in implicit solvent. These pairs are categorized into three interaction types: (a) Watson–Crick hydrogen-bonding base pairs, (b) intrastrand-stacking base pairs, and (c) interstrand-stacking base pairs, respectively.

Figure 9: Changes in the frontier molecular orbital energies on each fragment type of the 12-base-paired DNA duplex in vacuo and in implicit solvent. Each value is averaged for same type fragments excluding those in the terminal base pairs.

49



Figure 10: Frontier molecular orbitals on (a) base and (b) sugar–phosphate backbone fragments of dC9 compared in vacuo to those in implicit solvent.

Figure 11: The solvation free energy of the 12-base-paired DNA duplex in the explicit solvent model as a function of the shell thickness by comparison with those in implicit one. The number of counterions containing of explicit solvent is also plotted.

50


Page 50 of 51

Page 51 of 51 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


Graphical TOC Entry

51


Fragment Molecular Orbital Calculations with Implicit Solvent Based

Recommend Documents