Erix Wiliam Hernández-Rodríguez

Mechanisms of the Spectral Shifts for Retinitis Pigmentosa Mutants Explored by Quantum Mechanical/Molecular Mechanical Calculations.

Erix Wiliam Hernández-Rodríguez,^†,* Elsa Sánchez-García,^‡ Ana Lilian Montero-Alejo,^¶ Luis Alberto Montero-Cabrera,^¶ Walter Thiel^‡,*

^†Departamento de Bioquímica, Instituto de Ciencias Básicas y Preclínicas“Victoria de Girón”, 11600 Havana City, Cuba and Charité Centrum für Innere Medizin und Dermatologie, Biomedizinishes Forschungszentrum, Campus Virchow, Charité-Universitätsmedizin, 13353 Berlin, Germany.

^¶Laboratorio de Química Computacional y Teórica. Departamento de Química Física, Universidad de La Habana, 10400 Havana City, Cuba.

^‡Max-Planck-Institut für Kohlenforschung Mülheim an der Ruhr, 45470 Germany.
^*Corresponding author. Email: [e-mail] and [e-mail]

ABSTRACT: Retinitis Pigmentosa (RP) is a pathological condition associated to blindness due to a progressive retinal degeneration. RP-linked mutations leading to a less stable rhodopsin and changes at the retinal binding site (RBS) can also cause deviations of absorption spectra, affecting crucial functions; even when the stability could be reached natural or artificially. Mutation effects are largely unknown beyond stability; its solution sheds light onto the molecular mechanisms related to optical spectra and dark-state geometry. To evaluate the stability, geometries, electric influence and vertical excitation energies in the dark state of mutated human rhodopsins carrying the abnormal substitutions M207R or S186W at RBS, we mainly calculated rhodopsins within or out a lipidic bilayer using Molecular Dynamics, combined Quantum Mechanical/Molecular Mechanical approach, for ground state properties and the ab initio DFT/MRCI and TDDFT methods for excited state calculations. Both mutations diminished the rhodopsin stability, even when the folding and retinal binding could take place in these mutants. Substantial blue-shifted absorption spectra were observed as well as geometrical deviations and electric changes causing an unoptimal geometry for the photoisomerization and an increased energy in the dark state of mutated rhodopsins. The energy excess could lead to harmful reactions. The calculations explored in large regions of the conformational chromophore space were accurate, in very good agreement with available experimental data, providing a reliable atomistic insight on the mechanisms of these mutations near chromophore, which can open new therapeutic strategies with the purpose of minimizing the surplus energy and its consequences or/and using stabilizing molecules suggested in recent advances.

Vision is one of the most valued senses for humans. Ophthalmological pathologies as Retinitis Pigmentosa (RP), a condition associated to blindness, have been increasing since years ago as an important public health problem (1). RP refers to a heterogeneous group of inherited diseases characterized by progressive retinal degeneration due to death of the rod photoreceptor cells, the vertebrate photoreceptors dedicated to dim light vision (2, 3), and the outer field of vision is lost first when mutations affects the rhodopsin, the visual pigment in rods; central vision may be impaired at the beginning if the primary injury is in the visual pigments of cone photoreceptor cells, causing difficulty in tasks such as identifying colors and reading (4), although in some cases, other proteins and cells in retina some affected primarily (2). Even when a high genetic heterogeneity has been found, over 120 point mutations have been discovered in the rhodopsin gene. The large majority of these mutations cause autosomal dominant form (ADRP) of the pathology (3, 5). Patients suffering ADRP caused by rhodopsin mutations display night blindness, progressive loss of peripheral and, eventually, central vision and the characteristic accumulation of the intraretinal pigment deposits named lipofuscin, from which these dystrophies get its name (3-7), its pathogenesis and some mutation mechanisms are unclear; no a definitive treatment is available today.

Rhodopsin, the prototype of G-protein coupled receptor (GPCR) superfamily (8, 9), is a heptahelical transmembrane (TM) receptor protein expressed in retina, composed by four building blocks, (a) an apoprotein opsin with 348 amino acid residues, (b) a covalently bound chromophore (retinal), (c) two palmitoyl residues linked to Cys322 and Cys323 and (d) two sugars linked to Asn2 and Asn15 (10). In the dark state, the chromophore is the 11-cis-retinal, forming a retinal protonated Schiff´s base (RPSB) linkage with Lys296 at the rhodopsin binding pocket (RBP). Visible light absorption at ~500 nm (11) by the pigment triggers the isomerization of the 11-cis of the RPSB to the nonprotonated all-trans form of the retinal Schiff´s base (RBS) linked to the opsin upon a single photon absorption (10, 12), this reaction occurs with high efficiency characterized by a quantum yield in the range 0.65-0.67 (12, 13); providing rhodopsin with the energy to form the active state (14, 15); the primary photoproduct, photorhodopsin, is formed within a very short time (200 fs) and thermally relaxes within a few picoseconds to a distorted all-trans configuration, bathorhodopsin (10), which establishes the equilibrium with a blue-shifted intermediate (BSI) before forming lumirhodopsin. Metarhodopsin I intermediate is the next one followed by metarhodopsin II, the active conformation for G-protein coupling, Meta II state is formed by translocation of the Schiff’s base proton to residue Glu113, which plays the role of counterion for the RPSB in the dark state (12, 13). This sequence of events results in excitation of the visual nerve and perception of light in the brain (16).

RP-associated mutations can cause an impairment of protein folding or expression or retinal binding (17), affecting the rhodopsin function; but in other cases these processes can take place and/or could be favoured by drugs (18, 19), being the mutation influence on absorption and photoisomerization a more relevant aspect.

Extensive mutagenesis analyses on rhodopsin has shown that key residues as Pro267 are important, structurally; substitutions at that position can impair the transducin activation by affecting the structure of the G-protein interaction site and with implication for opsin folding, membrane insertion, assembly, and/or function (20). On the other hand, researches have been focused on RP-associated mutations for cysteines in rhodopsin, showing impossibility for regeneration with 11-cis retinal and transportation to plasma membrane (21). It is well documented by these studies that the disulfide bond forming Cys110 and Cys187 is necessary for an appropriate folding and receptor function (22).

Trp265 mutations have a large effect on the spectral properties of rhodopsin, changing the retinal binding site (20) as well as substitutions at Phe261 and Gly121, positions from the 6^th and 3^rd transmembrane helices, respectively, influences in the early steps of the photoisomerization, causing blue-shifted intermediates (23). In vitro, studies using recombinant rhodopsin for carrying amino acid substitutions associated with RP in different positions of the rhodopsin structure, including the RBP, show a modified spectral behavior (24, 25). Several other mutants studied related to, or not, with RP have been E113Q, E181Q, G90D, A292S, A269T, H211C and D83N displaying a diverse light-absorption pattern (16). A plentiful biochemical knowledge, molecular biology data on retinal proteins and high-resolution structural data by X-ray of bovine rhodopsin are available nowadays, for understanding the structure-function relationships at a pertinent molecular level for this prototypical GPCR (13, 26).

These advances with other experimental studies are a wealth used by theoretical studies to obtain new knowledge, besides to reproduce and to explain experimental data. The vast majority of researches on normal or wild-type (WT) and mutants structures have been performed for bovine rhodopsin and other retinal proteins with high-resolution structure data resolved experimentally; despite large experimental and computational studies in the last years as described and the plentiful knowledge reached on mechanism regulating the absorbance of 11-cis retinal chromophore in bovine rhodopsin (27), little data have been reported for human rhodopsin and its mutants, a not available crystal structure for this protein has been an important handicap.

Analyzing the mechanisms for RP-associated mutations needs structural information and hence it is crucial the human rhodopsin structure; although, its crystal structure is not solved until now, it is possible to achieve a reasonable three-dimensional (3D) structure, by means of Homology Modelling (28) combined with optimizations based on Molecular Mechanics (MM) and Molecular Dynamics (MD) methods (29), taking into account the high sequence similarity with bovine rhodopsin and enough biochemical data available. Subsequently, a reliable structure to calculate spectral properties and retinal geometry from human rhodopsin can be obtained, using QM/MM MD and QM/MM optimizations in different mediums (29, 30). On the other hand, robust methods for calculating the electronic spectra (31) can be used on an applicative optimized structure of human rhodopsin.

Hence, mechanisms of RP-associated mutations implicated with absorption shifts and photoisomerization, affecting the RBP in human rhodopsin, could be resolved appropriately, by theoretical methods at a very molecular level nowadays.

Obviously, rhodopsin mutations associated with RP act in many ways, from strong misfolding to severe impairments of folding/expression and/or binding retinal (32). However, other mutations can cause RP, changing amino acid residues near chromophore; even when the folding, binding retinal and light absorption take place somehow, its mechanisms are yet largely unknown.

M207R and S186W mutations cause ADRP (17) and affecting the RBP. A recent experimental study reports that the mutation M207R allows the rhodopsin folding, rhodopsin regeneration with 11-cis-retinal, increasing the experimental chromophore concentration and light absorption was detected in the dark state, outlining possible protonation-state changes for the SBL (33). A low misfolding at the second extracellular loop has been reported for the mutation S186W (34) and that thermally destabilize rhodopsin and to increase the rate of thermal isomerization, indirectly it is shown light absorption and folding in some grade (35); experimental studies about this mutation are not well documented. Previous calculation of the wavelength of absorption maxima (λ_max) for these RP-associated mutations were made on a very reduced RBP model; its geometry was optimized with MM methods, only, and no estimation of the protonation states was performed; the vertical excitation energies were calculated with a semiempirical method, giving some qualitative conclusions on light absorption pattern without reproducing a λ_max very near to experimental value (36).

On the other hand, a very recent computational study using FoldX analysis on disease-linked rhodopsin mutations included to M207R and S186W. The research found a change in protein stability, in a mutation more than in other, classifying these mutations within Class IIa; hence, no interference with retinal binding area can be possible for these mutants, since the mutations was not within Class IIb. Both substitutions could destabilize in some extent the protein according to FoldX predictions (the higher ∆∆G, the more severe folding effects). The study also reported that the disease mechanism could be especially valuable for the misfolding Class II mutants, taking into account the advances made using small molecules to stabilize a protein that otherwise tends to misfold (37). It is needful to deepen into the mechanisms of these mutations beyond the protein stability; since for these mutations the folding and retinal binding are possible naturally, or through drugs; no evaluation at a very molecular level on the mutation (M207R or S186W) effect concerning optical spectra and/or farther processes has been reported yet; studies in that sense are required for improving a deeper and necessary understanding for prospective therapies.

It is reported that the optical properties of the chromophores in retinal proteins are modulate by the protein moiety. Roughly speaking, three mechanisms are described for the spectral tuning: first, distortion of the retinal geometry as a result of steric interactions with the protein binding pocket; second, interaction of the RPBS with the counterion Glu113 balancing its positive charge; and third, interaction of the retinal with the rest polar amino acid residues lining the binding pocket (38).

Since 207 and 186 positions drop into the RBP and near chromophore in the 3D rhodopsin structure, M207R and S186W mutations could induce structural changes, affecting the dark state geometry of the RBP, including chromophore, and/or electric perturbations, impairing the appropriate mechanisms of photoisomerization and/or spectral tuning.

Calculating the spectra and absorption shifts with accuracy for human rhodopsin and retinal proteins as a very challenging problem, even using combined Quantum Mechanical/Molecular Mechanical (QM/MM) approaches; where the entire protein into account and thus are viable in understanding the visual event (16, 39). Several models, from the more simple to the more complexes or realistic and diverse methods have been employed for calculating the optical spectra (27, 36). At present, spectral properties of rhodopsin can be studied with reasonable exactness using combined QM/MM approaches (36).

Specific conditions are important for an accurate calculation of the absorption spectra: (a) highly accurate methods must be used, (b) geometrical parameters of the chromophore must be properly described because the spectrum is highly sensitive to the chromophore geometry, (c) chromophore polarization by the environment as well as environment polarization must be considered, (d) dispersion effects should be taken into account, (e) a sampling of the conformations is necessary to calculate absorption energies and spectra that are directly comparable to those from the experiment. Most approaches rely on geometry optimization, evaluating the spectrum at a single point in configuration space. However, this may not necessarily lead to representative structures (40).

Here, we apply combined QM/MM methods (30, 41), calculating the protonation states of all titratable amino acid residues of rhodopsins and with a large RBP as active region, using WT- and- mutants systems isolated or within a lipid membrane, obtained from previous homology modelling combined with MM optimization and MD; in order to evaluate the effect of M207R and S186W mutations on the absorption spectra, retinal geometry, electronic changes and its interaction with the protein environment in the dark state from mutants of human rhodopsin as well as to analyze possible mechanisms for these mutations related to energy balancing, photoisomerization and potential side or harmful reactions. Structural and stability characteristics are guided through classical MDs for protein-membrane complexes or QM/MM MD and QM/MM geometry optimization techniques using as QM component the self-consistent charge density functional tight-binding (SCC-DFTB) method and density functional theory (DFT, B3LYP functional), respectively. Vertical excitation energies are determined by time-dependent DFT (TDDFT) and a DFT-based multi-reference configuration interaction treatment (DFT/MRCI), it is not reported any study on these mutations with the methodology used and for reaching our aims.

Anticipating our results, we found a decreased stability due to these mutations; as well as, substantial blue-shifted absorption spectra in mutants carrying amino acid substitutions, caused by the RP-associated mutations M207R and S186W, with regard to WT rhodopsins. The light-absorption patterns were consistent with geometric and electronic changes leading to an unoptimal geometry for the photoisomerization and energy excess in the dark state.

METHODS

Starting Structure. Our models are based on the more recent and best-resolved X-ray structure of bovine rhodopsin [PROTEIN DATA BANK (PDB) (42), accession code: 1U19, resolution =2.2 Å], described in detail without missing amino acid residue (26). Only, the protein chain B (opsin, 11-cis retinal and two palmitic acid residues linked to Cys322 and Cys323) of bovine-WT rhodopsin was used with 29 crystal waters from this chain. Metal ions and other non-amino acid residues were omitted as well as sugar residues linked to Asn2 and Ans15, since these are oriented out from the structures and experimental studies report that it is not important for the correct opsin folding, palmitoylation, the formation of the characteristic rhodopsin chromophore (λ_max, ~500 nm) with 11-cis retinal and absorption spectra (43).


Homology Modelling. 3D-structure models of WT rhodopsins [both bovine-WT rhodopsin (bWT) and human-WT rhodopsin (hWT)], and RP-associated mutated rhodopsin, [both human-M207R mutant (hM207R) and human-S186W mutant (hS186W)] were generated by homology modelling using the starting structure as template; bovine rhodopsin was simulated (bWT) from its own experimental structure for a better homogeneity in later comparisons. M207R and S186W mutations were collected from International Retinitis Pigmentosa Association (IRPA) (17). The homologous proteins were identified with a BLAST (Basic Local Alignment Search Tool); the rhodopsin sequences, obtained from TrEMBL and SWISS-PROT databases of ExPASy Molecular Biology Server, reached a significant score of ~93 % with BLAST with the crystal structure of bovine rhodopsin, 1U19. Homology models were built using MODELLER (44) in version 9v7, the Cys110-Cys187 disulfide bridge was built and possible for all structures, the 29 crystal waters and palmitoyl moieties connected to Cys322 and Cys323 from the starting structure were conserved and five models with all hydrogen atoms was calculated for each system; crystal waters of WTs and mutants taken from the starting structure were placed appropriately. Internal waters, including two water molecules near chromophore within the RBP, are also in excellent agreement with the resolved water binding sites in the crystal structures for WT models and according to the expected for mutants. Retinal geometries were very similar between WT models and the crystalline structure; since mutants are analyzed, geometrical descriptions for chromophore are discussed later. Models were optimized by the variable target function method (VTFM, normal schedule) with conjugated gradient (CG) algorithm and refined using MD (fast level) with simulated annealing (SA). DOPE-HR method was used for selecting the model with lowest Discrete Optimized Protein Energy (DOPE) using the standard MODELLER energy function.

pK_a Calculations. pK_a values were calculated for all titratable amino acid residues of the structures bWT, hWT, hM207R and hS186W with all non-hydrogen atoms, both protein and ligands (11-cis retinal, palmitoyl moieties and waters) were included and PROPKA (45) in version 2.0 was used for the calculations, the protein stability in function of pH and residue interactions (i.e. hydrogen-bonding network at the retinal binding site) were also explored using this tool.

Protonation States. Visual inspection, previous obtained pk_a values, and the experimental pH values reported for vertebrate rod photoreceptors (46) were considered for assigning the protonation states to all titratable residues; the cytoplasmic pH (7.29±0.02) was taken into account for the transmembrane and cytoplasmic domains, it includes to residues from the RBP, and the intradiscal pH (6.5±0.07) for residues falling into intradiscal space, the limit for including residues into each region was established using our model with a lipidic membrane for every structure. Amino acid residues Asp83, Glu122 and Glu181 were protonated (neutralized) for the systems bWT, hWT and hS186W, in good agreement with UV–Vis spectral studies of site-directed mutants for the first two residues and experiments of FTIR in bovine-WT rhodopsin (47, 48). Asp83 and Glu181 were protonated and Glu122 was unprotonated (charged) for hM207R mutant. In vertebrate rhodopsin, the protonation state of Glu181 has not been established unequivocally, being a contentious issue to resolve. Two-photon spectroscopy studies (49), 3-D crystal structure and site-specific mutagenesis studies involving preresonance Raman vibrational spectra of the unphotolyzed E181Q mutant found Glu181 to be protonated in the dark state, paving the way for earlier and recent QM/MM studies on visual rhodopsin to assume Glu181 to be neutral (50). Arguments in favor of a unprotonated (charged) Glu181 have been put forward by NMR spectroscopic measurements (51, 52), FTIR, molecular dynamics simulations (13) and a very recent theoretical QM/MM study proposed a charged Glu181 and no significant implication of this residue on the bond length alternation (BLA) or optical spectra (53). Our study points to a neutral Glu181 for all studied systems, being homogeneous for this feature. His100 and His211 were protonated at δ nitrogen only, for WTs and mutants and additionally His195 for bWT. All systems had to His65 and His152 protonated at ε nitrogen only; all these histidine protonation states coincide with previous study (26); His278 (intradiscal) was protonated at both nitrogen atoms. Curiously, Cys316 was charged for mutants, this change is not significant for the absorption in the dark state, but for transducin activation (54), Cys110 and Cys187 were unprotonated, forming a disulfide bridge for all structures. Glu113 was negatively charged and compensated by the positive charge due to RPSB in bWT, hWT and hM207R; these states for both Glu113 and Schiff-base linkage (SBL) are well known for WTs. Experimental data hypothesized an unprotonated SBL for M207R mutation based on the absorption spectra, only (33). However, we found a similar spectral behavior and pK_a^SBL values point to a PRBS. hS186W showed to Glu113 protonated (pK_a^Glu113= 8.31) and an unprotonated SBL (pK_a^SBL= 5.88); so, the saline bridge between Glu113 and SBL is broken in hS186W mutant. Standard protonation states were used for all other amino acid residues. The protonation states assigned to significant residues for the optical spectra within the RBPs of WT systems, and for the closer to these, were fully agreement with experimental data.

MM Optimization. Once, protonation state were assigned using VMD (55) in version 1.8.7, the systems were then energy minimized to obtain a better homology model, using NAMD (56) in version 2.7 with CG algorithm, at 310 K and pressure (1 bar). These rhodopsin models were the starting point for other calculations.

Structure Validation. 3D-structure models were analyzed by three different structure validation tools, to test its reliability and internal consistency. Sequence–structure compatibility was evaluated by VERIFY-3D (57), backbone conformation was inspected by the Psi ⁄ Phi Ramachandran plot obtained from PROCHECK analysis (58) and the 3D-model packing quality was tested using the Z-score value calculation of the WHAT_CHECK (59). On the other hand, the root-mean-square deviation (RMSD) was calculated between each 3D-model and the crystalline structure for comparing the general structure of homology models with the crystal. RMSD was calculated, using VMD, for backbone only, because differences between systems for the number of atoms (Table 1). The programs were executed from different structure analysis and verification servers: http://www.doe-mbi.ucla.edu/Services/Verify_3D/ for VERIFY-3D, http://www.ebi.ac.uk/thornton-srv/software/PROCHECK/ for PROCHECK and http://swift.embl-heidelberg.de/whatcheck/ for WHAT_CHECK tool from WHAT IF program.

Membrane Embedding. Simulations of the complexes rhodopsin/membrane/water/ions were performed for determining if a crystal-like conformation of protein is conserved and stable after simulations adapting the membrane to the constrained protein, preventing large deviation, for both WTs and mutants; as well as, to pick snapshots from a final MD looking for extensive changes of conformation in order to evaluate its effects on the optical spectra, using a combined QM/MM approach. The stability was also analyzed during these classical MDs. Rhodopsin simulations within its native-like environment, a hydrated lipid bilayer, are possible as experimental data support that rhodopsins reconstituted in artificial membranes are functional (60) and several theoretical studies have reproduced very well bovine-visual-pigment properties with membranes (13). The rod outer segment disk membranes, in which rhodopsin is located are composed of unsaturated lipids (61) and phospholipids (10). Hence, the rhodopsins, oriented to Z axis, were embedded in a rectangular patch of palmitoyl-oleoyl-phosphatidyl choline (POPC) hydrated bilayer, using VMD. Additional solvation waters were placed using Helmut Grubmüller’s SOLVATE program (62), in version 1.0, for filling empty spaces inside the pores and surround these; water molecules within the rhodopsin-membrane interface were removed. The placement of rhodopsins in the POPC membrane was guided by the hydrophobicity of its residues; a very good membrane-protein alignment was obtained, bad lipids and water molecules overlapping with the rhodopsin or palmitoyl residues were eliminated. An appropriate room was made for every system WT or mutated in the POPC membrane; subsequently, complexes were placed in a box of water with a similar and slightly smaller size than the xy-plane, using Solvate plugin from VMD, any water molecule added inside lipid bilayer was removed. The systems were neutralized using Na⁺ and Cl^- ions at physiological concentration (154 mM); the complexes rhodopsin/membrane/water/ions with dimensions 83×84×87 ų contain the homology structures, 170 POPC molecules, waters and ions are included; mimicking rhodopsins, WTs or mutants, within a cell membrane (Figure 1).

Molecular Dynamic Calculations. A first energy minimization was performed on the complex, using CG algorithm at ambient temperature (300 K), pressure (1 bar) and constant volume followed by a MD for these conditions during 500 ps with the waters, ions, protein, 11-cis retinal, lipid headgroups fixed, except lipid tails and palmitoyl residues, using periodic boundary conditions (PBC). Electrostatic interactions for a full-system with PBC were calculated with the Ewald particle mesh method (63) and a grid spacing of about 1 Å, using a cutoff of 12 Å. The time integration step (timestep) was set to 1.0 fs. An appropriate disorder of a fluid-like bilayer was induced. Subsequently, the systems were guided to the nearest local energy minimum in configuration space with a second minimization at the same conditions before carrying out a new equilibration during 500 ps at 300 K, constant pressure (1 bar) with non-fixed atoms, keeping the conditions from the first one for PBC, PME (for full-system periodic electrostatics), grid spacing , cutoff and timestep; now, harmonic constraints were imposed on the protein and retinal (rhodopsin) and an exponent 2 was used for the harmonic constraint energy function. This simulation allowed that lipids, waters and ions to adapt to the protein in its crystal conformation, an appropriate environment relaxation and to avoid large deviations from the X-ray structure at the backbone and at the chromophore (64); in addition, the hydration of the membrane-protein interface was prevented. Consecutively, systems were heated to 310 K during a constant pressure (1 bar) MD run of 500 ps, with harmonic constraints released, PME and PBC, followed by a similar simulation; but now, at constant pressure (1 bar) and constant area during 1.5 ns, being possible to find some structural deviations regarding crystal in this final MD (FMD) in membrane. All minimizations and MDs were performed using NAMD with a total simulation time of 3 ns.

QM/MM Molecular Dynamics. All MD calculations were performed using the program CHARMM in version 33b1 (65, 66). The initial set of coordinates used for the MD simulations was taken from the structures bWT, hWT, hM207R and hS186W derived from homology modeling, which were optimized by MM. SHAKE constraints (67) were applied to all bonds involving hydrogen atoms. Due to the lack of suitable CHARMM force field parameters, the chromophore was treated using the semi-empirical SCC-DFTB method (68, 69) as implemented in CHARMM (70). The QM region for the systems bWT, hWT and hM207R included the chromophore 11-cis retinal, atoms N_ζ, H_ζ, C_ε, H_ε1, H_ε2, C_δ, H_δ1, H_δ2, C_γ, H_γ1 and H_γ2 of residue Lys296, atoms C_δ, O_ε1, O_ε2, C_γ, H_γ1 and H_γ2 of residue Glu113, and water (TIP3P) molecules 6 and 8; QM region for hS186W mutant included chromophore 11-cis retinal, atoms N_ζ, C_ε, H_ε1, H_ε2, C_δ, H_δ1, H_δ2, C_γ, H_γ1 and H_γ2 of residue Lys296, atoms C_δ, O_ε1, O_ε2, H_ε2, Cγ, H_γ1 and H_γ2 of residue Glu113, and water molecules 6 and 8 near to chromophore. The QM region was connected to the MM region using generalized hybrid orbital (GHO) boundary atoms (30, 71). Atoms C_γ of Lys296 and C_γ of residue Glu113 served as boundary atoms. The latter is a sp² carbon atom, while boundary atoms are generally assumed to be sp³ hybridized in the SCC-DFTB/CHARMM implementation. As a consequence, the average bond length to this boundary atom is 0.03 Å shorter than the experimental value, which is tolerable for our purposes. The bond involving the other boundary atom is only 0.01 Å too short.

An important issue in a QM/MM approach is the treatment of long-range electrostatic interactions. While it is recommended to use no cut-off for the interaction between the QM and the MM part of the system (41), it is usually necessary to truncate electrostatic interactions within the MM part to make the molecular dynamics feasible. To avoid an unbalanced treatment of the system, we used no cutoff for the QM-MM electrostatic interaction and applied the group-based extended electrostatics approach within the MM part (72). In this approach the electrostatic interactions between particles closer than a cutoff distance (14 Å in our case) are treated by the conventional pairwise additive scheme, while the interactions at larger distance are approximated by a computationally cheaper multipole approach.

After the initial setup the system was solvated in a pre-equilibrated sphere of TIP3P water with radius 30 Å, located at the origin of the coordinate system which was chosen to be at the position of atom C11 of the 11-cis retinal. In all MD calculations a spherical quartic boundary potential acting on the O atoms of the water molecules was used in order to retain the shape of the water droplet and to avoid evaporation of outer water molecules. All protein atoms more than 20 Å away from the center of the coordinate system were kept frozen in all subsequent calculations. After initial placement of the water sphere, overlapping water molecules (distance between Ow and any heavy atom < 2.8 Å) were deleted and all water molecules were energy-minimized in 250 cycles steepest-descent and 250 cycles ABNR (augment basis Newton-Raphson) minimizations. Positions of all non-water atoms were frozen at this stage. Thereafter the protein atoms were restrained by a strong harmonic potential, and the protein and water molecules were minimized using the same procedure as before. After this minimization protein and water atoms were subjected to a 15 ps constant-temperature MD with a time step of 1 fs, applying a starting temperature of 50 K, a final temperature of 300 K and a temperature step of 10 K every 100 MD steps. Subsequently, a new water sphere was superimposed on the system, overlapping atoms were deleted and the system was again subjected to a 15 ps MD, using weaker harmonic restraints on the protein atoms. This procedure was performed 12 times, each time lowering the harmonic restraints on the protein atoms. In the last three cycles the equilibration time was increased to 30 ps. Finally a 500 ps constant-temperature production MD with a time step of 1 fs was performed. Again the temperature was raised from 50 K to 300 K in steps of 10 K every 100 MD steps at the beginning. The VMD program was used for visualizing all the MD trajectories and resulting geometries (55).

QM/MM structure optimizations. QM/MM optimizations of picked snapshots from QM/MM MD (four) or FMD (three) were performed with the program ChemShell (73, 74) which is a modular package that allows the combination of several QM and MM codes. We used Turbomole, in version 5.7.1 (75), for the QM calculations and DL_POLY (76) as driver of the CHARMM22 force field. The QM part was described using the B3LYP density functional (77-81) and the TZVP basis set from the Turbomole basis set library.

Due to the low computational costs of the SCC-DFTB method we could use a relatively large QM region in the QM/MM MD calculations. This kept the number of QM boundary atoms low and avoided cutting of C-N bonds, for which the GHO method is currently not parametrized. Despite the higher-level QM/MM optimizations employing DFT, the QM regions selected were appropriated, even for the subsequent calculation of electronic spectra. We therefore chose the QM regions described for QM/MM MD calculations. Open valencies at the QM/MM border were saturated using hydrogen link atoms. We applied an electrostatic embedding scheme (82), i.e. the fixed MM charges were introduced into the one-electron part of the QM Hamiltonian and the QM/MM electrostatic interactions were evaluated from the QM electrostatic potential and the MM atomic charges. To avoid overpolarization of the QM region at the boundary, the charge on the first layer of MM atoms (those connected to the QM region) was moved to the second layer and a dipole (constructed from two point charges) was placed on the recipient atom of this charge shift, so that the charge and dipole of the MM system were conserved after the shift (83, 84). No electrostatic cutoffs were used.

Geometry optimizations were performed using the HDLC optimizer (85) implemented in ChemShell. The optimizations were done in hybrid delocalized internal coordinates. The coordinates of all atoms beyond 15 Å from the center (atom C11 of the 11-cis retinal) were fixed. Within the 15 Å active region, all atoms were allowed to move in each optimization step. The optimization was finished when the maximum gradient component was below 0.00045 a.u.

UV-Vis Spectra Calculations. We calculated vertical excitation energies and dipole moments in the protein environment applying both the combined density functional multi-reference configuration interaction (DFT/MRCI) method of Grimme and Waletzke (86) and the time-dependent-DFT (TDDFT) method based on the Runge–Gross theorem, which is the analogue of the Hohenberg–Kohn theorem in ground state DFT (87). The same QM regions and boundary scheme was used as in the preceding B3LYP/TZVP optimizations. Ten roots with singlet multiplicity were calculated for each snapshot and the TZVP basis set from the Turbomole basis set library was employed. All valence electrons (identified by their MO energies) were correlated.

Density functional theory (DFT) has emerged as a powerful method for calculating potential energy surfaces of ground electronic states. TDDFT determines excited state energies through a response formalism which avoids the direct computation of wavefunctions for electronically excited states. The basic principle in TDDFT is that a time-dependent system of interacting electrons can be mapped onto a model time-dependent system of non-interacting electrons with the same one-electron density as the true interacting system. Contradictions remain regarding the ability of method to describe states with double-excitation character, and it is now widely agreed that the method fails to appropriately describe excited states with significant charge transfer character (87). However, TDDFT reproduced very well the vertical excitation energies of WT rhodopsins studied in this work and additionally, the robust DFT/MRCI method was also used for studying the spectral behavior in our systems.

The basic idea of the DFT/MRCI method is to take major parts of the dynamic electron correlation into account by DFT, while static correlation effects are included by a short CI expansion. The configuration state functions in the MRCI expansion are built up from Kohn-Sham orbitals of a closed-shell reference state. Diagonal elements of the effective DFT/MRCI Hamiltonian are constructed from the corresponding Hartree-Fock based expression and a DFT-specific correction term. The effective DFT/MRCI Hamiltonian contains five empirical parameters, which depend only on the multiplicity of the excited state, the number of open shells of a configuration, and the density functional employed. Currently, optimized parameter sets are available for the DFT/MRCI Hamiltonian in combination with the BH-LYP functional (77-80) (88), to speed up our calculations we did not use the original parameterization from Grimme and Waletzke (86) but a modified set of parameters developed more recently in the Grimme group (89) (p₁=0.629, p₂=0.611, p_J=0.119, p[0]=8.000, α=0.503) which is currently only available for singlet electronic states. The required DFT wavefunctions were generated using Turbomole, version 5.71 (75).

As is common practice in the DFT/MRCI approach, the CI space was kept moderate by only selecting configurations whose estimated energy is beyond a threshold δE_sel above the highest reference space energy. While the usual choice for this parameter is 1.0 Hartree, we used a value of 0.8 Hartree in our calculations. This was possible due to the modified DFT/MRCI parameter set applied, which improves convergence of the electronic energies with respect to δE_sel. Using this combination of δE_sel and DFT/MRCI parameters, calculations are sped up by a factor of ~20, while excitation energies are only slightly changed (in our experience mostly higher by less than 0.1 eV) (90). Using standard parameters the DFT/MRCI method yields vertical excitation energies which are typically within 0.2 eV of the experimental values (86). Taking CASPT2 results as the reference, the standard DFT/MRCI approach typically underestimates excitation energies by 0.2 eV, as shown in a recent benchmark study (91).

The MRCI reference space was determined iteratively by first performing a DFT/MRCI calculation with a configuration selection threshold δE_sel=0.6 Hartree and a reference space that included all configurations which can be generated by up to doubly exciting ten electrons within the ten frontier MOs. In the next iteration δE_sel was increased to 0.8 Hartree and all configurations with a squared coefficient of at least 0.003 in the previous CI calculation were included in the reference space. Further iterations were performed until the energies of all calculated roots were converged to 10^-6 Hartree.

Parameters and Force Field. All MM calculations, be pure MM or combined QM/MM method, were carried out with CHARMM force field for all atoms and TIP3P water model. All-atom CHARMM force field was used for the protein; so, all force field charges, types and parameters of each atom for all amino acid residues were selected from CHARMM library (92). Retinal and Lys296 were tried as independent residues and linked between each other, through the appropriated patch for SBL in each system. All available force field parameters, types and charges of each atom for non-amino acid residues as retinal (Loccisano et al, 2005) or palmitoyl residues (Carlos F. López, 1999); as well as, for protonated or unprotonated SBL patched to corresponding system, and ions were also taken from CHARMM library (93-97).

Data Analysis. Results from every snapshot were summarized and expressed as mean ± standard deviation (SD) for each system. The assumptions for parametric test were explored and the significance was evaluated with F from one-way ANOVA (Post Hoc) for the difference between groups and correlation-coefficient r for relation between variables, inferring to structure population. Results were considered significant at p<0.05 and taking into account the information from 95 % confidence intervals (CI), using SPSS software (98), in version 18.0. WMD, NAMD or CHARMM were used for RMSD calculations, depending on each stage. VMD was used for drawing molecular structures.