Homology modeling using simulated annealing of restrained molecular dynamics and ...

15 Pages·2002·8.18 MB·English

Prorein Science (1997), 6956-970. Cambridge University Press. Printed in the USA.
Copyright 0 1997 The Protein Society
Homology modeling using simulated annealing of
restrained molecular dynamics and conformational
search calculations with CONGEN:
Application in predicting the three-dimensional
structure of murine homeodomain Msx-1
HAICHENG LI,’.’ ROBERTO TEJER0,1.3D ANIEL MONLEON? DONNA BASSOLINO-KLIMAS?
CORY ABATE-SHEN,’s5 ROBERT E. BRUCCOLERI? AND GAETANO T. MONTELIONE’z2
‘Center for Advanced Biotechnology and Medicine, Piscataway, New Jersey 08854-5638
*Department of Molecular Biology and Biochemistry, Rutgers University, Piscataway, New Jersey 08854-5638
’Departamento de Quimica Fisica. Universidad de Valencia, Dr. Moliner 50, 46100-Bujassot, Valencia, Spain
4Department of Macromolecular Structure, Bristol-Myers Squibb Pharmaceutical Research Institute, Princeton,
New Jersey 08543
’Department of Neuroscience and Cell Biology, UMDNJ-Robert Wood Johnson Medical School, Piscataway,
New Jersey 08854-5638
(RECEIVED Octobe1r5 , 1996; ACCEPTEDFe bruary 4, 1997)
Abstract
We have developed an automatic approach for homology modeling using restrained molecular dynamics and simulated
annealing procedures, together with conformational search algorithms available in the molecular mechanics program
CONGEN (Bruccoleri RE, Karplus M, 1987, Biopolymers 26137-168). The accuracy of the method is validated by
“predicting” structures of two homeodomain proteins with known three-dimensional structures, and then applied to
predict the three-dimensional structure of the homeodomain of the murine Msx-1 transcription factor. Regions of the
unknown protein structure that are highly homologous to the known template structure are constrained by “homology
distance constraints,” whereas the conformations of nonhomologous regions of the unknown protein are defined only
by the potential energy function. A full energy function (excluding explicit solvent) is employed to ensure that the
calculated structures have good conformational energies and are physically reasonable. As in NMR structure determi-
nations, information on the consistency of the structure prediction is obtained by superposition of the resulting family
of protein structures. In this paper, our homology modeling algorithm is described and compared with related homology
modeling methods using spatial constraints derived from the structures of homologous proteins. The software is then
used to predict the DNA-bound structures of three homeodomain proteins from the X-ray crystal structure of the
engrailed homeodomain protein (Kissinger CR et al., 1990, Cell 63:579-590). The resulting backbone and side-chain
conformations of the modeled yeast MatcY2 and D. melunogaster Antennapedia homeodomains are excellent matches to
the corresponding published X-ray crystal (Wolberger C et al., 1991, Cell 67517-528) and NMR (Billeter M et al.,
1993, J Mol Biol234:1084-1097) structures, respectively. Examination of these structures of Msx-1 reveals a network
of highly conserved surface salt bridges that are proposed to play a role in regulating protein-protein interactions of
homeodomains in transcription complexes.
Keywords: conformational energy calculations; conserved surface salt bridges; DNA-binding protein; transcription
factor
-
Reprinrt e quests to: Gaetano T. MontelioneC, ABMR, utgers University, 679 HoesL aneP, iscataway, New Jersey 08854-5638; e-mail:
[email protected].
Abbreviations: Antp, homeodomain of the D. melanogusrer Antennapedia protein; DG, distance geometry calculations using metric-mauix embedding
methods; Conf. E., total conformational energy including electrostatic effects, computed from the CHARMM potential function; Mata2h, omeodomain of the
yeast Matd protein; Msx-1, homeodomain of the murine Msx-1 protein; pdf, probability density function; RMSD, RMS deviation; SARMD, simulated
annealing with restrained molecular dynamics; VDW E., van der Waals energy computed from the Lennard-Jones portion of the CHARMM potential function.
956
Homology modeling of homeodomain Msx-1 with CONGEN 957
The homeodomain is a highly conserved sequence-specific DNA- Struhl, 1996; Rivera-Pomar et al., 1996). Thus, understanding the
binding domain that has been found in many transcription factors. atomic basis of specific homeodomaidnucleic acid interactions is
First discovered in Drosophila melanogaster, homeodomains have essential to the comprehension of the complex interactions that
been found in almost every organism from nematodes to humans lead to transcriptional (and perhaps translational) regulation in
(Kessel & Gruss, 1990; Wang et al., 1993), and have been found development.
to play a fundamental role in directing embryogenesis (Gehring, Msx-1 is a homeodomain protein encoded by a member of a
1987; Scott et al., 1989; for a recent review see Krumlauf, 1994). relatively small subfamily of homeobox genes expressed in cra-
The sequence of the homeodomain corresponds to 60 amino acid niofacial structures, the neural tube, and the limbs of the develop-
residuest hat assemblei ntot hreea -helicesa ndo ne flexible ing embryo (Davidson et al., 1991; Krumlauf, 1994). Expression
N-terminal arm (Scott et al., 1989;K issinger et al., 1990; Laughon, of Msx-1 in myoblasts inhibits terminal differentiation and induces
1991; Wolberger, 1996). Structures of several homeodomains in cell transformation (Song et al., 1992). Msx-1 interacts with DNA
both unliganded and DNA-bound states have been determined by sites that contain the consensus sequence (C/G)TAA?TG (Catron
NMR and X-ray crystallography (Otting et al., 1988; Qian et al., et al., 1993) and functions as a transcriptional repressor both in
1989, 1993, 1994; Kissinger et al., 1990; Phillips et al., 1991; vitro and in vivo (Catron et al., 1995, 1996). Like many other
Wolbergeer t al., 1991A ; ssa-Munt et al., 1993B ; illeter homeodomains, Msx-1 appears to regulate cellular proliferation by
et al., 1993, 1996; Ceska et al., 1993; Cox et al., 1993; Leiting its ability to repress differentiation-specific target genes.
et al., 1993; Klemm et al., 1994; Sivaraja et al., 1994; Li et al., Because the relative positions of certain structurally and/or func-
1995). These protein structures sharea common three-helical chain tionally crucial atoms should be similar among a family of homol-
fold in which the second and third helices are arranged in a helix- ogous proteins, the three-dimensional structure of a protein can
turn-helix motif. In thef ive homeodomain protein-DNA com- often be modeled reliably based on the known structures of ho-
plexes that are available (Kissinger et al., 1990; Wolberger et al., mologous proteins. Homology modeling has been applied success-
1991; Billeter et al., 1993; Klemm et al., 1994; Li et al., 1995), the fully to a number of proteins (see, for example, Greer, 1985, 1990;
second helix of the helix-turn-helix motif (helix 111 in Fig. 1) Chothia et al., 1986; Palmer et al., 1986; Havel & Snow, 1991;
contacts DNA in the major groove and the “N-terminal arm” (res- Bruccoleri & Novotny, 1992; Bajorath et al., 1993;B rocklehurst &
idues 1-9 in Fig. 1) contacts DNA in the minor groove. Perham, 1993; Fogolari et al., 1993; Havel, 1993; Srinivasan et al.,
Homeodomains have a high degree of conservation in primary 1993; for a comprehensive review see Sali, 1995). Such homology
sequence, tertiary structure, and their mode of interaction with models are very useful in certain protein engineering applications
DNA (Gehring et al., 1994; Wolberger, 1996). They all bind to that do not require an accurate high-resolution structure, and for
consensus sites that contain a TAAT core motif (Odenwald et al., accelerating experimental structure determinations by NMR and
1989; Laughon, 1991; Catron et al., 1993; Gehring et al., 1994; X-ray crystallography using molecular replacement methods. Sev-
Wolberger, 1996). Despite their high conservation, most homeo- eral successful approaches for homology modeling (reviewed in
domain proteins function by binding to specific duplex DNA se- Sali, 1995) have included knowledge-based interactive model build-
quences, and even small differences in their sequence specificity ing (Blundell et al., 1983; Claessens et al., 1989; Bazan, 1990;
and affinity can have biological significance (Laughon, 1991; Des- Bajorath et al., 1993), systematic conformational search (Brucco-
sain et al., 1992; Ekker et al., 1992; Catron et al., 1993; Kornberg, leri & Novotny, 1992), combinatorial side-chain conformational
1993; Gehring et al., 1994). Recently, homeodomain proteins have analysis (Ponder & Richards, 1987), polypeptide segment match-
also been found to recognize specific RNA molecules (Dubnau & ing (Levitt, 1992), conformational threading (Jones et al., 1992),
distance geometry and/or simulated annealing calculations using
homology constraints (Engh et al., 1990; Havel & Snow, 1991;
Fujiyoshi-Yoneda et al., 1991; Brocklehurst & Perham, 1993, Havel,
Helix I 1993; Snow, 1993; Srinivasan et al., 1993; Sudarsanam et al.,
5 IO 15 20
1994), and structure generation by satisfaction of spatial restraints
derived from sequence alignments and expressed as probability
density functions (Sali & Blundell, 1993). However, there is no
general agreement on what is the most reliable method for this
process. In this regard, it is important to test methods for homology
modeling on available NMR and X-ray crystal structures in order
1 Helix I1
to validate their precision and accuracy. It is also important to
23 B b c 25 30 35
develop a database of homology-modeled structures of proteins for
which no NMR or X-ray crystal structure is yet available, because
later these can be compared with experimentally determined struc-
tures to evaluate weaknesses and strengths of different methods.
Many of the existing methods for homology modeling men-
L Helix 111 1 tioned above result in structures with high-energy overlaps of non-
40 45 50 55 bonded atoms and other unreasonable energetic features; when
energy is simply minimized, these structures can become trapped
in local energy minima (Scheraga, 1984). Recently, distance ge-
ometry (Havel & Snow, 1991; Havel, 1993; Srinivasan et al.,
1993; Sudarsanam et al., 1994), simulated annealing with molec-
Fig. 1. Amino acid sequence alignments for four homeobox domains: en- ular dynamics (Engh et al., 1990; Fujiyoshi-Yoneda et al., 1991;
grailed (en), Msx- 1, antennapedia (antp). and Mata2 ((~2). Brocklehurst & Perham, 1993), and variable target function (Sali
958 H. Li et al.
& Blundell, 1993) minimization methods like those used for pro- Helix 111 of homeodomain Msx-1 from mouse
tein structure determination from NMR data have been found to be
Our approach for homology modeling by SARMD was first tested
adopted easily for generating families of homology-modeled struc-
by computing the 3D structure of the DNA recognition helix of
tures from “homology constraints.’’ Like NMR structure determi-
Msx-I. This polypeptide fragment contains 17 amino acid resi-
nations, the results of these structure generation calculations can be
dues, 11 of which are consensus residues that are conserved in
represented as a family of structures, each of which satisfies the
nearly all homeobox domains. Fourteen residues are identical in
“homology constraints.’’ Superposition of this family of solutions
this segment of Msx-1 and engrailed homeodomains (Fig. 1). A set
providesi nformationa boutt hec onsistency of theh omology-
of 25 structures of Msx-1 helix 111 was calculated using 323 “ho-
modeled structures given the input homology assumptions. Al-
mology constraints,” and the 10 with the smallest residual con-
thoughh omologym odelsg enerated with DG algorithms may
straint violations and lowest conformational energies (as explained
sometimes exhibit high energies and biases due to inadequate sam-
in the Materials and methods) were selected for structural analysis.
plings of solution space, these biases can be eliminated and the
The result of these calculations is a well defined a-helix. A sum-
conformational energies improved by further refinement of struc-
mary of structural statistics is reported in Table 1 and a summary
tures using restrained molecular dynamics (Havel, 1993; Sali &
of residual “homology constraint” violations and conformational
Blundell, 1993).
energies is presented in Table 2. The RMSDo f the backbone atoms
We have developed an automatic and objective approach for
within the bundle of structures is 0.3 8,. Most side chains are also
homology modeling using simulated annealing with restrained mo-
well-defined in the predicted structure, with RMSDs for all heavy
lecular dynamics and conformational search methods available in
atoms of 0.8 8, within the bundle. Van der Waals energies for the
the molecular mechanics program CONGEN (Bruccoleri & Kar-
predicted structure are very low, ranging from -4.2 to -4.7 kcall
plus, 1987, 1990; Bruccoleri, 1993; Bassolino-Klimas et al., 1996).
mol-residue, whereas deviations from ideal bond angles, bond
Our hybrid approach (SARMDkearch) avoids the metric matrix
lengths, and planar peptide bonds (Table 1) are within the range of
embedding step used in similar homology modeling methods (Havel
deviations typically seen in high-quality X-ray and NMR structure
& Snow, 1991; Havel, 1993; Srinivasan et al., 1993) or the use of
determinations. These results show that SARMD calculations with
variable-target function minimization (Sali & Blundell, 1993), and
CONGEN can be used to generate homology models of this poly-
generates structures directly from extended starting conformations
peptide segment that exhibit good conformational energies and
using restrained molecular dynamics calculations. These structures
satisfy the set of homology constraints.
are further refined using conformational search methods (Brucco-
leri & Karplus 1987, 1990; Bruccoleri et al., 1988; Bruccoleri &
Homeodomain Mata2 from yeast
Novotny, 1992) for improved sampling of low-energy conforma-
tions of poorly defined surface loops and side chains. Regions of In order to validate the predictive value of our SARMDlsearch
the unknown protein structure that are highly homologous to the method, we next homology modeled the three-dimensional struc-
known template structure are constrained by “homology distance ture of the yeast Mata2 homeodomain, whose crystal structure has
constraints,” whereas the conformations of nonhomologous re- been determined to 2.7 A resolution (Wolberger et al., 1991). The
gions of the unknown protein are defined only by the potential X-ray structure of the DNA-bound engrailed homeodomain protein
energy function. A full energy function (excluding explicit solvent) was again used as a template. Mate2 has 27% sequence identity
is employed to ensure that the calculated structures have good with the engrailed homeodomain. Relative to engrailed, the poly-
conformational energies and are physically reasonable. Informa- peptide sequence of Mata2 also includes a tripeptide insertion in
tion on the consistency of the structure prediction is obtained by the polypeptide loop segment between helices I and I1 (Fig. I),
superposition of the resulting family of protein structures. which is not present in the other homeodomains that we studied. In
In this paper, our homology modeling algorithms are described generating homology models of Mata2 from the structure of en-
and used to predict the three-dimensional structures of the yeast grailed, we assumed that the relative positions of the three helices
homeodomain Mata2 (Scott et al., 1989), the D. melanogaster are well conserved, because the overall homeodomain structure
Antennapedia homeodomain (Antp) (Schneuwly et al., 1986), and appears to be very well conserved in nature, but that the confor-
the murine homeobox domain Msx-1 (Catron et al., 1993) from the mation of this interhelical loop can adjust to accommodate the
X-ray crystal structure of the D. melanogaster engrailed homeo- tripeptide insert. Accordingly, homologous atoms in the Mata2
domain (Kissinger et al., 1990). The resulting backbone and side- homeodomain were defined relative to the engrailed template se-
chain conformations of the modeled homeodomains Antp and Mata2 quence, except in the vicinity of the inserted octapeptide segment
are in excellent agreement with the published NMR (Billeter eta l., of residues 23-27 (Fig. I). No homologous atoms are defined in
1993) and X-ray crystal (Wolberger et al., 1991; Li et al., 1995) this loop region (i.e., for polypeptide segment Asn 23-Ile 23a-
structures, respectively. Comparisons of these structures reveal a Glu 23b-Asn 23c-Pro 24-Tyr 25-Leu 26-Asp 27) because the
conserved network of surface saltb ridges common to thee n- insertion of three residues into the original five-residue-long se-
grailed, Mata2, Antp, and Msx-l homeodomains. This family of quence would likely change the conformation of the loop. As a
predicted structures of the Msx-1 homeodomain has been depos- result, no constraints are applied in the loop region during the
ited in the Brookhaven Protein Data Bank for future comparisons simulated annealing procedure and its conformation is determined
with NMR and X-ray structures once they are available. solely by energy considerations.
A set of 25 homology-modeled structures of Mata2 were cal-
culated from 919 “homology constraints,’’ and the 10 with the
Results
smallest residual constraint violations and lowest conformational
Summaries of our homology modeling algorithms and the proto- energies were selected for structural analysis. This family of 10
cols used for SARMDkearch calculations are presented in Mate- structures is shown in Figure 2A. The convergence within these IO
rials and methods. structures is quite good, except in the loop region. Atomic RMSDs
Homology modeling of homeodomain Msx-1 with CONGEN 959
Table 1. Statistics for homology-modeled structures
Msx-I Helix HIa Mata2 Antp MSX-1
Homology constraints
Total number of heavy atoms 177 532 555 507
Total number of homologous heavy atoms 142 347 444 375
Total number of homology constraintsb 919 323 1,000 1,000
Intraresidue [i = j] 10 2 4 4
Sequential [(i - j) = I] 37 28 34 30
Longer-range [(i - j) > I] 276 889 962 966
Structural statistics
Number of final structures 10 10 10 10
RMS constraint violation 0.0067 A 0.0088 A 0.0054 8, 0.0053 8,
RMSDs from ideal polypeptide geometries
Bond angles 2.39" 1.78" 2.23" 2.00"
Bond lengths 0.010 A 0.010 A 0.011 A 0.010 A
Peptide bond ws 2.53" 0.66" 0.63" 2.19"
Atomic RMSD' to average homology-modeled structure
Backbone (N, C", C') 0.3 Aa 0.5 Ad 0.4 & 0.3 8'
All heavy atoms 0.8 Aa 1.2 Ad 1.0 Ae 0.9 A'
Atomic RMSD' to X-ray crystal or average NMR structure
Backbone (N, C", C') - 0.9 Ad 0.8 Ae -
All heavy atoms - 2.0 Ad 1.7 k -
Atomic RMSD' of engrailed template structure to X-ray crystal or average NMR structure
Backbone (N, C", C') - 0.9 Ad 0.8 Ae -
aResidues 42-58 of Msx-1.
bEach homology constraint corresponds to one upper-bound and one lower-bound distance constraint.
'RMSDs are computed only for residues with backbone dihedral order parameters (Hyberts et al., 1992). S(q5) + S(@) > 1.5.
Residue numbers used here are defined in Figure 1.
dFor residue range 12-22, 32-60, excluding poorly defined N-terminal, loop, and C-terminal residues.
eFor residue range 4-57, excluding poorly defined N-terminal and C-terminal residues.
'For residue range 3-56, excluding poorly defined N-terminal and C-terminal residues.
for residues 7-58 (excluding thel oop region) are 0.5 8, for Homeodomain Antennapedia from D. melanogaster
backbone atoms and1 .2 8, for all heavy atoms (Table 1). Re-
sidual constraint violations and conformational energies for this We next modeled a second homeodomain, the Antennapedia ho-
family of homology-modeled Mata2 structures are presented in meodomain (Antp) from D. melanogaster. As before, we used the
Table 2. The VDW energies for the 10 calculated structures engrailed homeodomain crystal structure determined in complex
range from- 4.5 to -4.9 kcal/mol-residue, and are similar to with double-stranded DNA (Kissinger eta l., 1990) ast he structural
values obtained in CONGEN SARMD calculations on the small template. The sequence identity between the Antp and engrailed
protein crambin using distance constraints derived from the X-ray homeodomains is 51% (Fig. 1). A family of2 5 structures was
crystal structure (Bassolino-Klimas et al., 1996) and for the struc- calculated using 1,000 homology constraints and the 10 structures
tures of small proteins using experimental NMR data (Tejero with lowest constraint violations and conformational energies were
et al., 1996). selected for structural analysis. A summary of structural statistics
The backbone conformation of each homology-modeled Mata2 for this family is presented in Table 1 and a summary of the
structure is very similar to that of the X-ray structure (Fig. 2A). residual constraint violations and energies is showni n Table 2. The
Like the engrailed homeodomain template structure, the crystal convergence within the family of 10 calculated structures is very
structure of Mata2 (Wolberger et al., 1991) was determined in good, with RMSD values of 0.4 8, for the backbone and 1.0 8, for
complex with DNA. A statistical summary of the comparison be- all heavy atoms of residues 4-57. Residual violations and confor-
tween the homology-modeled and X-ray structures of Mata2 is mational energies for these 10 structures of Antp are summarized
also shown in Table 1. The average RMSD of the backbone con- in Table 2. The VDW energies of these homology-modeled Antp
formation (the loopr egion excluded) between the modeled and the structures range from -5.2 to -5.7 kcal/mol-residue; these values
X-ray structures is 0.9 A. When individual helices are compared, are typical of good-quality structures in the CONGEN force field
the corresponding RMSD values for backbone atoms range from (Bassolino-Klimas et al., 1996; Tejero et al., 1996).
0.3 to 0.7 8,. The loop region is not well-defined in the homology These 10 homology-modeled structures were superimposed on
models; instead, it is predicted to exist in two families of confor- the average NMR structure of DNA-bound Antp (Fig. 2B). The
mations, one of which includes the backbone conformation ob- backbone conformations of the homology-modeled structures are
served in the crystal structure (Fig. 2A). very similar to the average NMR-determined structure (Table l),
960 H. Li et al.
Table 2. Summary of residual homology-constraint violations and final energies
Numbers of violations and energies for each homology-modeled structure
Structure no. 1 2 3 4 5 6 7 8 9 10
Msx-1 Helix IIIa
>0.2 A 0 0 0 0 0 0 0 0 0 0
0.1-0.2 A 1 2 0 3 2 4 3 1 4 2
<0.1 A 13 13 IO 12 8 9 12 15 10 I1
Total 10 14 15 13 15 10 15 16 14 13
VDW E? -71.2 -74.2 -77.5 -73.8 -73.7 -78.1 -72.2 -80.1 -77.0 -75.0
Conf. E.b -317.0 -278.9 -287.2 -299.5 -283.2 -300.5 -282.9 -275.0 -285.0 -258.9
Mata2
>0.2 A 0 0 0 0 0 0 0 0 0 0
0.1-0.2 A 2 3 1 0 2 2 2 0 2 I
<0.1 A 66 65 56 43 70 49 58 40 65 43
Total 68 68 57 43 72 51 60 40 67 44
VDW E? -285.9 -279.7 -304.7 -297.1 -277.1 -277.4 -289.4 -304.4 -287.9 - 292. I
Conf. E? -1,216.6 - 1,269.1 - 1,298.6 - 1,338.2 - 1,156.9 - 1,204.9 - 1,207.4 - 1,335.6 -1,155.3 - 1,246.3
Antp
>0.2 A 0 0 0 0 0 0 0 0 0 0
0.1-0.2 A 0 0 0 0 0 1 0 0 0 0
<0.1 A 50 56 49 64 68 52 52 48 43 58
Total 50 56 49 64 68 53 52 48 43 58
VDW E? -325.9 -334.2 -323.9 -318.3 -317.2 -322.6 -312.7 -319.3 -314.3 -304.9
Conf. E? - 1,547.5 -1,511.3 -1,516.2 - 1,538.2 - 1,495.8 - 1,526.5 - 1,539.9 - 1,534.6 - 1,543.8 - 1,525.5
Msx-1
>0.2 A 0 0 0 0 0 0 0 0 0 0
0.1-0.2 A 0 0 1 0 0 0 0 0 0 0
<0.1 A 62 53 59 68 65 56 69 60 59 64
Total 62 53 60 68 58 56 69 60 59 64
VDW E? -292.3 - 292.4 - 290.0 -295.5 -291.1 -304.9 -280.4 -288.4 -295.6 -292.7
Conf. E? -1,186.5 -1,223.1 -1,230.4 - 1,201.5 -1,215.1 -1,223.6 - 1,197.5 - 1,281.0 - 1,248.7 - 1,208.0
aResidues 42-58 of Msx-I,
hVan der Waal (VDW E.) and conformational (Conf. E.) energies are defined in the text and are reported in units of kcal mol"
with backbone RMSD of 0.8 8, for residues 4-57; this value is dicted by SARMDlsearch with CONGEN; Msx-1 has 48% sequence
about the same as the backbone RMSD for residues 4-57 within identity with the engrailed homeodomain. Twenty-five conformers
the family of 16 NMR structures deposited in the Brookhaven were generated, and the1 0 with the lowest constraint violations and
Protein Data Bank (Billeter eta l., 1993). For each helix, as well as energies were selected to represent the predicted structure of Msx- 1.
for the helix-turn-helix motif in the structure, backbone RMSD A superposition of these 10 structures is shown in Figure 3 and sta-
values between the predicted structures and the average NMR tistics for the superpositions of backbone and heavy atoms for these
structure range from 0.3 to 0.8 8,. For all heavy atoms of the structures are alsop resented in Table 1. As before, the convergence
predicted Antp structure, the RMSD for residues 8-56 relative to of calculated structures was quite good; for residues 3-56, the
the average NMR structure is 1.78 ,, slightly higher than the RMSD RMSDs between each predicted structure and thea verage structure
of 1.3 8, observed for the heavy atoms of the 16 NMR structures are 0.3 8, for backbone atoms and 0.9 8, for all heavy atoms. VDW
relative to the average NMR coordinates. This deviation is under- energies for these homology-modeled structures of Msx- 1 range from
standable because some of the side chains of the NMR structure -4.7 to -5.1 kcal/mol-residue, which are comparable to the cor-
are constrained by interactions with DNA. In addition, some sur- responding values for Mata2 and Antp homeodomains and typical
face side chains interpreted as single conformers int he NMR of good-quality structures in the CONGEN forcef ield (Bassolino-
analysis may in fact adopt several isoenergetic conformations. Klimas et al., 1996; Tejero et al., 1996).
Homeodomain Msx-1 from mouse Discussion
From the predictions of atomic coordinates for Mata2a nd Antp, it
Homology modeling by SARMD and conformational
appears that homology modeling using SARMDhearch with CON-
search using CONGEN
GEN provides a robust and reliable method for predicting homol-
ogous homeodomain backbone chain folds. Having validated our This study provides a benchmark comparison of an automated
approach on two known homeodomain structures, the entire three- homology modeling technique using SARMD protocols (Bassolino-
dimensional structure of murine homeodomain Msx-1 was then pre- Klimas et al., 1996; Tejero et al., 1996) andc onformational search
Homology modeling of homeodomain Msx-1 with CONGEN 96 1
A
B
Fig. 2. A: Stereo diagram showing superpositions for backbone atomofs 10 homology-modeled structureosf Mata2 together with the
X-ray structure (red) from the Mata2-DNA complex (Wolberger a elt., 1991). B: Stereo diagram showing superpositionfso r backbone
atoms of 10 homology-modeled structures of Antp together with the average NMR structure (red) from the Antp-DNA complex
(Billeter et al., 1993).
algorithms (Bruccoleri & Karplus 1987, 1990; Bruccoleri et al., addition, in this work, the distance constraintsa re defined loosely
1988; Bruccoleri & Novotny, 1992) of the CONGEN computer in order to provide guidance to the information inherent in the
program.B y generating a family of structures, one obtains an CHAR" potentiale nergyf unction,w hilea llowing the mol-
estimate of the precision of the homology model given the as- ecule to adopt an energetically relaxed conformatioAn m. ore sim-
sumptions made regarding homology constraints and the validity plified approach excluding conformational energy would not provide
of the potential energy function. The accuraocf yo ur methodology one of the key features of our modeling method, i.e., low-energy
was determined by comparisons with experimental crystallo- side-chain packing that results from combination of homology
graphic and NMR structures. constraintsw itht hee nergyp otential.M oreover, the automated
For the homeodomains studied here, a small noufm cbareerf ully methodology outlined here is not restricted to a-helical proteins
selected a-helical constraints (e.g., enforcing all helical hydrogen and has also been used successfully to homology model &barrel
bonds) together with several properly chosen tertiary packing con- structures of Escherichia coli cold shock proteins (W. Feng, R.
straints would be enought o define the chain fold to -1 A. In our Tejero, & G.T. Montelione, unpubl. results).
approach, the bias of the modeler in selecting such constraints is Weh ave also compared backbone (N, C", C') RMS atomic
removed and replaced with a random selection ofm anym ore deviations between the experimental structuroefs Mata2 and Antp
homology constraints (in this case, -16 constraints per residue) and (1) the averaged coordinateso f the corresponding "predicted"
than the minimum number required to mimic the chain fold. In structures and (2)t he engrailed template structuref rom which the
962 H. Li et al.
Fig. 3. Stereo diagram showing superpositions of heavy atoms (backbone and side chain) for 10 homology-modeled structures of
murine Msx-1 .
homology constraints are derived (Table 1). ForM ata2, the back- culations. For this reason, theN -texminal polypeptide segmentso f
bone RMSDs of the X-ray crystal structure to the averaged coor- Matcu2, Antp, and Msxl are not packed against other portions of
dinates of the predicted structure (0.86 A) and to the engrailed the protein structure, and are not as tightly defined by the com-
template structure (0.85A ) are essentially identical. Similarly, for bined homology and energetic constrainatss other regionso f these
the Antp homeodomain, the backbone RMSDs of the averaged homology models (Fig. 2A,B).
NMR coordinates to the averaged coordinates of the predicted In the case of insertions and deletions, we also used the
structure (0.76 A) is only slightly lower than the corresponding conformational search algorithmosf CONGEN (Bruccoleri & Kar-
RMSD to the engrailed template structure( 0.79 A). In both cases, plus, 1987; Bruccoleri, 1993)to construct lowenergy polypeptide
the backbone conformations of the “predicted” structures are not loop conformations,a s described in the Materials and methodsI. n
significantly closer to the actual experimental structures than the Mata2, we defined residues 23, 23a, 23b, 23c2, 4,25,26, and 27
backbone conformationo f the engrailed template from which they as an unconstrained loop. Conformational searching (Bruccoleri,
are derived. This is because the extensive network of homology 1993) allowed us to explore the entire conformational space and
constraints used in this modeling ensures that the homology model find conformations with the lowest energies for this octapeptide
will be very similar to the template from which it is derived. loop. The search resultedi n two families of loop backbone struc-
On the other hand, the template structure is not an accurate tures (Fig. 2A), one of which includes the X-ray structure (Wol-
model of the target structures because is has different side-chain berger et al., 1991)T. his result indicates that there miginh tf act be
residues. Although potentially providing a reasonable predoifc tion more than one low-energy conformationf or this surface loop, al-
the backbone conformation, simple mutatoiof nth e engrailed struc- though the electron density apparently could be fit to a single
ture into the structures of Matcu2 or Antp (even if followed by backbone conformation.
energy minimization), would resiunl ts tructures that are physically
unreasonable and much higheri n energy than the structures gen-
Comparison with other approachesu sing satisfaction
erated by our procedure. By using simulated annealing methods,
of homology constraints
“predicted” structuresa re generated that satisfyb oth the homology
constraints and conformational energy constraints imposed byt he This section compares our hybrid SARMD search method using
CHAR” potential energy function. In addition, the SARMD CONGEN with related automated methodsfo r spatial satisfaction
procedure generatesa family of structures that provides a measure of homology constraints. The comparison focuses on the methods
of the uncertainty of the structure prediction, givetnh e assumption used for generating structures and the kinds ofh omology and
that “homologous atoms” will have similar relative positions in conformational energy information that are used. No efforts have
homologous structures. Although using the template structure as been made to datteo compare relative performance or reliabiliotfy
the predicted structure providesp roper orientation of the homeo- these various approaches.
domain helices, the challenge addressedb y the simulated anneal- Among the several general methods available for homology
ing methodology described here is to generate the solution space modeling, approaches that are most similar to ours include those
of low-energy conformations that are consistent with the homol- using DG calculations with homology constraint(sH avel & Snow,
ogy constraints. 1991;H avel,1 993;S rinivasan et al.,1 993;S udarsaname t al.,
Our approach generates the familoyf structures consistentw ith 1994), automated methods using restrained simulated annealing
a set of “homology constraints” and energy considerations, includ- with selected spatial constrainjtus dged to be important for the fold
ing steric packing effects. In the crystal structure of the engrailed and/or function (Brocklehurst& Perham, 1993), and the probabil-
homeodomain (Kissinger et al., 1990),t he N-terminal polypeptide ity density function approach employed by the program MOD-
segment interacts with the minor groove of duplex DNA. These ELLER( Sali & Blundell,1 993).L ikeo ura pproach,t heb asic
steric constraints are not included in our homology modeling cal- philosophy of these methods is to generate homology models au-
Homology modeling of homeodomain Msx-1 with CONGEN 963
tomatically with little or nou ser intervention. Significant dif- pdf s used in the method of Sali and Blundell (1993) are also
ferences among these various approaches and our method include derived from multiple homologous structures. In this study, we
(1) details of algorithmsu sedf or structure generation, (2)t he used a single template structure (i.e., the crystal structure of the
target functions and their relationships to homology structural in- engrailed protein), although it is also possible to define homology
formation, and (3) the useo f a set of superimposed conformers to constraints for SARMD from multiple template structures in the
interpret the consistency of the homology modeling prediction in future. In addition, other approaches (Brocklehurst & Perham, 1993;
different parts of the protein structure. Havel, 1993; Sali & Blundell, 1993) have included dihedral angle
The DG approaches (Havel & Snow, 1991; Havel, 1993; Srini- homology constraints, hydrogen bond homology constraints, andlor
vasan et al., 1993;S udarsaname t al., 1994)u se metric-matrix explicit water molecules in their homology modeling calculations,
methods to embed thes tructure into three-dimensional space, while all of which might be beneficial to include in our SARMDlsearch
minimizing homology and steric constraint violations. The method approach. Future homology modeling applications of SARMD/
of Sali and Blundell(l993) uses a variable-target function optimi- search with CONGEN could also benefit from more realistic de-
zation procedure with conjugate gradient minimization followed scriptions of electrostatics, solvent effects, counter-ion effects, and
by simulated annealing with molecular dynamics. It is unique in other interactions between the modeled protein and other mol-
that homology-constraint information is evaluated in the target ecules. In the case of these homedomains, the template engrailed
function in the form of a probability density function that describes structure comes from the structure of a complex between engrailed
various homologous features to be restrained. This continuous prob- and duplex DNA (Kissinger et al., 1990); modeling of the entire
ability density function differs from the flat-bottomed upperllower protein-DNA complex would provide better definition and more
bound constraints used by DG and restrained MD methods. Pre- correct predictions of conformations of interfacial side chains be-
viously described pure SARMD approaches have used a knowlege- cause the interactions with DNA would act to reduce the number
based selection of a relatively small number of homology constraints of isoenergetic side-chain conformations.
(Brocklehurst & Perham, 1993) or have used only homology con- Despite these possible shortcomings, the SARMDlsearch method
straints between backbone atoms and minimal conformational en- using CONGEN with homology constraints provides a robust and
ergy information (Engh eta l., 1990; Fujiyoshi-Yoneda et al., 1991). reliable approach for homology modeling of small single-chain
The procedures of Havel (1993) and Sali andB lundell(l993) also proteins. The predicted homeodomain structures described here
use molecular dynamics to refine structures that are first generated exhibit precision and accuracy similar to that available from medium-
by DG or variable-target function gradient minimization methods. resolution solution NMR structures. We have also had good suc-
Although the energies of homology models produced by these cess using this approach for homology modeling &barrel protein
methods are improved by the restrained MD calculations (Havel, structures (W. Feng, R. Tejero, & G.T. Montelione, unpubl. re-
1993; Sali & Blundell, 1993), these conformers may be trapped in sults). Overall, the results presented in this paper provide a good
local energy minima that cannot be overcome without a suffi- example of the utility of homology modeling using spatial homol-
ciently long high-temperature simulated annealing procedure. As a ogy constraints (Havel & Snow, 1991) applied to an important
result, the family of structural models may retain energetically class of transcription-regulating domains. The family of modeled
unfavorable features, and may not sample adequately the range of structures of the Msx-1 homeodomain has been deposited in the
conformations that are consistent with the assumptions of the ho- Brookhaven Protein Data Bank for use in future comparisons of
mology modeling. the predicted structure with NMR and X-ray structures once they
In our approach, a subset of homology constraints are selected are available or in molecular replacement approaches for deter-
randomly without user intervention from the “homologous dis- mining its structure by these experimental methods.
tances’’ in the template structure. SARMD with CONGEN is then
used right from the start to generate a family of structures that
Predicted structure of murine Msx-1
satisfies both homology and energetic constraints while maximiz-
ing the sampling of the solution space. The convergence rate of Figure 4 shows a ribbon model depicting the average, energy-
these SARMD calculations is enhanced by using the flexible dis- minimized predicted structure of Msx- 1. The structure consists of
tance restraint function of CONGEN, as described by Bassolino- three a-helices: helix I (residues 10-22), helix I1 (residues 28-37),
Klimas et al. (1996). Efficient sampling of the solution space of and helix I11 (residues 42-58). The helix content in the predicted
side-chain and loop conformations is provided both by the simu- structure is -67% and is consistent with estimates from analysis of
lated annealing protocol and by the subsequent directed confor- far-UV CD spectra of -70% helix content in aqueous solution at
mational search. Like the procedures of Brocklehurst and Perham pH 7.0 and temperature of 22 “C (Shang et al., 1994). Helices I and
(1993), Havel (1993), and Sali and Blundell (1993), our method I1 are antiparallel, whereas helix I11 is perpendicular to the other
generates families of homology-modeled structures that are used to two helices. The N-terminal polypeptide segment (residues 1-9) is
evaluate the consistency of the structure prediction in different in an extended conformation directed away from the three helical
parts of the model. bundle. The conformation of the segment, modeled from the bound-
state conformation of the engrailed-DNA complex, is not defined
uniquely in our simulated annealing calculations (Fig. 3) and may,
Further improvements in the CONGEN
in fact, be quite flexible in Msx-1 when it is not bound to DNA.
SARMDhearch approach
Dihedral angle order parameters S(4) and S($) (Hyberts et al.,
Each of the methods described above have their respective strengths 1992) and Ramachandranp lotsc omputed for thisf amily of
and weaknesses, and it may be advantageous to incorporate ideas homology-modeled Msx-1 structures are shown in Figures 5 and 6,
described by other workers into our CONGEN homology model- respectively. Only residues of well-defined regions of the struc-
ing calculations. The work described by Havel (1993) used mul- ture, defined as those with S(4) + S($) > 1.5, are shown on the
tiple template structures for generating homology constraints. The Ramachandran plots. All of these well-defined residues in all of the
964 H. Li et al.
models are in low-energy regions of the Ramachandran plot the Msx-1-DNA complex. The positively charged side chains of
(Fig. 6). The only residue on the right side (i.e., 9 > 0) of the residuesa rgininea ndl ysine( shown in bluei nF ig. 4) on the
Ramachandran map is residue Ser 39, with 939FJ 50" and 1,439 FJ N-terminal arm ando nt hes urface ofh elix IU ofM sx-1 are
60". Because residueL eu 38 has 9, I,4 values around -90" and O", available to interact with DNA in the minor and major grooves,
respectively, the polypeptides egmentS er 37-Leu3 8-Ser3 9- respectively. The highly charged N-terminal arm and third helix
Leu 40 forms a typNe /3 bend (Lewise t al., 1973) int he predicted can provide high-affinity electrostatic interactions with the poly-
structure of Msx-1. Position 39 is occupied by a glycine residuein phosphate backbone of DNA and the complex hydrophobic and
the template-engrailed protein structure (Kissingere t al., 1990). It ionic nature of this surface can provide specificity to dock the
will be interesting to learn if this prediction of a positive value protein at its correct binding position.
for residue Ser 39 is borne out by future NMR and/or crystal
structures of Msx-1 ori f it is an artifact of our homology modeling
A conserved network of salt bridges in the homeodomain
procedure.
chain fold
As illustrated in the stereodiagram of Figure 3, most of the
conformations of side chains in the interior of Msx-1 are tightly A network of ionic interactions preserves the interhelical packing
restrained by the combination of distance and energy constraints, of the homeodomain chain fold in Msx-1. In analyzing our pre-
whereas surface side-chain conformationarse more poorly defined dicted structureo f Msx- 1, we observae d n etwork osf everal surface
in the homology model, as would be expected from the energy salt bridges (Table 3) that appears to be conserved in other ho-
potential. Similar distributions of side-chain conformations were meodomain structures. These are illustrated on the three-dimensional
observed in the ensembles of homology-modeled Mata2 and Antp structure of Msx-1 inF igure 7. This network of surface salt bridges
structures. This kind of graphical information on the precisiono f in the homology-modeled structuroef Msx-1 results primarily from
side-chain conformations in these homology models is available constraints imposed by its predicted chain fold and in part as a
only because the modealsr e represented as ensembles of structures
and because MD and search methods are used to sample exten-
sively the solution space consistent wtithhe homology and energy
Table 3. Postulated salt bridges in homology-modeled
constraints.
structures of the Msx-1 homeodomain
Conserved residues, located mainly at the helical interfaces of
the three-helical bundle fold, form the hydrophobic core of the
Frequen cyC ation Anion of occurrencea
Msx-1 homeodomain. This core includes many highly conserved
hydrophobic residues (Scotte t al., 1989) with an invariant Trp 48 Glu 30 C'OO- Lys 19 NcHs+ 1. o
at the center (Fig. 4). Near-UV CD spectroscopy and tryptophan Glu 30 C'OO- LYS 23 N~H,+ 1. o
fluorescence quenching experiments on Msx-1 demonstrate that Glu 33 C'OO- LYS 19 N ~ H ~ + 0.7
residue Trp 48 is, in fact, buried in the hydrophobic core of the Glu 42 C'OO- 0.9
natives tructure of Msx-1( Shang et al., 1994), as iti s in our
Glu 17 C600- 0.9
predicted homology models.
Glu 17 C'OO- 0.7
In the engrailed, Antp, and the Mata2 homeodomain-DNA Glu 17 C'OO- 0.6
complexes, residues in helix IU contact the major groove, resi-
dues in the N-terminal arm contact the minor grove, and resi-
aCriteria used to identify postulated salt bridges are outlined in Mate-
dues inh elix I1 contact the phosphateb ackbone (Kissinger
rials and methods. The frequency of occurrence is the fraction of the 10
eta l.1, 990;W olbergere t al., 1991;B illetere t al., 1993;L i homology-modeled structures containing the postulatesda lt bridge. Poten-
et al., 1995). It is likely that these contacts will also be found in tial salt bridges with frequency of occurrence <0.6 are not reported.
Homology modeling of homeodomain Mss-I u‘ith CONGEN 965
1. 00 The side chains of residues Arg 31 and Glu 42 also form a salt
bridge between helix I1 and helix Ill in 9 of IO modeled Msx-I
0.90
structures (Fig. 7; Table 3). Homologous 31-42 salt bridges are
0.80 also present in the engrailed and Antp structures. Although Mata2
- 0.70 does not have a salt bridge between 31 and 42. it does have the
.- 0.60 “inverse” Glu 32-Arg 42 salt bridge. These compensating amino
s acid substitutions in Mata2 and the strong conservation of the
5* 0.50 31/32-42 ionic interaction suggest that this salt bridge between
0.40 helix I1 and helix 111 has an important structural andlor functional
role.
0.30
Another frequently observed salt bridge is Glu 17-Arg 52,
0.20 which connects helix I and helix 111; it is found in 7 of IO
0.10 homology-modeled structures (Table 3) and is conserved in the
Antp and Mata2 homeodomain structures. In the TTF-1 homeo-
0.00
5 10 15 20 25 30 35 40 45 50 55 domain, a salt bridge between Glu 17 and His 52 is also pre-
Residue Number dicted (Fogolari et al., 1993). However, a homologous 17-52
salt bridge is not present in engrailed homeodomain. where both
residues are lysines.
In a covariance analysis of 60 homeodomain protein sequences,
Clarke (1995) has observed that residue pairs 19/30. 3 1/42, and
1. oo 17/52 are among the most strongly correlated covariant residue
0.90 pairs. Moreover, the nature of these covariances generally func-
tions to preserve salt bridges at these three surface sites (Clarke,
0.80
1995). The energetic basis of the requirement for these surface salt
0.70 bridges is not yet certain. Residues Arg 3 1 of engrailed and Arg 42
0.60 of Mata2 each make contacts with corresponding phosphate atoms
in these homeodomain-DNA complexes, suggesting that there may
0.50
be a functional constraint requiring a basic residue at one (but not
0.40 both) of these two sites (Clarke, 1995). However, residue pairs
0.30 19/30 and 17/52 are distant from the DNA-binding site and do not
interact with DNA in the structures of complexes.
0.20
Considering that solvent-accessible surface ionic interactions
0. 10 will be suppressed by the high dielectric of bulk water and by
0.00 counterion salts, it is unlikely that these salt bridges contribute
s IO 1s x 2s 30 7s w 4s sn ss significantly to the thermodynamic stability of the free homeo-
domain. This suggests that the strong conservation of these salt
Residue Number
bridges reflects a requirement to balance charges because these
Fig. 5. Backbone dihedral angle order parameters (Hyhens et al.. 1992) faces of the homeodomain (Fig. 7) are buried in protein-protein
for the ensemble of hornolopy-modeled structures of Msx- I.
and protein-nucleic acid interactions that occur in the formation of
functional transcription complexes. For example, in the crystal
structure of the MataI/Mata2 homeodomain heterodimer bound to
consequence of the CHARMM force field, which drives oppositely DNA, the Lys 23-GIu 30 salt bridge of Matal is buried at the
charged atoms toward one another. interface between Matal and the carboxy-terminal tail of Mata2
In the Msx-l homeodomain structure, the ionizable side-chain (Li et al.. 1995). Because it is energetically very unfavorable to
atoms of residues Lys 19, Lys 23, and Glu 30 are spatially very bury an unbalanced charge in a protein-protein interface, the pres-
close together. forming salt bridges between helix I and helix I1 in ence or absence of these generally conserved salt bridges on the
all the IO modeled structures. Residues Lys 19 and Lys 23 are four surfaces of different homeodomains potentially can play an im-
residues apart in the C-terminus of helix-I and their side chains are portant role in modulating the energetics of interactions among
on the same side of the helix and adjacent to one another in space. homeodomains and between homeodomains and other proteins
Salt bridges between positions 19 andlor 23 in helix I and position andlor DNA in transcription complexes.
30 in helix I1 are also present in other homeodomains, including
the Glu 19-Arg 30 salt bridge in both the engrailed and Antp
Implications of the Mm-I structure for protein engineering
homeodomains. the Glu 22-Arg 30 salt bridge in engrailed, and
the “inverse” Arg 19-Glu 30 salt bridge in the predicted structure Using structural information derived from this homology model of
of thyroid transcription factor 1 (TTF-1) (Fogolari et al., 1993). Msx-I, we have designed and synthesized homeodomain Ala-Msx
Although a homologous salt bridge is not present in Mata2, in the (Shang et al., 1994). which contains most homeodomain consensus
structure of the related yeast, Mat-a1 homeodomain includes a salt residues (Scott et al., 1989), but replaces most nonconsensus res-
bridge between residues Lys 23 and Glu 30 (Li et al., 1995). The idues with alanine. Ala-Msx contains 46% Ala and binds to duplex
double mutant [Lys 19-GIu, Glu 30-Arg] Msx-I, with an “inverse” DNA with a sequence-specific dissociation constant K,, < 100 nM
salt bridge, exhibits a native-like DNA-binding affinity (Isaac (Shang et al., 1994). This alanine-substituted analogue of Msx-I,
et al.. 1995). which retains the hydrophobic core defined by the homology model