Virtual rigid body dynamics.

Virtual Rigid Body Dynamics TERESA HEAD-GORDON* and CHARLES 1. BROOKS 111

Department of Chemistry, Carnegie-Mellon University, Pittsburgh, Pennsylvania 1521 3

SYNOPSIS

An important direction in biological simulations is the development of methods that permit the study of larger systems and/or longer simulation time scales than is currently feasible by molecular dynamics. One such method designed with this objective in mind is stochastic boundary molecular dynamics (SBMD) . SBMD was developed for the characterization of spatially localized processes in proteins, and has been shown to successfully reproduce structural and dynamical properties of these macromolecules, as compared to a molecular dynamics control simulation, when concerted or global motions are not present. The virtual rigid body dynamics method presented in this work extends the range of applicability of the SBMD method, by providing a framework to include these important long time scale conformational transitions. In this paper we describe the two-step implementation of the virtual rigid body model: first, the reduction of the full atomic representation to a reduced particle (virtual bond) model, and second, the propagation of the dynamics of flexibly connected rigid bodies containing virtual atom sites.

INTRODUCTION Conventional simulation methods such as molecular dynamics have been very successful in application to proteins, but the evaluation of the nonbonded potential energy function and derivative is very expensive computationally due t o the large number of atoms in these systems ( t h e calculation scales as N 2 , where N is the number of atoms). Therefore, protein simulations are limited to relatively small proteins and short ( - 100 p s ) time scales. Novel methods such as stochastic boundary molecular dynamics'-' ( SBMD) , developed for studies of localized processes such as protein active sites, have successfully addressed the computational cost problem of molecular dynamics, thereby extending the size and time scale range of protein simulations. However, when global protein motions are present, methods like SBMD are not appropriate since large conformational changes are not adequately sampled. Furthermore, the truncation of long-ranged interBlopolyrners Vol 31, 77-100 (1991) (C 1991 JcBhn Wilry & Sons, Inc

CCC 00O6-3525/91/010077-24504 00

* Present address: AT&TBell Laboratories, Murray Hill, N J 07974.

actions in such reduced particle models may result in poor estimates of energetic quantities such a s relative, and certainly absolute, free energies of binding. The purpose of the virtual rigid body simulation method presented here is t o maintain the advantage of computational speed, while improving the accuracy and range of application beyond that which is currently possible with methods such as SBMD. In this work we have developed a protocol for generating a virtual potential energy function for proteins and polypeptides, which can be derived from the atomic and/or extended atom parameters of e ~ i s t i n g ~molecular -~ mechanics force fields. These coarse-grained potentials are employed in regions far removed from the spatially localized region of interest, thereby functioning as a classical bath. T h e existing full atomic force fields are used in regions of closer proximity to the active site. Because nonbonded interactions are the crucial component of protein simulations, due to their importance in the description of the distribution function and the computational cost associated with their evaluation, the virtual representation of protein systems is motivated by a reduction of the N atomic centers to ( fewer) virtual centers, and the desire to maintain a description of the original electrostatic anisotropy. Because a residue may contain 7-25 77

78

HEAD-GORDON AND BROOKS

atomic centers, a reduction to 1-5 virtual centers can be a substantial savings in computational expense. In addition, we must also develop a good approximation to the original potential energy surface. The virtual centers will be chosen so that the original full atomic nonbonded interactions are described; it will be shown that a n electrostatic multipole expansion truncated a t the dipole term, and an effective isotropic van der Waals radius about the virtual center, successfully mimic these original nonbonded interactions. A description of the connectivity between virtual centers, and the use of rigid body equations of motion to evolve the point dipoles, provide the correct anisotropic fluctuations of this nonbonded energy and force field. Previous studiesa have also attempted to reduce the complexity of the globular protein potential energy surface by reducing the original number of atomic centers. Our approach is quite different in a t least two important ways-the use of rigid body dynamics to propagate the time development of the system, and the use of a more conservative division in electron density. This paper is organized a s follows: In the next section we provide a qualitative discussion of the reduction of the original atomic centers to virtual centers for the 20 commonly occurring amino acids that make up proteins and polypeptides. In the third through fifth sections, we present the specific form of the potential energy function that will describe the interactions between these newly defined virtual centers. In the third section, the electrostatic multipole expansion model used to describe the chargecharge interactions is described, and is analyzed in terms of computational cost and its ability to reproduce the full atomic electrostatic potential energy surface. In the fourth section we develop the excluded volume and dispersion terms of the virtual potential that approximate the van der Waals envelope of the original centers. We complete the specific form of the virtual potential energy function in the fifth section, where the connectivity potential between virtual centers is described. In the sixth section, we explore the dynamical nature of the model, beginning in terms of the interplay of rigid body theory with the anisotropic potential functions, followed by the specifics on how protein and polypeptide fragments are reduced to rigid mass distributions. Finally in the seventh section, we test the reduced representation model to ascertain whether it provides the flexibility and accuracy necessary for describing local interactions and secondary structural elements of polypeptides or proteins.

THE CREATION OF VIRTUAL CENTERS The large number of degrees of freedom found in macromolecular systems may, in some cases, be rigorously reduced to “virtual degrees of freedom”,’ i.e., a smaller relevant subset that describes the configurational probability distribution function of the original system. These ideas were initially discussed by Floryg; a n example taken from his work will illustrate the basic strategy. Consider the polymer polyethylene, which is comprised of the repeating monomer unit - (CH,),-; the degrees of freedom that would describe a particular conformation include C,-C,+land C,-H, bond lengths, C,-C,+l-Cl+sand H,C,-H, bond angles, in addition to the Cr-CL+l-C1+2Cl+3torsion angles. Due t o the chemical nature of these degrees of freedom, bond lengths, and to a lesser extent, bond angles, may be considered as fixed by the high energy penalty on deformation. Hence the full configuration space distribution function of the polyethylene polymer may, in principle,” be replaced by its dihedral angle space distribution function. T o apply these ideas to globular macromolecules such as proteins however, we must account for the high degree of inhomogeneity of interaction, for which the simple dihedral angle potential description is insufficient; intermolecular interactions are required in addition t o the connectivity description. The polarity of the backbone and over half of the 20 or more different side chains make electrostatic considerations imperative. Furthermore, van der Waals nonbonded interactions impart the largest contribution to the three dimensional “shape” of proteins. Therefore, the potential energy description for use in simulations of proteins must account for the greater complexity of the configuration space distribution function of these biological macromolecule~.~ T h e reduced representation of a polypeptide or protein adopted in this work is motivated by electrostatic considerations based on the spatial division of the biomolecule in terms of its secondary structure. Secondary structures, such as a-helices and 0sheets, result when a hydrogen bond is formed between the polypeptide backbone amide of one residue and the carboxyl group on a different residue. The large net stabilizing electrostatic force due t o the alignment of the individual peptide dipoles is thought to be the reason why proteins so readily incorporate these structural elements. This feature of protein structure motivates a secondary structure electrostatic model in which the amino acid residue

VIRTUAL RIGID BODY DYNAMICS

is divided into a backbone region and a side-chain region. In this section, we concern ourselves with reducing the backbone and side chains of all 20 commonly occurring amino acids into a smaller number of virtual centers, which should provide significant computational savings compared to the full atomic potential. T h e reduction will be performed so that important connectivity degrees of freedom are maintained, i.e., those that are expected to contribute significantly to motion, and that the nonbonded interactions of a virtual center mimic the original centers it replaces. The quality of the electrostatic multipole model considered in the next section, truncated a t second order, is dependent on a localized spatial distribution of the original charges. Describing the van der Waals envelope of the original centers will also require careful consideration when centers are reduced. In this section we will keep in mind the limitations of a one-center description of the original anisotropic groups, but leave the discussion of the specific form of the virtual potential function describing the interaction between virtual centers to later sections. In Figure 1( a ) we schematically show a molecular mechanics representation of the glycine dipeptide, from which the backbone pattern that is characteristic of all of the amino acids except proline is evident. The a-carbon and terminating methyl blocking groups are represented as extended atoms ( CH2and CH:?),in which the aliphatic hydrogens are absorbed into the carbon van der Waals radius. The absorption of aliphatic hydrogens into carbon or sulfur

Cn3

c Y

Vm

VCO

VNH

Figure 1. ( a ) Schematic representation of the extended atom model of the polypeptide backbone. In the extended atom representation: all atomic centers are represented explicitly, except for the aliphatic centers involving carbon and sulfur. In this case, the hydrogens are absorbed into the carbon or sulfur center to which they are attached, i.e., CH, CH2, CH3, and SH are each represented as one effective center. ( b ) Schematic representation of the virtual atom model of the polypeptide backbone. Those centers that are virtual centers are labeled so that the first letter is always V. The subscript denotes the original centers present.

79

centers is a commonly employed molecular mechanics potential representation5 and is the “full atomic” potential we refer to in the remainder of this paper. The peptide groups that lie between the extended atom carbons are chemically similar to formamide; like formamide, the partial double-bond character about the C-N bond allows these moieties t o be considered as rigid. As mentioned above, the peptide group oxygen and hydrogen atoms are involved in the formation and stabilization of secondary structure. Because secondary structural elements are ubiquitous in proteins and are considered to be possible seeds for initiating protein folding, a good electrostatic description of the backbone is imperative. Furthermore, the alignment of dipoles in secondary structures requires the 6 ( Cl-l-Nl-Cci-Ci)and $ ( N,-C,,-C,-N,+I) dihedrals of the backbone to adopt particular values; the a-helix, for example, takes on 6 and $ values of -60” and -60°, respectively. We therefore have chosen the electron density division as follows: the extended a-carbon atom is kept explicitly in order to preserve the original connectivity of the backbone, i.e., 6 and #. For the peptide moiety itself, electrostatic interactions for stabilizing secondary structure dictate that the carbonyl group ( C - 0 ) be redefined as one virtual center (VCO), and the amide ( N - H ) group as one virtual center ( V N H ). This division of the backbone will be shown later to preserve a-helical secondary structures. The extended atom methyl blocking groups are retained in the virtual description. In Figure 1( b ) we depict the virtual representation of the glycine dipeptide. We note here that the virtual atom labels given in the text of this section will be important for defining their associated virtual molecular mechanics parameters in later sections (see Table I). Figures 2-20 show the virtual representation of the remaining 20 commonly occurring amino acid side chains, where each amino acid is depicted in its blocked dipeptide form. T h e role of side chains in the formation and stabilization of secondary structural elements may play a role equal to that of the backbone atoms; however, they may also detract from stability. Proline (Figure 2 0 ) , for example, prevents the formation of the secondary structure hydrogen bond. Adjacent bulky side chains in a n ahelix or a P-sheet may require the structure to adopt nonoptimal +,$ values in order to overcome unfavorable steric interactions between the pair of side chains. Favorable electrostatic interactions between side chains, such as an adjacent lysine and aspartic acid, may result in exceptionally stable secondary structures. As a general rule, side chains will be re-

-

-

80


CH3

Figure 2. A schematic representation of the virtual model of the alanine side chain. See Figure 1( b ) for details.

duced t o a proportionately small number of virtual centers than the backbone due t o the fact that there are significant computational advantages t o including fewer centers. However, the reduction will be based on preserving gross electrostatic and excluded volume interactions, and the connectivity degrees of freedom t h a t contribute to the largest amplitude motion possible for the side chain. It is intuitive that motion involving the dihedral angle, X1 (torsion about the C,-C, bond) ,will result in the largest motion of the side chain as a whole. In cases involving side chains of significant length (three or more nonbranching atoms) the x1 dihedral will be preserved by retaining the 0-carbon explicitly. T h e side chain glycine has no functional group that branches off the extended a-carbon, and therefore its division has already been discussed in terms of the backbone. For the aliphatic side chains alanine, valine, leucine, and isoleucine, electrostatic considerations are unnecessary, and we will only consider the excluded volume and allowed conformation space of the original number of centers. T h e alanine side chain (Figure 2 ) is represented as a n extended P-carbon atom ( C H 3 ) 5t h a t is responsible for limiting the accessible conformation space of this residue compared t o glycine. Therefore to keep it distinct from glycine, we keep the extended P-carbon center of alanine, so that this side chain remains unchanged. For valine (Figure 3 ) , the side chain is not of sufficient length t o keep the extended P-carbon atom explicitly; in fact, if a geometric origin

Figure 4. A schematic representation of the virtual model of the isoleucine side chain. See Figure 1( b ) for

details.

is chosen for the C,H,, CT1H3,and CT2H3extended atoms, one virtual center with an effective excluded volume radius serves a s a good approximation to the original centers. We therefore represent this side chain with one virtual center (VCCC ) . Isoleucine and leucine (Figures 4 and 5 ) , on the other hand, are sufficiently extended so that motion about XI has significant amplitude. We retain the extended P-carbon center for leucine; the remaining atoms of leucine are equivalent to the valine side chain, and hence are reduced to a similarly defined virtual center (VCCC). T h e isoleucine side chain is less well described by the valine / isoleucine procedure due to the asymmetry of the excluded volume envelope of the side-chain atoms when the extended @-carbon atom is kept explicitly. Instead, we can reduce this side chain to two centers by representing the P- and y2-methyl extended atoms with one virtual center (VCC) , and condensing the y1and 6 extended atom carbons into a different virtual center (VCC) . Next we consider amino acid side chains that have electrostatic functional groups. Although serine (Figure 6 ) is not a side chain with large-amplitude motion, we will preserve the extended @-carboncen-

Ct.3 CH3

Figure 3. A schematic representation of the virtual model of the valine side chain. See Figure 1( b ) for details.

Figure 5. A schematic representation of the virtual model of the leucine side chain. See Figure 1( b ) for details.


81

CH2

Figure 6. A schematic representation of the virtual model of the serine side chain. See Figure 1( b ) for details.

ter in order t o retain a good description of the electrostatic dipole of the hydroxyl group, which itself will be reduced to one virtual center (VOH) . Similarly, the hydroxyl group of threonine (Figure 7 ) will be a single virtual center ( V O H ) , but we condense the extended @-carbonand y-carbon atoms to one virtual center due to the aliphatic nature of these centers (VCC) . Aspartic acid (Figure 8) is similar to serine in that the extended @-carbon is represented explicitly and the carboxyl moiety defines one virtual center ( V C 0 2 ) . This provides a good electrostatic representation of the carboxylate moiety. The amide form of aspartic acid, asparagine (Figure 9 ) , is also divided to accommodate a reasonable electrostatic representation of many atoms with one center. Therefore, the extended P-carbon atom remains explicitly, the carbonyl group is reduced to a virtual center (VCOS), and the amide group is also reduced to one center (VNH2 ) . Glutamic acid (Figure l o ) , on the other hand, is of sufficient length that the amplitude of the x1torsion is important; hence we wish t o keep the extended P-carbon center. In order to maintain a good electrostatic description of the carboxylate group, we retain the C, extended atom as well. The acid functional group is represented as one virtual atom, which preserves the electrostatic description of the carboxylate group ( V C 0 2 ) . Similarly, for glutamine

VC33

Figure 8. A schematic representation of the virtual model of the aspartic acid side chain. See Figure 1( b ) for details.

(Figure 11) we maintain the extended @- and y-carbon centers and again reduce the carbonyl and amide groups to (two) distinct virtual centers ( VCOS and VNH2 ) . For the amino acids lysine (Figure 1 2 ) and arginine (Figure 1 3 ) , there is considerably more latitude in the number of centers that can be reduced. In the case of lysine, we reduce the C, and C, extended atoms to one center (VCC) and the Clgand C atoms to a n additional virtual center (VCC). T h e amino group will be reduced to one center in order to preserve the electrostatic nature of this basic residue ( V N H 3 ) . For arginine, the two aliphatic ( y and 6 ) extended atom carbons are reduced to one virtual atom (VCC) . T h e electrostatic interactions of the remaining atoms of the side chain dictate that we preserve the integrity of 6-guanido group. Therefore we have defined the functional groups N,H ( V N H E ) and C,(C), and the two amide N,H2 (VNH2 ) , as virtual centers. Next we consider the side chains that contain sulfur. Cysteine (Figure 1 4 ) is a unique residue that often reacts with other cysteines to form disulfide bonds in proteins. We shall preserve this possible

.. .,......... ..... ...-, .... .... ..,:... . I .

.I

v~tlcl:

Figure 7. A schematic representation of the virtual model of the threonine side chain. See Figure 1( b ) for details.

$L-

Figure 9. A schematic representation of the virtual model of the asparagine side chain. See Figure 1( b ) for details.

82


Figure 10. A schematic representation of the virtual model of the glutamic acid side chain. See Figure 1( b ) for details.

covalent interaction by leaving the cysteine residue unchanged. The moment arm of the methionine side chain (Figure 15) is of sufficient length to warrant retaining the extended @-carbonexplicitly. However, the methionine sulfur forms no disulfide linkages and the side chain is largely hydrophobic. Therefore in the interest of computational savings we choose to replace the C,, S6,and C, extended atoms with one virtual center (VCSC) . The aromatic side chains phenylalanine (Figure 16), tyrosine (Figure 1 7 ) , and tryptophan (Figure 18) can be classified as hydrophobic amino acids with long moment arms; hence for these residues we will define the extended @-carbonatom as one center, and look for reductions in the remaining atoms so that a gross approximation to the original van der Waals envelope is achieved. If the aromatic rings of these side chains were able to spin freely about the C,-C, bond, a very reasonable approxi-

it

I

'

j

Figure 12. A schematic representation of the virtual model of the lysine side chain. See Figure 1( b ) for details.

mation to the excluded volume of this functional group would be to equate the van der Waals radius with the volume radius swept out b the ring? However, the packing forces in proteins are quite strong so that a large concerted motion of the protein matrix surrounding these groups is required in order for the ring flip transition to occur; therefore these motions are rather infrequent, and the aromatic group should be viewed as anisotropic on the molecular dynamics time scale. Any model should capture this anisotropy, at least in a gross way, by defining a larger van der Waals envelope in the plane of the ring than in the direction perpendicular to this plane. We accomplish this by breaking up the

....

it, .... .....',

..... .... ........ .... .......

....../ 4 3 c - 3 +5?:......... ......... .... I

i

I

VC?

G

Figure 11. A schematic representation of the virtual model of the glutamine side chain. See Figure 1( b ) for

details.

..*-

Figure 13. A schematic representation of the virtual model of the arginine side chain. See Figure 1( b ) for details.


C

83

c.3 CHJ

Figure 14. A schematic representation of the virtual model of the cysteine side chain. See Figure 1( b ) for details.

six-membered ring into two virtual centers, one defined by the extended C,, Cal, Ctl plane of atoms (VC3R). and the second defined by the extended Cr2, Cf2,C, plane (VC3R). T h e new plane defined by the two virtual centers and the P-carbon atom will accomplish the desired result of describing a larger excluded volume in the plane of the original aromatic ring. This reduction describes phenylalanine, and tyrosine only requires a n additional virtual center to describe the hydroxyl dipole (VOH ) . For tryptophan we must also consider the reduction of the indole functional group. Electrostatically, we wish to preserve the integrity of the NC1Hgroup ( VNH ) ; hence we treat it as one virtual center; the extended C, and Csl will comprise a different virtual center (VCCR), thereby completing our reduction of the aromatic sidechains. Histidine is another side chain whose reduction to a virtual description would benefit from retaining the C, extended atom explicitly. For the remainder of the side chain, the imidazole group, several protonation states are possible. We feel it is important to maintain these in the virtual description. We therefore use a n explicit description of all unprotonated nitrogens ( N R ) , and reduce to different virtual centers the protonated nitrogens ( V N H ) and the extended C,-Cs2 group (VCCH) . The extended Cr2atom cannot be conveniently reduced with any other atoms; hence it is explicitly retained ( C R l E ) .

Figure 16. A schematic representation of the virtual model of the phenylalanine side chain. See Figure 1( b ) for details.

Finally, we consider the imino acid residue proline. This residue is largely hydrophobic and does not undergo large amplitude motion about x1 due to the fact that the side chain C, bonds covalently to the backbone nitrogen. The importance of this residue is t h a t the rigid constraints of the ring restrict the rotation about the N-C, bond of the backbone, so that the adjacent peptide bond is more likely to adopt a cis conformation." Therefore, we have reduced the extended C, and C, atoms to one virtual center (VCCP) , and have retained extended Cs explicitly in order to maintain the N-Cs bond.

ELECTROSTATIC MULTIPOLE EXPANSION MODELS In this section we present a n electrostatic multipole expansion model that will be shown to provide a good static description of the secondary structural elements of the polypeptide or protein while reducing the number of atomic centers. In the following discussion, the multipole model is assessed for its computational efficiency, i.e., the number of virtual centers and the complexity of the potential energy

ct

vc-2c,-

i

... .......

Mi.3 ......., : . ....

i

Figure 15. A schematic representation of the virtual model of the methionine side chain. See Figure 1( b ) for details.

Figure 17. A schematic representation of the virtual model of the tyrosine side chain. See Figure 1( b ) for details.

84


CHJ CH3

VNC

Figure 18. A schematic representation of the virtual model of the tryptophan side chain. See Figure 1 ( b ) for

Figure 20. A schematic representation of the virtual model of the proline side chain. See Figure 1( b ) for details.

details. The coefficients Q M , Q u , and QQ are, respectively, scalar, vector, and tensor quantities defined a s function, and the accuracy with which it can reproduce the original potential energy hypersurface. T h e electrostatic multipole expansion of the original multisite potential has been derived previand is a n expansion of the form

(where we denote all vector quantities by a circumflex and all matrices by open boldface). The ria = 11 ti - to I( -'is the distance between the test charge qi and the expansion origin, @( r W') is the potential, l?( ?io) is the electric field, and F is the field gradient, due t o a point charge.

c-3

H3

Figure 19. A schematic representation of the virtual model of the histidine side chain. See Figure 1( b ) for details.

QM

=

C ql I

QD =

;= C I q l ( t l - F0)

(3b)

where LY and 0are the x , y , and z components of the atomic charge position vector ijodefined from the expansion origin ?o, and 6 is the Kronecker 6 function. It is evident that the multipole expansion of the electrostatic potential energy function is a more complicated form of the electrostatic potential than the simple Coulomb's law sum. In addition to charge-charge interactions, one must consider the dipole interaction with the field due to a point charge, a s well as the quadrupole interaction with the field gradient due to a point charge. Furthermore, multipole-multipole terms must also be considered. Although we can replace a large number of atomic centers with only a few multipole centers, thereby reducing the N 2 problem, the coefficient of what is now the M 2 loop (where M < N ) may be very large, especially when quadrupole terms are included. We will therefore only consider the multipole expansion case when the series is truncated a t the dipole term. The division of amino acid residues into virtual centers has been outlined in the second section. However, the quality of the fit for the multipole series will depend on where the (arbitrary) expansion origin is placed. We have chosen the virtual origins to lie on atoms that maintain a good description of some of the original connectivity, and/or on centers for which a spherical van der Waals radius is a good approximation to the original excluded volume envelope. For example, we will place the two backbone multipole centers on the carbon and the nitrogen in

85


order to maintain the dihedral angle space of the backbone. For side chains, a good excluded volume description is of greater concern; hence all virtual side-chain origins will be determined from the van der Waals weighted center of geometry of the original group of atoms.

where (T is the CHARMM5 van der Waals radius parameter. In neither case is the origin optimized with respect to the optimal electrostatic description. T h e following test therefore concerns the electrostatic performance of these origin definitions. The crucial test of the virtual center model is how well it can reproduce a static potential energy surface. We take a s our test case ribonuclease A, and perform a typical SBMD 1 ~ type 4 division; all residues with a t least one atom outside a 14-A sphere centered a t the a carbon of His-12 are labeled as “heat-bath region residues.” We consider the full atomic electrostatic potential beyond this region as the control. Each residue in the heat bath is further divided into a backbone region and a side-chain region, which in turn are further divided into the virtual centers discussed in the second section. The input needed to generate either the full or electrostatic multipole model is a protein crystal structure and a set of atom point charges. The crystal coordinates of ribonuclease A are taken from Richards.’* The point charges we have used are taken from the topology file used in the molecular mechanics program CHARMM.5 In Figures 21-23 we compare the full atomic and multipole expansion (truncated a t the dipole term) electrostatic surfaces, which were generated by placing a test charge a t discrete values of 19,$ on the spherical surfaces defined by radial distances of 7.5, 10.5, and 13.5 A, respectively. T h e origin is defined as the position of the a-carbon atom of His-12. The test charge interacts with the electrostatic centers outside the 14-A boundary only. We have chosen test distances that correspond t o typical interaction distances in a molecular mechanics simulation. T h e interaction of the test charge on the 7.5-A sphere with multipole or atomic centers on and beyond the boundary a t 14.0 8, is typically in the range of the “cutoff” criteria used in most simulations. The 13.5A test sphere was chosen as the most sensitive probe of the convergence for the multipole expansion series. As is evident from these figures, the dipole model reproduces the full atomic surface a t the three radial

distances quite well, despite the nonoptimized placement of the multipole origins. This is due to the rather conservative reduction of atomic centers to virtual centers, so that the expansion series converges quickly. The good agreement a t a distance of 7.5 A is due to the rapid convergence of the multipole series. There seems to be negligible degradation as one moves outward radially from the origin (Figure 22), which is especially striking a t the 13.5-A boundary (Figure 23), where effects arising from the limited multipole series should be expected. Overall the dipole model seems to be very adequate for reproducing the static electrostatic field of the

120 60

s

0 -60

-120 -180

-180 -120

-60

0

60

120

180

-60

0

60

120

180

120 60

3.

0

-60

-120

-180 -180 -120

9 Figure 2 1. The full atomic and electrostatic multipole expansion (monopole, dipole) potential energy 4, surface a t 7.5 A. A test charge of 1 eu interacts with the above representations outside of 14 A. Each contour represents 0.015 kcallmole. ( a ) Full atomic and ( b ) electrostatic multipole expansion.

+

86


180

These range from a choice of dielectric responses, such as a constant dielectric

120

60

-+

where c is the dielectric constant, or a distance-dependent dielectric

0 -60

-120 -180 -180

-120

-60

0

60

120

180

@

which takes into account electrostatic screening. This is not a rigorous treatment of screening, because such rigor for microscopic descriptions is especially difficult (and expensive) €or inhomogeneous 180

120 120

60 60

0

0 -60 -60 -120 -120 -180 -180 -180

180 -120 -120

-60

0

60

120

-60

0

60

120

180

60

120

180

180

9 9

180

Figure 22. The full atomic and electrostatic multipole expansion (monopole, dipole) potential energy I$,+ surface at 10.5 A. A test charge of 1 eu interacts with the above representations outside of 14 A. Each contour represents 0.033 kcal/mole. ( a ) Full atomic and ( b ) electrostatic multipole expansion.

120

60

0

full atomic surface. We will return t o further testing of the local nonbonded field, by examining the potential as a function of conformation in the seventh section. Finally, we discuss the implementation of the electrostatic multipole expansion model in the molecular mechanics program CHARMM.5The actual atomic parameters (charges, for example) used in the definition of the multipole coefficients are simply those of the CHARMM force field.5 In addition, one often has a choice of several nonbonded options in a molecular mechanics package when performing energy evaluations and/or dynamic simulations.

-60

- 120

'

-180 -180

I

-120

-60

0

9 Figure 23. The full atomic and electrostatic multipole expansion (monopole, dipole) potential energy 4, surface a t 13.5 A. A test charge of 1 eu interacts with the above representations outside of 14 A. Each contour represents 0.15 kcal/mole. ( a ) Full atomic and ( b ) electrostatic multipole expansion.

+


biological systems. Indeed, the functional form of the screening in Eq. ( 6 ) is ad hoc, and simply performs as another "parameter" in the empirical force field (for describing DNA interactions in particul a r ) . We therefore limit the implementation of the multipole model to a constant dielectric form only. A practice usually employed by those who run simulations is t o truncate the number of N2 interactions by defining a cutoff radius around a given atomic center and evaluating only those interactions found within the cutoff. The main problem with such a technique is that energy and derivatives are discontinuous a t the boundary, resulting in unstable minimizations and dynamics, particularly when highly mobile solvent molecules are included.'' A number of correcting functions have been applied to the energy and its derivatives to address this problem. One such method is known as the atom-atom shifting potential,'

87

action. In order to quantify the savings in this case, we employ ribonuclease A as a n illustrative example. We divide all residues into a backbone and sidechain region, and further reduce these regions so that the backbone is comprised of two multipole centers and the side chains are represented by two multipole centers. We then evaluate the cost of the N2 loop as follows: For a given residue i, we count the number of dipole-dipole interactions within 5 A; outside this region we consider the multipole interaction series to have converged a t the first term, and count the number of monopole-monopole interactions. We then repeat the process for residue i 1, and so on for all M residues. For the division of ribonuclease A, a given residue sees an average of K = 10 other residues within 5 A,and L = 114 residues outside this sphere (note that K L = M ) . Residue-residue interactions within 5 A (where each residue has four dipole centers) require 16 X 51 = 816 multiplies, where 51 multiplies are required to evaluate the energy and force interactions between a pair of dipole centers. Coulomb energy and forces require 17 multiplies, so that the total cost of an electrostatic evaluation for ribonuclease is 8 1 6 ( K ) 2 17(L)'. This should be compared to 1 7 ( N 2 ) for the conventional method, where N = 1192 for the ribonuclease example. Although this is a n artificial case, the computational savings of the multipole method implemented in this fashion is over 75 times that of the full atomic case, and offers the attractive advantage that the savings will scale with increasing system size. For applications that genuinely require the incorporation of all longrange interactions (such as absolute free energies, or large systems such as lipid bilayers) , the virtual method should find broad application.

+

+

+

which is multiplicatively applied to the energy and derivatives

The shifting function, in conjunction with the constant dielectric, multipole expansion model, has been incorporated into the molecular mechanics program CHARMM.5 It is interesting to note that for the shifted potential, the long-range energy and forces go to zero smoothly a t the cutoff radius, but that the short-range electrostatic field is modified as well. The result of this (unphysical) modification is that probes such as rms deviations (structure) and rms fluctuations (dynamics) can be significantly larger than those found with no cutoffs. One possible solution for smoothing the nonbonded interactions with these truncation schemes is to reparameterize the empirical force fields to compensate for the modified structure and dynamics." An advantage of the multipole expansion used in the virtual potential function, however, is that the interaction between multipole centers will decay as r ,I' at large enough distances, in which case one only needs to evaluate the monopole-monopole inter-

EXCLUDED VOLUME AND DISPERSION POTENTIALS All virtual atoms, or multipole expansion centers, must have a van der Waals potential energy contribution that mimics the excluded volume effect of the original set of atoms it replaces. The van der Waals function implemented in CHARMM5 has the Lennard-Jones form, which describes the repulsive wall and attractive dispersion well as

88


where r y n is the van der Waals radius parameter for the interaction of centers i and j , ti] is the well depth parameter for interactions i and j , and ri, is the distance between the centers i and j . T h e energetic well depth mixing rule is multiplicative ti; =

[

112

and the van der Waals radius of a pair is determined by a n additive mixing rule

and

The simplicity and computational advantages of this isotropic potential function will be exploited by the virtual potential function as well. The objective of this section is to define an effective van der Waals radius and well depth for each virtual center i, which can then be used with the mixing rules defined in Eqs. ( 11)and ( 1 2 ) for atomic-virtual and virtualvirtual interactions. We choose to redefine the original van der Waals radii and well depth parameters for virtual centers, and simply use the CHARMM parameters for atoms retained explicitly (such as the a- and &carbon). For the virtual van der Waals parameterization, we view the effective radius as being the more important parameter since it manifests itself in a power law (6-12), whereas the well depth is only a scaling factor. Therefore for both the backbone and side-chain virtual centers, we define the virtual well depth to be a sum over the number of centers:

In the following discussion we consider the virtual backbone and side chains separately since their virtual origins are defined differently. As discussed previously, the electrostatic origins of the backbone virtual centers were chosen t o maintain the full dihedral angle distribution accessible to the main chain, and not necessarily for a good description of the excluded volume. However, the conservative reduction of these centers, and the nature of the peptide group itself, permit these centers to be reasonably well described by a n isotropic envelope. In general, we will define the backbone virtual van der Waals radii as

“where bil is the bond length between atomic centers i and j.” If we consider the amide group, the CHARMM van der Waals radii of nitrogen and hydrogen are 1.6 and 0.8 A,respectively, and the N-H bond length is 1.0 8.As discussed previously, the virtual center origin is defined as the original nitrogen center. T h e excluded volume in the plane of this functional group is 1.8 8, in the direction of the bond, and 1.6 A in the diametrically opposed direction. Perpendicular to the plane of this group, and perpendicular to the bond in the plane, there is a n excluded volume of 1.6 A from nitrogen, which diminishes t o 0.8 A about hydrogen. An effective radius of 1.7 A should represent the excluded volume about the nitrogen end extremely well. Such a radius definition in the plane of the group, where the direction of approach is end on toward hydrogen, should also be a n excellent approximation. This is a n especially important direction for linear hydrogen bonds, such as those found for secondary structural elements. In order t o maintain the correct excluded volume in this direction, we will define the effective radius of the amide group as 1.7 A,a t the price of a larger excluded volume about the remaining surface of the hydrogen atom. For the carbonyl group, the CHARMM5 van der Waals radii of carbon and oxygen are 2.2 and 1.6 A, respectively, and the C - 0 bond length is 1.2 A. As discussed, we have chosen the origin to be the original carbon atom. Therefore the excluded volume from the carbon origin in the plane of this linear group will be 2.8 A in the direction of oxygen and 2.2 A in the opposite direction. Perpendicular to the plane, we have a n excluded volume of 2.2 A above the carbon center that diminishes to 1.6 A above the oxygen center. An effective radius of 2.5 8, in the plane of the carbonyl group should be a reasonable approximation to the interaction of the oxygen end with a hydrogen for forming a hydrogen bond, and for the volume about carbon, although the “quality of fit” is not quite as good as the amide description. Again, we will have a larger excluded volume in all other directions about oxygen. It should be noted that the CHARMM5molecular mechanics potential has separate 1-4 van der Waals interactions (i.e., pairs of atoms involved together in a connectivity element) for all carbons, in which case the carbonyl carbon has a van der Waals radius of 1.9 A.No important hydrogen bonds are formed by 1-4 interactions; hence we define the effective vander Waalsradiusfor 1-4interactionstobesmaller, 2.35 A.


The virtual center origins of side-chain functional groups were chosen to optimize an isotropic excluded volume representation of the original centers. The definition of the virtual origin for side chains (third section) is a van der Waals weighted average of the original centers. Given this definition of the virtual origin, we may now discuss the effective van der Waals radius for all of the side-chain functional groups. For virtual centers representing aliphatic groups, we first define the average of the maximum and minimum radius in the plane of the original set of atoms,

as well as the average of the maximum and minimum radii perpendicular to this plane

i(J-) = f timax(1)

+ ;min(J-)I

(17)

where the maximum and minimum radii are with respect to the origin defined in Eq. ( 4 ) .T o illustrate this procedure, we consider the effective radius for the virtual center VCCC, which represents the excluded volume of the original group CH3-CH-CH3 (see Figure 3 ) . Given the origin defined in Eq. ( 4 ) , typical bond lengths and angles for this group, and the van der Waals radii of 2.35 and 2.15 A for CH3 and CH, respectively, we determine a maximum extent radius of 3.57 A in the plane of this group, and 2.37 A perpendicular to this plane. Given the aliphatic nature of this group, we therefore simply average these two effective radii directions to give a value of 2.97 A. For 1-4 interactions the CHARMM van der Waals parameters for these same atoms are 1.9 A, so that the maximum radius in the plane and perpendicular to the plane is 3.30 and 1.9 A, respectively; the effective radius of the VCCC center for 1-4 interactions is therefore 2.60 A. When the virtual center represents an original fragment with significant electrostatic interactions, we would like to optimize the excluded volume for the direction of approach that is favored electrostatically. We do this so that strong electrostatic interactions do not dominate, or even become catastrophic, due to a van der Waals radius that is too small. Since all of the electrostatic virtual centers involve hydrogen and oxygen, both of which are important for hydrogen-bond formation, we define the effective van der Waals radius as the distance between the origin defined in Eq. ( 4 ) and the excluded volume envelope of these atomic centers. We again illustrate these principles with an example. For the virtual center VOH, which represents all hydroxyl

89

groups, the direction of approach of interest is along the 0 - H bond coming in toward the hydrogen, an optimal approach for hydrogen-bond formation as well. Because the oxygen van der Waals radius is the same as that of nitrogen, the effective radius of 1.7 A found for the amide is also optimal for the hydroxyl; in fact, this effective radius is exact, given the origin defined in Eq. (4),for this line of approach. At this point we have laid out the protocol for reparameterizing the original nonbonded full atomic force field in terms of virtual centers. In Table 11, we have listed the VRBD van der Waal parameters. In the final section we test how well the isotropic excluded volume description of the original anisotropic envelope reproduces local interactions, where the van der Waals parameters play a crucial role. Before proceeding with these tests, however, we require a description of the interatom connectivity between virtual atoms.

VIRTUAL CONNECTIVITY MODELS The next goal is to develop potential terms that are functions of the virtual degrees of freedom in order to mimic the covalent interactions of the full system. A reduced representation of protein systems requires a description of the connectivity involving virtualvirtual and virtual-atomic centers. In the previous sections we have developed a reduced representation of the nonbonded interactions that provides significant computational savings. In this section we define the connectivity necessary to describe the re-

Table I

Virtual Atom Types

Virtual Atom Type

VCO VNHl VNH2 VNH3 VNHE

vcos vcsc vccc VC3R vc02 VOH

vcc

VCCP VCCR VCCH

Original Atomic Fragments Backbone carbonyl (CO) Backbone amide (NH1) Side chain amide (NH2) Side chain amino (NH3) Imidazole amide (NH) Side chain carbonyl (CO) Methionine (CH2-S-CH3) Propane (CH2-CH2-CH3) Ring (CH-CH-CH) Side chain carboxyl ( 0 - C - 0 ) Side chain hydroxyl ( 0 - H ) Ethane (CH2-CH3) Ethyl (CH2-CH2) Ring (CH-CH) Ring (CH2-CH)

90


mainder of the potential surface on which we may then perform dynamics. We first consider the backbone virtual degrees of freedom, since the virtual origins have been chosen so that their connectivity maps directly onto to the original atomic degrees of freedom. We then describe the procedure for determining the connectivity of the 20 amino acid side chains, where a direct mapping procedure is not generally possible. The representation of the backbone degrees of freedom in terms of two virtual centers will give rise to the virtual connectivity degrees of freedom displayed in Figure 1 ( b ) , where Ca, VCO, and VNH, denote an a-carbon, virtual backbone carbonyl group, and virtual amide group, respectively. The benefit of picking virtual origins that preserve the original connectivity is that the virtual bonds, angles, and dihedrals completely map onto the fullsystem connectivity. The virtual connectivity potential will be of the form bonds

V

=

2

angles

kbi(bi

- bio)'

+ C

k,i(Bi

-

Bio)'

i

1

impropers

+ 2

kui(ai- Wio)

2

i torsions

+ C

cosine Fourier series

($i)

( 18)

i

which is the same form used in most molecular mechanics potentials. Therefore, the harmonic bond and angle force constants and equilibrium geometries of the backbone connectivity degrees of freedom defined by CHARMM5 will also be used for the virtual potential function. In general, the virtual side-chain degrees of freedom do not map directly onto the original full atomic center connectivity due to the more drastic reduction in the number of side-chain atoms. However, all sidechain virtual center origins have been defined to optimize the representation of the excluded volume of the original set of atoms. Therefore equilibrium geometries for the potential function defined in Eq. (18) are easily defined. In Figures 2-20, we have presented the geometry of all the virtual side-chain representations. Bond lengths are simply determined by evaluating distances between virtual centers using the equilibrium geometry of the full atomic system and Eq. ( 4 ) . Defining force constants for the connectivity degrees of freedom, on the other hand, is less straightforward. In some cases, the virtual center origin is very close spatially to the original center involved in the bond or angle definition. For example, the force constant for the bond between

the @-carbonand the oxygen of serine can be used for the bond between the @-carbonand virtual hydroxyl center as well, since the virtual origin lies very close to the original oxygen center. For other virtual connectivity degrees of freedom, the mapping is much less obvious, especially when original soft torsions become hard degrees of freedom; if we consider lysine, for example, the C,-C,-C, -C, dihedral is now a virtual bond C,-VCCC. Because most of the CHARMM bond force constants lie between 400-600 kcal mole-' k', and angle force constants between 20-85 kcal mole-' rad-', we choose the lower bound for virtual hard degrees of freedom when torsions were present originally. The force constants for dihedrals are very close to zero in general, except for cases involving torsions about bonds with significant double-bond character, namely carbon-nitrogen bonds. For virtual dihedrals in which the torsion is about a bond where one center has incorporated a nitrogen and the other center a carbon, we use this larger force constant. In Tables IIIVI, we give the VRBD parameters (equilibrium geometries and force constants) for all virtual connectivity degrees of freedom. Although often the choice of force constant is not rigorously defined, we will show in the next section that the interaction of the virtual potential function with the rigid body equations of motion make the specific values of the force constant unimportant. It is appropriate to note at this point that certain original degrees of freedom have been lost on reduction to virtual centers, such as the bonds and angles involving the original carbonyl and amide group of the backbone. Evolving the electrostatic dipoles of the virtual backbone and side-chain groups in a physically correct way requires that we do not permit the dipole orientation to violate the original geometric constraints. In the following section, we briefly discuss the recovery of these degrees of freedom by defining proteins and polypeptides to be composed of rigid body fragments.

THE VIRTUAL POTENTIAL FUNCTION A N D RIGID BODY M O T I O N To complete a virtual description, we must consider the dynamical orientation of the anisotropic electrostatic interactions. The correct dynamical evolution of the electrostatic origin and dipole orientation in the reduced system depends on the mass, connectivity, and the nonbonded forces of the original system. For example, due to the significant torsional barrier about the C-N bond of the peptide


91

Table I1 Virtual Rigid Body Dynamics van der Waals Parameters"

Atom Type

vco VNHl VNH2 VNH3 VNHE

vcos vcsc vccc VCSR vc02 VOH

vcc

VCCP VCCR VCCH ~

2.5000 1.7000 1.7000 1.7000 1.7000 2.3000 2.6370 2.7260 2.5180 2.7000 1.5000 2.6450 2.6150 2.6450 2.6450

-0.1396 -0.1441 -0.1127 -0.09 70 -0.1142 -0.1396 -0.1128 -0.1396 -0.1200 -0.4713 -0.1045 -0.1142 -0.1142 -0.1142 -0.1142

-0.1000

2.3500

-0.1000 -0.1000 -0.1000 -0.1000 -0.1000 -0.1000 -0.1000 -0.1000 -0.1000

2.0600 2.4550 2.3610 2.3200 2.1100 2.2800 2.2800 2.2800 2.2800

~~~~~

a

Well depths

(0in units of kcal/mole and van der Waals radii

unit, an electrostatic multipole expansion about the carbon and nitrogen centers of the peptide link should always result in a dipole that lies in the plane of the original peptide group; the small displacements of the remaining hard degrees of freedom (bond lengths and angles) ensure that there is very little torque on the dipoles within the plane. The

(rmIn/2)in units of

A.

orientational effects of the softer electrostatic interactions of the dipole with the surrounding field will be insignificant compared to the connectivity forces in describing the orientation of the dipole. The question to be answered is how to describe these original hard degrees of freedom after the reduction to vii%ualcenters.

Table I11 Virtual Rigid Body Dynamics Bond Parameters* Bond Type VCO-VNH VCO-N VCO-CH1E VCO-CH2E VCO-CHSE VNH-CH1E VNH-CH2E VNH-CH3E VNH-C VNH-VCCH VNH-CR1E VNH-VC3R VNH-VCCR VNH2-VCOS VNH2-C VNH3-CH1E VCC-CH1E VCC-CH2E VCC-VNHE a

Bond Type 1.33 1.33 1.52 1.52 1.52 1.45 1.45 1.49 1.33 1.73 1.80 3.72 2.45 1.95 1.50 1.45 1.92 1.92 2.10

Equilihrium bond lengths (b,) are in units of

471.0 471.0 405.0 405.0 405.0 422.0 422.0 422.0 471.0 300.0 300.0 300.0 300.0 471.0 471.0 422.0 300.0

300.0 300.0

VC3R-VCCR VC3R-CH2E VC3R-VC3R VC3R-VOH VCSC-CH2E VCCH-NR VCCH-CH2E VCCP-CH1E VCCP-CH2E VNHE-C VCCC-CHlE VCCC-CH2E VCO2-CH1E VC02-CH2E VOH-CH2E VCCR-CH1E VCOS-CH2E VCC-VOH

vcc-vcc

8,and force constants (kb) in units of kcal/(mole-AZ).

3.00 2.42 1.81 2.29 2.75 2.30 1.87 1.90 1.75 1.84 1.94 1.94 1.52 1.87 1.55 1.87 1.85 1.89 2.14

300.0 300.0 300.0 300.0 300.0 300.0 300.0 300.0 300.0 471.0 300.0 300.0 405.0 300.0 300.0 300.0 300.0

300.0 300.0

92


Table IV Virtual Rigid Body Dynamics Angle Parameters" Angle Type

00

ks

Angle Type

00

ks

VCO-VNH-CH1E VCO-VNH-CH2E VCO-N-CH1E VCO-CHlE-VNH VCO-CH2E-VNH VCO-CHlE-VNH3 VCO-CHZE-VNH3 VCO-CH2E-NH VCO-CHlE-N VCO-CHlE-NH VCO-VNH-CH3E VCO-CHlE-CH2E VCO-CHlE-CH1E VCO-CHlE-CH3E VCO-N-CH2E VCO-CHlE-VCCP VCO-CHlE-VCC VCO-CHlE-VCCC N-VCO-CH1E N-VCO-CH2E N-CHZE-VCCP N-CHlE-VCCP VNH3-CHlE-CH1E VNH3-CHlE-CH2E VNH3-CCCL-CH2E VNH2-VCOS-CH2E VNH2-C-VNHE VNHE-VCC-CH2E VCCC-CHZE-CH1E VCCC-CHlE-NH VCCC-CHlE-C VC3R-VC3R-CHlE VC3R-VC3R-CHZE VC3R-VC3R-VCCR VC3R-VCCR-CH2E VC3R-VNH-VCCR

120.0 120.0 120.0 111.6 111.6 111.6 111.6 111.6 111.6 111.6 120.0 109.5 109.5 106.5 120.0 120.9 99.3 115.9 117.5 117.5 87.2 81.6 110.0 110.0 153.0 91.6 98.9 160.0 123.0 114.5 115.9 111.6 108.7 62.0 128.5 46.0

77.5 77.5 77.5 45.0 70.0 45.0 70.0 70.0 45.0 45.0 77.5 70.0 70.0 70.0 80.0 20.0 45.0 45.0 20.0 20.0 20.0 20.0 65.0 65.0 45.0 45.0 5.0 45.0 45.0 45.0 45.0 45.0 45.0 45.0 45.0 45.0

VCO2-CHlE-VNH VCO2-CHZE-VNH VC02-CHlE-NH VC02-CH2E-NH VCO2-CHlE-CH2E VCO2-CHZE-CH1E VC02-CH2E-CH2E VCO2-CHlE-VCCC VNH-CHlE-C VNH-CHZE-C VNH-VCO-CH1E VNH-VCO-CH2E VNH-CHLE-CH3E VNH-CHlE-CH1E VNH-CHlE-CH2E VNH-VCCH-NR VNH-CHlE-VCCC VNH-CHlE-VCC VNH-VCCR-CH2E VNH-VC3R-VC3R VNH-VCCR-VC3R VOH-CHZE-CH1E VOH-VC3R-VC3R VOH-VCC-CH1E VCOS-CHZE-CH1E VCOS-CHZE-CH2E VNH2-C-VNH2 VCC-CHZE-CH1E VCC-VNHE-C VCC-VCC-CH1E VCC-CHlE-C

111.6 111.6 111.6 111.6 112.5 112.5 113.0 120.0 111.6 111.6 117.5 117.5 108.5 110.0 110.0 104.3 114.5 96.9 175.0 40.5 56.7 121.8 121.9 86.2 106.0 127.0 120.0 134.4 100.0 103.7 99.3

45.0 70.0 45.0 70.0 45.0 45.0 45.0 45.0 45.0 70.0 20.0 20.0 65.0 65.0 65.0 45.0 45.0 45.0 45.0 45.0 45.0 50.0 45.0 45.0 45.0 45.0 45.0 20.0 45.0 45.0 45.0

VCCH-CHZE-CH1E VCCH-NR-CR1E VCCH-VNH-CR1E CHZE-VCCH-NR

111.7 60.9 75.3 73.1

45.0 45.0 45.0 45.0

Equilibrium angles (0,) are in units of degrees, and force constants (he)in units of kcal/(mole d e g ) .

For the case of the backbone, we can define the atoms C , , C , 0, N, and H as a rigid body, and use rigid body equations of motion to represent the torque on the nonbonded dipoles. Therefore these atoms will comprise one rigid body and the carbonyl and amide dipole vectors will lie in the plane of the peptide group. The definition of a rigid body dictates that the distances between these atoms (and the virtual atoms that replace them) remain fixed in time, with the result that the geometric constraints of the dipole orientation will be satisfied. The appropriate rigid body definition for side chains will

also allow their electrostatic dipoles to evolve in a meaningful way by changes in the orientation of the rigid body. In addition to the advantages of evolving electrostatic dipoles in a physically correct way, we may also incorporate information concerning the original mass distribution, instead of defining virtual centers as point masses. Because the translational Verlet equation of motion for evolving atomic positions is essentially a statement of Newton's law, the motion of the original atomic centers depends on their masses. We can represent this dependence on mass


Table V Virtual Rigid Body Dynamics Torsion Parameters" Torsion Type

40

k,

Phase

CHlE-VCO-N-CH1E CH2E-VCO-N-CH1E X-VCO-VNH-X X-VCO-NH-X X-VCO-N-X X-C-VNH-X X-CHlE-VNH-X X-CH2E-VNH-X X-VCO-CHlE-X X-VCO-CH2E-X X-VCOS-CH2E-X X-VCC-CH2E-X X-VCC-VNHE-X X-VNHE-C-X X-CH2E-VCCH-X X-CHlE-VCC-X X-CHlE-VCCP-X X-VCCP-CH2E-X X-CH2E-CH2E-VCCC X-CH2E-VCCR-X X-VC3R-VC3R-X X-VC3R-CH2E-X

180.0 180.0 180.0 180.0 180.0 180.0 0.0 0.0 0.0 0.0 0.0 0.0 180.0 180.0 180.0 180.0 0.0 0.0 180.0 180.0 180.0 180.0

10.0 10.0 8.2 8.2 8.2 8.2 0.3 0.3 0.0 0.0 0.0 0.0 0.3 8.2 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3

2 2 2

2 2 2

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

93

about the carbon-nitrogen bond. A rigid body implementation involving a full atomic representation of the backbone would therefore require the reduction of the atoms C, 0, N, and H to a constant mass distribution, with four atomic force sites per rigid peptide body. For the virtual potential function description we may still represent the original mass distribution, but the peptide rigid body will contain only two virtual force sites: the VCO and VNH centers. T h e VCO and VNH centers are the sites of the electrostatic multipole expansion of the original carbonyl and amide groups; hence these dipoles will be defined as rigidly embedded in the plane of the peptide group, with constant orientations relative to each other (figure 24). However, we must ensure that there are no unphysical deformations of the rigid body dipoles. For example, forces that act on the body may give rise to a torque or a couple that would violate the integrity of the original hard degrees of freedom. Improper torsions are CHARMM5 potential terms that maintain, in the case of the backbone, planarity of the peptide and a-carbon atoms. The form of the potential is harmonic,

Equilibrium torsions (&) in units of degrees, and force constants ( k , ) in units of kcal/mole * deg2).

by incorporating the original mass distribution through the inertia tensor in the rigid body equations of motion. Therefore, for functional groups that are now represented by virtual centers, we define a rigid body derived from the original point masses of the system, with the virtual center as a potential interaction site in the body. The specifics of rigid body theory and the discretized algorithm used is relegated to Appendix I. Given a means for evolving rigid protein fragments, we must specify the actual mass division of the backbone and the 20 amino acid side chains. We would like to fulfill two criteria when reducing portions of a protein or polypeptide to rigid bodies. First, we wish to retain a good approximation of the original flexibility, when possible. Second, we will define groups as rigid in order t o approximate the original electrostatic response of the system, i.e., attaching electrostatic dipoles to rigid bodies so that the timedependent orientation of the dipole is approximately correct. As discussed, the peptide group of the protein or polypeptide backbone may be considered rigid to a good approximation due to the large torsional barrier

where the angle I9 defined for the carbonyl improper is determined by the inclination of the plane C,,-C,N,+l t o the plane CaI-Ol-Nl+l,where d,, = 0 ( t h e corresponding planes for the amide group are those defined by the atoms C,-N,+l-C,,+l and CI-HI+lC,,,) . The significant force constant for the C-N proper torsion along with the backbone improper torsions provide the molecular mechanics5 description of the backbone rigidity. The creation of virtual centers results in a loss of the 0 and H atoms, which are necessary for the improper torsions and therefore a description of the rigidity; a loss of rigidity may result in unphysical electrostatic interactions of the now-present dipoles. In order t o recover these missing degrees of freedom we include the a-carbon, in addition to the VCO and VNH virtual atoms of the peptide group, as one rigid body. This is a reasonable approximation to the original system since the atoms Cnr-CIOl-N,+lH,+l-C,,+l always define a plane due t o the diminished torsion about the C,-N,+, bond. We still retain a n approximate description of the two significant flexible degrees of freedom, which are the torsion around 4, $ (where we now refer to the conventional names applied to the dihedral angles C,-l -N,-C,,-C, and N,-C,,-C,-N,+l). The approximation of peptide rigidity modifies the energy

94


and forces due to these potential energy terms in a simulation, since they are a function of the atoms now comprising the rigid body. For example, the COiC,-Ni+lbond angle is determined from a dot product of the C,i-Ci and Ci-N,+lbond vectors. Because the C-N bond is in reality deformable, the values of 0 sampled in the rigid body approximation will be different from a nonrigid calculation; however, the already large force constants of the nonrigid harmonic bond and, t o a lesser extent, angle degrees of freedom, will ensure that a n approximate distribution of the original 0’s will be reasonably represented in the rigid body approximation. Although the correct distribution of internal degrees of freedom is not sampled when rigid constraints are present in a simulation, lo we emphasize that the constraints are present only in the “bath” region, while the “reaction region” of atomic dynamics experiences no constraints other than perhaps hydrogen-bond constraints. Similar considerations of flexibility and electrostatic response will be required for side-chain virtual atoms as well; again, the original mass distribution of each side-chain fragment (since some side chains will have several rigid bodies) will be maintained. In some cases the 0-carbon has been retained explicitly in the virtual potential in order to describe the original large-amplitude torsion about the CcyC, bond. Therefore in general the 0-carbons will be treated as a n atomic site, and not as part of a mass distribution involving other side-chain atoms unless specified otherwise. The side-chain atoms of cystiene and alanine will be treated as point masses. Valine, isoleucine, leucine, methionine, proline, and phenylalanine are side chains whose electrostatic centers are completely isotropic. Therefore we use a reduction that incorporates a s much of the original flexibility as possible. The @-carbonsof alanine, isoleucine, methionine, and phenylalanine are retained explicitly in the virtual description. For valine, the entire side chain will be treated as one rigid body, but with only one virtual atom force site (see the second section); the remainder of the leucine side-chain atoms will be reduced in the same way as valine. The isoleucine side chain will be composed of two rigid bodies, one containing the original 0and yz-carbons, the other the y1and 6 extended carbon atoms; each body is composed of one virtual force site, due t o the reduction of sites discussed in the second section. T h e remaining side-chain atoms of phenylalanine can be considered rigid to a very good approximation; hence the entire ring is reduced t o one rigid body, with two virtual sites. T h e me-

thionine plane of atoms C,-S6-C, are now represented by a mass distribution with one atomic force site. Proline retains the &carbon explicitly, and the remaining ring 0-and y-carbons are treated as a linear body with one potential site. For side chains such as tryptophan and histidine, the mass reduction is unambiguous, since all torsions are about bonds with double-bond character. Hence we reduce each of these side chains to a 0-carbon and one rigid body with four force sites. The remaining side chains contain one to several anisotropic electrostatic sites. In addition to permitting flexibility, we must also maintain the integrity of the original electrostatic response by associating dipoles with rigid bodies that permit them to evolve realistically. For example, rotation around the C,-O, bond for tyrosine is certainly very permissible. However, the reduction of the hydroxyl group to an electrostatic dipole results in a loss of the C,-0,-H, angle; therefore treating the hydroxyl group as a rigid body may result in motion that would violate this original hard degree of freedom. Therefore we include both the hydroxyl group and the ring as one rigid body in order to retain the good electrostatic description, with a loss of some flexibility. For similar reasons involving the hydroxyl group, we reduce the entire side chains of serine and threonine to one rigid body. The side chains of asparagine and aspartic acid are very similar to the main chain to which they are attached; improper torsions impose planarity on the &carbon, carbonyl, and amide groups for asparagine, and for the 0-carbon and carboxyl group for the aspartic acid side chain as well. We therefore treat the entirety of these side chains as rigid. The only difference between the related side chains glutamine and glutamic acid is that the explicit 0-carbon has been replaced by a virtual 0-carbon; we therefore also treat the side chains of these amino acids as rigid bodies. For the arginine side chain, the guano group itself is rigid; however, the t-amide group requires a n improper torsion to make it coplanar with the surrounding C6 and C, extended atoms. We treat the guano and C6 as one rigid body, and reduce the pand y-carbon t o a second rigid body. Finally, the large reduction of atomic sites to three virtual sites for lysine has reduced the conformational space of this amino acid to one dihedral (C,-C,-VCCCVNH3) and one angle ((2,-VCCC-VNHS). We therefore treat each virtual site as one rigid body since the amino dipole will evolve reasonably due to the hard angle degree of freedom.

95


TESTS OF THE VIRTUAL RIGID BODY METHOD In this section we wish t o provide evidence that the virtual potential function is performing adequately, before attempting large scale calculations with the virtual rigid body dynamics ( VRBD ) method. We have shown in the third section that a t both short and long ranges the electrostatic multipole expansion reproduces the full atomic Coulomb’s potential very well. In this section we will address three points concerning the performance of the entire virtual potential function and its interaction with the associated rigid body definitions by performing additional tests. The first test concerns understanding the effect of the rigid body approximation on allowed peptide conformations. Once one separates the rigid body approximation from the approximations inherent in the virtual potential function, we are in a position to see how well local interactions are reproduced by the virtual potential function. This is a n especially stringent test since a poor performance relative to long-range behavior is expected due t o the loss of degrees of freedom and the isotropic van der Waals description. Finally, we provide evidence that the VRBD method is capable of reproducing secondary structural elements in model systems. In Figures 25-27 we present a fully relaxed $,I) map for the simplest dipeptide without a side chain, glycine dipeptide, for the CHARMM potential function (Figure 2 5 ) , the CHARMM potential function with backbone rigid bodies as defined in the sixth section composed of five atomic force sites (C, C, 0, N, H; Figure 26), and the virtual potential function with three force sites (C, VCO, V N H ) per rigid body (Figure 2 1 ) . Points on this surface correspond to restraining the molecule a t particular $,I) values, and allowing all other degrees of freedom to relax; this surface should be largely representative

150 100

50 0

3

-50 100 150

150-100 -50

0

50

100

150

@ Figure 25. A fully relaxed CHARMM 4, 4 surface for glycine dipeptide. Dashed contours denote low-energy regions ranging from 0.0 to 5.0 kcallmole, each contour corresponding to 0.5 kcallmole. Solid contours range from 5.5 to 23.5 kcal/mole, with each contour representing 1 kcal/mole.

of the full conformation space for this backbone fragment. As is evident from Figures 25 and 26, the treatment of the a-carbon and peptide group as a rigid

150 100

50 3

0

-

50

~100 -

150 -150-100

-50

0

50

100

150

@

Figure 24. The orientation of the carbonyl and amide electrostatic dipoles at a conformational minimum. See Figure 1( b ) for further details.

Figure 26. A fully relaxed CHARMMlrigid body 4, $ surface for glycine dipeptide. Dashed contours denote lowenergy regions ranging from 0.0 to 5.0 kcallmole, each contour corresponding to 0.5 kcallmole. Solid contours range from 5.5 to 26.5 kcal/mole, with each contour representing 1 kcallmole.

96


Table VII The 4 and $ Values for a Heptapeptide a-Helix

150

CHARMM/ Powell

100

Creighton" 50

41 =

-47 -47 4 3 = -47 44=-47 4 5 = -47 42 =

0

3

-50 -

-57 -57 $3 = -57 $4= -57 $5 = -57

-75

-71 -71

$1

-71 -73

-87 -85 -85 -85 -86

-57 -54

-56

-52

-51 -47

-54 -58

-50 -50

100

- 150

-150-100

-50

50

0

100

150

@ Figure 27. A fully relaxed VRB 4, $ surface for glycine dipeptide. Dashed contours denote low-energy regions ranging from 0.0 to 5.0 kcallmole, each contour corresponding to 0.5 kcal/mole. Solid contours range from 5.5 to 21.5 kcal/mole, with each contour representing 1kcall mole.

body has not significantly altered the available conformation of glycine dipeptide, although barriers between conformations are higher in the rigid body approximation. The consequence of this is that the polypeptides and proteins composed of rigid bodies can adopt conformations that give rise to secondary structure, although transition between such structures will occur on slightly longer time scales, or will require higher temperatures for a given time scale. We conclude that the rigid body treatment of backbone fragments still provides a very reasonable approximation to the original flexibility. The independent question of the performance of the virtual potential function can be addressed by comparing Figures 25 and 27. Most of the gross fea-

Table VI Virtual Rigid Body Dynamics Improper Torsion Parameters" $0

ko

CHlE-X-X-VCCC CHlE-X-X-VCC N-CHlE-CH2E-VCO CHlE-N-VCO-VCCP

35.5 35.5 0.0 48.5

55.0 55.0 45.0 45.0

Equilibrium angles in units of degrees, and force constants,

( 4 )in units of kcal/(mole-de2).

tures of Figure 25 are reproduced by the VRBD method, in that there are similar regions of conformational accessibility and exclusion. In addition, most stationary points are present and positioned in the same region of space as the full atomic surface, and all barriers are very reasonable approximations to the full atomic case. Therefore, we conclude that the virtual potential function is adequately reproducing the local interactions of the original surface. Finally, we provide evidence that the VRBD potential surface has a well-defined minimum that corresponds to a typical secondary structural element. An all-alanine heptapeptide was placed in a conformation where all 6 and rc/ dihedral angles were initialized to -47" and -57", respectively, based on accepted values by Creighton l1 for a right-handed a-helix. We then performed Powell minimization using the fully flexible CHARMM potential and the virtual potential with rigid degrees of freedom, and examined the resulting minimized structures. In Table VII we tabulate the five pairs of 6 and rc/ values for the right-handed a-helix of the alanine heptapeptide, for the starting structure," and the final minimized molecule for both the CHARMM5 and VRBD potentials. Clearly the VRBD potential has accomplished one of the original goals set out in this paper, which was to describe important secondary structural elements such as a-helices.

CONCLUSION

Torsion Type

a

=

$2 =

VRB/ Powell

In this paper we have presented a simulation method that begins to address the problems of computational efficiency, accuracy, and suitability of current simulation implementations (such as molecular dynamics and SBMD ) for molecular-level descriptions of polypeptide and protein systems. The approach we have developed involves a reduced representation


of the system in terms of electrostatic multipole expansions, excluded volumes, virtual connectivity degrees of freedom, and rigid body motion of individual amino acids in terms of either minimization or dynamical equations of motion. The computational efficiency and accuracy of the VRBD method will permit applications to probing structural, dynamical, and energetic properties of proteins with ligand binding sites, and the simulation of larger protein systems and longer times than has been reasonably possible until now. In an article to follow shortly, we explore the use of the VRBD method for active site studies. The influence of portions of the protein far removed from the active site should be neglible in the region of binding, and a full atomic molecular dynamics study may be unnecessarily expensive computationally for understanding local molecular interactions at the active site. The SBMD method has successfullydealt with the computational efficiency and accuracy for understanding local active site interactions, when global motions are unimportant. In this forthcoming work we will compare the simulation methods VRBD, SBMD, and molecular dynamics by evaluating the gas-phase trajectory of ribonuclease -A, which is shown to exhibit such large-scale conformational changes.

97

having double-bond character, are present. In these cases the body or bodies present may be treated to a good approximation as a collection of flexibly connected rigid bodies. For rigid systems, the original 3N translational degrees of freedom of the N atomic sites can be replaced by 3 translational degrees of freedom for the center of mass of M bodies, and 3 rotational degrees of freedom about the center of mass of M bodies (minus, in each case, the translation and rotation of the molecule as a whole). The original distribution of point masses comprising each body J are then replaced by a constant mass distribution in terms of moments and products of inertia. Any rigid body motion implementation requires a definition of the center of mass coordinates,

the mass at this site

and the inertia matrix i J , whose elements are the moments and products of inertia

APPENDIX 1. RIGID BODY EQUATIONS OF MOTION In the body of this work we have introduced the description of a protein or a peptide as a collection of flexibly connected rigid bodies. Consider as an example the N-methyl, N-acetyl blocked dipeptide Phe-Phe. One could divide this system into a collection of nine rigid bodies-two phenyl rings, three peptide units, two a-carbons, and two methyl capping groups. The phenyl rings would be flexibly connected by a bond to an a-carbon; the rigid peptide groups would be flexibly attached to the a-carbons and/or the capping methyl groups. In order to evolve such a system in time, we must develop a dynamical scheme based on combining atomic and rigid body equations of motion. The purpose of this appendix is to outline the discretized algorithms used for the time evolution of a system of rigid bodies. In general, the distances between all atoms in a system change due to the forces acting on them; however, in some cases the elasticity of the body may be neglible, especially when degrees of freedom that give rise to bond and angle forces, or torsions

Rigid body theory requires that we define a frame of reference, the body frame, where all vectors are at the same relative orientation to other vectors and the principal (body) axes system. This is accomplished by diagonalizing G to form the diagonal matrix IJ’

where the matrix Ut is used to initialize the rotational variables and the rotation matrix that transforms lab frame vectors (denoted by superscript L ) to body frame vectors (superscript B ) . The initial orientation of all moment arms (the vector between the center of mass and a force site in the body) is obtained by the following transformation:

98


Similarly, electrostatic dipoles also have a body frame reference:

Once all of these initialization conditions are met, we have completely defined the fragment or molecule as a rigid body. In the course of a dynamics or minimization calculation, we will update positions in the lab frame using the relation

We next turn to the consideration of the equations of motion for the translation and rotation of rigid bodies in the microcanonical (NVE) ensemble, which requires that the total energy of the system be a constant of motion. A commonly employed numerical algorithm for evolving the position vector due to translation through time is the Verlet algorithm'?

i (t

+ at) = -i(t

-

at)

+ 2i(t)

+ m - ' f ( t ) ( 6 t ) 2 +O ( 6 t ) 4 where the changes in the lab frame positions are due to changes in U only. Obviously, the convenient definition of the original 3N independent variables for minimization or dynamical motion are the 3 Cartesian directions for each atomic center. However, care must be taken in choosing the rotational variables. The 3 rotational variables known as Euler angles--, 4, and $which are related to the direction cosines-cos 0, cos 4, and cos$-are unstable dynamical variables, due to singularities in the Euler equations of motion, as 0 + 0 or K. An elegant solution to the singularities in the Euler equations of motion is to redefine the three independent variables in terms of four dependent variables known as quaternions, ''*I7 which are related to the Euler angles through the relations

Clearly there is a dependency among the four quaternion variables due to the fact that there are only 3 rotational degrees of freedom; this dependency is expressed as

The quaternion rotation matrix A (which is U at t = 0 ) , is given as

C(t)=

i(t

+ 6t) - i(t

-

26t

at)

+

o(w3

(A12) (A 1 3 )

This algorithm is often used for full atomic molecular dynamic^,^ and shows quite reasonable energy conservation when compared to higher order predictor-corrector methods"; it is also appropriate for the translation of the center of mass of a rigid body, where the force on the center of mass is determined by summing over all of the site forces of the rigid body NEJ I

It is also the appropriate algorithm for evolving the motion of atomic sites in the bath region, as well as for those explicitly represented atoms in the reaction region. For the case of rotational motion (rigid bodies in the bath region), we must update the quaternions using an algorithm that is numerically stable and incurs an error of high order in 6 t so that the total system energy is conserved. The quaternion numerical integration scheme due to Finchem" is based on a two-step Lax-Wendrof algorithm.21One simply has to recognize the analogy between Newton's equation of motion

and the quaternion equations of motion:


=

-41w:

41 =

qowx

=

q3wx

43 =

-42w:

90

42

B

-

B B

99

-

q2w; - q3wz

934

+ qow;

+ qzw? - QlO?

+ q10; + qow:

(AX)

& = [-q-rY3 I - 3(cj.?i,)r;5];ij ?I

Both are of the same form

where q and 2 are 4-vectors, ax, wyt %) and

where the forces in Eq. (A22) are defined in the lab frame as

( 40,

q l , q2, q3), and ( 0 ,

+

(A23)

The purpose of the auxiliary part of the algorithm (A19) is to form the quaternions a t the half step in order t o complete the actual updating of the quaternions a t the full time step 6 t :

2

qo

Lq3

q3 -q2

-q3

40

41

q2

-41

90

I

(AM)

J

The Lax-Wendrof algorithm,21 like the Verlet" method, is based on a central difference expansion of these equations of motion. The Lax algorithm is comprised of the following auxiliary steps:

We prepare for the next energy evaluation by updating the lab frame site positions using the new center of mass positions and quaternions: (:;

t

+ s t ) = ?J( t + s t ) + A-'( t + 6 t ) ?a

(A25)

where i y is the body frame vector from the center of mass t o the site i defined in Eq. ( A 6 ) . Similarly, the dipole vectors are advanced: where J1,and jB are the angular momentum in the lab frame and body frame, respectively, A is the rotation matrix defined in Eq. ( A l l ) , I, is the principal moment of inertia in the x, y, or z direction in the body frame, and 6 t is the full time step. At this point one can determine the kinetic energy a t time t from

T h e lab frame torques are determined from all cross products between orientational vectors and site forces of a rigid body, i.e.,

There may also be contributions t o the torque due t o couples, such as that due to electrostatic dipoles for example,

c;

cc(t

+ s t ) = A-'(t + 6.t);:

(A26)

This completes the dynamical description of polyatomics treated a s rigid bodies. Linear molecules, however, require special treatment due to the fact that one principal moment of inertia is zero. Special algebraic solutions have been proposed'"'' that reduce the number of rotational degrees of freedom by one, and may offer a n implementation benefit due to the reduced amount of storage; this is only true if one is performing a simulation on a homogeneous mixture of linear molecules, otherwise the procedure is cumbersome. Due to the inhomogeneity of protein systems, we require a rigid body implementation for both linear and polyatomic protein fragments. In practice this is handled quite simply by setting up the initial conditions (see the beginning of this section) so that it is always the x component of the principal moment of inertia that is zero. The N loop of dynamics will then be performed in the order of atomic, linear molecules, and finally

100


polyatomic molecules, to avoid if tests as to which angular velocity components should be evaluated in Eqs. (A19) and (A24). This work was supported by a grant from the National Institutes of Health (GM37554-02). THG would like to thank Dr. Martin Head-Gordon for many interesting discussions, and William Young for help with the figures.

REFERENCES 1. Brooks, C. L. & Karplus, M. (1983) J. Chem. Phys. 78,6312-6325. 2. Brunger, A., Brooks, C. L. & Karplus, M. ( 1984) Chem. Phys. Lett. 105,495-500. 3. Brooks, C. L., Brunger, A. & Karplus, M. ( 1985) Biopolymers 24,843-865. 4. Brunger, A., Brooks, C. L. & Karplus, M. (1985) Proc. Natl. Acad. Sci. USA 82,8458-8462. 5. Brooks, B. R., Bruccoleri, R. E., Olafson, B. D., States, D. J., Swaminathan, S. & Karplus, M. (1983) J. Comput. Chem. 4,187-217. 6. Weiner, S . J., Kollman, P. A., Case, D. A., Singh, U. C., Ghio, C., Alagona, G., Profeta, S. & Weiner, P. (1984) J. Am. Chem. SOC.106, 765-784. 7. Jorgensen, W. L. & Tirado-Rives, J. (1988) J. Am. Chem. SOC.110,1657-1666. 8. Levitt, M. L. (1976) J. Mol. Biol. 1 0 4 , 59.

9. Flory, P. J. (1969) Statistical Mechanics of Chain Molecules, Interscience Publishers, New York, chap. 9. 10. Fixman, M. (1974) Proc. Nut. Acad. Sci. USA 7 1 , 3050-3053. 11. Creighton, T. E. (1984) Proteins: Structures and Molecular Principles, w. H. Freeman and Company, New York, chaps. 5 and 6. 12. Jackson, J. D. ( 1975) Classical Electrodynamics, 2nd ed., Wiley, New York, chaps. 3 and 4. 13. Price, S. L. & Stone, A. J. (1988) J. Phys. Chem. 92, 3325-3335. 14. Richards, F. M. & Wyckhoff, H. W., (1971) The Enzymes, 3rd ed., Academic Press, New York, p. 647. 15. Loncharich, R. J. & Brooks, B. R. (1989) PROTEINS: Struct. Funct. Genet. 6, 32-45. 16. Evans, D. (1977) Mol. Phys. 34,317-325. 17. Goldstein, H. ( 1971) Classical Mechanics, AddisonWesley, New York, chap. 4. 18. Verlet, L. (1967) Phys. Rev. 159,98-103. 19. Allan, M. P. & Tildesley, D. J. ( 1987) Computer Simulation of Liquids, Clarendon Press, London, chap. 3. 20. Finchem, D. (1981) CCP5 Quart. 2 , 6-9. 21. Press, W. H., Flannery, B. P., Teukolsky, S. A. & Vetterling, W. T. (1986) Numerical Recipes, The Art of Scientific Computing, Cambridge Univeristy Press,

New York. Received May 2, 1990 Accepted August 31, 1990

New Langevin and gradient thermostats for rigid body dynamics.

Computing the Free Energy along a Reaction Coordinate Using Rigid Body Dynamics.

In silico single-molecule manipulation of DNA with rigid body dynamics.

First-Person Perspective Virtual Body Posture Influences Stress: A Virtual Reality Body Ownership Study.

Ring polymer dynamics for rigid tops with an improved integrator.

Collapsed fat navigators for brain 3D rigid body motion.

The Efficacy of Rigid Bronchoscopy for Foreign Body Aspiration.

Rapid determination of RMSDs corresponding to macromolecular rigid body motions.

Quantitative Validation of a Human Body Finite Element Model Using Rigid Body Impacts.

Virtual Worlds versus Real Body: Virtual Reality Meets Eating and Weight Disorders.

Modeling behavior dynamics using computational psychometrics within virtual worlds.

A Virtual Out-of-Body Experience Reduces Fear of Death.

Drift and ownership toward a distant virtual body.

Developing a Novel Measure of Body Satisfaction Using Virtual Reality.

Modulation of pain threshold by virtual body ownership.

A Virtual Reality Full Body Illusion Improves Body Image Disturbance in Anorexia Nervosa.

Owning an overweight or underweight body: distinguishing the physical, experienced and virtual body.

The body fades away: investigating the effects of transparency of an embodied virtual body on pain threshold and body ownership.

Diversity in neural firing dynamics supports both rigid and learned hippocampal sequences.

A simple molecular mechanics integrator in mixed rigid body and dihedral angle space.

mtsslSuite: In silico spin labelling, trilateration and distance-constrained rigid body docking in PyMOL.

friction tensors of assemblies of spheres.

Bubble dynamics in a compressible liquid in contact with a rigid boundary.

Conformational dynamics of nucleic acid molecules studied by PELDOR spectroscopy with rigid spin labels.