Subsystem Analysis for the Fragment Molecular Orbital Method and Its

Subscriber access provided by University of South Dakota

Article

Subsystem Analysis for the Fragment Molecular Orbital Method and Its Application to Protein-Ligand Binding in Solution Dmitri G. Fedorov, and Kazuo Kitaura J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b00163 • Publication Date (Web): 07 Mar 2016 Downloaded from http://pubs.acs.org on March 8, 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.

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

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

Subsystem Analysis for the Fragment Molecular Orbital Method and its Application to ProteinLigand Binding in Solution Dmitri G. Fedorov†*, Kazuo Kitaura‡ Research Center for Computational Design of Advanced Functional Materials (CDFMat), National Institute of Advanced Industrial Science and Technology (AIST), Central 2, Umezono 1-1-1, Tsukuba, 305-8568, Japan; Graduate School of System Informatics, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan †AIST ‡Kobe University AUTHOR EMAIL ADDRES S: [email protected], Phone: +81-29-861-3170

ABSTRACT A subsystem analysis is derived incorporating interfragment interactions into the fragment properties, such as energies or charges. The relative stabilities of three alanine isomers, the α-helix, the β-turn and the extended form are studied and the differences in fragment properties are elucidated. The analysis is further elaborated for studies of binding energies. The binding of the Trp-cage protein (PDB: 1L2Y) to two ligands is studied in detail. Binding energies defined for each fragment can be used as a convenient descriptor for analyzing contributions to binding in solution. 1 ACS Paragon Plus Environment


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 58

1 Introduction During recent years the field of computational chemistry has shown tremendous progress in developing faster and more efficient algorithms.1,2 For large molecular systems,

linear

scaling

3

,

4

7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25

,

5

,

6

and

fragmentation

approaches

can be used to do calculations in a feasible time.

Fragment-based methods have an additional benefit of providing properties of fragments. In particular, interfragment interactions are attractive means of understanding the physical picture in a large and complex system. The importance of polarization and charge transfer in biochemical systems has been stressed,26,27 prompting applications of quantum-mechanical methods.28 The fragment molecular orbital method (FMO) 29,30,31,32,33 is a fragmentation approach, in which a system is subdivided into fragments, and the total properties of large molecular systems are obtained in a many-body expansion,

34 , 35 , 36 , 37

combining the

properties of fragments. In the two-body method, FMO2, fragments and their pairs (dimers) are calculated. The interfragment pair interaction energies are abbreviated as PIE or IFIE and are used to analyze the molecular recognition, in particular, of ligands by proteins.

38 , 39 , 40 , 41 , 42 , 43 , 44 , 45 , 46 , 47

PIEs have been used extensively in FMO-based

quantitative structure-activity relationship (QSAR) studies. 48 , 49 , 50 , 51 , 52 Improving predictive power of FMO is of importance to biochemical applications, for example, in drug discovery. 53,54,55 The assumption that PIEs determine the observed binding energies or activities is flawed because PIEs do not include the destabilization polarization and desolvation energies of fragments. Moreover, PIEs are often obtained from gas phase calculations

2 ACS Paragon Plus Environment

Page 3 of 58

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


whereas the goal is to model binding in solution. While the importance of solvent to binding energies has been discussed,56,43 the need to analyze binding energies directly within the FMO formalism has not been answered, and it is addressed here, using a subsystem analysis. There are many kinds of analyses developed for quantum-mechanical methods. As Mulliken wrote, “'the more accurate the calculations become, the more the concepts tend to vanish into thin air”,57 and it is an elaborate task to condense them back out of it. In the realm of density functional theory (DFT), there is a variety of so called subsystem DFT methods, 58 but the concept of subsystems there corresponds to fragments in FMO, while a subsystem in this work is a set of fragments. Clearly, for molecules made of a dozen atoms, an atom (or a functional group in chemistry) is a useful building block for an analysis, and for small proteins consisting of a few dozens of residues, these residues are an appropriate unit. When hundreds or thousands of residues comprise a protein, it is useful to consider such units as helices or random coils, that can be called subsystems in the context of the present work. Some attempts to study the change in the internal energies of a protein upon solvation59 and to separate the properties of the protein in a protein-ligand complex have been reported earlier,40 and here these ideas are systematized and applied to the study of isomer stabilities and binding in solution. To describe solvent, polarizable continuum model (PCM) 60 is used in the FMO/PCM method, 61,62,63 and the pair interaction energy decomposition analysis (PIEDA) 64,65 is employed to obtain more physical details.

2 Methodology



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 58

2.1 Subsystem analysis The subsystem analysis (SA) for FMO is formulated for the decomposition of the total energy into subsystem values. The decomposition of other properties is described below separately. The equations are written for FMO2, and a three-body (FMO3) formulation can also be derived. The total system in FMO is divided into N fragments denoted as I, and fragments are grouped into M subsystems denoted as i. As an example, in a proteinligand complex, one can define two subsystems, the protein and the ligand. The total energy E in FMO2 is given as the sum of internal energies of fragments EI′ and fragment pair interaction energies ∆E IJint . For subsystems it becomes N

N

M

M

E = ∑ E I′ + ∑ ∆E IJint = ∑ Ei′ + ∑ ∆Eijint I >J

I

i

(1)

i> j

The fragment internal energy E ′ in case of PCM can be further divided into the internal solute energy of fragment I ( E I′′ ) and its solvation energy ( E Isolv ) as E I′ = E I′′ + E Isolv

(2)

For subsystems, one defines the internal Ei′ and interaction ∆Eijint energies as E i′ = ∑ E I′ + I ∈i

∆Eijint =

∑ ∆E ( )

int IJ

(3)

J > I ∈i

∑ ∆E

int IJ

(4)

I ∈i , J ∈ j

In these expressions and below, for simplicity the upper sum limit N is omitted. Subsystem energies Ei′ can be decomposed into partial energies of fragments.

Ei′ = ∑ E Ipart ,i

(5)

I ∈i


Page 5 of 58

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 partial fragment energies are ′ E Ipart ,i = E I +

1 ∑ ∆E IJint 2 ( J ≠ I )∈i

(6)

Here, 1/2 is used to avoid double counting ( ∆E IJint is added with a factor of 1/2 to both part int int int E Ipart ,i and E J ,i ). ∆E IJ are defined for I ≠ J , and they are symmetric ∆E IJ = ∆E JI .

Some comment should be made regarding the distinction between the internal and partial energies. The notion of “internal” is used to refer to the energy of fragment I (subsystem i) excluding its interactions with other fragments (subsystems); while the notion of “partial” includes these interactions partitioned into fragment (subsystem) contributions. The whole (total energy) is decomposed either as the sum of parts (partial energies) or as the sum of internal energies plus interactions. The internal energies in FMO are a direct measure of the destabilization component of the polarization.64 The partial energies, however, also include the stabilization component of the polarization (via the electrostatic component of the interaction energy). 64 Analogously to partial fragment energies, partial subsystem energies can be defined:

Eipart = Ei′ +

1 M ∑ ∆Eijint 2 j ≠i

(7)

M

E = ∑ Eipart

(8)

i

For molecular clusters the above analysis is well suited as formulated. However, there is a complication in the application of it to fragments connected by covalent bonds (for example, in proteins). In this case the interaction energy can be large, because of the treatment of fragment boundaries (for fragments connected by C-C bonds, the interaction energy is about -15 hartree). 66 If one is interested in energy differences (as below in 5 ACS Paragon Plus Environment


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 58

protein-ligand binding) then in the relative sense the above components can be used even for covalently connected fragments, but if one considers absolute values, then the analysis above can be applied, but in order to make it more useful, one more step is necessary. The previously suggested recipe64 is applied to the subsystem analysis. There are two types of fragment pairs, those where two fragments are connected by a covalent bond (C) and those unconnected (U). Then all pair interactions can be divided into these two types, and the former type can be added to the internal energies forming the backbone energies. For systems with covalently connected fragments it is more convenient to use the backbone energies rather than the internal energies. ′ E Ipart ,i = E I +

1 1 1 1 ∆E IJint = E I′ + ∑ ∆E IJC,int + ∑ ∆E IJU,int = E IBB + ∑ ∆E IJU,int ∑ 2 ( J ≠ I )∈i 2 ( J ≠ I )∈i 2 ( J ≠ I )∈i 2 ( J ≠ I )∈i

(9) where the backbone energy of fragment I in subsystem i is, E IBB = E I′ +

1 ∆E IJC,int ∑ 2 ( J ≠ I )∈i

(10)

The total unconnected interaction energy of fragment I is, ∆E IU,int =

1 ∆E IJU,int ∑ 2 ( J ≠ I )∈i

(11)

Likewise, for subsystems the backbone and partial energies are E iBB = Ei′ +

1 ∑ ∆EijC,int 2 j ≠i

(12)

E ipart = Ei′ +

1 1 ∆Eijint = EiBB + ∑ ∆E ijU,int ∑ 2 j ≠i 2 j ≠i

(13)

where


Page 7 of 58

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


E i′ = ∑ E IBB + I ∈i

1 ∆E IJU,int = ∑ E IBB + ∑ ∆E IJU,int ∑ ∑ 2 I∈i ( J ≠ I )∈i I ∈i ( I > J )∈i

(14)

The following equation is then obtained from eq 1, M

M

E = ∑ EiBB + ∑ ∆EijU,int i

(15)

i> j

The internal energy of a subsystem Ei′ or the total energy E is given by the sum of the backbone energies and unconnected interactions. These expressions (eqs 14-15) are useful in practice. For molecular clusters (or, in general, if subsystems are not connected by covalent bonds), there are no connected interactions and the related equations (e.g., eqs 1 and 15) become identical. part E Ipart can also be called the site energies of a fragment and a subsystem, ,i and E i

respectively. What is the purpose of introducing these quantities? For example, consider performing MD simulations of a protein solvated in explicit solvent with counterions. E i′ can be used to define the energy of the protein, solvent and counterions. One can plot E i′ as a function of time and discuss the fluctuations of the site energies. Another possible use is a development of force field models, when subsystems can be used to generate potentials.67,68 One can also use partial energies to identify the sources of the structure instability. For example, one can look at the site energies of individual residues and discuss their role to the overall stability of structures. 2.2 Other properties As an example of some property other than the energy, the partial fragment charge is defined,



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 58

N

Q Ipart = Q I + ∑ ∆Q J → I

(16)

J ≠I

where QI is the formal charge and ∆Q I → J is the charge transfer from fragment I to fragment J. The normalization of charges is as follows, N

N

N

N

N

∑ QIpart = ∑ QI + ∑ ∑ ∆Q J →I = ∑ QI = Q I =1

I =1

I =1 J ≠ I

(17)

I =1

where Q is the total charge of the system, and the relation ∆Q I → J = − ∆Q J → I is used. In addition to the energy components described above, three properties are also discussed below: the dipole moment of each fragment, d I , and the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) I I energies of each fragment, ε HOMO and ε LUMO , respectively. Here, the most simple and

easy to use way of obtaining molecular orbitals in FMO is used, by taking the orbitals from fragment calculations.69 IP can be useful in predicting what fragment is most prone to be affected by radiation damage, for instance. HOMO and LUMO are also directly related to electronic excitations.70 The orbital energies are discussed here, because they can be used in various contexts related to the concepts of hardness, chemical reactivity and conductivity. 71,72 The relations between the ionization potential (IP), electron affinity (EA) and HOMOLUMO gap of each fragment according to the Koopmans’ theorem, adopted here, are: I IP = −ε HOMO I EA = ε LUMO

(18)

I I gap = ε LUMO − ε HOMO = EA + IP

Other useful quantities to look at when discussing solvent effects are induced solvent charges qI, also called apparent solvent charges, summed over atoms in a fragment, and 8 ACS Paragon Plus Environment

Page 9 of 58

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


solvent coverages sI (computed for the cavity used for the electrostatics in PCM) showing how much of the total surface is exposed to solvent, for each fragment (for a single standalone atom, 100%; for an atom buried inside, 0%).64 One should also note that using the idea of the partial energy it is possible to add pair corrections to extensive fragment properties. For example, the solvation energy E Isolv of fragments above is computed for the electronic states of individual fragments. There is also a two-body (pair) correction to the solvation energy65 and using it one can define partial solvation energies (this, however, is not done here, because the two-body corrections to solvation energies are included in the solvent screening,65 forming a part of pair interactions, and adding it to fragment solvation energies will result in double counting and confuse the analysis). Likewise, the dipole moment d I is computed for individual fragments but one can sum up pair corrections and define a partial fragment dipole moment. The general expression for a property A is (for which pair corrections are symmetric, ∆AIJ = ∆AJI ), 36

AIpart ,i = AI +

1 N ∑ ∆AIJ 2 ( J ≠ I )∈i

(19)

For antisymmetric properties ∆AIJ = − ∆AJI one can use (compare to eq 16) N

AIpart ,i = AI +

∑ ∆A

(20)

JI

( J ≠ I )∈i

However, some properties cannot be obtained simply following this expression. An example is orbital energies. Although one can use many-body expansions in FMO to



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 58

obtain the Fock matrix, 73 , 74 , 75 the resultant orbital energies are not decomposable following eq 19, because these properties are not extensive. 2.3 Discussion of definitions Some thought should be given to the splitting of interactions in the partial energies. A general discussion of reducing many-body effects into a smaller order (for example, here a two-body quantity, interaction energy, is reduced to a one-body quantity, partial energy) is found elsewhere, 36 and this idea finds use for three and four-body terms as well.

35

However, in fact two-body properties PIEs are already a contraction of three-

body effects. It is clear by looking at the explicit definition of PIE,65 ∆E IJint = (E IJ′′ − E I′′ − E ′J′ ) + Tr ( ∆D IJ V IJ ) + ∆E IJsolv

(21)

where the second term is Tr (∆D IJ V IJ ) =

∑ Tr (∆D

IJ

V IJ ( K ) )

(22)

K ≠I ,J

and ∆D IJ is the density transfer matrix for dimer IJ and V IJ ( K ) is the matrix of the contribution of fragment K to the electrostatic potential acting upon dimer IJ. This is the reason why FMO2, although on the surface looks like a 2-body method as judged from is energy expression, is in fact a many-body method, with many-body polarization effects64 included in fragment energies EI′′ and three-body effects (two-body charge transfer coupled to one-body electrostatics) in Tr (∆D IJ V IJ ) . One can also consider the conceptually related case of atomic charges. There is no unique way of dividing molecular density into atomic contributions. A very simple recipe by Mulliken 76 is to rely on atomic orbitals, a formal division without regard for their spacial distribution, and to divide contributions from two atoms equally between them.


Page 11 of 58

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


There are many alternative formulations, many of which use spacial distribution of density. Here, molecular interactions are also divided among individual fragments. As argued elsewhere,36 one can do a division with equal factors, or use some other weights. Here, two simple cases are considered as food for thought encouraging further future studies on the subject. The first is a hydrogen bond between two fragments; as an example, it can be written as O-H…O. Clearly, it is asymmetric and there is no specific reason to assign 1/2 of its energy to each fragment in the partial energy definition. Which side contributes more? Some analyses developed for atoms, may be able to give some quantitative answers.77,78 The partial energies are in general size-dependent. The larger the subsystem, the more interactions it would have under the same conditions, and, consequently, the partial energies change as well. In a water cluster, the partial energy of a water fragment in a subsystem can become more negative if one adds more water molecules to the same subsystem. Interestingly, there is an attempt to eliminate the size dependence of interactions in the FMO context, where the interaction is normalized by the number of heavy atoms, called fragment efficiency.53 This is an example of a concept that can be called energy density. 2.4 Binding in solution In PIEDA/PCM, one can obtain PIE for protein-ligand binding. However, these PIEs in the complex by themselves are very misleading. The reason is that they do not reflect the desolvation penalty.56 If one is to ask the question of what the driving factors for the ligand binding in solution are, an answer cannot be given only by PIEs in the complex.



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 58

Physically, one should look at the change in the total energy between the isolated state and the complex, and decompose it into physically meaningful components. Before a protein binds a ligand, the surface of the whole protein is open to interactions with solvent. When a ligand approaches, it removes some solvent from the protein surface (and vice versa for the ligand), resulting in a desolvation penalty. Therefore, strongly attractive interactions in the complex can be accompanied by a large desolvation penalty, and the actual contribution of each residue to the binding can be different. In order to increase the descriptive power of FMO for studying binding or chemical reactions in solution, the following formulation is proposed. Let us consider the following general process of a complex formation between A and B, of which protein-ligand binding in water is one important example and there A and B denote the protein and ligand, respectively.

A + B → AB

(23)

~ The energy of the complex formation ∆E (the binding energy) is given by ~ ~ ~ ∆E = E AB − E A − E B

(24)

where E AB is the energy of the complex AB, and the energies of A and B in their

~ ~ isolated energy-minimum state (denoted by tilde) are E A and E B , respectively. When the complex is formed, usually the geometries of A and B change (a counterexample is an Ar2 complex formation from two argon atoms). The energies of A and B at their geometries in the complex are denoted as E A and E B , respectively. They are in general different from the isolated state, and the difference is called the deformation energy.38

~ ~ ~ ~ ~ ∆E = E AB − E A − E B = E AB + E A − E A − E A + E B − E B − E B ~ ~ = ∆E + ∆E def, A + ∆E def,B


(25)

Page 13 of 58

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


where ∆E = E AB − E A − E B ~ ~ ∆E def, A = E A − E A ~ ~ ∆E def, B = E B − E B

(26)

~ The complex binding energy ∆E is decomposed into the direct binding energy ∆E ~ (relative to A and B in the complex geometry), and the deformation energies ∆E def, A and ~ ∆ E def, B . The direct formation energy is somewhat analogous to the vertical excitation

energy for electronic excited states, whereas the total binding energy corresponds to the adiabatic transition between the respective minima. Next, ∆E is decomposed by applying the subsystem analysis for the case of two subsystems, M=2, and the deformation energies are dealt with later. In FMO/PCM, the following equation is used for the total energy of a system (A, B or AB): E = ∑ E I′′ + ∑ ∆E Isolv + ∑ ∆E IJint I

I

(27)

I >J

These terms can be decomposed into contributions (capital letter superscripts (ES etc) are used for solute-solute interactions, while small letter superscripts (es etc) are used for solute-solvent terms),65 ∆E Isolv = ∆E Icav + ∆E Ies + ∆E Idisp + ∆E Irep

(28)

∆E IJint = ∆E IJES + ∆E IJEX + ∆E IJCT + mix + ∆E IJDI + ∆E IJsolv

(29)

∆E IJsolv = ∆E IJes + ∆E IJCT⋅es + ∆E IJdisp + ∆E IJrep

(30)

In eq 28, the total solvation energy of a fragment representing solute-solvent interactions is decomposed into the electrostatic (es), dispersion (disp) and nonelectrostatic repulsion (rep) contributions (following earlier work65 but in the notation the subscript I(I) was simplified to I here). In eq 29, the total pair interaction (int) between 13 ACS Paragon Plus Environment


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 58

fragments representing solute-solute interactions is decomposed into electrostatic (ES), exchange-repulsion (EX), charge transfer and mixed (higher order) terms (CT+mix), dispersion (DI) and solvent screening (solv). The latter representing solute-solvent interactions is decomposed in eq 30 into the electrostatic (es, labeled as es2 earlier65), coupling of the intersolute charge transfer to the solvent electrostatics (CT·es, labeled as es3 earlier65), dispersion (disp) and non-electrostatic repulsion (rep) contributions. Substituting the total energies (eq 27) into the binding energy (eq 26), one gets: ∆E =

∑ (E ′′ + ∆E ) + ∑ ∆E solv I

I

I ∈AB

int IJ

I > J ∈ AB

(

) ∑ ∆∆E

= ∑ ∆E I′′ + ∆∆E Isolv + I ∈A

= ∆E ′A + ∆E B′ + ∆E

(

) ∑ ∆E − ∑ (E ′′ + ∆E ) − ∑ ∆E + ∑ (∆E ′′ + ∆∆E ) + ∑ ∆∆E + ∑ ∆E

− ∑ E I′′ + ∆E Isolv − I∈A

int IJ

I > J ∈A

I

int IJ

I > J ∈A

I ∈B

solv I

I ∈B

solv I

I

I > J ∈B

int IJ

I > J ∈B

int IJ

I ∈A, J ∈B

int AB

(31) where for A (and likewise for B):

(

) ∑ ∆∆E

∆E ′A = ∑ ∆E I′′ + ∆∆E Isolv + I ∈A

int ∆E AB =

∑ ∆E

int IJ

(32)

I > J ∈A

int IJ

(33)

I ∈ A , J ∈B

The binding energy ∆E in eq 31 is decomposed into the total destabilization polarization energies64 in solution ∆E ′A and ∆E B′ (representing the polarization of A by int B and vice versa, in solution) and the interaction energy between A and B, ∆E AB (which

includes the stabilization polarization). The definition of quantities in eqs 32-33 is (the superscripts A, B and AB denote the system for which the property is calculated) ∆E I′′ = E I′′ AB − E I′′ A if I ∈ A or ∆E ′′ = E I′′ AB − E I′′ B if I ∈ B ∆∆E Isolv = ∆E Isolv , AB − ∆E Isolv , A if I ∈ A or ∆∆E Isolv = ∆E Isolv , AB − ∆E Isolv , B if I ∈ B ∆∆E

int IJ

= ∆E

int, AB IJ

− ∆E

int, A IJ

if I , J ∈ A or ∆∆E

int IJ

= ∆E

int , AB IJ

− ∆E


int , B IJ

if I , J ∈ B

(34)

int IJ

=

Page 15 of 58

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


int The PIE contribution ∆E AB is accompanied by the changes in the internal energy of A

and B, ∆E ′A + ∆E B′ containing the desolvation penalty. By dividing the internal energy of A (and, likewise, of B) into partial fragment contributions, ∆E ′A = ∑ ∆E Ipart ,A

(35)

I ∈A

solv ′′ ∆E Ipart + , A = ∆E I + ∆∆E I

1 ∆∆E IJint ∑ 2 ( J ≠ I )∈A

(36)

one can rewrite the binding energy as

(

)

int part ∆E = ∆E ′A + ∆E AB + ∆E B′ = ∑ ∆E Ibind , A + ∑ ∆E I , B I ∈A

∆E

bind I,A

= ∆E

part I,A

+ ∑ ∆E

I ∈B

int IJ

(37)

J ∈B

int Some bias (asymmetry) is introduced by merging the ∆E AB term with the partial

energy of A, to focus on the binding contribution of fragments in A (protein) to the fragments in B (ligand). To discuss the binding contribution of fragment I in solution, one has to consider not just

∑ ∆E

int IJ

(residue-ligand PIEs) but rather the residue binding

J ∈B

energy ∆E Ibind , A , which includes the desolvation energy. This is a very important message of this work. If B (e.g., ligand) consists of just one fragment L, and one wants to know what is the role of a residue I ∈ A in the total binding, then the equations become part ∆E = ∑ ∆E Ibind , A + ∆E L , B

(38)

I∈A

part int solv ′′ ∆E Ibind + , A = ∆E I , A + ∆E IL = ∆E I + ∆∆E I

1 ∑ ∆∆E IJint + ∆E ILint 2 ( J ≠ I )∈A


(39)


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 58

Now the physical meaning of terms can be discussed. The direct protein-ligand binding energy ∆E in eq 38 is given as the sum of the residue binding energies ∆E Ibind , A plus the ligand partial energy (showing how much the ligand energy within the complex with a protein differs from the energy of L by itself; this includes the desolvation energy of L). The binding energies of individual residues ∆E Ibind , A are decomposed in eq 39 into the residue-ligand interaction energies ∆E ILint and the changes due to complexation in (a) the internal energies ∆EI′′ of the amino acid residues, (b) the desolvation energies of residues ∆∆E Isolv , and (c) the residue-residue pair interactions ∆∆E IJint . To summarize, it is suggested that instead of using ∆E ILint , one should look at ∆E Ibind , A as a descriptor of protein-ligand binding. It should be noted that here the desolvation energies are calculated from the solvation energies of monomers, i.e., they do not reflect pair contributions. Finally, the deformation energies are addressed. They are always endothermic (positive) because the energy of the isolated state is a minimum energy structure (if the minimum energy structure is not the global minimum, and the structure in the complex is lower in energy, then the deformation energy can turn out to be negative but that would imply that one should find the global minimum for the isolated state). One can decompose the deformation energy and find out which part of the system contributes more to the deformation. For A (and likewise for B): ~ ~ ∆E def, A = E A − E A = ∑ E I′′ + ∆E Isolv +

(

) ∑ ∆E

I∈A

~ ~ = ∑ ∆E I′′ + ∆∆E Isolv +

(

I ∈A

) ∑ ∆∆E~ I > J∈A

int IJ

I > J ∈A

int IJ

~ ~ − ∑ E I′′ + ∆E Isolv −

(

I ∈A

~ part = ∑ ∆E Idef, ,A I ∈A

where 16 ACS Paragon Plus Environment

) ∑ ∆E~ I > J ∈A

int IJ

=

(40)

Page 17 of 58

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


~ ~ ∆E I′′ = E I′′ − E I′′ ~ ~ ∆∆E Isolv = ∆E Isolv − ∆E Isolv ~ ~ ∆∆E IJint = ∆E IJint − ∆E IJint

(41)

1 ~ part ~ ~ ~ ∆E Idef, = ∆E I′′ + ∆∆E Isolv + ∆∆E IJint ∑ ,A 2 ( J ≠ I )∈A

(42)

~ ~ ~ The tilde in E I′′ , ∆E Isolv and ∆E IJint indicates that these quantities are computed for the geometry of the isolated state. In other words, the deformation energy is decomposed into

~ the change due to the structure deformation in the fragment internal energies ∆EI′′ , ~ ~ fragment solvation energies ∆∆EIsolv , and interaction energies ∆∆E IJint . It is possible to do the analysis without separating the deformation energy. Instead, one can take the difference between the energy of AB and isolated states directly. The equations are the same as above (eqs 31-33), only the meaning of the reference states of A and B changes to assume that their geometries are for the isolated energy-minimum state rather than the same ones as in the complex. There are many various quantities defined, which allow a magnifying glass like investigation of details upon need. In order not to get lost in the relations between them, the tree of components is shown in Figure 1. In order to get more details one can go one level down and trace where a certain change comes from.

3 Computational details The methods were implemented79 into GAMESS, 80 and this program was used for all calculations. C-PCM with FIXPVA and 60 tesserae per sphere was employed with the van-der Waals radii and the cavitation and dispersion models (ICAV=1 IDISP=1).



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

Spherical AOs (5d) were used. The ES-DIM threshold in FMO was set to 2.0, and the analytic FMO/PCM gradient was used.81 In the FMO/PCM method, the induced solvent charges on the cavity are determined using the potential calculated from the electronic state of monomers (individual fragments) neglecting pair contributions. Polypeptides and proteins were fragmented as 1 residue per fragment at Cα atoms using Facio,82 the ligands formed separate fragments. The hybrid orbital projection treatment of fragment boundaries was used.29 Three isomers of capped polyalanines containing 10 residues each, were chosen: αhelix, β-turn and the extended form. The initial geometry was taken from the previous study and the coordinates of all atoms were optimized83 with RHF/D/PCM and the 6-31G* basis set. The optimized structures are shown in Figure 2. The protein-ligand binding in solution is examined for the Trp-cage protein (PDB: 1L2Y), 84 which, due to its small size, has been the object of many computational studies.85,86,87,88,89 Two ligands are studied, the neutral (HOPhCOOH) and deprotonated (HOPhCOO–) forms of the p-phenolic acid. The docking poses were chosen following the earlier study. 90 The geometry was constructed as follows. For the protein, the initial structure was taken from the previous study63 and it was reoptimized at the FMORHF/D/PCM level (D is the empiric dispersion 91 ). Isolated ligands were also optimized at the same level (for N=1, FMO is not used, and the method is RHF/D/PCM). For all calculations, 6-31(+)G* basis set was used, i.e., 6-31+G* for carboxylic groups COO– and 6-31G* for all other groups. All geometries were optimized with the threshold (OPTTOL) of 10-4 hartree/bohr.


Page 18 of 58

Page 19 of 58

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


For the protein ligand complexes a preliminary single point calculation was performed to identify the important residue fragments in the binding pocket based on the interaction energies and the distance between fragments. For the neutral ligand the important residues were found to be Tyr-3 and Pro-18, for the anionic form, Gln-5 and Arg-16. The geometry of each protein-ligand complex was optimized for these two residues and the ligand, all other atoms were kept frozen.

4 Results and Discussion 4.1 Comparative analysis of polyalanine isomers. The subsystem analysis applied to the polyalanine isomers (Figure 2), presents an energy decomposition different from similar attempts done in the past.62, 92 The two terminal fragments include caps, and these fragments have a different number of atoms than others; thus, they are excluded from some comparisons. The fragments are numbered from the C-terminus, capped with NHCH3 to the N-terminus, capped with COCH3, so that fragments 2-9 correspond to identical alanine residue fragments, whereas fragments 1 and 10 are the C and N-terminus, respectively. Only one subsystem is defined in polyalanine. The ionization potential of fragments (Figure 3) in the α-helix shows the trend to decrease from the C to N-terminus (thus the tip is the easiest to ionize or to donate electron density). On the other hand, the β-turn and extended forms resemble each other except at the turn (i.e., fragments 5 and 6). The N-terminus has similar ionization potentials in all three forms. The electron affinity in the α-helix shows a uniform increase from the C to N-terminus. Overall, the β-turn and extended forms resemble each other



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

except at the turn. The HOMO-LUMO gap has an interesting feature that for the middle fragments 2-9 it is almost constant for the α-helix and the extended form, while for the βturn it features 3 peaks (the highest is at the turn). The properties shown in Figure 4 are absolute fragment energies, thus they are large. The internal energies EI′′ reflect the destabilization polarization of the solute, and the backbone energies E IBB add to it the solvation energies and the connected interactions (covalent boundary corrections). E Ipart also include unconnected interactions (in polyalanine there is only one subsystem so that the subsystem index i is omitted in E Ipart ,i ). Because for fragments 2-9 the fragment definitions in all three isomers are the same, these quantities can be compared. For each isomer, the shapes of the curves for all three properties E I′′ , EIBB or E Ipart in general are similar, but the curves shift up and down relative to each other. The destabilization polarization energy, determined by the sum of the internal solute energies E I′′ over all fragments, is the largest for the β-turn, followed by the extended form and the α-helix. To elucidate it further, gas phase calculations were also performed. The order was: the β-turn, α-helix and extended form, and the difference implies a relative reduction of the polarization of the α-helix by the solvent, in other words, in gas phase the extended form is the least polarized, but in solution it is the α-helix. For all three properties in Figure 4, the extended form and the β-turn feature oscillations, with a period of two, and for the extended form is shape is sinusoidal. The partial energy of the extended form is shifted up relative to the backbone energy by the unconnected interactions, which are weak in that isomer. The fragment partial energy answers the


Page 20 of 58

Page 21 of 58

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


question: what fragment is most stable? The answer is that fragment 9 in the α-helix is the most stable alanine fragment among all three isomers; the least stable fragment is number 4 in the β-turn. For the turning point, fragments 5 and 6 show quite a difference in the stability (about 12 kcal/mol). The solvation energy for the middle part (fragments 2-9) is quite uniform for the extended form, as expected due to its linear structure (Figure 5). The β-turn shows two peaks in the energy for fragments 3 and 8. There is also an approximate mirror symmetry in the distribution so that the property for fragment I is similar to that of fragment 11-I, observed for the β-turn and extended form. The termini are special in that they have more solvent accessible surface. The solvent coverage has an approximate mirror symmetry, and is nearly constant for the middle fragments except the β-turn, where the turn has a larger coverage. The induced solvent charge shows a nearly flat profile for the middle fragments in the extended form (actually, very small oscillations occur). In the β-turn a sinusoidal behavior is observed, disrupted into a plateau at the turn (fragments 5 and 6). For the α-helix there is a steady increase in the induced charge from the C to N-terminus. The fragment dipole moments (Figure 6) in FMO1 include many-body polarization effects due to the self-consistent treatment of the electrostatic potential. The α-helix has the smallest dipole moment for each fragment than the other two isomers, except for the N-terminus and fragment 5. The direction of the dipole moment vector (not shown) in the extended form alternates between up and down along the chain (apparently mainly because of the C=O groups). It is interesting that the presence of the electrostatic embedding does not lead to a unidirectional increase in the fragment dipole moments in the α-helix. The β-turn shows some variety at the turn (fragments 5-6), featuring a drop. 21 ACS Paragon Plus Environment


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 58

The dipole moment of fragments are similar in the β-turn and extended form, excluding the turn. The fragment charge (Figure 6) is uniformly close to 0 for the extended form (with very small oscillations), because the main source of charge is hydrogen bonding absent in the extended form. For the β-turn it shows a sinusoidal shape disrupted at the turn, where charges flatten out to nearly zero. The α-helix shows a uniform decrease of the charge from the C to N-terminus. The sign is flipped in the middle, and the C and N terminus shoulders are charged positively and negatively, respectively. The solute charge Q should correlate to the induced solvent charge q, because the two properties are related in PCM (approximately because of the escaped charge93 issues) as: q≈−

ε −1 Q ε

(43)

where ε is the dielectric constant. Fragment charges q I in FMO/PCM are determined by summing the charges on tesserae assigned to atoms in fragment I. These charges are not related to formal charges QI according to eq 43, because the induced charges q I are determined not by the solute potential of fragment I only, but of the total system. The many-body effects (the polarization of solvent charges by all fragments in the calculation of q I and interfragment charge transfer in QIpart ) make the dependence much more complicated.65 The unconnected interaction energies for each fragment decrease in the order: α-helix, β-turn and the extended form. The interactions are mainly due to hydrogen bonding as well as steric repulsions in the β-turn, as further discussed below. The α-helix and extended form have an almost flat profile for fragments 3-8, but fragments 1, 2 and 9 in 22 ACS Paragon Plus Environment

Page 23 of 58

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 α-helix have less attractive values. The turn in the β-turn has a nearly zero value, and the stabilization is mainly to the terminal groups, as well as fragments 3, 4 and 8. The interactions for each pair I,J are shown in Figure 7. For the α-helix the trend is that I,I+2 pairs for the C-terminus shoulder and I,I+3 pairs for the N-terminus part of the helix have an interaction of -7…-10 kcal/mol (because fragmentation in FMO differs from conventional residues by the CO group, I,I+3 for fragments corresponds to the I,I+4 pattern for residues; whereas for 4 fragments near the C-terminus the structure is somewhat distorted relative to the canonical α-helix). Other interactions are weaker and are uniformly attractive. The extended form shows an interesting chess-board pattern, with the far interactions alternating between weakly repulsive and attractive. This is related to the fragment dipole – for pairs in which dipoles moments are aligned (J-I=2) and anti-aligned (J-I=1) repulsion and attraction is observed, respectively, affected by the solvent screening. The pairs I,I+2 have a repulsion of about 2 kcal/mol (in two cases the values slightly exceed the threshold of 2 kcal/mol and are shown in a different color). The β-turn shows more features than the rest: it has strong repulsions mixed with very strong attractions, due to its compact form. The largest attraction, -21.3 kcal/mol, is between the two termini, 10,1. Pairs of fragments directly opposite to each other tend to form hydrogen bonds and show strong attraction due to hydrogen bonding: pairs 9,3 and 8,4, -13.6 and -15.6 kcal/mol, respectively. Some pairs involving a fragment at the turn show a strong repulsion (8,5 and 10,3), 6.0 and 5.1 kcal/mol, respectively. The origin appears to be the two C=O groups (see Figure 2), with the oxygens separated by 3.0 and 3.3 Å, respectively. The ES, EX, CT+mix, DI and solv components in PIEDA for the 8,5 pair are: 7.5, 0.8, -0.2, -0.8 and -1.2 kcal/mol, respectively. This means that the nature of



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 repulsion is electrostatical (namely, of the dipole-dipole kind; the solvent screens it to some extent, by -1.2 kcal/mol). Finally, the total stabilities (the differences of total energies E) are as follows: the most stable is the α-helix, relative to which the β-turn and extended form have the energies of 31.5 and 48.0 kcal/mol, respectively. These values are very strongly affected by the dispersion interaction. Without it, the α-helix is still the most stable isomer, but the energy differences to the β-turn and the extended form are only 10.6 and 6.2 kcal/mol, respectively so that the order of these two isomers is also reverted. It is not surprising at all that the dispersion interaction is strong in the compact α-helix and β-turn, and the relative stabilities are larger than those found in the previous study, using structures optimized without dispersion,62 which is attributed to the role of the dispersion on the structures.

4.2 Analysis of protein-ligand binding in solution The optimized structures of the protein-ligand complexes are shown in Figure 8 and some numeric results are also shown in Figure 1. For the neutral ligand, the total binding energy of -8.0 (here and below the units are kcal/mol) is divided into the direct binding energy of -13.7, and the deformation energies of the protein and the ligand of 0.5 and 5.2, respectively. The protein is not deformed much, whereas the ligand undergoes a rotation of the OH in the COOH group (Figure 8c). The deformation energy of 5.2 is further ~ decomposed into the change in the internal energy ∆E I′′ of 10.5 and the solvation energy ~ ∆∆E Isolv of -5.3 so that the solvent stabilizes the ligand isomerization, mainly due to the ~ electrostatic contribution ∆∆E Ies of -4.9 (see eq 28), because of the more favorable for

solvation redistribution of the charges in the solute. 24 ACS Paragon Plus Environment

Page 24 of 58

Page 25 of 58

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 sum of all protein-ligand pair interactions is -39.6, that can be compared with the total binding energy of -8.0. The difference is the change in the internal energy of the protein and ligand. For the ligand the change of 16.7 is mainly due to the desolvation penalty (16.8) plus the destabilization polarization of the ligand by the protein. The reason for the large desolvation penalty is the loss of the solvent accessible surface. The protein has a change of 9.2 in its internal energy. Using partial energies, this change is found to be mainly due to fragments Leu-2 (2.1), Tyr-3 (-2.4), Gln-5 (2.7), Trp-6 (-1.7), Arg-16 (3.4), Pro-17 (4.3), Pro-18 (-2.5), Pro-19 (-1.5) and Ser-20 (3.1). Some of these changes are exothermic, implying that the fragment is stabilized upon complexation due to more favorable interactions or its internal energy). The total desolvation penalty of the protein i.e.,

∑ ∆∆E

solv I

, is 10.7 (can be compared with the internal protein energy change

I ∈A

of 9.2). The residues with large desolvation contributions are Gln-5 (1.4), Arg-16 (6.0), Pro-18 (1.4) and Ser-20 (0.9). Some residues have an exothermic desolvation penalty implying they gained in the solvation energy upon complexation (the largest is Asn-1, 0.3). Thus, for the neutral ligand the sum of the protein-ligand interactions (-39.6 kcal/mol) is a large overestimate of the total binding energy (-8.0 kcal/mol), and the difference is mainly due to the deformation energy of the ligand (5.2 kcal/mol), and the desolvation penalties of the ligand and protein, 16.8 and 10.7 kcal/mol, respectively. Turning now to the anionic ligand, its total binding energy of -9.5 (here and below the units are kcal/mol) is divided into the direct binding energy of -11.1 and the deformation energies of the protein and ligand of 0.3 and 1.3, respectively. The total protein-ligand interaction energy is -66.9, counterbalanced by the changes in the internal energies of the 25 ACS Paragon Plus Environment


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

protein and ligand of 23.2 and 32.6, respectively. The rather large change in the internal energy of the ligand is mainly due to its desolvation penalty (33.5). The protein has a change of 23.2 in its internal energy. The main contributors to it determined using the partial energies of fragments are: Asn-1 (3.1), Leu-2 (6.0), Gln-5 (7.8), Lys-8 (4.4), Asp9 (-3.3), Arg-16 (8.2) and Ser-20 (-2.0). The total desolvation penalty of the protein is 11.8 (the main contributors are Leu-2 (0.7), Tyr-3 (0.7), Gln-5 (2.2), Arg-16 (6.2) and Pro18 (1.2)). This is only about one half of the internal energy change of 23.2. The rest comes from the polarization of the protein and interactions inside the protein, the sum of which and the desolvation penalty is the internal energy contribution of the protein. The anionic ligand destabilizes the protein more than the neutral ligand. The total sum of the protein-ligand interactions (-66.9 kcal/mol) is a large overestimate of the binding energy -9.5 kcal/mol. The desolvation penalty is a large contributor to the difference, 11.8 and 33.5 kcal/mol, for the protein and ligand, respectively, offsetting the large gain due to the protein-ligand interactions. Another way of looking at it is that the favorable protein-solvent interactions are replaced by the protein-ligand interactions (in agreement with the cluster hydration model56), and it is this difference that forms the binding energy. The partial charges on the protein and ligand for the neutral (anionic) ligand are: 1.016 (1.002) and -0.016 (-1.002), respectively. Thus the total protein-ligand charge transfer is relatively small, and the total charge on the ligand is close to its formal charge of 0 and -1 for the neutral and anionic forms, respectively. However, this is the total amount and there is some cancellation due to the charge being transferred between the protein and ligand in both directions.


Page 26 of 58

Page 27 of 58

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 protein-ligand energies are shown in Figure 9. Three properties are plotted, ∆E ILint , bind int ∆E Ipart , A and ∆E I , A related by eq 39. The PIE ∆E IL has been used formerly in many

applications of FMO discussing protein-ligand binding. The fragment binding energy part ∆E Ibind , A is proposed in this work, and the sum of all such values and ∆E L , B gives the

direct binding energy ∆E (eq 38). PIEs are useful to analyze a complex, whereas fragment binding energies are appropriate to understand a complex formation from components. For the neutral ligand, the two quantities ∆E ILint and ∆E Ibind , A are similar to each other in most cases; the differences are local to the binding site, and for both Tyr-3 and Pro-18 int ∆E Ipart , A enhances PIE, ∆E IL . However, for Arg-16, Pro-17 and Ser-20 they have an

opposite sign and attractive PIEs turn into repulsive binding energies, due to the desolvation penalty. For the anionic ligand, the partial energies of residues and PIEs are of the opposite sign for the main residues, Asn-1, Gln-5, Trp-6, Lys-8 and Arg-16. For Arg-16, the PIE of -31.4 kcal/mol turns into the binding energy of -23.2 kcal/mol, whereas for Gln-5 the change is from -11.1 to -3.3 kcal/mol, making Gln-5 relatively insignificant. Summarizing, focusing on PIEs misses the desolvation penalties, especially significant for charged ligands. These penalties are in some sense correspond to concept of the solvent screening65 reducing interactions in solution, applied to the binding. In addition to the solvent screening, the destabilization of the protein energy (as seen in the partial energies of fragments) is another factor. Likewise, the destabilization and the desolvation



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 28 of 58

of the ligand are also important, and the protein and ligand deformation energies also contribute. If one is to compare several ligands by analyzing the binding energies, the very important contribution of ∆E Lpart , B should not be omitted. A suggestion is made here to include it in the plot at the last contribution (Figure 10). One can consider ∆E Lpart , B as the part binding energy contribution of the ligand, and denote ∆E Lbind , B = ∆E L , B . Then the set

{∆E } = {∆E bind I

bind I,A

}

, ∆E Ibind , B , where I runs over all fragments in the complex (both in the

protein and ligand), is the quantity proposed here to analyze protein-ligand binding energies. N

∆E = ∑ ∆E Ibind

(44)

I =1

In this definition protein-ligand interactions ∆E IJint , I ∈ A, J ∈ L are not split equally between the ligand and protein and, instead, all of it is added to the protein contributions ∆E Ibind , contrary to what is done for the partial energies (eq 5). Such a definition may be more suitable for such systems as protein-ligand complexes. Even if a ligand is split into several fragments, eq 44 can still be applied, and the corresponding plot (Figure 10) will have several ligand energies rather than 1. In drug discovery many ligands are relatively large, and it may be advantageous to divide ligands into several fragments in order to obtain more information. The binding energies in Figure 10 now show quite clearly all contributions to the complex binding energy; for both ligands, the ligand energy change (mainly, its desolvation penalty) is one of the main destabilizing factors; for the charged ligand the


Page 29 of 58

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


repulsion with the likewise charged residue (Asp-9) is also large, as well as the contribution of Leu-2. For the charged ligand, oppositely charged fragments (Asn-1, Lys8 and Arg-16) are the main contributors to the binding. The solvent screening ∆E IJsolv makes a large contribution to ∆E IJint (not shown), reducing it drastically especially for the two termini, which have a large solvent accessible surface and large induced solvent charges.65 If the deformation energies of the protein and ligand are not separated, as discussed at the end of section 2.4, one can obtain binding energies with respect not to the structures of protein and ligand in the complex (as in eq 25), but rather with respect to the minimum ~ geometry of their isolated states. In this case, the sum of these binding energies ∆E Ibind

will give the total binding energy.

~ N ~ ∆E = ∑ ∆E Ibind

(45)

I =1

~ An explicit definition of ∆E Ibind (without a subsystem index the total set is meant,

{∆E~ } = {∆E~ bind I

bind I,A

~ , ∆E Ibind ) is not given here as it is the same as ∆E Ibind except that the ,B

}

properties of the protein and ligand by itself, which are subtracted in eq 34 are computed at the optimized geometries of the isolated protein and ligand rather than for their ~ bind geometry in the complex, as is done in eq 34. ∆E Ibind , A is related to ∆E I , A as ~ ~ def,part bind ∆E Ibind , A = ∆E I , A + ∆E I , A

(46)

and likewise for B (see also eq 42).

5 Conclusions



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 subsystem analysis has been formulated for the fragment molecular orbital method and applied to study energy changes during complex formation. The analysis is presented for the polarizable continuum model of solvation but with some modifications it can be adapted to gas phase studies or other solvent models. The binding energy is represented as a sum of the destabilization polarization and interaction energies, and the former includes the desolvation penalty. The analysis has been applied to study the stabilities of polyalanine isomers and protein-ligand binding. The α-helix shows a uniform increase or decrease in some properties as a function of the position., which is related to the length of helices necessary to accomplish certain biochemical functions. The β-turn has interesting properties related to the turn, as described in detail above. For the protein ligand complexes a novel concept of the fragment binding energies has been proposed with the important relation that their sum is equal to the direct or total binding energy, which is expected to be a better descriptor for QSAR studies. The importance the internal energy, composed of the desolvation penalty as well as the mutual polarization effects of the protein and ligand have been elucidated. In biological systems entropy is often important.94 This analysis may be extended in future to include some averaging over structures obtained from MD simulations 95 or using classical models96 for statistical effects.

ACKNOWLEDGMENT


Page 30 of 58

Page 31 of 58

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


We thank the Next Generation Super Computing Project, Nanoscience Program (MEXT, Japan) and Computational Materials Science Initiative (CMSI, Japan) for financial support.



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 Captions Figure 1. Tree of the components in the binding energy analysis. Dashed boxes imply summation (over I or J). Blue boxes show solvent components. Broad arrows not terminated with boxes mean that a further decomposition is possible similar to other likewise components. Branches of the tree on the same level can be summed yielding the value of the component one level above. The components with and without a tilde (e.g., ~ ∆E I′′ and ∆E I′′ ) differ in the reference (compare eqs 34 and 41). The numbers in the grey

boxes are the energies (kcal/mol) for the complexes of Trp-cage with the ligands, HOPhCOOH (top) and HOPhCOO– (bottom). Figure 2. Fragmentation of polyalanine isomers. Fragments are numbered from 1 to 10. Bonds connecting fragments are shown in grey. Figure 3. Fragment properties of polyalanine isomers α-helix (alpha), β-turn (beta) and extended form (ext): IP, EA and HOMO-LUMO gaps, (FMO-RHF/D/ PCM/6-31G*). Figure 4. Fragment properties of polyalanine isomers: internal solute E I′′ , backbone EIBB and partial E Ipart energies (FMO-RHF/D/ PCM/6-31G*). Figure 5. Fragment properties of polyalanine isomers: ∆E Isolv , solvent coverage s I and induced solvent charges q I (FMO-RHF/D/ PCM/6-31G*).

Figure 6. Fragment properties of polyalanine isomers: dipole moments d I , partial charges QIpart , and unconnected interaction energies ∆E IU,int (FMO-RHF/D/ PCM/631G*). 32 ACS Paragon Plus Environment

Page 32 of 58

Page 33 of 58

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 7. Unconnected pair interaction energies ∆E IJU,int . The grey area in the middle corresponds to excluding terms for connected pairs and self terms I=J. The color legend at the bottom shows brackets for each color, for example, values between -25 and -20 kcal/mol are shown in navy blue (FMO-RHF/D/ PCM/6-31G*). Figure 8. Optimized complexes of Trp-cage with (a) HOPhCOOH and (b) HOPhCOO–. The residue fragments included in the geometry optimization are colored by element, otherwise protein atoms are shown in purple. The ligand atoms are shown as balls and sticks. The change in the ligand (HOPhCOOH) geometry (c) due to complex formation: the optimized isolated and complex geometries are shown as red and green, respectively. Figure 9. Pair interaction energies ∆E ILint between residue fragments I and ligand L, and bind changes in the partial ∆E Ipart , A (Epart) and binding ∆E I , A (Ebind) energies of fragment I in

the complexes of Trp-cage with the ligands (FMO-RHF/D/ PCM/6-31G*). Figure 10. Binding energies ∆E Ibind of fragments I for the complexes of the Trp-cage protein with the neutral and anionic ligands (FMO-RHF/D/ PCM/6-31G*).



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 34 of 58

REFERENCES

(1) Ufimtsev, I. S.; Luehr, N.; Martinez, T. J. Charge Transfer and Polarization in Solvated Proteins from Ab Initio Molecular Dynamics. J. Phys. Chem. Lett. 2011, 2, 1789–1793. (2) Akimov, A. V. ; Prezhdo, O. V. Large-Scale Computations in Chemistry: A Bird’s Eye View of a Vibrant Field. Chem. Rev. 2015, 115, 5797–5890. (3) Scuseria, G. E. Linear Scaling Density Functional Calculations with Gaussian Orbitals. J. Phys. Chem. A 1999, 103, 4782–4790. (4) Zalesny, R.; Papadopoulos, M.G.; Mezey, P.G.; Leszczynski, J. (Eds.), Linear-Scaling Techniques in Computational Chemistry and Physics, Springer, Berlin, 2011. ( 5 ) Reimers, J. R., Ed. Computational Methods for Large Systems: Electronic Structure Approaches for Biotechnology and Nanotechnology, Wiley, New York, 2011. (6) Orimoto, Y.; Yamamoto, R.; Xie, P.; Liu, K.; Imamura, A.; Aoki, Y. Ab Initio O(N) Elongation-Counterpoise Method for BSSE-Corrected Interaction Energy Analyses in Biosystems. J. Chem. Phys. 2015, 142, 104111. (7) Gordon, M. S.; Pruitt, S. R.; Fedorov, D. G.; Slipchenko, L. V. Fragmentation Methods: a Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632–672. ( 8 ) Otto, P.; Ladik, J. Investigation of the Interaction Between Molecules at Medium Distances: I. SCF LCAO MO Supermolecule, Perturbational and Mutually Consistent Calculations for two Interacting HF and CH2O Molecules. Chem. Phys. 1975, 8, 192–200.

ACS Paragon Plus Environment

34

Page 35 of 58

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


(9) Yang, W. Direct Calculation of Electron Density in Density-Functional Theory. Phys. Rev. Lett. 1991, 66, 1438–1441. (10) Gao, J. Toward a Molecular Orbital Derived Empirical Potential for Liquid Simulations. J. Phys. Chem. B 1997, 101, 657–663. (11) Hua, W. J.; Fang, T.; Li, W.; Yu, J. G.; Li, S. H. Geometry Optimizations and Vibrational Spectra of Large Molecules From a Generalized Energy-Based Fragmentation Approach. J. Phys. Chem. A 2008, 112, 10864–10872. (12) Kobayashi, M.; Nakai, H. How Does it Become Possible to Treat Delocalized and/or Open-Shell Systems in Fragmentation-Based Linear-Scaling Electronic Structure Calculations? The Case of the Divide-and-Conquer Method. Phys. Chem. Chem. Phys. 2012, 14, 7629–7639 (13) Söderhjelm, P.; Ryde, U. How Accurate Can a Force Field Become? A Polarizable Multipole Model Combined with Fragment-Wise Quantum-Mechanical Calculations. J. Phys. Chem. A 2009, 113, 617–627. (14) Kiewisch, K.; Jacob, C. R.; Visscher, J. Quantum-Chemical Electron Densities of Proteins and of Selected Protein Sites From Subsystem Density Functional Theory. J. Chem. Theory Comput. 2013, 9, 2425–2440. (15) Frank, A.; Möller, H. M.; Exner, T. E. Toward the Quantum Chemical Calculation of NMR Chemical Shifts of Proteins. 2. Level of Theory, Basis Set, and Solvents Model Dependence. J. Chem. Theory Comput. 2012, 8, 1480–1492.


35


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 36 of 58

(16) Collins, M. A. Systematic Fragmentation of Large Molecules by Annihilation. Phys. Chem. Chem. Phys. 2012, 14, 7744–7751. (17) Yu, H.; Leverentz, H. R.; Bai, P.; Siepmann, J. I.; Truhlar, D. G. Water 26-Mers Drawn From Bulk Simulations: Benchmark Binding Energies for Unprecedentedly Large Water Clusters and Assessment of the Electrostatically Embedded Three-Body and Pairwise Additive Approximations. J. Phys. Chem. Lett. 2014, 5, 660–670. (18) He, X.; Merz, K. M. Jr. Divide and Conquer Hartree−Fock Calculations on Proteins. J. Chem. Theory Comput. 2010, 6, 405–411. (19) Steinmann, C.; Fedorov, D. G.; Jensen, J. H. Effective Fragment Molecular Orbital Method: a Merger of The Effective Fragment Potential and Fragment Molecular Orbital Methods. J. Phys. Chem. A 2010, 114, 8705–8712. (20) Sahu, N.; Gadre, S. R. Molecular Tailoring Approach: a Route for Ab Initio Treatment of Large Clusters. Acc. Chem. Res. 2014, 47, 2739–2747. (21) Huang, L.; Massa, L. Quantum Kernel Applications in Medicinal Chemistry. Future Medicinal Chemistry 2012, 4, 1479–1494. (22) Gao, J.; Truhlar, D. G.; Wang, Y.; Mazack, M. J. M.; Löffler, P.; Provorse, M. R.; Rehak, P. Explicit Polarization: A Quantum Mechanical Framework for Developing Next Generation Force Fields. Acc. Chem. Res. 2014, 47, 2837–2845.


36

Page 37 of 58

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


(23) Jose, K. V. J.; Raghavachari, K. Evaluation of Energy Gradients and Infrared Vibrational Spectra through Molecules-in-Molecules Fragment-Based Approach. J. Chem. Theory Comput.

2015, 11, 950–961. ( 24 ) Liu, J.; Herbert, J. M. Pair–Pair Approximation to the Generalized Many-Body Expansion: An Alternative to the Four-Body Expansion for ab Initio Prediction of Protein Energetics via Molecular Fragmentation. J. Chem. Theory Comput. 2016, 12, 572–584. (25) Liu, J.; Zhang, J. Z. H.; He, X. Fragment Quantum Chemical Approach to Geometry Optimization and Vibrational Spectrum Calculation of Proteins. Phys. Chem. Chem. Phys. 2016, 18, 1864-1875. (26) Vaart, A. van der; Merz, K. M. Jr.The Role of Polarization and Charge Transfer in the Solvation of Biomolecules. J. Am. Chem. Soc. 1999, 121, 9182–9190. (27) Duan, L. L.; Mei, Y.; Zhang, Q. G.; Tang, B.; Zhang, J. Z. H. Protein's Native Structure is Dynamically Stabilized by Electronic Polarization. J. Theor. Comput. Chem. 2014, 13, 1440005. (28) Merz, K. M. Jr. Using Quantum Mechanical Approaches to Study Biological Systems. Acc. Chem. Res. 2014, 47, 2804–2811. (29) 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. ( 30 ) Fedorov, D. G.; Kitaura, K. Extending the Power of Quantum Chemistry to Large Systems with the Fragment Molecular Orbital Method. J. Phys. Chem. A 2007, 111, 6904–6914.


37


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 38 of 58

(31) The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems. Fedorov, D. G.; Kitaura, K., Eds., CRC Press, Boca Raton, FL, 2009. ( 32 ) Fedorov, D. G.; Nagata, T.; Kitaura, K. Exploring Chemistry with the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 2012, 14, 7562–7577. (33) Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. Electron-Correlated Fragment-Molecular-Orbital Calculations for Biomolecular and Nano Systems. Phys. Chem. Chem. Phys. 2014, 16, 10310–10344. (34) Fedorov, D. G.; Kitaura, K. The Importance of Three-Body Terms in the Fragment Molecular Orbital Method. J. Chem. Phys. 2004, 120, 6832–6840. (35 ) 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 Structure-Based Drug Design. J. Mol. Graph. Modell. 2013, 41, 31–42. ( 36 )

Fedorov, D. G.; Asada, N.; Nakanishi, I.;

Kitaura, K. The Use of Many-Body

Expansions and Geometry Optimizations in Fragment-Based Methods. Acc. Chem. Res. 2014, 47, 2846–2856. (37) Richard, R. M.; Lao, K. U.; Herbert, J. M. Aiming for Benchmark Accuracy with the Many-Body Expansion. Acc. Chem. Res. 2014, 47, 2828–2836. (38) Fukuzawa, K.; Kitaura, K.; Uebayasi, M.; Nakata, K.; Kaminuma, T.; Nakano, T. Ab Initio Quantum Mechanical Study of the Binding Energies of Human Estrogen Receptor with its


38

Page 39 of 58

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


Ligands: an Application of Fragment Molecular Orbital Method. J. Comput. Chem. 2005, 26, 1– 10. (39) Amari, S.; Aizawa, M.; Zhang, J.; Fukuzawa, K.; Mochizuki, Y.; Iwasawa, Y.; Nakata, K.; Chuman, H.; Nakano, T.

VISCANA: Visualized Cluster Analysis of Protein-Ligand

Interaction Based on the Ab Initio Fragment Molecular Orbital Method for Virtual Ligand Screening. J. Chem. Inf. Comput. Sci. 2006, 46, 221–230. (40) Nakanishi, I.; Fedorov, D. G.; Kitaura, K. Molecular Recognition Mechanism of FK506 Binding Protein: an all-Electron Fragment Molecular Orbital Study. Proteins: Struct., Funct., Bioinf. 2007, 68, 145–158. (41) 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. (42) Ohno, K.; Mori, K.; Orita, M.; Takeuchi, M. Computational Insights into Binding of Bisphosphates to Farnesyl Pyrophosphate Synthase. Curr. Med. Chem. 2011, 18, 220–233. (43) 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. (44) Ozawa, T.; Okazaki, K.; Kitaura. K. CH/π Hydrogen Bonds Play a Role in Ligand Recognition and Equilibrium between Active and Inactive States of the b2 Adrenergic Receptor: an Ab Initio Fragment Molecular Orbital (FMO) Study. Bioorg. Med. Chem. 2011, 19, 5231– 5237.


39


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 40 of 58

(45) Devarajan, A.; Markutsya, S.; Lamm, M. H.; Cheng, X.; Smith, J. C.; Baluyut, J. Y.; Kholod, Y.; Gordon, M. S.; Windus, T. L. Ab Initio Study of Molecular Interactions in Cellulose Iα. J. Phys. Chem. B 2013, 117, 10430–10443. (46) Kurauchi, R.; Watanabe, C.; Fukuzawa, K.; Tanaka, S. Novel Type of Virtual Ligand Screening on the Basis of Quantum-Chemical Calculations for Protein-Ligand Complexes and Extended Clustering Techniques. Comput. Theor. Chem. 2015, 1061, 12–22. (47) Otsuka, T.; Okimoto, N.; Taiji, M. Assessment and Acceleration of Binding Energy Calculations for Protein-Ligand Complexes by the Fragment Molecular Orbital Method. J. Comput. Chem. 2015, 36, 2209–2218. (48) Yoshida, T., Yamagishi, K., Chuman, H. QSAR Study of Cyclic Urea Type HIV‐1 PR Inhibitors Using Ab Initio MO Calculation of Their Complex Structures with HIV‐1 PR. QSAR Comb. Sci. 2008, 27, 694–703. (49) Mazanetz, M. P.; Ichihara, O.; Law, R. J.; Whittaker, M. Prediction of Cyclin-Dependent Kinase 2 Inhibitor Potency Using the Fragment Molecular Orbital Method. J. Cheminf. 2011, 3, 2. (50) Munei, Y.; Shimamoto, K.; Harada, M.; Yoshida, T.; Chuman, H. Correlation Analyses on Binding Affinity of Substituted Benzenesulfonamides with Carbonic Anhydrase Using Ab Initio MO Calculations on their Complex Structures (II). Bioorg. Med. Chem. Lett. 2011, 21, 141–144.


40

Page 41 of 58

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


(51) Zhang, Q.; Yu, C.; Min, J.; Wang, Y.; He, J.; Yu, Z. Rational Questing For Potential Novel Inhibitors of FabK From Streptococcus Pneumoniae by Combining FMO Calculation, CoMFA 3D-QSAR Modeling and Virtual Screening. J. Mol. Model. 2011, 17, 1483–1492. (52) Hitaoka S.; Chuman H.; Yoshizawa K. A QSAR Study on the Inhibition Mechanism of Matrix Metalloproteinase-12 by Arylsulfone Analogs Based on Molecular Orbital Calculations. Org. Biomol. Chem. 2015, 13, 793–806. ( 53 ) Alexeev, Y.; Mazanetz, M. P.;

Ichihara, O.; Fedorov, D. G. GAMESS as a Free

Quantum-Mechanical Platform for Drug Research. Curr. Top. Med. Chem. 2012, 12, 2013–2033. (54) Phipps, M. J. S.; Fox, T.; Tautermann, C. S.; Skylaris, C.-K. Energy Decomposition Analysis Approaches and their Evaluation on Prototypical Protein–Drug Interaction Patterns. Chem. Soc. Rev. 2015, 44, 3177–3211. (55) Heifetz, A., Chudyk, E. I.; Gleave, L.; Aldeghi, M.; Cherezov, V.; Fedorov, D. G.; Biggin, P. C.; Bodkin, M. J. The Fragment Molecular Orbital Method Reveals New Insight into the Chemical Nature of GPCR-Ligand Interactions. J. Chem. Inf. Model. 2016, 56, 159–172. (56) Murata, K.; Fedorov, D. G.; Nakanishi, I.; Kitaura, K. Cluster Hydration Model for Binding Energy Calculations of Protein-Ligand Complexes. J. Phys. Chem. B 2009, 113, 809– 817. (57) Mulliken, R. S. Molecular Scientists and Molecular Science: Some Reminiscences. J. Chem. Phys. 1965, 43, S2.


41


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 42 of 58

(58) Jacob, C. R.; Neugebauer, J. Subsystem Density-Functional Theory. WIREs Comput. Mol. Sci. 2014, 4, 325–362. (59) 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. (60) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999–3093. ( 61 ) 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. (62) 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. (63) 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. (64) Fedorov, D. G.; Kitaura, K. Pair Interaction Energy Decomposition Analysis. J. Comput. Chem. 2007, 28, 222–237. (65) 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.


42

Page 43 of 58

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


(66) Nakata, H.; Fedorov, D. G. ; Yokojima, S.; Kitaura, K.; Nakamura, S. Derivatives of the Approximated Electrostatic Potentials in Unrestricted Hartree-Fock Based on the Fragment Molecular Orbital Method and an Application to Polymer Radicals. Theor. Chem. Acc. 2014, 133, 1477. (67) 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. (68) Chang, L.; Ishikawa, T.; Kuwata, K.; Takada, S. Protein-Specific Force Field Derived from the Fragment Molecular Orbital Method Can Improve Protein-Ligand Binding Interactions. J. Comput. Chem. 2013, 34, 1251–1257. (69) Sekino, H.; Sengoku, Y.; Sugiki, S.-I.; Kurita, N. Molecular Orbital Analysis Based on Fragment Molecular Orbital Scheme. Chem. Phys. Lett. 2003, 378, 589–597. (70) Nakata, H.; Fedorov, D. G.; Yokojima, S.; Kitaura, K.; Sakurai, M.; Nakamura, S. Unrestricted Density Functional Theory Based on the Fragment Molecular Orbital Method for the Ground and Excited State Calculations of

Large Systems. J. Chem. Phys. 2014, 140,

144101. (71) Fukui, K.; Yonezawa, T.; Shingu, H. A Molecular Orbital Theory of Reactivity in Aromatic Hydrocarbons. J. Chem. Phys. 1952, 20, 722.


43


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 44 of 58

(72) Hoffmann, R. Interaction of Orbitals through Space and through Bonds. Acc. Chem. Res.

1971, 4, 1–9. (73) Fedorov, D. G.; Kitaura, K. The Role of the Exchange in the Embedding Electrostatic Potential for the Fragment Molecular Orbital Method. J. Chem. Phys. 2009, 131, 171106. ( 74 ) Tsuneyuki, S.; Kobori, T.; Akagi, K.; Sodeyama, K.; Terakura, K.; Fukuyama, H. Molecular Orbital Calculation of Biomolecules with Fragment Molecular Orbitals. Chem. Phys. Lett. 2009, 476, 104–108. (75) Kobori, T.; Sodeyama, K.; Otsuka, T.; Tateyama, Y.; Tsuneyuki, S. Trimer Effects in Fragment Molecular Orbital-Linear Combination of Molecular Orbitals Calculation of OneElectron Orbitals for Biomolecules. J. Chem. Phys. 2013, 139, 094113. (76) Mulliken, R. S. Electronic Population Analysis on LCAO-MO Molecular Wave Functions. J. Chem. Phys. 1955, 23, 1833–1840. (77) Nakai, H. Energy Density Analysis with Kohn–Sham Orbitals. Chem. Phys. Lett. 2002, 363, 73–79. (78) Parrish, R. M.; Sherrill, C. D. Spatial Assignment of Symmetry Adapted Perturbation Theory Interaction Energy Components: the Atomic SAPT Partition. J. Chem. Phys. 2014, 141, 044115. (79) Fedorov, D. G.; Olson, R. M.; Kitaura, K.; Gordon, M. S.; Koseki, S. A new Hierarchical Parallelization Scheme: Generalized Distributed Data Interface (GDDI), and an Application to the Fragment Molecular Orbital Method (FMO). J. Comput. Chem. 2004, 25, 872–880.


44

Page 45 of 58

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


(80) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., Jr. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347–1363. (81) 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. (82) Suenaga, M. Development of GUI for GAMESS/FMO Calculation. J. Comput. Chem. Jpn. 2008, 7, 33 (in Japanese). (83) Fedorov, D. G. ; Ishida, T.; Uebayasi, M.; Kitaura, K. The Fragment Molecular Orbital Method for Geometry Optimizations of Polypeptides and Proteins. J. Phys. Chem. A 2007, 111, 2722–2732. (84) Neidigh, J. W.; Fesinmeyer, R. M.; Andersen, N. H. Designing a 20-Residue Protein. Nat. Struct. Biol. 2002, 9, 425–430. ( 85 )

Pitera, J. W.; Swope, W. Understanding Folding and Design: Replica-Exchange

Simulations of “Trp-Cage'' Miniproteins. Proc. Nat. Acad. Sc. U.S.A. 2003, 100, 7587–7592. (86) Pascheka, D.; Hempela, S.; Garcíab, A. E. Computing the Stability Diagram of the TrpCage Miniprotein. Proc. Nat. Acad. Sc. U.S.A. 2008, 105, 17754–17759.


45


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 46 of 58

(87) Zhu, T.; He, X.; Zhang, J. Z. H. Fragment Density Functional Theory Calculation of NMR Chemical Shifts for Proteins with Implicit Solvation. Phys. Chem. Chem. Phys. 2012, 14, 7837– 7845. ( 88 ) Tsukamoto, T.; Mochizuki, Y.; Watanabe, N.; Fukuzawa, K.; Nakano, T. Partial Geometry Optimization with FMO-MP2 Gradient: Application to TrpCage. Chem. Phys. Lett.

2012, 535, 157–162. (89) Brorsen, K. R.; Zahariev, F.; Nakata, H.; Fedorov, D. G.; Gordon, M. S. Analytic Gradient for Density Functional Theory Based on the Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2014, 10, 5297–5307. (90) Fedorov, D. G.; Alexeev, Y.; Kitaura, K. Geometry Optimization of the Active Site of a Large System with the Fragment Molecular Orbital Method. J. Phys. Chem. Lett. 2011, 2, 282– 288. (91) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (92) Nakata, H.; Nagata, T.; Fedorov, D. G.; Yokojima, S.; Kitaura, K.; Nakamura, S. Analytic Second Derivatives of the Energy in the Fragment Molecular Orbital Method. J. Chem. Phys. 2013, 138, 164103.


46

Page 47 of 58

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


(93) Cammi, R.; Tomasi J. Remarks on the Use of the Apparent Surface Charges (ASC) Methods in Solvation Problems: Iterative Versus Matrix-Inversion Procedures and the Renormalization of the Apparent Charges. J. Comput. Chem. 1995, 16, 1449–1458. ( 94 ) Jensen, J. H. Predicting Accurate Absolute Binding Energies in Aqueous Solution: Thermodynamic Considerations for Electronic Structure Methods. Phys. Chem. Chem. Phys.

2015, 17, 12441–12451. (95) Ishikawa, T.; Ishikura, T.; Kuwata, K. Theoretical Study of the Prion Protein Based on the Fragment Molecular Orbital Method. J. Comput. Chem. 2009, 30, 2594–2601. (96) Tanaka, S.; Watanabe, C.; Okiyama, Y. Statistical Correction to Effective Interactions in the Fragment Molecular Orbital Method. Chem. Phys. Lett. 2013, 556, 272–277.


47


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

Total binding energy

-8.0

~ of A and B, E

Direct binding

-13.7

energy E

-11.1

Page 48 of 58

-9.5

Deformation energy

0.5

~ def, A

Deformation energy

~ def, B

0.3

of A, E

5.2 1.3

of B, E

Change in partial energy Change in energy of A, E A

9.2

Change in energy

23.2

of B, E B

16.7

interaction energy,

E

32.6

Change in partial energy

Change in partial energy

of I in A, E I , A

of I in B, E I , B

part

int AB

-39.6

~ def, part

of I in A, E I , A

-66.9

part

E

int IJ

Change in

internal energy

solvation energy

interactions of I,

~

Change in

Change in

internal energy

solvation energy

interactions of I,

of I, E

Change in

solv I

~ solv

of I, E I

Change in of I, E I

Change in

interactions

1 E IJint  2 J A

of I, E I

Change in individual interactions

Change in

Change in

Change in

Change in

cavitation

electrostatic

dispersion

repulsion

energy of

energy of I,

energy of

energy of

I, E Icav

E Ies

I, E Idisp

I, E Irep

Figure 1. Fedorov et al.

1 ~ E IJint  2 J A

1 E IJint 2

Change in

Change in

Change in

Change in

Change in

electrostatic,

exchange-

charge transfer

dispersion,

screening,

1 E IJES 2

repulsion,

+ mixed terms

1 E IJDI 2

1 E IJsolv 2

1 / 2E

ES IJ

1 / 2E

CT  mix IJ

Change in

Change in

Change in

Change in

ES screening,

CT solvation,

DI screening,

EX screening,

1 E IJCTes 2

1 E IJdisp 2

1 E IJrep 2

1 2

ACS Paragon Plus Environment E IJ

es

Page 49 of 58

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





Ionization potential, eV

Fragment ionization potential 12 12 11 11

alpha

10 1

2

3

4

beta

5 6 7 fragment

ext 8

9 10

Electron affinity, eV

Fragment electron affinity 6.0 5.5 5.0 4.5 1

2

3

4

5 6 7 fragment

8

9 10

HOMO-LUMO Gap 17.0 Gap, eV

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

16.5 16.0 15.5 1

2

3

4

5 6 7 fragment

8

9 10



Page 50 of 58

Backbone energy, kcal/mol

Partial energy, kcal/mol

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


Internal energy, kcal/mol

Page 51 of 58

Internal fragment solute energy alpha

-144,730.0

beta

ext

-144,740.0 -144,750.0 -144,760.0 2

3

4

5 6 7 fragment

8

9

8

9

8

9

Fragment backbone energy -154,270.0 -154,280.0 -154,290.0 -154,300.0 2

3

4

5 6 7 fragment

Fragment partial energy -154,270.0 -154,280.0 -154,290.0 -154,300.0 2

3

4

5 6 7 fragment



Solvent coverage, %

Fragment solvation energy 5 0 -5 alpha

-10 1

2

3

4

beta

5 6 7 fragment

ext 8

9

10

9

10

Fragment solvent coverage 25 20 15 10 5 1

2

3

4

5 6 7 fragment

8

Induced solvent charge Solvent charge, a.u.

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

Solvation energy, kcal/mol


0.5 0.3 0.1 -0.1 -0.3 -0.5 1

2

3

4

5 6 7 fragment

8

9

10



Page 52 of 58

Page 53 of 58

Dipole moment, Debye

Fragment dipole moment 12.0

alpha

beta

ext

7.0

2.0 1

2

3

4

5 6 7 fragment

8

9 10

Fragment charge, a.u.

Fragment partial charge 0.10 0.05 0.00 -0.05 -0.10 1

2

3

4

5 6 7 fragment

8

9 10

Unconnected interaction energy 5.0 Energy, kcal/mol

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


0.0 -5.0 -10.0 -15.0 -20.0 1

2

3

4 5 6 7 fragment

8

9 10




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 54 of 58

Page 55 of 58

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


(a)

(b)

(c)




Neutral ligand HOPhCOOH energy, kcal/mol

5 0 -5 -10 -15 -20

Epart

PIE

Ebind

ASN-1 LEU-2 TYR-3 ILE-4 GLN-5 TRP-6 LEU-7 LYS-8 ASP-9 GLY-10 GLY-11 PRO-12 SER-13 SER-14 GLY-15 ARG-16 PRO-17 PRO-18 PRO-19 SER-20

-25

fragment

Anionic ligand HOPhCOO– 20 energy, kcal/mol

0 -20 -40 ASN-1 LEU-2 TYR-3 ILE-4 GLN-5 TRP-6 LEU-7 LYS-8 ASP-9 GLY-10 GLY-11 PRO-12 SER-13 SER-14 GLY-15 ARG-16 PRO-17 PRO-18 PRO-19 SER-20

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

fragment Figure 9. Fedorov et al.


Page 56 of 58

Page 57 of 58

Ligand binding energy contributions 35 energy, kcal/mol

25

neutral

15

anionic

5 -5 -15 -25 ASN-1 LEU-2 TYR-3 ILE-4 GLN-5 TRP-6 LEU-7 LYS-8 ASP-9 GLY-10 GLY-11 PRO-12 SER-13 SER-14 GLY-15 ARG-16 PRO-17 PRO-18 PRO-19 SER-20 ligand

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


fragment Figure 10. Fedorov et al.



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

TOC image. Fedorov et al.


Page 58 of 58

Subsystem Analysis for the Fragment Molecular Orbital Method and Its

Recommend Documents