1Makineni Theoretical Laboratories, Department of Chemistry, University of Pennsylvania, 231 South 34th Street, Philadelphia, PA 19104, USA and 2Neutron Science Center and Center for Promotion of Computational Science and Engineering, Japan Atomic Energy Research Institiute, 8-1, Umemidai, Kizu-cho, Souraku-gun, Kyoto 619-0215, Japan
3 To whom correspondence should be addressed. e-mail: saven{at}sas.upenn.edu
![]() |
Abstract |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Keywords: computational protein design/protein oligomers/p53 tetramerization domain/quaternary structure
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Most efforts in computational design apply directed methods. Here directed protein design refers to the search for a sequence (or a small set of sequences) consistent with a predetermined backbone structure. The search is directed toward those sequences having low energy (or a favorable score) when they assume the target structure. Pioneering efforts in design identified proteins having substantial ordering but not necessarily well defined tertiary structures (Bryson et al., 1995). Since exponentially large numbers of sequences are possible even for small proteins, computational methods have accelerated successful design. Exhaustive searching of all mN possible sequences is usually intractable, where N is the number of variable residues and m is the number of residue degrees of freedom, e.g. the number of allowed amino acids and side chain conformations. Genetic algorithms (Jones, 1994
) and simulated annealing methods (Shakhnovich and Gutin, 1993
; Hellinga and Richards, 1994
) search sequence space in a partially random fashion, on average progressing toward lower energy sequences. Such stochastic methods have been used to redesign a variety of natural proteins (Desjarlais and Handel, 1995
, 1999; Desjarlais and Clarke, 1998
; Johnson et al., 1999
; Kraemer-Pecore et al., 2001
) and novel helical proteins (Bryson et al., 1998
). Elimination methods such as dead end elimination provide better estimates of global optima and have also been used to automate both the redesign of sub-regions within natural proteins (Malakauskas and Mayo, 1998
; Strop and Mayo, 1999
; Bolon and Mayo, 2001
; Looger et al., 2003
), as well as full sequence design (Dahiyat and Mayo, 1997
). Directed approaches are sensitive to the energy functions used, which may be problematic given that such energy functions are necessarily approximate. Uncertainties in the energy function may not merit the search for global optima. For cases where information about a problem is incomplete, probabilistic rather than directed approaches to design may be appropriate.
Statistical methods complementary to directed protein design have been developed (Zou and Saven, 2000; Kono and Saven, 2001
). Such methods reveal the features of sequences likely to fold to a particular structure but which may not be thermodynamically optimal. Such probabilistic methods provide site-specific information about the range of allowed amino acid substitutions. Within this methodology, much of the formalism of statistical thermodynamics is recast so as to reveal the properties of sequences likely to fold to a target three-dimensional structure. The site-specific probabilities of the amino acids are determined by maximizing an effective entropy function, subject to constraints on the sequences. Such constraints can be physically based, such as the energy of sequences after they acquire the target backbone structure, or functionally based, such as the patterning of amino acids at predetermined positions to confer metal binding. The theory takes as input (i) a given target structure, (ii) energy functions for quantifying sequencestructure compatibility and (iii) a set of constraints on the sequences. For some forms of the constraints, the approach reduces to a form of heterogeneous mean field theory (Koehl and Delarue, 1996
). The theory yields estimates of the number of sequences and, most importantly, the site-specific probabilities of the amino acids and their side chain conformations. The computation time of the calculations scale is (Nm)2. With statistical methods, in a much shorter time larger numbers N of variable residues can be examined using a larger diversity of residue states m than with other computational methods. Sequences are not explicitly sampled; the calculations yield the amino acid probabilities directly, which are useful in protein design for identifying allowed amino acids at each position in a target structure. The probabilities may also be used both to identify sites tolerant of mutations, where functional properties may be engineered, and to guide the construction of a combinatorial library of protein sequences.
Most computational efforts in protein design have focused on the design of tertiary structures, but some of the earliest design targets were oligomeric structures. The DeGrado group has crafted a variety of dimeric, trimeric and tetrameric helical bundles (Bryson et al., 1995; DeGrado et al., 1999
). Directed computational methods have been used in the design of coiled coils, where the symmetry of the structure facilitates the design (Harbury et al., 1998
; Keating et al., 2001
). Just as in nature, quaternary structure provides a facile route to large, well structured protein systems. It is of interest to extend the probabilistic methods for designing protein and protein libraries to include such an intermolecular structure, which is presented herein. For symmetric arrangements of protomers, a symmetry assumption greatly reduces the computational overhead so that arbitrary oligomers and multidimensional arrays may be included as design targets.
![]() |
Theory |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Overview of statistical theory
A statistical, entropy-based formalism has been developed to identify the features of sequences that are likely to fold to a given backbone structure (Saven and Wolynes, 1997; Zou and Saven, 2000
; Kono and Saven, 2001
). Generally, the method takes as input a given template backbone structure and a function for characterizing sequencestructure compatibility. The calculation yields the site-specific probabilities of the amino acids. Since the method does not involve the explicit generation of sequences, information about exponentially large numbers of sequences may be obtained. This is particularly useful given the large numbers of possible protein sequences (20N for an N-residue protein) and the wide natural diversity of sequences sharing common folds. The generality of the method allows a number of requirements to be prescribed upon sequences as constraints. Such constraints can include both the patterning of residues, e.g. so as to confer solubility, or global energetic constraints, so as to identify the properties of sequences likely to have a particular target structure as their folded state.
The fundamentals of the method have been presented in previous work (Saven and Wolynes, 1997; Zou and Saven, 2000
; Kono and Saven, 2001
), and are briefly reviewed here. The most probable set of site-specific amino acid probabilities at each position is determined by maximizing an effective entropy function subject to imposed constraints using the method of Lagrange multipliers (McQuarrie, 1976
), thus specifying a variational functional V:
V = S 1f1
2f2 ...(1)
The sequence entropy S quantifies the number of sequences likely to fold to a particular structure. The fk are functions that impose the constraint conditions and the k are their corresponding Lagrange multipliers. The variational function V and the constraint function fk depend upon the site-specific probabilities wi(
, r(
)). Each residue position (site) in the structure is labeled by the index i (i = 1, ..., N), where N is the total number of residues. Here wi(
, r(
)) is the probability that residue position i is in a particular state, where the state of the residue is specified by both the identity of the amino acid
and the conformation of its side chain r(
). Generally,
labels any one of the possible amino acids (natural or non-natural), and r(
) labels the members of a discrete set of side chain conformations for each amino acid, so-called rotamer states (Dunbrack, 2002
). Here S is written solely in terms of the one-body probabilities wi(
, r(
)):
In identifying the state probabilities consistent with particular values of the constraints, the k-th constraint function fk is constrained to have a particular value fok:
fok = fk({wi(, r(
))})(3)
For example, such constraints may be used to specify the normalization of the probabilities, to pattern residue types, or to specify the energy sequences acquire when in the folded structure. As a result, the properties of sequences sharing these common constraints can be examined. The site-specific probabilities wi(, r(
)) are determined as those that maximize V subject to given values of fo1, fo2, ... (see Equation 3). It is important to note that although the form of S in Equation 2 would seem to imply that the wi(
, r(
)) are independent, the constraint conditions in Equation 3 will cause them to be coupled to one another. The probability of an amino acid at a particular position may be obtained by summing over the probabilities of its rotamers: wi(
) =
r(
)wi(
, r(
)). This formalism may be applied to different representations of protein structure. In many cases, it is useful to simplify the representation of each residue so as to eliminate explicit side chain degrees of freedom, thus yielding effective united-residue representations and energy functions often used in protein science (Miyazawa and Jernigan, 1996
; Liwo et al., 1998
). For such representations, each amino acid has effectively only one conformational state: wi(
) =
r(
)wi(
, r(
)) = wi(
, r(
)), and for such cases we may omit the variable indicating the rotamer state.
Constraints and protein sequence energetics
The theory may accommodate a wide variety of constraints of the form in Equation 3. Each residue site i must be occupied, i.e. the probabilities wi(, r(
)) are normalized, which leads to:
The hydrophobic effect may be included as an effective one-residue energy that is dependent upon the local density of ß-carbons (Kono and Saven, 2001). The sum of these solvation or environment scores may also be constrained by specifying the average environmental energy as summed over the positions in a particular target structure:
Other multi-body representations of the hydrophobic effect and effective pair interactions between residues may also be incorporated in the theory (Miyazawa and Jernigan, 1985).
For realistic, atom-based representations of proteins, design algorithms that focus primarily on optimizing inter-atomic interactions within the folded state have had substantial success (Hellinga and Richards, 1994; Dahiyat et al., 1997
; Kortemme et al., 1998
; Baker and DeGrado, 1999
; Hellinga, 1999
; Strop and Mayo, 1999
; Bolon and Mayo, 2001
; Looger et al., 2003
). Such methods use atom based potentials to account for both covalent and non-covalent interactions, e.g. van der Waals forces, hydrogen bonds, and electrostatic interactions. The statistical theory of sequence ensembles may be formulated so as to include such atom-based descriptions of both intra- and inter-molecular interactions.
For a protomer structure in an oligomeric complex, the intra-molecular energy Ef of a single chain depends upon the set of amino acid identities {1, ...,
N} and the rotameric state of each of these amino acids r(
i). Ef may be written as a sum of effective one-residue and two-residue interactions:
The indices i and j refer to residue positions in the structure, and the second term sums only over unique interactions between pairs of residues. The one-residue energy i(1)(
i, r(
i)) is the energy associated with locating the amino acid
with conformation r(
) at site i within a single protomer structure. This energy is determined by side-chain backbone interactions or intrinsic structural tendencies of the amino acids, such as preferences for solvent exposure or secondary structure. Similarly, the pair energy
i(2)(
i, r(
i);
j, r(
j)) is the sequence dependent interaction energy between residue i and residue j, each of which is a member of the same protein chain. Such interaction energies may be inferred from a database using reduced descriptions of amino acids (Miyazawa and Jernigan, 1996
; Onuchic et al., 1997
) or, using an atomically detailed model, as a sum over inter-atomic interactions of the identity-rotamer states at sites i and j, as specified by (
i, r(
i)) and (
j, r(
j)).
Similarly, the intermolecular association energy Ea is that due to the intermolecular interactions among sites in a complex of M protein chains, where chains are labeled by the indices m and n. The association energy of the complex Ea may be written in terms of the identity-rotamer states of each site in the complex:
where i(1)(
im, r(
im)) is the additional effective energy associated with locating the amino acid
im with conformation r(
im) at site i on chain m in the complex. For example, a residue that is exposed to solvent in the folded structure of an isolated protomer may become sequestered from solvent or interact with the backbone of an adjacent protomer upon complex formation;
i(1)(
im, r(
im)) quantifies the effective change in energy upon association. For given values of the identities and rotamers at sites i and j on two different protomers, the pair interaction
ij(2)(
im, r(
im);
jn, r(
jn)) is the sequence-dependent intermolecular interaction energy between residue i on chain m and residue j on chain n.
Equation 7 may be applied to any given complex, including hetero-dimers and other proteins having asymmetric quaternary structure. For symmetric complexes comprising identical chains, however, symmetry implies that equivalent positions on each chain have the same identity and conformational state. This symmetry assumption may be expressed as follows:
1. The amino acid identities at equivalent sites on each chain are the same: im =
i for each protomer m = 1, ..., M.
2. Side chains at equivalent sites on different chains take on equivalent rotamer states: r(im) = r(
i) for each m = 1, ..., M.
This assumption is introduced to expedite the design process, though it is perhaps overly stringent, particularly with regard to side chain conformation. For example, distant, exposed residues of a protein need not always have precisely the same conformation in solution. For systems with only one state per site, e.g. the lattice model discussed in the Applications and results (section Lattice model), this symmetry assumption is exact.
Within this symmetry assumption, Equation 7 may be written in a form analogous to Equation 6, involving effective interactions between sites on a single chain:
where
We note that Ea now has a form exactly analogous to Ef. Each residue sees other residues on the same chain and images of itself and other residues via intermolecular interactions with other chains.
For a given set of constraint conditions, the values of Ef and Ea are assumed to be well represented by their sequence averages (Zou and Saven, 2000; Kono and Saven, 2001
). With this approximation and the factorization implicit in Equation 2, Ef and Ea may be expressed as functions of the site probabilities wi(
, r(
)):
Note that we sum over amino acids at each site and no longer consider just a single sequence as would be indicated by the set {i}. The Ef and Ea are now each functions of the site-specific probabilities wi(
, r(
)) and may now appear as constraints of the form suggested in Equation 3. With the symmetry assumption, the number of variables decreases dramatically relative to a calculation that treats each site i in the complex separately, i.e. a calculation where N is the total number of residues in the complex rather than in a single protomer. For symmetric structures, only the wi(
, r(
)) of residues on a single chain need be determined. There are typically Nxm such independent variables, where m is the number of allowed identity rotamer states at each site and N is the number of residues per protomer.
Other energetic constraints may be imposed on the sequences as well. For simplified representations of the amino acids, in designing protein sequences it is often important to account for non-target structures that a sequence may also acquire. For a foldable sequence, the target state should be energetically removed from these other possible collapsed structures. There are many ways to achieve this (Saven, 2001), but perhaps the simplest is to optimize a stability gap
f = Ef
Ef
u, where
Ef
u is an average over the an appropriate ensemble of unfolded structures. Similarly, we may also define a binding stability gap
b = Ea
Ea
b, where
Ea
b represents an energy average over configurations of the complex other than the target quaternary structure. Herein, these alternate configurations involve different relative orientations of the individually folded protomers. The nature of the configurational averaging
...
b depends on the particular model and may be impractical for atomically detailed models due to both the large number of unfolded conformations and the large numbers of identity-rotamer states. For contact type energy functions, however, this averaging is straightforward (Zou and Saven, 2000
), and an example is presented in the next section.
![]() |
Applications and results |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Lattice model
The model consists of a 27 residue self-avoiding polymer whose monomers occupy sites on a cubic lattice (Lau and Dill, 1990; Shakhnovich and Gutin, 1990
; Leopold et al., 1992
). Two types of monomers, hydrophobic (H) and polar (P), are used to construct sequences. Each of these residue types has a single conformational state: wi(
, r, (
)) = wi(
), where
= H or P. The 227 possible sequences may be exhaustively enumerated. As a result, this model may be solved exactly, and the results may be compared with theoretical results. A simplified potential developed for HP type monomers (Li et al., 1996
; Zou and Saven, 2000
) is used that contains only two-body interactions, those interactions make the calculations non-trivial. The two body interactions are non-zero only if the amino acids are nearest neighbors on the lattice but are not bonded to one another. Let rij be the distance between residues i and j and r0 is the distance between nearest neighbors on the lattice. Then the contact parameter
ij is non-zero only if rij
r0, for which
ij = 1. The contact variable for intermolecular interactions appears as
ijmn, which is unity only if residues i and j of chains m and n, respectively, are in contact with one another and
ijmn = 0 otherwise. For this model, the folded state energy and the association energy then take the following form (see Equations 12 and 13):
The same energy function is used for both folding and for association: (2)(
;
') =
(2)(
;
'). The contact energies are chosen as (Li et al., 1996
; Zou and Saven, 2000
):
(2)(Hi, Hj) = 3
,
(2)(Hi, Pj) =
(2)(Pi, Hj) =
, and
(2)(Pi, Pj) = 0. Equating
(2)(
,
') =
(2)(
,
') corresponds to treating inter-residue interactions in the same fashion for both intra-molecular (folding) and inter-molecular (specific oligomerization) organization. As necessary, distinct potentials for intra-molecular folding and for inter-molecular association may be used.
The energetic constraints in this lattice model study are the folding stability gap (f) and binding stability gap (
b). For such folding criteria, ensembles of unfolded and mis-associated states are necessary. The choice of folded structure is arbitrary, and a structure that is the conformational ground state for a large number of sequences is chosen in this study (Li et al., 1996
). The remaining 103 345 compact, cubic structures of the 27-mer are chosen as the ensemble of unfolded states in the calculation of
Ef
u. The target structures of the complex are those depicted in Figure 1. As an ensemble of mis-associated states for the lattice proteins, we consider arrangements of the oligomer for which each chain takes on the target folded structure and the interface between two protomers has nine contacts, i.e. the faces of the protomers are in registry with one another (Figure 1). This loosely mimics the expectation that the mis-associated states of a particular complex most likely to compete with the target structure are those involving large numbers of residues in inter-molecular contact. For the dimer, there are 84 unique configurations of the two protomers.
|
|
|
For realistic representations of proteins, an all-atom model is used that includes side chain conformational degrees of freedom. Side chains assume discrete conformations, rotamers (Ponder and Richards, 1987; Dunbrack and Cohen, 1997
; Tuffery et al., 1997
). Here we use the backbone-dependent rotamer library of Dunbrack and coworkers (Dunbrack and Cohen, 1997
) The Amber force field (Weiner et al., 1984
) is used to calculate non-bonded interactions. This same potential is used to quantify both the intra- and inter-molecular energies. To address the hydrophobic effect, the environmental energy (Equation 5) for the complex is constrained to the value of the wild-type sequence (Kono and Saven, 2001
). The site-specific probabilities are determined for multiple (constrained) values of E = Ef + Ea, and conjugate to this energy is an effective inverse temperature ß. In what follows, the probabilities are presented at ß = 0.5 mol/kcal, which is also the effective temperature at which (unfolded) reference energies of the amino acids are determined, as described in Calhoun et al. (2003
).
We select the tetramerization domain of p53 tumor suppressor (PDB code: 1C26) as the target structure (Figure 4). This transcription factor is involved in a number of important physiological roles, including regulation of the cell cycle, apoptosis, DNA repair and angiogenesis, and mutations of p53 have been linked to cancers (Vogelstein et al., 2000). The domain considered here is a tetramer with four identical 32 residue chains. In this study, all 20 amino acid residues are allowed at each position. Thus, a total of 2032
4x1041 sequences are possible for a single chain. Taking side chain conformations into account, there are a total of 320 states per residue position. In the absence of the symmetry assumption, there are 128 variable positions and 320128 possible states. The imposition of the symmetry constraint reduces this complexity to 32 variable positions and 32032 possible states; the number of independent site probabilities to be determined is reduced from 128x320 to 32x320.
|
|
![]() |
Summary |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
Acknowledgements |
---|
![]() |
References |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Bolon,D.N. and Mayo,S.L. (2001) Proc. Natl Acad. Sci. USA, 98, 1427414279.
Bryson,J.W., Betz,S.F., Lu,H.S., Suich,D.J., Zhou,H.X., ONeil,K.T. and DeGrado,W.F. (1995) Science, 270, 935941.[Abstract]
Bryson,J.W., Desjarlais,J.R., Handel,T.M. and DeGrado,W.F. (1998) Protein Sci., 7, 14041414.
Calhoun,J., Kono,H., Lahr,S., Wang,W., DeGrado,W.F. and Saven,J.G. (2003) J. Mol. Biol., 334, 11011115.[CrossRef][ISI][Medline]
Dahiyat,B.I. and Mayo,S.L. (1997) Science, 278, 8287.
Dahiyat,B.I., Sarisky,C.A. and Mayo,S.L. (1997) J. Mol. Biol., 273, 789796.[CrossRef][ISI][Medline]
DeGrado,W.F., Summa,C.M., Pavone,V., Nastri,F. and Lombardi,A. (1999) Annu. Rev. Biochem., 68, 779819.[CrossRef][ISI][Medline]
Desjarlais,J.R. and Clarke,N.D. (1998) Curr. Opin. Struct. Biol., 8, 471475.[CrossRef][ISI][Medline]
Desjarlais,J.R. and Handel,T.M. (1995) Protein Sci., 4, 20062018.
Desjarlais,J.R. and Handel,T.M. (1999) J. Mol. Biol., 290, 305318.[CrossRef][ISI][Medline]
Dunbrack,R.L. (2002) Curr. Opin. Struct. Biol., 12, 431440.[CrossRef][ISI][Medline]
Dunbrack,R.L.J. and Cohen,F.E. (1997) Protein Sci., 6, 16611681.
Harbury,P.B., Plecs,J.J., Tidor,B., Alber,T. and Kim,P.S. (1998) Science, 282, 14621467.
Hellinga,H. (1999) FASEB J., 13, A1430.
Hellinga,H.W. and Richards,F.M. (1994) Proc. Natl Acad. Sci. USA, 91, 58035807.[Abstract]
Johnson,E.C., Lazar,G.A., Desjarlais,J.R. and Handel,T.M. (1999) Struct. Fold. Design, 7, 967976.[ISI]
Jones,D.T. (1994) Protein Sci., 3, 567574.
Keating,A.E., Malashkevich,V.N., Tidor,B. and Kim,P.S. (2001) Proc. Natl Acad. Sci. USA, 98, 1482514830.
Koehl,P. and Delarue,M. (1996) Curr. Opin. Struct. Biol., 6, 222226.[CrossRef][ISI][Medline]
Kono,H. and Saven,J.G. (2001) J. Mol. Biol., 306, 607628.[CrossRef][ISI][Medline]
Kortemme,T., Ramirez-Alvarado,M. and Serrano,L. (1998) Science, 281, 253256.
Kraemer-Pecore,C.M., Wollacott,A.M. and Desjarlais,J.R. (2001) Curr. Opin. Chem. Biol., 5, 690695.[CrossRef][ISI][Medline]
Kuhlman,B.K. and Baker,D. (2000) Proc. Natl Acad. Sci. USA, 97, 1038310388.
Lau,K.F. and Dill,K.A. (1990) Proc. Natl Acad. Sci. USA, 87, 638642.[Abstract]
Leopold,P.E., Montal,M. and Onuchic,J.N. (1992) Proc. Natl Acad. Sci. USA, 89, 87218725.[Abstract]
Li,H., Helling,R., Tang,C. and Wingreen,N. (1996) Science, 273, 666669.[Abstract]
Liwo,A., Kazmierkiewicz,R., Czaplewski,C., Groth,M., Oldziej,S., Wawak,R.J., Rackovsky,S., Pincus,M.R. and Scheraga,H.A. (1998) J. Comput. Chem., 19, 259276.[CrossRef][ISI]
Looger,L.L., Dwyer,M.A., Smith,J.J. and Hellinga,H.W. (2003) Nature, 423, 185190.[CrossRef][ISI][Medline]
Malakauskas,S.M. and Mayo,S.L. (1998) Nat. Struct. Biol., 5, 470475.[ISI][Medline]
McQuarrie,D.A. (1976) Statistical Mechanics. Harper and Row, New York.
Miyazawa,S. and Jernigan,R.L. (1985) Macromolecules, 218, 534552.
Miyazawa,S. and Jernigan,R.L. (1996) J. Mol. Biol., 256, 623644.[CrossRef][ISI][Medline]
Onuchic,J.N., Luthey-Schulten,Z. and Wolynes,P.G. (1997) Annu. Rev. Phys. Chem., 48, 539594.
Ponder,J.W. and Richards,F.M. (1987) J. Mol. Biol., 193, 775791.[ISI][Medline]
Raha,K., Wollacott,A.M., Italia,M.J. and Desjarlais,J.R. (2000) Protein Sci., 9, 11061119.[Abstract]
Saven,J.G. (2001) Chem. Rev., 101, 31133130.[CrossRef][ISI][Medline]
Saven,J.G. and Wolynes,P.G. (1997) J. Phys. Chem. B, 101, 83758389.[CrossRef][ISI]
Shakhnovich,E. and Gutin,A. (1990) J. Chem. Phys., 93, 59675971.[CrossRef][ISI]
Shakhnovich,E.I. and Gutin,A.M. (1993) Protein Eng., 6, 793800.[ISI][Medline]
Strop,P. and Mayo,S.L. (1999) J. Am. Chem. Soc., 121, 23412345.[CrossRef][ISI]
Tuffery,P., Etchebest,C. and Hazout,S. (1997) Protein Eng., 10, 361373.[CrossRef][ISI][Medline]
Vogelstein,B., Lane,D. and Levine,A.J. (2000) Nature, 408, 307310.[CrossRef][ISI][Medline]
Weiner,J.S., Kollman,P.A., Case,D.A., Singh,U.C., Ghio,C., Alagona,S.J.,G. Profeta and Weiner,P. (1984) J. Am. Chem. Soc., 106, 765784.[ISI]
Zou,J. and Saven,J.G. (2000) J. Mol. Biol., 296, 281294.[CrossRef][ISI][Medline]
Zou,J. and Saven,J.G. (2003) J. Chem. Phys., 118, 38433854.[CrossRef][ISI]
Received August 20, 2003; revised October 24, 2003; accepted October 28, 2003