J. Mol. Biol.
(1979) 132, 243-368
Dynamics
of Ligand D. A.
CASE?
Binding M.
AND
to Heme Proteins KARPLUS
Department of Chemistry Harvard University Cambridge, Mass. 02138, U.S.A. (Received
7 December 1978)
The dynamics of CO or 0, binding to myoglobin are examined theoretically with particular emphasis on the effect of the globin on the barriers along possible reaction paths. By use of a diabatic representation, in which a model ligand moves in the potential of a rigid protein, two paths in or out of the heme pocket are found in myoglobin ; one is by the E helix near the distal histidine and the other near the intersection of the B, D and E helices. Classical trajectories calculated for a photoligand in the diabatic approximation yield complicated motion; the dissociated ligand tends to spend considerable time undergoing multiple collisions with the protein matrix in the wells between barriers. To investigate the height of the barriers, the adiabatic limit is used ; that, is, the protein is allowed to relax in the presence of a perturbation due to the ligand-protein interaction. A detailed description is obtained for the displacements of residues required for the ligand to enter or leave the globin. The barrier heights calculated in this manner are found to be in qualitative agreement with observed values. Some differences between myoglobin and hemoglobin are described. Although the present calculations make use of a simplified model of the ligandmyoglobin system, they provide a useful first step in the analysis of the effect of the globin matrix on the binding process.
1. Introduction A molecular description of substrate and ligand binding by enzymes and other proteins can for convenience be divided into two parts. The first deals with the initial and final geometries of the species involved, while the second is concerned with the reaction path and the dynamics of the binding process. Transport and storage proteins, such as hemoglobin and myoglobin, are favorable cases for study since both the ligated and unligated species are stable and there exists a wealth of structural, spectroscopic, and kinetic data concerning them (Rifkind, 1973; Perutz, 1976). Further, the simplicity of the diatomic ligands (e.g., 0, and CO) permits one to focus on the role of the protein and its conformational changes. Nevertheless, it has not been possible as yet to give a complete molecular interpretation of the binding process, due to the difficulty of obtaining direct information concerning intermediates and the inherent complexity of the reaction (Austin et al., 1975). It is essential, therefore, to supplement the available experimental data by a theoretical analysis of the problem. In the present paper we limit ourselves to an investigation of the motion of the ligand through the protein matrix; binding to the heme group is not considered. Such a study begins with the known protein co-ordinates, locates possible t Present U.S.A.
address:
Department
of Chemistry,
University
of California,
Davis,
Calif.
95616,
343 0022-2836/79/230343-26
$02.00/O
0 1976 Academic
Press Inc. (London)
Ltd.
344
i).
.I.
CASE
r\Sl)
31.
ii:~R~'Li.~S
erkance and exit paths, and dot~erminc~s 11~3 tluct2ua,tions in thr,a tibruct8rii’4x rll’ t iit, protein t’hat can affect the barriers along these pa.th:~. Thr lat,tcr iti ~IFsrly c# importance for hemoglobin and myoglobin, since exa,minations of t,hr X-ray structure (Perutz & Mat#thews, 1966; Nobbs, 1966; Takano, I977a,0) have ihdicatecl that is blocked by cert,ain amino acid entrance of the ligand into the heme “pocket” side-chains (particularly HisE7). Magnuason (1971), in a useful survey of ligantl binding to myoglobin, has pointed out, that the residues HisE'I: VaEXl and PheCDl would all have short van der Waals’ contacts with one or another bound liga,nd. To facilitate the theoretica,l calculations, it is convenient to distinguish t,wo limiting regimes of dynamical behavior in prot,ein-substrat,e interactions. In the first, called the “adiabatic” limit (Gelin & Karplus, 1975), the protein conformational degrees of freedom relax on a time-scale short,er tfhan that, of substrate motion. For this case it is appropriate to speak of induced conformational changes a,long t’he binding path. The opposite regime is the “diabatic” limit in which substrate motion is faster than protein rela,xation, so that, the protein acts on the substrate as a rigid force field. Neither limit is likely to be an accurate approximabion to the true dynamics in most, situations, but each has its uses. In particular, the greater simplicity of calculations in the diabatic limit make them highly efficient for identifying possible reaction pa’ths, while conformation Auctua,tions along such pat,hs can be probed most easily in the adiabat)ic limit. In this paper, we examine theoretically the effect of t’he protein on t.he dynamics of the binding of oxygen and carbon monoxide to myoglobin. Empirical energy functions (Gibson & Sheraga, 1967; Levitt &, Lifson, 1969: Gelin & Karplus, 1975,1979) are used to describe the protein and the protein-ligand interactions. Both diabatic and adiabatic calculations are made. The former, supplemented by ligand trajectories, serve to determine possible paths for the ligand in ent’ering and leaving the heme pocket. The case of t,he ligand starting at its position bound to the heme and escaping through the protein into the external solution is the one investigated by trajectory calculations. To examine the effect, of protein relaxation on the magnitude of the barriers along the binding paths, calculations are made of the energy required for the rotation of various side-chains that, hinder the ligand. Finally, t’he ligand is positioned at the important barrier maxima and adiabatic calculations are made to determine the contributions of protein relaxation effects other than side-chain rotations. The empirical energy functions and the methodology used for t,he various calculations are described in Details of the Calculations. In Traject’ory Search for Binding Pathways, ligand traject,ories in the diabatic limit are used to search out and study possible binding paths. These a,re examined in more detfail in Conformational Relaxation, which is concerned with the consequences of protein relaxation. The results of the va,rious ca,lculations are discussed and compared with experiment in the Discussion.
2. Details of the Calcuiatians In this section, we describe the empirical used for the various calculations. (a) Potential
energy
function
and
the methodology
enwyy
This work employs an empirical energy function of t,he “‘extended atom” type (Gibson & Sheraga, 1967) for the protein conformationa. energy and for the prot,ein-
LIGAND
BINDING
TO HEME
PROTEINS
345
ligand interactions ; that is, only heavy atoms are trea’ted explicitly and the hydrogen atoms are subsumed into the heavy atoms to which they are attached. The total conformational energy is written as dihedrals
+
nonM2ded pairs
gz
-
;
+
y
+
h,,~g,, bond
g
-
g5>
(1)
pairs
where A, B, C and D are appropriate parameters (Gelin & Karplus, 1979). In equation (l), b is a bond length, 0 and + are bond and torsional angles, respectively, and r is the distance between non-bonded atoms. The parameters appearing in equation (1) are constants chosen to match the known conformational and vibrational properties of a series of molecules (Gelin & Karplus, 1979). The non-bonded interaction parameters reflect the fact that hydrogens have not been included ; e.g., the radius of a methyl carbon is set at 1.9 A, instead of 1.7 A for a “bare” carbon. Tests on small polypeptides have shown that the errors imposed by neglect of hydrogens are within an acceptable range (Gelin, 1976) and that the considerable increase in computational efficiency achieved by their neglect is worthwhile. The additional terms required to introduce the heme group have been described recently in connection with studies on the effects of ligand binding on the cc-chains of human hemoglobin (Gelin, 1976; Gelin & Karplus, 1977). Further details and values of the parameters are given in these references. The protein-ligand energy is introduced by assuming that the diatomic ligand can be approximated as a neutral sphere so that it interacts with the atoms of the protein through a Lennard-Jones (Hirschfelder et al., 1954) potential,
I;=&;
(2)
as in equation (1). The pa,rameters A and B are chosen in the same manner as the protein parameters; that is, the attractive term is determined by the Sla,ter-Kirkwood (Hirschfelder et al., 1954) formula and the repulsive term is adjusted to give the desired radius. The test pa,rticle has the characteristics of an inert oxygen atom, i.e. it is a sphere with a non-bonded diameter of 3.2 A. This value was deliberately chosen to be slightly smaller than the effective dia,meter of a diatomic ligand like 0, or CO( ~4.4 A), so that the potential might simulate slight rearrangements in the protein that are not included explicitly. For the trajectories, which were done in the diabatic. limit, the ligand was modelled with a still smaller radius (2.8 A). Table 1 lists the values of the Lennard-Jones parameters for the protein atom-ligand interactions obtained for the two assumed radii. (b) Diabatic
calculations
The starting conformation is the X-ray structure of deoxymyoglobin (154 amino acids) determined by Takano (1977a,b) by the refinement of 2 A dat,a. In this highspin derivative the iron atom is 0.55 A from the mean plane of the heme and 0.42 A from the plane of the four pyrrole nitrogen atoms; that is, the porphyrin group is slightly domed with pyrroles tilted toward the iron. A total of 141 hydrogen bond pairs are included in the energy function-the large number reflects the high helical
346
II Protein
_Il--,l---l-,l~~~-,~-,-,,-,,~atom
typet
ES 1’m 9: 1___-----~
-~Carbonyl or carboxy oxygen Alcoholic -OH Nitrogen (His,Pro) -N= Peptide -NH-NH, or -NH; Aliphatic -CHAlipha.tic -CH,or -C& Carbonyl or aromatic carbon Aromatic CH Sulfur
0.16 0.18 0.19 0.20 0.20 0.16 0.17 0~17 @I.8 0.05 0.20 0.19
IlYOll
Porphyrin
nitrogen
3.2 3..3 3.2 3.3 34 3.4 3.6 3.4 3.5 3.6 3.6 3.2
-t At,om t,ypes are from in terms of E and r,: $ Well-depth § Potent.ial 11Potential
Table
18 of Gelin
(1976).
The parameters
1’rnl‘I 2.8 2.9 2.8 2.9 3.0 3.0 3.1 3.0 3.1 3.2
3.2
2.8 -
ofeqn
(2) in the text
are defined
,4 = #‘f. B 2 &r” rn’ In. in kcal/mol. minimum (8) for Figs 1 to 7 and 10. minimum (A) for trajectories.
content of myoglobin. The list includes all those given by Takano for metmyoglobin (1977a, Table 10): w
(1979) 132, 243-368
Dynamics
of Ligand D. A.
CASE?
Binding M.
AND
to Heme Proteins KARPLUS
Department of Chemistry Harvard University Cambridge, Mass. 02138, U.S.A. (Received
7 December 1978)
The dynamics of CO or 0, binding to myoglobin are examined theoretically with particular emphasis on the effect of the globin on the barriers along possible reaction paths. By use of a diabatic representation, in which a model ligand moves in the potential of a rigid protein, two paths in or out of the heme pocket are found in myoglobin ; one is by the E helix near the distal histidine and the other near the intersection of the B, D and E helices. Classical trajectories calculated for a photoligand in the diabatic approximation yield complicated motion; the dissociated ligand tends to spend considerable time undergoing multiple collisions with the protein matrix in the wells between barriers. To investigate the height of the barriers, the adiabatic limit is used ; that, is, the protein is allowed to relax in the presence of a perturbation due to the ligand-protein interaction. A detailed description is obtained for the displacements of residues required for the ligand to enter or leave the globin. The barrier heights calculated in this manner are found to be in qualitative agreement with observed values. Some differences between myoglobin and hemoglobin are described. Although the present calculations make use of a simplified model of the ligandmyoglobin system, they provide a useful first step in the analysis of the effect of the globin matrix on the binding process.
1. Introduction A molecular description of substrate and ligand binding by enzymes and other proteins can for convenience be divided into two parts. The first deals with the initial and final geometries of the species involved, while the second is concerned with the reaction path and the dynamics of the binding process. Transport and storage proteins, such as hemoglobin and myoglobin, are favorable cases for study since both the ligated and unligated species are stable and there exists a wealth of structural, spectroscopic, and kinetic data concerning them (Rifkind, 1973; Perutz, 1976). Further, the simplicity of the diatomic ligands (e.g., 0, and CO) permits one to focus on the role of the protein and its conformational changes. Nevertheless, it has not been possible as yet to give a complete molecular interpretation of the binding process, due to the difficulty of obtaining direct information concerning intermediates and the inherent complexity of the reaction (Austin et al., 1975). It is essential, therefore, to supplement the available experimental data by a theoretical analysis of the problem. In the present paper we limit ourselves to an investigation of the motion of the ligand through the protein matrix; binding to the heme group is not considered. Such a study begins with the known protein co-ordinates, locates possible t Present U.S.A.
address:
Department
of Chemistry,
University
of California,
Davis,
Calif.
95616,
343 0022-2836/79/230343-26
$02.00/O
0 1976 Academic
Press Inc. (London)
Ltd.
344
i).
.I.
CASE
r\Sl)
31.
ii:~R~'Li.~S
erkance and exit paths, and dot~erminc~s 11~3 tluct2ua,tions in thr,a tibruct8rii’4x rll’ t iit, protein t’hat can affect the barriers along these pa.th:~. Thr lat,tcr iti ~IFsrly c# importance for hemoglobin and myoglobin, since exa,minations of t,hr X-ray structure (Perutz & Mat#thews, 1966; Nobbs, 1966; Takano, I977a,0) have ihdicatecl that is blocked by cert,ain amino acid entrance of the ligand into the heme “pocket” side-chains (particularly HisE7). Magnuason (1971), in a useful survey of ligantl binding to myoglobin, has pointed out, that the residues HisE'I: VaEXl and PheCDl would all have short van der Waals’ contacts with one or another bound liga,nd. To facilitate the theoretica,l calculations, it is convenient to distinguish t,wo limiting regimes of dynamical behavior in prot,ein-substrat,e interactions. In the first, called the “adiabatic” limit (Gelin & Karplus, 1975), the protein conformational degrees of freedom relax on a time-scale short,er tfhan that, of substrate motion. For this case it is appropriate to speak of induced conformational changes a,long t’he binding path. The opposite regime is the “diabatic” limit in which substrate motion is faster than protein rela,xation, so that, the protein acts on the substrate as a rigid force field. Neither limit is likely to be an accurate approximabion to the true dynamics in most, situations, but each has its uses. In particular, the greater simplicity of calculations in the diabatic limit make them highly efficient for identifying possible reaction pa’ths, while conformation Auctua,tions along such pat,hs can be probed most easily in the adiabat)ic limit. In this paper, we examine theoretically the effect of t’he protein on t.he dynamics of the binding of oxygen and carbon monoxide to myoglobin. Empirical energy functions (Gibson & Sheraga, 1967; Levitt &, Lifson, 1969: Gelin & Karplus, 1975,1979) are used to describe the protein and the protein-ligand interactions. Both diabatic and adiabatic calculations are made. The former, supplemented by ligand trajectories, serve to determine possible paths for the ligand in ent’ering and leaving the heme pocket. The case of t,he ligand starting at its position bound to the heme and escaping through the protein into the external solution is the one investigated by trajectory calculations. To examine the effect, of protein relaxation on the magnitude of the barriers along the binding paths, calculations are made of the energy required for the rotation of various side-chains that, hinder the ligand. Finally, t’he ligand is positioned at the important barrier maxima and adiabatic calculations are made to determine the contributions of protein relaxation effects other than side-chain rotations. The empirical energy functions and the methodology used for t,he various calculations are described in Details of the Calculations. In Traject’ory Search for Binding Pathways, ligand traject,ories in the diabatic limit are used to search out and study possible binding paths. These a,re examined in more detfail in Conformational Relaxation, which is concerned with the consequences of protein relaxation. The results of the va,rious ca,lculations are discussed and compared with experiment in the Discussion.
2. Details of the Calcuiatians In this section, we describe the empirical used for the various calculations. (a) Potential
energy
function
and
the methodology
enwyy
This work employs an empirical energy function of t,he “‘extended atom” type (Gibson & Sheraga, 1967) for the protein conformationa. energy and for the prot,ein-
LIGAND
BINDING
TO HEME
PROTEINS
345
ligand interactions ; that is, only heavy atoms are trea’ted explicitly and the hydrogen atoms are subsumed into the heavy atoms to which they are attached. The total conformational energy is written as dihedrals
+
nonM2ded pairs
gz
-
;
+
y
+
h,,~g,, bond
g
-
g5>
(1)
pairs
where A, B, C and D are appropriate parameters (Gelin & Karplus, 1979). In equation (l), b is a bond length, 0 and + are bond and torsional angles, respectively, and r is the distance between non-bonded atoms. The parameters appearing in equation (1) are constants chosen to match the known conformational and vibrational properties of a series of molecules (Gelin & Karplus, 1979). The non-bonded interaction parameters reflect the fact that hydrogens have not been included ; e.g., the radius of a methyl carbon is set at 1.9 A, instead of 1.7 A for a “bare” carbon. Tests on small polypeptides have shown that the errors imposed by neglect of hydrogens are within an acceptable range (Gelin, 1976) and that the considerable increase in computational efficiency achieved by their neglect is worthwhile. The additional terms required to introduce the heme group have been described recently in connection with studies on the effects of ligand binding on the cc-chains of human hemoglobin (Gelin, 1976; Gelin & Karplus, 1977). Further details and values of the parameters are given in these references. The protein-ligand energy is introduced by assuming that the diatomic ligand can be approximated as a neutral sphere so that it interacts with the atoms of the protein through a Lennard-Jones (Hirschfelder et al., 1954) potential,
I;=&;
(2)
as in equation (1). The pa,rameters A and B are chosen in the same manner as the protein parameters; that is, the attractive term is determined by the Sla,ter-Kirkwood (Hirschfelder et al., 1954) formula and the repulsive term is adjusted to give the desired radius. The test pa,rticle has the characteristics of an inert oxygen atom, i.e. it is a sphere with a non-bonded diameter of 3.2 A. This value was deliberately chosen to be slightly smaller than the effective dia,meter of a diatomic ligand like 0, or CO( ~4.4 A), so that the potential might simulate slight rearrangements in the protein that are not included explicitly. For the trajectories, which were done in the diabatic. limit, the ligand was modelled with a still smaller radius (2.8 A). Table 1 lists the values of the Lennard-Jones parameters for the protein atom-ligand interactions obtained for the two assumed radii. (b) Diabatic
calculations
The starting conformation is the X-ray structure of deoxymyoglobin (154 amino acids) determined by Takano (1977a,b) by the refinement of 2 A dat,a. In this highspin derivative the iron atom is 0.55 A from the mean plane of the heme and 0.42 A from the plane of the four pyrrole nitrogen atoms; that is, the porphyrin group is slightly domed with pyrroles tilted toward the iron. A total of 141 hydrogen bond pairs are included in the energy function-the large number reflects the high helical
346
II Protein
_Il--,l---l-,l~~~-,~-,-,,-,,~atom
typet
ES 1’m 9: 1___-----~
-~Carbonyl or carboxy oxygen Alcoholic -OH Nitrogen (His,Pro) -N= Peptide -NH-NH, or -NH; Aliphatic -CHAlipha.tic -CH,or -C& Carbonyl or aromatic carbon Aromatic CH Sulfur
0.16 0.18 0.19 0.20 0.20 0.16 0.17 0~17 @I.8 0.05 0.20 0.19
IlYOll
Porphyrin
nitrogen
3.2 3..3 3.2 3.3 34 3.4 3.6 3.4 3.5 3.6 3.6 3.2
-t At,om t,ypes are from in terms of E and r,: $ Well-depth § Potent.ial 11Potential
Table
18 of Gelin
(1976).
The parameters
1’rnl‘I 2.8 2.9 2.8 2.9 3.0 3.0 3.1 3.0 3.1 3.2
3.2
2.8 -
ofeqn
(2) in the text
are defined
,4 = #‘f. B 2 &r” rn’ In. in kcal/mol. minimum (8) for Figs 1 to 7 and 10. minimum (A) for trajectories.
content of myoglobin. The list includes all those given by Takano for metmyoglobin (1977a, Table 10): w
