Modelling and simulation of light-activated membrane proteins Dynamical transitions in bacteriorhodopsin Christian Simon,a Malika Aalouacha and Jeremy C. Smith*ab a L aboratoire de Simulation Moleculaire Section de Biophysique des Proteç ines et des Membranes DBCM CEA-Saclay 91191 Gif-sur-Y vette CEDEX France b L ehrstuhl f ué r Biocomputing IW R Universitaé t Heidelberg Im Neuenheimer Feld 368 D-69120 Heidelberg Germany Receiøed 2nd September 1998 Many of the functions of membranes are carried out by proteins associated with them. A knowledge of atomic-detail membrane protein structures and dynamics is required for a full understanding of these functions. We brie—y discuss recent progress in this –eld using modelling and simulation. One of the best characterised membrane proteins bacteriorhodopsin undergoes dynamical transitions with temperature.Here we present preliminary results of molecular dynamics simulation of this protein as a function of temperature indicating the presence of dynamical transitions at approximately the temperatures seen experimentally. Introduction Many of the functions of biological membranes are performed by proteins bound to them. Among the roles of these proteins are the reception/transmission of messages and/or the transport of materials. However due to difficulty in their crystallization only a small number of atomic-detail three-dimensional structures exist for membrane-spanning proteins. Among the few known structures are those of the light-transducing proteins the photosynthetic reaction centre and bacteriorhodopsin.1h4 In the present paper we brie—y review some recent progress in the modelling and simulation of light-activated membrane proteins before presenting new results on dynamical transitions in bacteriorhodopsin as a function of temperature. The paucity of crystallographic structures has led to a bottleneck in structural membrane protein research and adds impetus to the development of computer modelling techniques for determining their structures. One of these techniques is homology modelling it can be possible to determine an unknown protein structure by using a known X-ray structure as a three-dimensional template if there is sequence homology between the two. The higher the sequence homology the higher the probability of obtaining a reliable model structure.An example of homology modelling at high sequence identity is the recent work on the photosynthetic reaction centre protein from the bacterium Rhodobacter capsulatus.5 This protein has been the subject of a considerable amount of molecular biological and spectroscopic work aimed at improving our understanding of the primary steps of photosynthesis. A structural model was derived by combining information from the experimental structure of the highly homologous (54% sequence identity) reaction centre from Rhodopseudomonas viridis1 with molecular mechanics 95 Faraday Discuss. 1998 111 95»102 and simulated annealing calculations. In the Rb. capsulatus model the orientations of the bacteriochlorophyll monomer and bacteriopheophytin cofactors on the pathway inactive in electron transfer diÜer signi–cantly from those in the reaction centre of Rps.viridis. The orientational diÜerence was found to be in agreement with linear dichroism measurements.6 Moreover the pattern of cofactor hydrogen-bonding to the protein was found to be in agreement with optical spectroscopic experiment.7 The Rb. capsulatus model was used to provide an explanation as to why a partially symmetrized mutant Rb. capsulatus which has been of particular interest for experiments on primary excited states in photosynthesis lacks an electron acceptor bacteriopheophytin (BPhL).8h10 Conformational energy calculations on the partially symmetrised mutant and several BPh -binding revertants also provided an explanation for the relative BPh -binding properties of L L the proteins in terms of interactions involving two residues in the binding pocket these being a tryptophan and a methionine.10 Modelling at lower sequence homology although less reliable can be useful for suggesting experiments as part of an iterative procedure to obtain structural information on a membrane protein of particular interest.An example of this is the recent modelling of the photosystem II reaction centre core in plants for which a model was constructed by exploiting homology existing with the bacterial reaction centre proteins.11 In the rare cases where high-resolution experimental structures do exist modelling and simulation can be undertaken so as to re–ne structural detail and to understand physically how structure leads to function.A good example of such a system is bacteriorhodopsin (bR) a membrane protein that functions as a light-driven proton pump in the purple membrane of the bacterium Halobacterium halobium.12 The light-absorbing chromophore in bR is a retinal molecule that is covalently bonded via its SchiÜ base to the e-amino group of Lys 216.13 The characteristic purple colour of bacteriorhodopsin is due to absorption by the chromophore. The absorption is redshifted with respect to that of related model compounds in solution an eÜect that has been proposed to originate from interactions between the retinal and its polar environment in the protein.14 The retinal interactions may include hydrogen bonds with the SchiÜ base.Structures for bR at high resolution have been obtained.2,3 These revealed a channel through the protein that includes the SchiÜ base. Site-directed mutagenesis experiments suggest that the channel contains the pathway for proton transfer through bR.15h18 A considerable amount of data exist that suggest that the proton transfer channel is at least partially hydrated. Low resolution neutron diÜraction using contrast variation has indicated that about four water molecules are present in the neighborhood of the SchiÜ base although their positions in the direction perpendicular to the membrane plane could not be accurately determined.19 There is however considerable other evidence that water molecules are directly associated with the SchiÜ base. A resonance Raman study suggests that a negatively charged counterion located near the SchiÜ base group is stabilized by water molecules.20 Solid state 13C and 15N NMR experiments led to a model being proposed in which a water molecule is directly hydrogen-bonded to the SchiÜ base.21 Other solid state 1H and 15N NMR experiments suggest that there is a direct exchange of the SchiÜ base NH hydrogen with bulk water.22 A recent resonance Raman study of the SchiÜ base hydrogen»deuterium exchange also led to the conclusion that a water molecule is directly hydrogen bonded to the SchiÜ base NH proton.23 Finally the recent crystallographic structure of Pebay-Peyroula et al.has directly identi- –ed some water molecules associated with the SchiÜ base.3 Clearly a detailed understanding of SchiÜ base hydrogen bonding in the various stages of the photocycle will be required for a complete description of bR function.Computational chemistry has an important role to play in resolving such questions by identifying and quantifying hydrogen-bonding geometries and energies of pertinent model systems. For example quantum chemistry and molecular mechanics techniques have been combined to determine the geometries and energetics of retinal»water interactions.24,25 Ab initio molecular orbital calculations were used to determine potential surfaces for water»SchiÜ base hydrogen bonding and to characterize the energetics of rotation of the CwC single bond distal and adjacent to the SchiÜ base NH group. The ab initio results were combined with semiempirical quantum chemistry calculations to produce a data set used for the parameterization of a molecular mechanics energy function for retinal.Using the resulting molecular mechanics force –eld the hydrated retinal and associated bR protein environment were energy minimized and the resulting geometries examined. Two distinct Faraday Discuss. 1998 111 95»102 96 sites were found in which water molecules can make hydrogen-bonding interactions one near the NH group of the SchiÜ base in a polar hydrophilic region directed towards the extracellular side and the other near a retinal CH group in a relatively hydrophobic region directed towards the cytoplasmic side. To enable further investigations of internal hydration in bR and other systems a statistical mechanical formulation was derived that can be employed using molecular dynamics (MD) simulation to calculate the free energy of transfer of a small molecule from one environment to a speci–c site in another using molecular dynamics simulation.26 The method was used to calculate the free energy of transfer of water molecules from the bulk to individual sites in the proton transfer channel of bR.The channel contains a region lined primarily by nonpolar side-chains. The results obtained indicate that the transfer of water molecules from bulk water to this apparently hydrophobic region is thermodynamically favorable. The presence of two water molecules in direct hydrogen-bonding association with the SchiÜ base was also found to be thermodynamically allowed. Once a complete structural model of bR is obtained theoretical investigations into the photocycles of this protein can be envisaged.One interesting aspect of this in bR is the phenomenon of dark-adaptation in which retinal is found to exist in both all-trans and (13,15) syn conformations in approximately equal proportions. A theoretical investigation into dark-adaptation has been initiated. Initial free energy molecular dynamics calculations on a model of the isolated retinal suggested that the all-trans form is strongly favoured in vacuo.27 Calculations of factors in—uencing the conformational free energy diÜerence in the protein are now in progress. To fully understand bR function structural and thermodynamic examination must necessarily be complemented by dynamical investigations.Several analyses using molecular dynamics have been reported.28h33 In this respect it is of considerable interest that dynamical transitions have been found in bR as a function of temperature and have been correlated with function.34,35 In what follows we present preliminary results on the transitions investigated with molecular dynamics. The atomic position mean-square displacements are computed from a number of simulations at temperatures between 20 K and 300 K. Dynamical transitions are observed in the simulations at D150 K and D240 K. Methods Molecular dynamics The model system consists of 3544 atoms of bR of which 1806 are hydrogens. Four internal water molecules were included placed according to crystallographic data A 3 and each within 1 é of the positions derived in ref.27. The model system was subjected to molecular dynamics simulation using version 25 of the CHARMM program36 with the potential function described in ref. 37 26 and 27. The function includes bonded interactions (bond stretches bond angle bendings and dihedral and improper torsions) and nonbonded pairwise interactions represented by 12»6 Lennard- Jones and Coulombic electrostatic potentials which were cut-oÜ at 12 Aé . The Coulomb term was smoothed by multiplying by a switching cubic function between 8 and 12 Aé . The bR model molecule was simulated without an explicit environment. To approximately mimic the eÜect of the environment the relative permittivity was set to 1 and the a carbon atoms of the residues most surface exposed were harmonically restrained :28 the force constant used for this has a small value of 0.2 kcal mol~1 Aé ~2.This value was chosen so as to prevent gross deviation from the experimental structure while allowing internal —exibility. The system was energy minimized using 500 steps of Steepest Descent minimization followed by 2500 steps of Adopted Basis Newton»Raphson.36 The –nal RMS gradient was 0.15 kcal mol~1 Aé ~1. The energy-minimized structure was used as a starting point for the MD simulations. The equations of motion were solved using the Verlet algorithm. The SHAKE algorithm was applied to –x the lengths of the bonds involving hydrogen atoms. A 2 fs time step was used for integration of the equations of motion. MD simulations were performed at thirty temperatures from 10 to 300 K at 10 K intervals as follows.The minimized structure was heated to 10 K over 1 ps equilibrated over 5 ps then 10 ps of production was performed in microcanonical ensemble. The –nal frame of the production was then used for 1 ps heating to 20 K and the procedure 97 Faraday Discuss. 1998 111 95»102 Fig. 1 Temperature»time series calculated from the MD simulations. repeated. Subsequently the production runs were each extended by 100 ps yielding thirty 110 ps-long production trajectories at each temperature. The con–gurations were dumped to disk every 100 fs (50 steps) i.e. 1100 conformations per trajectory. Experimental connection Dynamical transitions of bR have been observed by incoherent neutron scattering experiments.The measurable quantity in incoherent neutron scattering the dynamical structure factor is the time Fourier transform of the intermediate scattering function I(q t) where q is the scattering wavevector and t is the time. I(q t) is the sum of the I(q t) of the individual atoms.38 In the case of scattering by a single atom in a harmonic potential I(q t) is Gaussian in q:39,40 (1) I(q t)\expS[q2c(t)V where and Sd2(t)T is the mean-square displacement de–ned as c(t)\16Sd2(t)T (2) Sd2(t)T\S[R(t)[R(0)]2T where R(t) is the atomic position vector and S… … …T indicates an ensemble average. The in–nite time limit of Sd2(t)T is (3) Su2T\ limSd2(t)T t?= Assuming limSR(t)R(0)T\0 then t?= (4) Su2T\2SdR2T where dR(t)\R(t)[SRT the displacement from the mean at instant t.The quantities Sd2(t)T and SdR2(t)T can be extracted from molecular dynamics trajectories by replacing the ensemble average with a time average. However the in–nite time Su2T can be obtained only from simulations long enough to sample all the motions involved. Correspondingly the experimental Su2T can be obtained only when the instrumental energy resolution is sufficiently good to resolve all the contributing motions. The neutron cross-section of a protein is dominated by the hydrogen atoms. Faraday Discuss. 1998 111 95»102 98 Results Stability of the simulation Fig. 1 shows time series of the temperature for each simulation. No signi–cant drift is seen at any temperature. The total energy was also found to be stable.The root mean square positional deviation (RMSD) between each bR conformation of the trajectory at each temperature and the initial energy-minimized structure is plotted in Fig. 2. The RMSD increases with temperature but remains approximately constant with time for most temperatures. Existence of dynamical transitions In Fig. 3 the simulation-derived hydrogen-atom Su2T is plotted against temperature. The increase in Su2T between 20 and 300 K is approximately 0.6 Aé 2 in accord with the neutron results obtained for dry purple membrane (PM).41 An in—ection is seen at D150 K a temperature at which a dynamical transition in bR has been experimentally reported.35,41 Fig. 2 Time series of the RMSD from the initial energy minimized structure at each temperature averaged over all the bR atoms.Fig. 3 Simulation-derived mean square displacement averaged over the hydrogen atoms vs. temperature. 99 Faraday Discuss. 1998 111 95»102 Fig. 4 Temperature dependance of the normalized variational contribution of the individual residues to the mean square displacement. The residues displacements were averaged over the hydrogen atoms. We introduce the variation with temperature of Sui2T the mean-square displacement of residue i (5) *Sui2T(T )\Sui2T(T ]dT )[Sui2T(T ) where dT is 10 K in this case. The normalized variational contribution of residue i is (6) **Sui2T(T )\ *Sui2T(T ) *Su2T(T ) where *Su2T\;i *Sui2T. **Sui2T is plotted against T and the residue number i in Fig. 4. Below 150 K all residues have **Su approximately equal i2T close to zero.At D140 K a dynamical transition is revealed by **Su variation of i2T for some residues. A second transition occurs at D240 K and above involving **Su larger variation of i2T than the D150 K transition and concerning a larger number of residues. A dynamical transition at D240 K has also been reported in neutron scattering work and has been correlated to the activation of bR function.34 Conclusions The modelling and simulation of membrane protein structures and dynamics is still in its infancy but will be of growing importance as more and more sequences of membrane proteins are determined. As structural research progresses the investigation of associated dynamical properties can also be expected to gain importance. The presence of a dynamical transition with temperature in water-soluble proteins has been recognized for some while.The recent neutron results on bacteriorhodopsin have demonstrated the presence of transitions also in a membrane protein.35,41 The present preliminary results indicate that transitions may also be apparent in molecular dynamics simulation. More research will be required to test further the present –ndings. In particular simulations of bR with an explicit membrane environment i.e. with the protein in trimer form with lipid and water surroundings rather than with harmonic constraints can be envisaged although they will be computationally demanding. Signi–cant variation in the dynamical transition properties with environmental changes has been documented.Tests of the temperature hysteresis would also be interesting to make. The question also arises as to how long a simulation would have to be performed at any given temperature to obtain converged mean-square displacements. Finally an Faraday Discuss. 1998 111 95»102 100 examination using MD of the eÜect of the application of the Gaussian approximation [Eqn. (1)] to the intermediate scattering function would be of interest. Two transitions are seen in Fig. 4 at D150 and D240 K. These are close to the temperatures at which transitions were seen experimentally. The higher-temperature transition is at about the water-melting temperature. But the present simulation was performed in the absence of water indicating that water is not required for it. However the D240 K transition is not apparent in the Su2T data in Fig.3 and may not therefore correspond to that seen experimentally. Work on this and other related questions is in progress. Acknowledgements We thank E. Pebay-Peroula for providing the crystallographic structure and D. Mihailescu J. Baudry and B. Costescu for useful discussions and preliminary calculations. References 1 J. Deisenhofer O. Epp K. Miki R. Huber and H. Michel Nature (L ondon) 1985 318 618. 2 R. Henderson J. M. Baldwin T. A. Ceska. F. Zemlin E. Beckmann and K. H. Downing J. Mol. Biol. 1990 213 899. 3 E. Pebay-Peyrouola G. Rummel J. P. Rosenbusch and E. M. Landau Science 1997 277 1676. 4 Y. Kimura D. G. Vassylyev A. Miyazawa A. Kidera M. Matsushima K. Mitsuoka K. Murata T. Hirai and Y.Fujiyoshi Nature (L ondon) 1997 389 206. 5 N. Foloppe M. Ferrand J. Breton and J. C. Smith Proteins Structure Function Genetics 1995 22(3) 226. 6 J. Breton E. J. Bylina and D. C. Youvan Biochemistry 1989 28 6423. 7 T. Mattioli personal communication. 8 M. H. Vos F. Rappaport J. C. Lambry J. Breton and J. L. Martin Nature (L ondon) 1993 363 320. 9 S. J. Robles J. Breton and D. C. Youvan Science 1990 248 1402. 10 N. Foloppe M. Ferrand and J. C. Smith Chem. Phys. L ett. 1995 242 238. 11 B. Svensson C. Etchebest P. TuÜery P. van Kan J. C. Smith and S. Styring Biochemistry 1996 35 14486. 12 D. Oesterhelt and W. Stoeckenius Nature (L ondon) New Biol. 1971 233 149. 13 K. J. Rothschild P. V. Argade T. N. Earnest K-S. Huang E. London M-J. Liao H. Bayley H.G. Khorana and J. Herzfeld J. Biol. Chem. 1982 257 8592. 14 R. A. Mathies S. W. Lin J. B. Ames and W. T. Pollard Annu. Rev. Biophys. Biophys. Chem. 1991 20 491. 15 T. Mogi L. J. Stern N. R. Hackett and H. G. Khorana Proc. Natl. Acad. Sci. USA 1987 85 5595. 16 T. Mogi L. J. Stern T. Marti B. H. Chao and H. G. Khorana Proc. Natl. Acad. Sci. USA 1988 84 5595. 17 L. J. Stern and H. G. Khorana J. Biol. Chem. 1989 264 14202. 18 T. Marti H. Otto T. Mogi S. J. Roé sselet M. P. Heyn and H. G. Khorana J. Biol. Chem. 1991 266 6919. 19 G. Papadopoulos N. Dencher G. Zaccaïé and G. Bué ldt J. Mol. Biol. 1990 214 15. 20 P. Hildebrandt and M. Stockburger Biochemistry 1984 23 5539. 21 H. J. M. De Groot S. O. Smith J. Courtin E. van der Berg C. Winkel J. Lugtenburg R.G. Griffin and J. Herzfeld Biochemistry 1990 29 6873. 22 G. S. Harbison J. E. Roberts J. Herzfeld and R. G. Griffin J. Am. Chem. Soc. 1988 110 7221. 23 H. Deng L. Huang R. Callender and T. Ebrey Biophys. J. 1994 66 1129. 24 M. Nina B. Roux and J. C. Smith in Structures and Functions of Retinal Proteins ed. J. L. Rigaud Colloque INSERM/John Libbey Eurotext Ltd. 1992 vol. 221 pp. 17»20. 25 M. Nina J. C. Smith and B. Roux J. Mol. Struc. (T HEOCHEM) 1993 286 231. 26 B. Roux M. Nina R. Pomes and J. C. Smith Biophys. J. 1996 670. 27 J. Baudry S. Crouzy B. Roux and J. C. Smith J. Chem. Inf. Comp. Sci. 1997 37(6) 1018. 28 M. Ferrand G. Zaccaïé M. Nina J. C. Smith C. Etchebest and B. Roux FEBS L ett. 1993 327 256. 29 W. Humphrey I. Logunov K. Schulten and M.Sheves Biochemistry 1994 33 3668. 30 W. Humphrey E. Bamberg and K. Schulten Biophys. J. 1997 72(3) 1347. 31 I. Logunov and K. Schulten J. Am. Chem. Soc. 1996 118(40) 9727. 32 D. Xu C. Martin and K. Schulten Biophys. J. 1996 70(1) 453. 33 D. Xu M. Sheves and K. Schulten Biophys. J. 1995 69(6) 2745 34 M. Ferrand A. J. Dianoux W. Petry and G. Zaccai Proc. Natl. Acad. Sci. USA 1993 90 9668. 35 V. Reç at H. Patzelt M. Ferrand C. P–ster D. Oesterhelt and G. Zaccai Proc. Natl. Acad. Sci. USA 1998 95 4970. 36 B. R. Brooks R. E. Bruccoleri B. D. Olafson D. J. States S. Swaminathan and M. Karplus. J. Comput. Chem. 1983 4 187. 37 A. D. MacKerell Jr. D. Bashford M. Bellot R. L. Dunbrack Jr. J. D. Evanseck M. J. Field S. Fischer J. Gao H. Guo S. Ha D. Joseph-McCarthy L. Kuchnir K. Kuczera F. T. K. Lau C. Mattos S. 101 Faraday Discuss. 1998 111 95»102 Michnick T. Ngo D. T. Nguyen B. Prodhom W. E. Reiher III B. Roux M. Schlenkrich J. C. Smith R. Stote J. Straub M. Watanabe J. Wioç rkiewicz-Kuczera D. Yin and M. Karplus J. Phys. Chem. B 1998 102 3586. 38 M. Beç e Quasi Elastic Neutron Scattering Principles and Applications in Solid State Chemistry Biology and Materials Science Adam Hilger Bristol 1998 UK. 39 A. Rahman K. Singwi and A. Sjoé lander Phys. Rev. 1962 126(3) 986. 40 A. Rahman Phys. Rev. 1963 140(4) 1334. 41 V. Reç at G. Zaccai M. Ferrand and C. P–ster in Biological Macromolecular Dynamics. Proceedings of a W orkshop on Inelastic and Quasielastic Neutron Scattering in Biology ed. S. Cusack H Bué ttner M. Ferrand P. Langan and P. Timmins Adenine Pres New York 1996 p. 117. Paper 8/06840B Faraday Discuss. 1998 111 95»102 102