| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
Biophys J, June 2002, p. 3063-3071, Vol. 82, No. 6
and
*Cambridge Centre for Molecular Recognition, Department of
Biochemistry, University of Cambridge, Cambridge CB2 1GA, United
Kingdom;
Department of Biochemistry, University of
Oxford, Oxford OX1 3QU, United Kingdom; and
The
Alexander Silberman Institute of Life Sciences, Department of
Biological Chemistry, The Hebrew University of Jerusalem,
Givat-Ram, Jerusalem 91904, Israel
| |
ABSTRACT |
|---|
|
|
|---|
Molecular interactions between transmembrane
-helices
can be explored using global searching molecular dynamics simulations (GSMDS), a method that produces a group of probable low energy structures. We have shown previously that the correct model in various
homooligomers is always located at the bottom of one of various
possible energy basins. Unfortunately, the correct model is not
necessarily the one with the lowest energy according to the
computational protocol, which has resulted in overlooking of this
parameter in favor of experimental data. In an attempt to use energetic
considerations in the aforementioned analysis, we used global searching
molecular dynamics simulations on three homooligomers of different
sizes, the structures of which are known. As expected, our results show
that even when the conformational space searched includes the correct
structure, taking together simulations using both left and right
handednesses, the correct model does not necessarily have the lowest
energy. However, for the models derived from the simulation that uses
the correct handedness, the lowest energy model is always at, or very
close to, the correct orientation. We hypothesize that this should also
be true when simulations are performed using homologous sequences, and
consequently lowest energy models with the right handedness should
produce a cluster around a certain orientation. In contrast, using the wrong handedness the lowest energy structures for each sequence should
appear at many different orientations. The rationale behind this is
that, although more than one energy basin may exist, basins that do not
contain the correct model will shift or disappear because they will be
destabilized by at least one conservative (i.e. silent) mutation,
whereas the basin containing the correct model will remain. This not
only allows one to point to the possible handedness of the bundle, but
can be used to overcome ambiguities arising from the use of homologous
sequences in the analysis of global searching molecular dynamics
simulations. In addition, because clustering of lowest energy models
arising from homologous sequences only happens when the estimation of
the helix tilt is correct, it may provide a validation for the helix
tilt estimate.
| |
INTRODUCTION |
|---|
|
|
|---|
Prediction of transmembrane domain structure is
facilitated by its tendency to adopt (in most cases) an
-helical
fold, limiting the number of possible conformations. Brunger and
co-workers (Treutlein et al., 1992
; Adams et al., 1995
) have developed
a procedure to explore transmembrane helix interactions based on global
searching molecular dynamics simulations. In this method, multiple
symmetric bundles of helices are constructed, each differing from the
other by the rotation of the helices about their axes (Fig.
1). These are then used as starting
positions for molecular dynamics and energy minimization protocols. The
output structures from these simulations are compared and grouped into
clusters that contain similar structures. An average of the structures
forming a cluster represents a model with characteristic interhelical
interactions and helix tilt.
|
The correct model is selected amongst the several different clusters,
based on existing experimental data, either from mutagenesis (Lemmon et
al., 1992a
,b
; Arkin et al., 1994
) or orientational data from site
specific infrared dichroism used as spatial restraints (Kukol et al.,
1999
; Torres et al., 2000
). Alternatively, we have also used a purely
computational approach where simulations are performed on close
sequence variants that are likely to share the same structure (Briggs
et al., 2001
).
The value of energy as a discriminating tool is usually overlooked
during the analysis due to the fact that the correct model is not
necessarily the one with lowest energy according to the computational
protocol. Recently however, in an exhaustive evaluation of the energy
of energy-minimized helical bundles as a function of their helix tilt
and rotational orientation, we obtained a picture of the energy
landscape of
-helical bundles as a function of their interhelical
interactions (Torres et al., 2001
). The conclusion was that, although
more than one energy basin was present, the model that was in agreement
with experimental data, (i.e., the correct model), was always located
at the bottom of one of these basins. In some cases, only two basins
were found, one for each handedness, i.e., left (
> 0) or
right (
< 0). In other cases, more basins were found, although
the search for the correct model is limited to a slice of the energy
surface due to the helix tilt constraint.
An example of this has been shown previously for phospholamban (Torres
et al., 2000
) and in the case of one particular variant of the M2
H+ of Influenza A (Kukol et al., 1999
). When the helix tilt
was restrained to the experimentally determined value by Fourier
transform infrared, site-specific infrared dichroism, the
correct model was also the lowest energy model when the simulation was
performed using the correct handedness (known a priori). This is
despite that the lowest energy model with the opposite, wrong
handedness, could have had even lower energy and could have been
equally compatible with experimental data.
Herein, we have attempted to discriminate between the energy basin that
contains the correct model and those that do not by using homologous
sequences. We have shown previously (Briggs et al., 2001
) that silent
amino-acid substitutions, present in close homologues or unveiled by
mutagenesis studies, do not affect the stability of the native
structure during the simulation, although they destabilize some
of the non-native structures also found. Global searches carried
out on enough variants produce a set of similar structures that appear
in all of the searches, forming a "complete set." We argued that
the structure represented by this "complete set" should be the
native one, because it had not been destabilized by any of the
conservative/silent mutations.
Similarly, we expect that a particular basin that does not contain the correct model, regardless of how low its energy is, will shift or disappear in at least one of the simulations (of the different homologues), whereas the basin containing the correct model will remain in all simulations.
When helix tilt and handedness are correct, therefore, upon use of
different homologous sequences, the lowest energy structures should
cluster around a particular orientation. Conversely, lowest energy
structures should spread over many different orientations if either of
these two parameters is incorrect. We show this to be true for the
glycophorin A homodimer for which a structure has been reported
(MacKenzie et al., 1997
).
The use of the model's energy values using close homologues is an
essential tool to resolve ambiguities that arise in the new approach
reported recently (Briggs et al., 2001
) that also uses close variants.
For example, if the number of close homologues available is not
sufficient, more than one model may form a "complete set."
Alternatively, if the correct model is not found amongst the candidate
structures and the number of sequences available is too small, a wrong
model may still be found to persist in all of the simulations tested.
Finally, when mutagenesis data are used, results are not always clear
cut, because function is not always the parameter under monitoring,
thereby creating ambiguity when evaluating these results. This can have
the effect of precluding the formation of any "complete set," even
when the models arising from global searching molecular dynamics
contain the correct structure. In all of these cases, it should be
possible to select the correct model by looking at the energy of each
model and determining where clustering of lowest energy models occurs
around certain orientations.
We provide examples of, and overcome, this ambiguity with a
homotetramer and a homopentamer. In one case, we use the
-helical (Duff et al., 1992
; Kukol et al., 1999
) transmembrane domain of M2, a
protein from the Influenza A virus that forms homotetrameric H+-selective ion channels. For M2 global searching
molecular dynamics simulations (GSMDS), we have used in the simulations
close homologues that are likely to have the same native structure.
However, we have deliberately included some variants that are known to
confer resistance to the antiviral drug amantadine (Hay et al., 1985
). Since it has been shown that the transmembrane part of the protein is
necessary and sufficient for the functional and structural properties
of the protein (Duff and Ashley, 1992
; Duff et al., 1994
), it is
possible that the variants of M2 from amantadine-resistant viruses
possess structural properties that are somewhat different from the
amantadine-sensitive variants. This is incompatible with the main
premise in the analysis of GSMDS by using close variants in which all
sequences should share the same structure.
The other system used, the homopentameric phospholamban (Plb), a
52-residue protein located in the cardiac sarcoplasmic reticulum, contains a hydrophobic domain (residues 26-52) that forms a
left-handed
-helical bundle (Arkin et al., 1994
; for recent reviews
see Arkin et al., 1997
; Simmerman and Jones, 1998
). No evolutionary
data is available for Plb, although a wealth of information exists concerning the residues that are important in maintaining the pentameric structure, mainly based on mutagenesis data (Fujii et al.,
1989
; Arkin et al., 1994
; Simmerman et al., 1996
; Kimura et al., 1997
).
Ambiguity here is introduced by the results of mutagenesis studies,
which did not define conservative mutations as those preserving
function but those that did not prevent pentamerization. Thus, some of
the conservative changes reported may had affected function.
The correct orientations and handedness for M2 and Plb are known. For
M2, the rotational orientation and the helix tilt have been
experimentally determined by two independent methods: solid-state nuclear magnetic resonance (NMR) (Kovacs and Cross, 1997
) and site-directed infrared dichroism (Kukol et al., 1999
). For Plb, two
models were capable of "accommodating" the mutagenesis data (Arkin
et al., 1994
; Simmerman et al., 1996
), although the leucine zipper
model (Simmerman et al., 1996
) has been confirmed by labeling experiments performed in sodium dodecyl sulphate (SDS) (Karim et al.,
1998
) and a site-directed infrared dichroism study (Torres et al.,
2000
). Herein, new ways in which examination of energy values can
contribute to the analysis of results of GSMDS are presented, using the
aforementioned examples.
| |
MATERIALS AND METHODS |
|---|
|
|
|---|
Global search molecular dynamics
Glycophorin
The protocols for the simulations and the homologues used for glycophorin homodimers were as described previously (Briggs et al., 2001M2
Simulations with M2 were performed using the transmembrane sequence from residues 26 to 43. M2 variants (Fig. 2) with at least 60% identity were obtained searching the OWL data base (Bleasby et al., 1994
|
Plb
For Plb, conservative mutations in the wild-type transmembrane sequence (residues 35-52) FCLILICLLLICIIVMLL (see text) were used from either Engelman and co-workers (Arkin et al., 1994GSMDS protocol
All calculations were performed using parallel crystallography and NMR system, the parallel-processing version of the Crystallography and NMR System (CNS Version 0.3) (Brunger et al., 1998
-carbon root mean square
deviation (RMSD) from any other structure within that cluster.
Some clusters therefore overlap, and output structures may be members
of more than one cluster. The structures belonging to each cluster were
averaged and subjected to a further simulated annealing protocol. This
final structure was taken as the representative of the cluster.
Finally, the side chain rotamer angles are free to rotate during the
course of the molecular dynamics run as well as the energy minimization
protocol. The initial values are those found most prevalent in the
rotamer library (Ponder and Richards, 1987Analysis of the simulations
The results from the global searching molecular dynamics
simulations were represented graphically by plotting each cluster representative as a function of two parameters, the helix tilt
and
the rotational orientation
, as described previously (Briggs et al.,
2001
) and shown schematically in Fig. 1. The tilt angle of the model,
, was taken as the average of the angles between each helix axis and
the bundle axis, which is coincident with the normal to the bilayer and
was calculated by CHI (Adams et al., 1995
, 1996
). The helical axis is a
vector with starting and end points above and below a defined residue,
where the points correspond to the geometric mean of the coordinates of
the five
carbons N-terminal and the five
carbons C-terminal to
the defined residue.
The rotational orientation
, is defined relative to an arbitrarily
chosen, specified residue. The angle
is defined by the angle
between a vector perpendicular to the helix axis, oriented towards the
middle of the peptidic C==O bond of the residue, and a plane that
contains both the helical axis and the normal to the bilayer (Fig. 1).
This angle is 0° when the residue is located in the direction of the
tilt (Arkin et al., 1997
). The structures identified were plotted
against the
angle of residue 83 for glycophorin A, 35 for M2, and
42 for Plb.
Precise comparisons between similar clusters obtained from different
variants were made by calculating the RMSD between their
-carbon
backbones. In the simulations, the handedness of the bundle is
indicated by the helix tilt sign, positive or negative, which
corresponds to left and right handed bundles, respectively. Note that
in a helix dimer, the sum of both tilt angles is equal to the commonly
used helix crossing angle,
. Such a relationship does not
necessarily exist for higher order oligomers.
M2 simulations in which the helix tilt was restrained were performed as
described previously (Torres et al., 2000
). Specifically, the angle
between the vectors connecting every C
of residue i and C
of residue i + 7 and
the z axis was restrained to 31.5°, which is the
experimental value obtained for the helix tilt using site-directed
dichroism (Kukol et al., 1999
).
Energy calculation
The energies calculated correspond to the total energy of the
system, including both bonded (e.g. bond, angle, dihedral, improper) and nonbonded (i.e. Van der Waals, electrostatic) terms (Brünger et al., 1998
). No account is made for contributions from the lipids to
the energetics of the system (Torres et al., 2001
).
| |
RESULTS AND DISCUSSION |
|---|
|
|
|---|
Glycophorin homodimer
The results of the simulations for various homologous sequences of
GpA are shown in Fig. 3. The lowest
energy structures for each of the sequences at left handed
configurations (Fig. 3 a) appear at various orientations
(Fig. 3 b). In contrast, the right handed structures (Fig.
3 c) show a clustering of the lowest energy structures
around
~
80° (Fig. 3 d), although the lowest
energy structures from pig and gibbon are 25° and 50° apart from
this orientation. These results are consistent with a previous approach in which a structure at
=
80° and
=
10°
(right handed) was identified at the bottom of an energy basin (Torres
et al., 2001
). This model is also indistinguishable to that obtained
computationally using the effect of conservative substitutions
(Briggs et al., 2001
) and experimentally by solution NMR in detergent
micelles (MacKenzie et al., 1997
).
|
M2
No restraints
Fig. 4 depicts the result of the simulations using both left- and right-handed configurations of TM-M2, without restraining the helix tilt. Fig. 4 a shows that no structure could be found with a helix tilt
of 31.5° or
35 = ~
100° compatible with that obtained
experimentally (see vertical broken line (Kukol et al., 1999
35 = 25° (see gray squares) is present
in all the simulations carried out with the amantadine sensitive
variants (s1-s5) and also with two of the amantadine resistant
variants, r1 and r3. However, even relaxing the clustering parameters
so that clusters should contain a minimum of six structures instead of
10, r2 failed to produce a model at
35 = ~25°
(not shown). This is a clear example of how ambiguity can easily arise when using the approach reported previously (Briggs et al., 2001
35 ~ 25° might have been taken as correct. Alternatively, one may have thought that r2 possesses a structure that is somewhat different and
should not have been included in the simulation in the first place,
which again would justify the model at
35 ~ 25°
as correct for the other sequences. Thus, due to a limited number of
sequences or to uncertainty regarding the structural identity of the
variants used, a wrong model is obtained. We show here that this can be overcome by simply examining the relative energy of the models generated.
|
35 = 25° (see gray square)
are not the lowest energy structures within a particular sequence for
any of the simulations. On the contrary, it has almost the highest
energy in each of the simulations, suggesting that this model is wrong. In addition, as no other model forms a "complete set," we can infer
that the correct structure is not present amongst the low energy
clusters generated by the global searching molecular dynamics protocol.
When the helix tilt is incorrect, it is not possible to generate a
"complete set" where the components are all lowest energy structures for a certain handedness. Therefore, the lack of results described above may indirectly help to identify when a helix tilt is
incorrect and can be used as a feedback on the validity of experimental measurements.
Use of homologues restraining the helix tilt
Fig. 5 shows the results of the global search molecular dynamics protocol for different M2 variants when the helix was restrained to 31.5°, the experimentally determined value (Kukol et al., 1999
~ 31°,
35 ~
100°, and
~ 31°,
35 ~ 100° (both left handed) clusters are found
for all the homologues that are amantadine sensitive (first five
variants). No structure within either of the "complete sets," at
35 ~
100° and
35 ~ 100°, differs from any other in the set by more than 0.8 Å C
RMSD. Therefore, these two orientations are in principle possible according to this computational approach. One of these models, however
(see gray rectangles), is equivalent to the model previously determined
using site-directed infrared dichroism (Kukol et al., 1999
|
125°, see arrow) that is also the lowest energy
cluster. In fact, the
-carbon RMSD between this cluster and any
other model belonging to the "complete set" at
35 ~
100° is less than 1.4 Å, which shows
that this structure is in fact very similar to those corresponding to
the "complete set."
|
|
Phospholamban
For Plb, as no sequences are available from evolutionary data,
variants resulting from mutagenesis data (Arkin et al., 1994
; Simmerman
et al., 1996
) were used. Fig. 8 shows
that when no restraints were used in the left-handed configuration no
"complete set" was found. In contrast to the simulations on M2 in
the absence of restraints, however, a number of structures are
identified in the vicinity of the experimentally determined helix tilt
(11°, Torres et al., 2000
), so it is unlikely that the correct
structure falls outside the conformational space searched.
|
Alternatively, some of the sequence variants used may contain some ambiguity regarding their ability to preserve the native structure, because the mutants were assayed for their ability to form pentamers and may not necessarily preserve function. Or in other words, they may pentamerize in a slightly different way than in the native sequence.
Clearly, not all mutations described in the literature can be
conservative/silent. For example, the now accepted leucine zipper model
(Simmerman et al., 1996
), equivalent to
35 ~
55° in Fig. 8, displays C41 facing the lipidic environment, whereas
C36 is located in the interhelical region and C46 is facing the lumen of the pore. The other alternative model, also in agreement with mutagenesis data (Arkin et al., 1994
), at
35 ~
10°, displays C41 and C46 in the interhelical contacts and C36
facing the pore. Only these two models, (see gray squares in Fig.
8) were present in all of the simulations but one. It is perhaps not
surprising, therefore, that the results that involve C41V and C36F are
similar. Whereas, the C41V simulation does not display a model at
~
10°, C36F did not produce a model at
~
55°.
We note that as with other systems simulated, although the
value
for a given residue may change somewhat, the structure can still be
very similar. For example, the model for L43F at
~
35°
belongs in fact to the set at
~
55° because is at less
than 1 Å C
RMSD from any structure within that set.
As in M2, however, the model with lowest energy for each of the
simulations (see lower panel in Fig. 8 at ~
55°) is also in
agreement with previous experimental data (Torres et al., 2000
). Therefore, our data suggest that C36F may preserve the pentameric properties of the protein but may not be functional due to slight changes in the pentameric structure.
| |
CONCLUSION |
|---|
|
|
|---|
In summary, the approach previously described (Briggs et al.,
2001
) is limited by the inability of the search protocol to include the
correct structure among its output and by ambiguities present in the
input sequences, either from evolutionary or mutagenesis data. We have
shown that the incorporation of a simple restraint such as the helix
tilt is sufficient to allow the native structure to be found. Once the
helix tilt is correct, this in turn allows one to discriminate between
the models according to their energy.
Examples of ambiguities have been provided and overcome by the use of energy values as an analysis tool. The main problems facing these methods are the imperfections in the simulation protocol, i.e., the force field does not reflect accurately the properties of the system and the simulations are performed in vacuo. We stress that it is because of the imperfections of the system referred to above, which may distort somewhat the energy landscape, that lowest energy models form a "loose" cluster, where orientations for some sequences can be as far as 50° from the main cluster.
It is likely that the contribution of lipid-protein interactions in the different models is similar, and therefore the correct model is also the one with lowest energy, even in the absence of lipids. This is important in three ways. It allows us to know when the correct model has not been found amongst the candidate structures, making necessary the incorporation of experimental restraints such as helix tilt. Also, it allows for the selection of the correct model when 1) more than one model form a "complete set," in a situation where the number of sequences is limited, or 2) when no "complete sets" are found due to ambiguous results from site-directed mutagenesis. Restraints such as the helix tilt can be relatively easy to obtain for homooligomers using techniques such as Fourier transfer infrared, and solid-state NMR or from low resolution structures from electron microscopy, although in the latter case, information on relative helix register, which can be obtained from electron paramagnetic resonance, is also necessary.
| |
ACKNOWLEDGMENTS |
|---|
This work was supported by grants from the Wellcome Trust and the Biotechnology and Biological Sciences Research Council to I.T.A.
| |
FOOTNOTES |
|---|
.
Address reprint requests to Isaiah T. Arkin, The Alexander Silberman Institute of Life Sciences, Department of Biological Chemistry, The Hebrew University of Jerusalem, Givat-Ram, Jerusalem 91904, Israel. Tel.: 972-0-2-658-4329; Fax: 972-0-2-658-4329; E-mail: arkin{at}cc.huji.ac.il.
Submitted April 25, 2001, and accepted for publication February 12, 2002.
| |
REFERENCES |
|---|
|
|
|---|
Biophys J, June 2002, p. 3063-3071, Vol. 82, No. 6
© 2002 by the Biophysical Society 0006-3495/02/06/3063/09 $2.00
This article has been cited by other articles:
![]() |
K. Parthasarathy, L. Ng, X. Lin, D. X. Liu, K. Pervushin, X. Gong, and J. Torres Structural Flexibility of the Pentameric SARS Coronavirus Envelope Protein Ion Channel Biophys. J., September 15, 2008; 95(6): L39 - L41. [Abstract] [Full Text] [PDF] |
||||
![]() |
S. W. Gan, L. Ng, X. Lin, X. Gong, and J. Torres Structure and ion channel activity of the human respiratory syncytial virus (hRSV) small hydrophobic protein transmembrane domain Protein Sci., May 1, 2008; 17(5): 813 - 820. [Abstract] [Full Text] [PDF] |
||||
![]() |
L. Bu, W. Im, and C. L. Brooks III Membrane Assembly of Simple Helix Homo-Oligomers Studied via Molecular Dynamics Simulations Biophys. J., February 1, 2007; 92(3): 854 - 863. [Abstract] [Full Text] [PDF] |
||||
![]() |
J. Torres, J. Wang, K. Parthasarathy, and D. X. Liu The Transmembrane Oligomers of Coronavirus Protein E Biophys. J., February 1, 2005; 88(2): 1283 - 1290. [Abstract] [Full Text] [PDF] |
||||
![]() |
W. Im, M. Feig, and C. L. Brooks III An Implicit Membrane Generalized Born Theory for the Study of Structure, Stability, and Interactions of Membrane Proteins Biophys. J., November 1, 2003; 85(5): 2900 - 2918. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |