The quantum topological electrostatic potential as a probe for functional group transferability

The electrostatic potential can be used as an appropriate and convenient indicator of how transferable an atom or functional group is between two molecules. Quantum-chemical topology (QCT) is used to define the electron density of a molecular fragment and the electrostatic potential it generates. The potential generated on a grid by the terminal aldehyde group of the biomolecule retinal is compared with the corresponding aldehyde group in smaller molecules derived from retinal. The terminal amino group in the free amino acid lysine was treated in a similar fashion. Each molecule is geometry-optimized by an ab initio calculation at B3LYP/6-311G+(2d,p)//HF/6-31G(d) level. The amino group in lysine is very little influenced by any part of the molecule further than two C atoms away. However, the aldehyde group in retinal is influenced by molecular fragments six C atoms away. This dramatic disparity is ascribed to the difference in saturation in the carbon chains; retinal contains a conjugated hydrocarbon chain but lysine an aliphatic one.

1. Introduction

The transferability of atomic or functional properties is closely related to the issue of partitioning of a molecular system into subgroups, which has received considerable attention in the literature (see §2). If it were necessary to express each atom in a molecule as a perfectly unique entity then force-field design would be futile. Moreover, the concept of a functional group, an important cornerstone of chemistry, would evaporate. A well documented approach to the issue of transferability is given by quantum-chemical topology (QCT) (Bader, 1990, 1991; Popelier, 2000; Popelier & Smith, 2002). The theorem of Hohenberg & Kohn (1964) was extended to subsystems by Riess & Munch (1981) and subsequently adopted by Bader & Becker (1988). Their work laid the basis (Bader, 1995) for the application of QCT in the context of transferability. Several studies focused on the related question of the reconstruction of molecular properties from the constituent atomic fragments. They include work (Bader et al., 1987; Laidig, 1992) on molecular multipole moments, group additivity (Bader & Bayles, 2000), the library of transferable atomic equivalents (TAE) (Breneman et al., 1995; Breneman & Rhem, 1997), the transferability of atomic and bond properties of aldehydes and ketones (Grana & Mosquera, 2000), the construction of a polypeptide (Chang & Bader, 1992), the construction of pripavine PEO, enkephalins and morphine (Matta, 2001) from a combination of fragments, the effect of conformation on the atomic properties of amino acids (Matta & Bader, 2000) and oligopeptides (Popelier & Bader, 1994).

Recently, we have extracted atom types (Popelier & Aicken, 2003a,b) from the topologically partitioned electron densities of all 20 natural free amino acids and smaller derived mol­ecules. Cluster analysis and a statistical separation criterion enabled us to arrange the 760 (unique) QCT atoms into classes with given degrees of transferability in atomic volume, energy and multipole moments up to hexadecapole. The fact that a computational scheme generated the atom types makes them less arbitrary than the largely guessed atom types of typical force fields, which may therefore benefit from suggested recommendations (Popelier & Aicken, 2003c). Amino acids with long side chains, such as lysine and arginine, are ideal to study (Aicken, 2000) how the multipole moments of terminal atoms vary when appearing in derived molecules with systematically shortened chains. In view of its general importance (Popelier & Rafat, 2003, and references therein) as a physical observable, we deemed it more useful to study the electrostatic potential, which can be obtained as a series expansion from the multipole moments. Since its inception two and a half decades ago (Scrocco & Tomasi, 1978), the molecular electrostatic potential has been recognized as an informative physical observable. Amongst its wide range of applications, we quote its role in zeolite catalysis (White & Hess, 1993), in the prediction of adsorption sites on crystalline oxide surfaces (McCarthy & Hess, 1992), in hydration models for proteins (Wodak et al., 1990) and nucleic acids (Hausheer et al., 1990), and as an indicator for evolutionary conservativeness near enzyme active sites (Desideri et al., 1992), to name a few.

We have chosen two biological systems: retinal and lysine. Retinal bears an aldehyde incorporated into an extended conjugated chain. Lysine exhibits an amino group at the terminus of a long aliphatic chain. Retinal, the prosthetic group of the membrane protein bacteriorhodopsin, forms a Schiff base with Lys 216 of this protein (Voet & Voet, 1995). Cleavage of the covalent CN bond of the Schiff base leads to a free amino group in Lys and an aldehyde group in retinal. We will report on each functional group (NH₂ and CHO) in turn, after a review of the various approaches to transferability and a section on computational methods and technical details.

2. Review of approaches to transferability

It is useful to review alternative partitioning schemes and their potential to detect transferability. Subsystem properties may be defined by partitioning the orbital space, the electron density or the Hamiltonian. Within these main categories, many different techniques have been applied for a variety of purposes. Where insight is sought, this may take the form of population assignment, analysis of bonding effects or the decomposition of intra- and intermolecular energies.

2.1. Partitioning of orbital space

Assigning orbitals to particular atoms is a very popular means of acquiring atomic properties from the molecular wavefunction. Mulliken population analysis is a well documented precursor to many other localized orbital techniques. Basis functions are assigned to particular atoms and hence the properties calculated for a particular function are associated with the atom on which the function is centred. Mulliken populations are very sensitive to basis set, particularly as the basis set is improved to higher accuracy. They are inadequate when modelling ionic systems (Reed et al., 1985). Löwdin recommended more realistic partitioning (Lowdin, 1953). Natural population analysis (NPA) (Reed et al., 1985) constructs an orthonormal set of natural atomic orbitals (NAOs) for a given molecule in an arbitrary atomic orbital basis set. NPA precludes some of the problems inherent in Mulliken population analysis while preserving computational economy.

Delocalized molecular orbitals do not always recover the particular types of bonding patterns that are detected experimentally. Localized molecular orbital (LMO) methods reconcile MO theory with chemical intuition of bonding. LMOs are defined as those orthonormal MOs that minimize the total interorbital charge-cloud repulsion energy. This method, originally defined by Lennard-Jones and Pople, was implemented for many molecular systems by Edmiston & Ruedenberg (1963) and Boys (Foster & Boys, 1960).

Another method that employs LMOs, called natural bond orbital (NBO) (Reed et al., 1988) analysis was originally designed to study hybridization and covalency in polyatomic wavefunctions. NBOs were formulated to correspond closely with chemists' ideas of localized bonds and lone pairs as basic units of molecular structure. NBOs are obtained by a transformation of the original basis set via NAOs and are themselves transformed into natural LMOs (NLMOs). The NLMO expansions assess the transferability of NBOs and LMOs from one molecule to another with the result that NBOs are intrinsically more transferable than LMOs. Attempts to transfer LMOs from one system to another have failed because of the sensitivity of the properties to changes in the vicinal groups (Christoffersen, 1972; Shipman & Christoffersen, 1973; Kapuy et al., 1980). The LAOs are found not to be entirely localized and, in fact, contain the basis functions of atoms in vicinal positions with coefficients of appreciable magnitude. The overlapping tails are found to be characteristic of the neighbouring atoms and any changes incurred in the tails are magnified in the exchange interaction (Kapuy et al., 1979).

Stone defined a distributed multipole analysis (DMA) (Stone, 1981) method, enabling the calculation of the electrostatic potential and interaction energy in regions close to the molecule. Faerman & Price (1990) used DMA in the construction of a transferable distributed multipole model for peptides and amides. They investigated the construction of distributed multipole models of polypeptide molecules from data generated from ab initio calculations performed on smaller molecules. They found that atoms exhibiting similar charges also exhibit similar higher multipole moments in the local axis system. Saturated hydrocarbon fragments were found to be more sensitive to effects induced by more distant parts of the molecule. DMAs of molecular fragments could only be transferred between molecules where the functional groups bonded to the fragment are the same. Dipeptides were recommended as the smallest realistic unit to model the peptide residue.

2.2. Partitioning of the Hamiltonian: perturbation techniques

Other than determining the correlation energy, the most common application of the partitioning of the Hamiltonian is to achieve insight into intermolecular interactions. The energy can be partitioned into its constituent parts: electrostatic, polarization, exchange repulsion and charge transfer (Kitaura & Morokuma, 1976). Korchowiec & Uchimaru (2000) considered the energy breakdown in the same terms as Kitaura & Morokuma except that they include a term to account for the deformation. This method improves on the Kitaura & Morokuma model by being less basis-set dependent and it also avoids the basis-set superposition-error correction.

Some of these techniques are applied in deriving transferable entities. Nicolas & Durand (1979) proposed a theoretical methodology for deriving transferable atomic potentials in molecules. The transferable quantities are derived from the ab initio wavefunction. Their methodology is not basis-set dependent and able to generate effective potentials from an approximate wavefunction.

Unger (1997) applied density functional theory in the local density approximation of Becke's scheme to decompose the total potential of a system into its atomic contributions. The partition function depends on the atomic configuration of the system to which it is applied. Becke (1988) defined a scheme for generating partition functions that comply with these conditions. This partitioning method produces geometrical cells that coincide qualitatively with the atoms defined using QCT. Becke's method also produces fuzzy cells by smoothing out discontinuities at the boundaries. For homonuclear entities, the cells take the form of Wigner–Seitz cells and, for disordered systems, diffuse Voronoi polyhedra are generated.

This method of dividing molecules into their substituent fragments is very similar to the divide-and-conquer technique prescribed by Yang (1991). The divide-and-conquer method has been defined in terms of the partitioning of the Hamil­tonian and also as a projection of the electron and energy density. Zhou's method (Zhou, 1995) operates by the manipu­lation of basis sets.

2.3. Partitioning of the electron density

Another means of fragmenting a molecule is to assign portions of the electron density to different nuclei. Wigner & Seitz (1933) explored the constitution of metallic sodium with regard to its chemical properties. They simplified the energy calculation by applying the symmetry of the crystal to define subsystems within the lattice. As the derivative of the wavefunction is zero perpendicular to every crystallographic symmetry plane, they employ this as their boundary condition. Becke's partitioning scheme (Becke, 1988) recovers the same geometric entities for a crystalline homonuclear system.

Politzer & Harris (1970) defined atomic fragments that recover the charge on each atom as an unambiguous property of the molecular electron distribution, independent of the formalism used to derive the distribution. The fragments are somewhat artificially bounded by an array of partition planes. Bonded atoms with plane faces are not intuitively compatible with quantum mechanics and are also awkward to deal with mathematically. The atoms were defined within the context of diatomic molecules and the extension of this model to non-linear molecules is dubious.

Hirshfeld (1977) proposed a remedy to these and other problems. As mentioned above, the partitioning of orbital space is plagued by basis-set dependence. Stone & Alderton (1985) suggested following Hirshfeld's stockholder partitioning method to lessen this effect. The generality of this methodology permits application to any system, unlike the very restrictive Politzer–Harris technique. As the atomic deformation densities vanish for the promolecule, they are well suited to describing the molecular properties.

Crystallographers treat bonded atoms as point nuclei with a spherically symmetrical electron distribution. Stewart and co-workers used a spherical approximation of the H atom in its bonded state to deduce its X-ray form factors (Stewart et al., 1965). The atomic density formalism is based around the definition of the scattering amplitude. The systematic error by which the model breaks down is minimized when the atomic density function best describes the electron density of the atom in a least-squares sense (Gill, 1996).

Walker & Mezey have been involved in the construction of molecular density fragments using the molecular electron density lego assembler (MEDLA) (Walker & Mezey, 1995). The MEDLA program uses either calculated or experimentally determined densities to construct the electron density for any conformation of the molecule. The fragments are generated according to a simple electron-density fragment additivity principle (Walker & Mezey, 1994). Once assembled, the analysis undergone is strictly visual. The molecular isodensity contour (MIDCO) surfaces are compared at different threshold levels. Exhaustive visual shape analysis of the surfaces is however impractical, and the study of the MIDCOs should be undertaken by non-subjective non-visual computer-based techniques for a more reliable outcome.

In summary, the partitioning of orbital space is problematic because of the dependence on the basis-set level and also the presence of orbital tails resulting in high sensitivity of the properties to changes in the vicinal groups. The decomposition of the Hamiltonian is most commonly used to gain insight into the stability of complexes and their intermolecular interactions. Where the method is utilized for the construction of transferable fragment properties, we find that the introduction of some necessary approximations is equivalent to the errors introduced in the orbital space partitioning upon basis-set truncation. Other than for QCT, the partitioning of the electron density is problematic because of the arbitrariness of defining the boundaries.

3. Computational methods

Fig. 1 shows how five smaller molecules were generated from lysine by retaining its terminal NH₂ group whilst systematically shortening the side chain, starting with its Cα end. For example, by cleaving the Cα—Cβ bond in molecule (6) (lysine) and substituting with an H atom, we obtain molecule (5). Similarly, molecule (4) is obtained from molecule (5) by cleaving and capping the bond to the terminal methyl group (i.e. corresponding to the Cβ—Cγ bond in lysine). Both conformers of retinal are shown schematically in Fig. 2 and also in a three-dimensional representation in Fig. 3, labelled as molecules (15) and (16). All `derived molecules' stem from cis-retinal. Each of the methyl groups attached to the ring in retinal can be either left intact or omitted, leading to a set of 2³ = 8 possible molecules, labelled from (8) to (15). Molecule (7) is obtained from (8) by chopping off the ring from the chain. Molecule (6) is obtained from (7) by omitting two chain C atoms or a methyl group, a process that is repeated until molecule (1) (methanal) is reached. Each molecule, from either the retinal or the lysine set, was geometry optimized at HF/6-31G* level using GAUSSIAN98 (GAUSSIAN98, 1998). This allows for small variations in the geometries across the molecules but keeps their essential conformation intact. Single-point wavefunctions were generated at the B3LYP/6-311+G(2d,p) level (Foresman & Frisch, 1996) and passed on to the computer program MORPHY98 (1998) for the computation of atomic multipole moments according to spherical tensor formalism (Stone, 1996), which is more compact than the Cartesian formalism.

Figure 1 Lysine (6) and its set of smaller derived molecules (1)–(5).

Figure 2 Schematic representation of the (a) trans and (b) cis conformers of retinal.

Figure 3 cis-Retinal (15) and trans-retinal (16) and `derived molecules' (1)–(14) of cis-retinal.

A local (unreleased) version of MORPHY calculated the electrostatic potential due to the CHO group (in the retinal set) or the NH₂ group (in the lysine set) on a grid of points lying on the ρ = 10⁻⁶ a.u. electron-density contour surface. These grid points lie at about 1.6–1.7 Å outside the ρ = 10⁻³ a.u. contour surface. Fig. 4 shows an example of such a grid. The `mean absolute difference', denoted by Δ, served as a simple statistical measure to gauge the discrepancy between two grids. It is defined by

where N is the number of grid points (equal in both grids) and V_i^j is the potential at the ith grid point of the jth grid, that is, the jth molecule. We found it useful to compare all grids with the one of the `reference' molecule, which is cis-retinal (15) or lysine (6). An alternative measure in which the absolute difference is replaced by a square does not alter our conclusions.

Figure 4 Example of a grid constructed around the terminal CHO group, consisting of points at the ρ = 10^−6 a.u. electron-density contour surface of molecule (3) in the retinal set. The electrostatic potential (kJ mol^−1), solely generated by the CHO group, is less than −50 (purple) or varies from −45 to −30 (orange), −30 to −25 (pink) or −25 to −10 (green).

How can one map a given grid point in one molecule to a corresponding one in another molecule, given that the contour surfaces may differ from one molecule to another? We explain our solution in detail for CHO. A local axis system with origin on the C nucleus for CHO enables the identification of a grid point relative to the given molecule and hence a 1:1 mapping across molecules. More specifically, a local x axis connects the local origin on C and the nearest C nucleus; so it is mounted along the C—C bond. The y axis was put in the plane determined by the three nuclei constituting the respective functional group, with the z axis orthogonal to this plane. Spherical coordinates were installed with respect to the local axis system in the conventional way, labelling a grid point as (r,θ,φ). A grid is then set up in molecule A by varying the azimuthal angle θ in equal steps. In order to make the grid points fairly equally spaced (see Fig. 4), the number of φ angles is linearly varied so as to set it to zero when θ is 0 or 180° and maximize it to twice the number of θ values, when θ = 90°. As such, each grid point is labelled and uniquely defined by two angular coordinates (θ,φ). The 1:1 mapping between grid points in mol­ecules A and B is then automatically established by realizing that, in principle, there is a grid point with identical (θ,φ) values in molecule B for every given point (θ,φ) in molecule A. We ensured that each grid point lies in any of the atomic basins that constitute the functional group generating the potential. Current tests and conclusions from earlier work (Kosov & Popelier, 2000a,b) ensured that the multipole expansions converged to the correct exact electrostatic potential.

How many points should be included on a grid? Tests with fine and coarse grids demonstrated that changing the resolution of the grid from 40 points per 180° to 20 points has very little effect on the results, obtained by comparing the two grids point by point. This change in resolution on retinal's aldehyde group represents a change from ~800 points to ~200 points.

4. Results and discussion

4.1. Transferability of retinal's aldehyde group

Fig. 5 shows how the mean absolute difference Δ decreases from molecule (1) to (15) of the retinal set. It will not come as a surprise that the average discrepancy in the electrostatic potential between the CHO fragment in methanal (1) and in cis-retinal (15) is as much as 29 kJ mol⁻¹. This discrepancy immediately falls off to 13 kJ mol⁻¹ when the H atom next to the CHO group (in methanal) is replaced by CH₂C=C (in molecule 2). As retinal's hydrocarbon chain grows incrementally, Δ decreases steadily. It is only when retinal's full conjugated chain, with the exception of the ring, is attached to the CHO group (in molecule 7) that Δ reaches a value as low as 2 kJ mol⁻¹. Another way of viewing the influence of the conjugated carbon chain on CHO is to focus on molecule (5). Its terminal C atom (left in Fig. 3) is six C atoms from the CHO group. Yet, if both H atoms attached to this terminal C atom are replaced by carbon-containing fragments (e.g. methyl), then Δ varies by 3.6 kJ mol⁻¹. Although still under the typical 1 kcal mol⁻¹ (or 4.2 kJ mol⁻¹) threshold, it is remarkably large given the distance. However, variations in the presence of methyl groups attached to the ring – as embodied by mol­ecules (8) to (14) – all have effects on the potential of less than 1 kJ mol⁻¹. Interestingly, the discrepancy between cis- (15) and trans-retinal (16) is again below 1 kJ mol⁻¹. In a wider context, it was found (Devereux & Popelier, 2004) that conformational change had less impact on local atomic and group properties than changes in covalently bonded partners.

Figure 5 The mean absolute difference in electrostatic potential, Δ (kJ mol^−1), between points on a grid for each molecule of the retinal set. The reference grid is that of molecule (15), which is cis-retinal.

In Fig. 6, we examine the change in the population and the magnitude of the dipole and quadrupole moment of the O atom in the CHO group throughout the set of retinal mol­ecules. Starting with the smallest molecule, methanal (1), the O atom gradually loses its population, by 0.03 e in total, as the conjugated C-atom chain grows. The profile in Fig. 6(a) looks roughly similar to that of Fig. 5, indicating a prominent role of this term in the multipole expansion of the electrostatic potential. The profiles of the dipole and quadrupole moment (Figs. 6b and c) drop off as slowly as that of the population. However, they are expected to contribute less to the potential because of the higher inverse power of their distance dependence. The significant effect of remote groups on the O-atom charge must be due to conjugation allowing electron density to be transferred to the aldehyde O atom through several bonds.

Figure 6 Variation of the magnitude of the multipole moments of the O atom in the aldehyde group in the retinal set. (a) Population, (b) dipole moment, (c) quadrupole moment.

4.2. Transferability of lysine amino group

Fig. 7 shows how the mean absolute difference Δ decreases from molecule (1) to (6) of the lysine set. Contrasting this figure with Fig. 5 makes it immediately clear that the influence of remote groups on the potential of NH₂ is quickly dampened. Only when the C∊H₂—CδH₂—CγH₂—CβH₂—CαH(COOH)—NH₂ is replaced by a methyl group does Δ become significant for the first time and is raised above 2 kJ mol⁻¹. It is obvious that non-conjugated methylene groups insulate the terminal amino group from chemical change elsewhere in the molecule. Overall, Figs. 5 and 7 suggest that non-conjugated systems have more transferable functional groups than conjugated systems.

Figure 7 The mean absolute difference in electrostatic potential, Δ (kJ mol^−1), between points on a grid for each molecule of the lysine set. The reference grid is that of molecule (6), which is lysine.

Fig. 8 shows the change in the population of the N atom in the amino group throughout the set of lysine molecules. Going from ammonia (molecule 1) to lysine, the N atom loses almost 0.04 e to the growing alkyl chain of molecules (2) to (6). This amount is remarkably similar to that lost by the O atom in the CHO group. Overall, the population (and the magnitude of the dipole and quadrupole moment, not shown) drop off more quickly than the corresponding values for CHO's O atom. This is in line with the different behaviour of Δ between the lysine and retinal sets.

Figure 8 Variation of the population of the N atom in the amino group in the lysine set.

5. Conclusions

The transferability of the electrostatic potential of two functional groups, aldehyde and amino, has been investigated in a set of molecules derived from retinal and lysine, respectively. The functional groups consisted of topological atoms as defined by quantum-chemical topology (QCT). In the lysine set of molecules, alterations in fragments bonded to the C atom that is bonded to the NH₂ cause a change of ~2 kJ mol⁻¹ in the mean absolute difference of the electrostatic potential generated by the NH₂ group. However, in the retinal set of molecules, the mean absolute difference in electrostatic potential of the CHO group differs by 1–2 kJ mol⁻¹ for chain modifications as far as six C atoms away. Our findings suggest that saturated alkyl chains insulate groups from remote change far more effectively than conjugated systems. This conclusion benefits the design of a peptide force field because each amino acid residue is shielded from the backbone by a methylene group. Long conjugated chains, and most likely also aromatic fragments, are best modelled with data drawn from the full fragment.