Molecular mechanism of long-range synergetic color tuning between multiple amino acid residues in conger rhodopsin.

The synergetic effects of multiple rhodopsin mutations on color tuning need to be completely elucidated. Systematic genetic studies and spectroscopy have demonstrated an interesting example of synergetic color tuning between two amino acid residues in conger rhodopsin's ancestral pigment (p501): -a double mutation at one nearby and one distant residue led to a significant λ(max) blue shift of 13 nm, whereas neither of the single mutations at these two sites led to meaningful shifts.To analyze the molecular mechanisms of this synergetic color tuning, we performed homology modeling, molecular simulations, and electronic state calculations. For the double mutant, N195A/A292S, in silico mutation analysis demonstrated conspicuous structural changes in the retinal chromophore, whereas that of the single mutant, A292S, was almost unchanged. Using statistical ensembles of QM/MM optimized structures, the excitation energy of retinal chromophore was evaluated for the three visual pigments. As a result, the λ(max) shift of double mutant (DM) from p501 was -8 nm, while that of single mutant (SM) from p501 was +1 nm. Molecular dynamics simulation for DM demonstrated frequent isomerization between 6-s-cis and 6-s-trans conformers. Unexpectedly, however, the two conformers exhibited almost identical excitation energy, whereas principal component analysis (PCA) identified the retinal-counterion cooperative change of BLA (bond length alternation) and retinal-counterion interaction lead to the shift.

Adaptation of living organisms to various photic environments has been enabled by a broad range of color tuning of visual pigments such as rhodopsins. Since the transition dipole moment of the protonated Schiff base of retinal (PSBR), responsible for the light absorption of rhodopsin, is sensitive to the interaction between the chromophore and the surrounding protein environment, visual pigments often exhibit dramatic shifts in the maximal absorption wavelength (λmax) with amino acid replacements. In fact, λmax values of the PSBR vary with the environment: 610 nm in vacuo1, 442 nm in solution2, and 480–520 nm in rhodopsin3.

Rhodopsins of most terrestrial vertebrates and those living in shallow water are tuned to around 500 nm, which exactly corresponds to the peak wavelength of the solar irradiance spectrum at the Earth’s surface. In contrast, in deep water, the available wavelength is limited to around 480 nm. During evolution, deep sea fishes’ ancestors, which lived in shallow water, adapted to poor-light environments by shifting their rhodopsin λmax to the blue range. Rhodopsins of the Japanese conger (Conger myriaster) and its ancestor, subjects of the present study, are a typical example of color tuning, with λmax values of 486 and 501 nm, respectively. The molecular basis for this adaptive evolution is an important aspect to be studied. Yokoyama et al. recently reported that among the 27 amino acid sites mutated in the two pigments, only 3 sites, 194, 195 and 292, are primarily responsible for the λmax shift. Pro-194, Asn-195, and Ala-292 in the ancestral pigment have been replaced by Arg, Ala, and Ser in conger rhodopsin4 Interestingly, although single mutations at each site induce little or no λmax shift, λmax of double mutants, P194R/A292S and N195A/A292S, were dramatically shifted from 501 to 488 nm. Notably, the amino acid sites 194 and 195 are 17 Å away from the retinal in bovine rhodopsin, in contrast to site 292, which is located closer to the retinal (Fig. 1). The question is how these residues control λmax of the visual pigment rhodopsins.

Figure 1.

Protonated Schiff-base of retinal and amino acid residues at the sites 195 and 292 in bovine rhodopsin.

The wavelength selectivity of rhodopsin is due to the variability in the electronic transition energy of the retinal from the ground to excited state; Intrinsic retinal structure and retinal-protein interactions play essential roles determining the electronic state. Because a previous mutagenesis study demonstrated that all 25 residues responsible for noticeable λmax shifts are located within 11 Å of the chromophore5, the short-range retinal-protein interactions were implicitly assumed. Compared to short-range interactions, the effects of long-range interactions are more difficult to understand coherently. Although there is no direct interaction with PSBR, a mutation of a remote residue may affect the local structure of PSBR as implied in a previous study6. Considering this evidence, the double mutations at sites 292 and either 194 or 195 are quite likely to induce conformational change in the protein in a synergetic manner, resulting in the dramatic λmax shift, which is not observed for respective single mutations.

Until date, many computational approaches have been used to investigate the color tuning mechanism of rhodopsins. However, achieving the predictable accuracy of an experimental λmax value from computations has caveats: methodology of excitation calculations, appropriate structure optimization before excited-state calculations, and spectral broadening due to static disorder.

Various theoretical approaches (such as CASPT27–8, SACCI9, TD-DFT10, and SORCI11) have been applied to the excited state of the retinal chromophore. Unfortunately, even CASPT2, which is one of the most accurate methods, has an error of up to 0.1 eV for the excitation energy of benzene12. This value corresponds to a wavelength of around 10 nm in the visible range. Thus, present computational techniques are obviously inadequate for discussing the absolute absorption wavelength value of more complex retinal molecule, with an accuracy of 1.0 nm. However, this does not necessarily dismiss qualitative discussion. In addition, an excitation energy difference between calculations for the same system with different configurations can still be more accurate because of error cancellation. The present study also focuses on energy differences between the mutants of conger rhodopsin, as with many previous studies have succeeded in reproducing mutational wavelength shifts in rhodopsin.

Second, consistency between structure optimization and excitation energy calculation guarantees a reliable theoretical analysis. Static and dynamics correlation should be considered for both excitation energy calculations and structure optimization. This study adopted the CASPT2 method for excitation energy calculations and the B3LYP method for structure optimization, since structure optimization with CASPT2 is not practical for PSBR. Although B3LYP is a density functional theory (DFT) method with no systematic way to improve accuracy, Page et al. showed that PSBR structure optimized with B3LYP agree with that with CASPT2, which incorporates the static and dynamic correlations in a balanced manner13. This study analyzed the consistency between the B3LYP and CASPT2 calculations.

Third, structural disorder of the retinal accommodated or “solvated” in protein environments causes spectral broadening. The retinal in vacuo has double minimum, namely, cis and trans, isomers, on its potential energy surface. Thus, the structure optimization results in one of these two isomers irrespective of the initial conformation used for optimization. In contrast, the potential energy surface of the retinal in an explicit protein environment exhibits a number of local minima, leading to variation of the final structure of energy optimization, namely, static disorder. Accordingly, the distribution of excitation energies causes inhomogeneous broadening of the optical absorption spectrum14–15.

To overcome these problems, we formulated a coherent strategy: statistical evaluation of excitation energies based on quality optimized structures using thermodynamics ensemble of the visual pigments. Unfortunately, we did not have the experimentally derived structural information regarding the ancestral pigment of conger rhodopsin (p501), the single mutant A292S (SM), and the double mutant N195A/A292S (DM). A homology modeling method was performed before nano-second scale molecular dynamics simulations of the proteins in explicit membrane environments to obtain initial models for these visual pigments. Then multireference perturbation theory calculations were used to obtain the statistical ensembles of excitation energies for the visual pigments. To provide consistent rationale for the color tuning mechanism, we used principal component analysis (PCA) and identified the structural determinants for λmax shifts.

Material and methods

Homology modeling

The program Modeller v8.116 provided the 3-D models for p501, A292S and N195A/A292S, using the bovine rhodopsin [PDB code 1U19;17] structure as a template. The CLUSTALW18 program aligned the template sequence and each of the three sequences. The top five models out of an ensemble of 100 models were selected on the basis of best fit to the restraints, and an assessment of stereochemical quality was then performed using PROCHECK19 to determine the overall best model. Each model was then inserted into a lipid bilayer environment using NAMD20.

Structural model

Since the experimental structures of conger rhodopsin have not been available, we substituted the structures modeled by sequential homology with bovine rhodopsin as initial coordinates for the following MD simulations. In the atomic model, the N-terminal ends of the three rhodopsins were blocked by acetyl groups. All ionizable amino acid residues, except for Asp-83 and Glu-122, were assumed to be in charged states. Protonation of Glu-181 in the dark state remains an ongoing debate21–24. A recent QM/MM study24 indicates that the switching of the protonation state of Glu-181 cause no detectable change in the rhodopsin structure. On the other hand, a MM study21 demonstrated that the local structure at the Glu-181 site was stable only for the deptoronated model. It is likely that the rhodopsin structure does not depend on the protonation state of Glu-181 in a real molecule, and the instability of the protonated model might come from the inaccuracy of the MM parameters. Unfortunately conformational sampling of the rhodopsin with the QM/MM method is not feasible with the current computational power. In the present study, we adopted the protonated model for Glu-181, because the stability of the molecular trajectory around the experimental X-ray structure is of primary importance for conformational sampling. Cys-110 and Cys-187 were connected via a disulfide bridge. The topology of PSBR was derived from a previous study25 in which the equilibrium bond lengths and angles were determined by ab initio calculations26–27, parameters for dihedral angles were based on DFT calculation28, and non-bonded interactions suitable for bacteriorhodopsin were taken from the literature29. The partial atomic charges of the retinal chromophore were taken from QM/MM studies30–31.

Molecular dynamics simulation

The following MD simulations were performed using NAMD 2.220, which used CHARMM22 and 27 as the potential energy function for proteins and lipid molecules, respectively32–33. Periodic boundary conditions were applied at every side of the rectangular simulation box. Electrostatic interactions were calculated using the particle mesh Ewald method. Berendsen pressure coupling was used with a relaxation time parameter, τp=0.5 ps, and compressibility parameter, β=4.57×10−5 bar−1. For the initial model obtained by the homology modeling method, energy minimization was performed for the structures modeled by homology with harmonic restraints imposed on the backbone atoms and followed by simulated heating for 100 ps. We generated the initial velocities of atoms assuming a Maxwell-Boltzmann distribution at 50 K. The temperature was increased by 10 K every 2.0 ps until it reached 300 K. The system was allowed to equilibrate at 300 K for 500 ps during which the atomic velocities were updated every 10 ps. Subsequently, a 300 ps Langevin dynamics simulation was performed using the Berendsen constant pressure method, and harmonic restraints on each peptide backbone atom were linearly decreased, until they totally vanished in the first 100 ps. We then continued an 8.0-ns simulation with the above system as the initial coordinate, and four coordinate sets were extracted from the 8.0-ns simulation at every 2.0 ns. These four coordinate sets were used as the initial coordinate sets of subsequent production runs of 12.0 ns each. As a result of the production run, we effectively obtained a 48.0-ns trajectory for each of the three visual pigments.

Constraint optimization

Constraint optimizations at the B3LYP/6-31G* level were performed to illustrate the potential energy surface along the C6-C7 dihedral angle. All inner coordinates were relaxed in these optimizations, except for the C6-C7 dihedral angle. The system comprised of the retinal and Schiff base capped by a hydrogen atom instead of the methyl group of Lys-296 (Fig. 2A). All electronic state calculations were assumed in vacuo. These optimizations were performed using Gaussian 0334.

Figure 2.

Protonated Schiff-base of retinal (PSBR). (A) QM region of QM/MM optimization. (B) The model structure for excited state calculations.

QM/MM optimization

From MD trajectories, 100 different snapshots were extracted at every 0.5-ns time point, and these structures were subjected to a two-stage energy minimization: CHARMM and B3LYP/CHARMM. MM energy minimizations using the CHARMM22 force-field parameters32–33 were followed by QM/MM energy minimizations in a two-layer ONIOM-ME (mechanical embedding) scheme35. The electronic state of the QM layer was evaluated at the B3LYP/6-31G* level, and the energy of the MM layer was calculated using the CHARMM22 force field functions.

Absorption wavelength calculation

To compute excitation energies for the QM/MM optimized structures, multi-state CASPT2 (hereafter simply noted as CASPT2) calculations were performed based on the state-average four-root wave functions obtained by CASSCF (12e/12o) using the ANO-L (3s2p1d/2s) basis set. To deal with the problem of intruder states, we adopt level shift parameter of 0.3 for all the CASPT2 calculation. The excitation energy calculation treats the electrons of system made of PSBR and Glu-181 as shown Figure 2B, where the Nitrogen atom of Schiff-base and Cε atom of Glu-181 were capped by hydrogen atom respectively.

Results and Discussion

Excitation energy

For the calculation, we adopted three different models: vacuo model, transparent Glu-181 model, and deprotonated Glu-181 model. In the vacuo model, the external charges of surrounding environment were neglected, while they are treated as point charges derived from atomic partial charges of the CHARMM force field for the transparent Glu-181 and the deprotonated Glu-181 models. Note that the atomic partial charges in Glu-181 were neglected for the transparent Glu-181 model. Figure 3 illustrates statistical distributions of excitation energies for the three models. The distribution of DM markedly shifted to higher energy than p501 and SM for the transparent Glu-181 model (Fig. 3 DEF), while this is not the case for the deprotonated Glu-181 model (Fig. 3 GHI). The values of λmax shift from p501 to SM (DM) were +1 nm (−8 nm) for the Glu-181 transparent model (Table 1), being consistent with the experimental observations of λmax shifts for SM (+1 nm) and DM (−13 nm). In contrast, notable shift of λmax were not observed among p501, SM and DM for the deprotonated Glu-181 model. It should be noted that the characteristic blue shift of DM was also recognized for the vacuum model (Fig. 3 ABC), indicating that the electrostatic interaction between retinal and surrounding environment is excluded from the color tuning mechanism of conger rhodopsin, which is one of the possible determinants for the excitation energy as described below. The excitation energies were larger by 0.4 eV from the experimental values. The overestimates of exciation energy are discussed in another section.

Figure 3.

Distributions of excitation energies of p501 (A, D, G), SM (B, E, H) and DM (G, H, I), repectively. The model consists of PSBR and counterion placed in (A, B, C) vacuo (D, E, F) the protein environment without Glu-181 and (G, H, I) the environment of the entire protein with deprotonated Glu-181. The histrograms were constructed using 100 structures for each of the three pigments (red). λmax values were determined by Gaussian fitting (black).

Table 1.

Peak excitation energies (eV) of the spectra in Fig. 3 and estimated λmax (nm)/full width at half-maximum (eV)

	In vacuo	Transparent E181	Deprotonated E181	Experiment
p501	2.96 (419)/0.12	2.91 (426)/0.12	2.96 (419)/0.14	2.47 (501)/(0.53)*
SM	2.94 (421)/0.14	2.90 (427)/0.15	3.00 (413)/0.21	2.47 (502)/ —
DM	3.03 (409)/0.05	2.96 (418)/0.11	2.97 (417)/0.31	2.54 (488)/ —

Experimental value of FWMH of p501 was taken from that of bovine rhodopsin for reference60.

A recent QM/MM study on the protonation state of Glu-181 demonstrated that the protonation made no difference in the NMR spectra24. Therefore, we assumed that the rhodopsin structure is not affected by the protonation state of Glu-181, and the present calculations support the charge-neutral, protonated model for Glu-181. According to the literature36, most mutations incorporated in the site do not affect the λmax significantly, except for E181Q and E181K, indicating that the electronic charge at the site has little influence on the excitation energy.

During the last three decades, at least three factors have been proposed to be responsible for the color tuning mechanisms of visual pigments37–39. The first factor originated from the intrinsic features of the chromophore structure such as the distortion of the conjugated polyene backbone chain. This distortion causes π-electron localization, which leads to stronger bond length alternation along the polyene chain and a λmax blue shift. The equilibrium PSBR structure varies with the surrounding environment, as seen in the X-ray structures of various retinal proteins.

The second is a charge separation between the positively charged PSBR and its counterion, such as Glu-113 for bovine rhodopsin. Upon excitation, the positive charge localized around the Schiff base proton relocates toward the β-ionone ring40. Thus, the electrostatic interaction between PSBR and its counterion stabilizes the ground state more than the excited state. The proximity of the counterion to the Schiff base proton increases the excitation energy, thereby resulting in a blue shift.

Electrostatic interactions between the retinal and its surrounding protein environment is the third factor, which results in a 33 kJ/mol (0.34 eV) reduction in the excitation energy41.

In the following sections, possible origins of the reproduced blue shift are examined from the above mentioned viewpoints, after which PCA is used to identify the determinants for the color tuning mechanism of conger rhodopsin.

Mutation effects on chromophore structure

For DM, MD simulation demonstrated a skewed β-ionone ring of PSBR that was remarkably different from those for SM and p501. The β-ionone ring of PSBR observed for SM and p501 were similar to each other (Fig. 4). This characteristic difference in the β-ionone ring was attributed to the C6-C7 dihedral angle (Fig. 5). p501 and SM have single peaks around −50° (6-s-cis), whereas DM exhibits double peaks around −150° (6-s-trans) and −50° (6-s-cis). Both of these double peaks were well-conserved through the QM/MM optimization (Fig. 5). Although the highest resolution X-ray structure17 exhibits a 6-s-cis retinal conformer, the existence of the 6-s-trans conformer for the dark state of rhodopsin has been implied from both experimental42–44 and theoretical45–46 studies.

Figure 4.

Average structures of the protonated Schiff-base of retinal and the surrounding amino acid residues are shown. The thermal average structure was obtained from MD simulations for p501 (red), SM(orange), and DM (green).

Figure 5.

(A) Distributions of the C6–C7 dihedral angle. The statistical populations were obtained from MD simulations for p501 (black), SM (orange), and DM (green). (B) Distributions of the C6–C7 dihedral angle for p501 (black)), SM (orange), and DM (green). The statistical populations were analyzed after the QM/MM optimizations performed for each 100 structures.

In PSBR, the π-electron population is delocalized and extended from the Schiff base to the C5 atom in the β-ionone ring. Thus, the distortion in the C6–C7 dihedral angle could affect λmax. Notably, squid rhodopsin has a twisted retinal β-ionone ring, with its λmax value of 480 nm, which is blue-shifted from that of the major vertebrate rhodopsins47–48.

A brief discussion regarding the characteristic features of the C6–C7 dihedral angle may be required. Until date, various theoretical methods have been used for the retinal structure in vacuo 11,27,40,49–50. Notably, all of these supported non-planar conformers for both 6-s-trans and 6-s-cis forms, where steric hindrance occurred between methyl groups bound to the C1 or C5 atom and polyene branch. Although both positive and negative values are permitted for the dihedral angle in vacuo, the asymmetric protein environment allows only negative values (Fig. 5). In addition, the distortions of the β-ionone ring to the polyene chain are not symmetric between the 6-s-trans and 6-s-cis conformers. The distortion degree depends on the QM theory and remains controversial, with values of −60.4° reported for the HF geometry49, −68.2° for CASSCF40, and −33.9° for B3LYP50.

The potential energy surface along the C6–C7 dihedral angle is illustrated in Figure 6B. The CASPT2 calculations were based on structures optimized using B3LYP. The overall potential shapes of B3LYP are similar to those of CASPT2//B3LYP. It agrees with the previous work by Page et al.13 showing that in retinal model B3LYP reproduced the result of CASPT2, which deals with static and dynamic correlations in a balanced manner. If dynamic correlations were neglected, such as by structure optimization by the CASSCF/6-31G* method, the overall shape of the potential surface would be different from that of CASPT2. In fact, the CASSCF relaxes into a 6-s-cis minimum even when starting from the 6-s-trans form, as reported by Wanko et al.11. Although the overall shapes of the potential energy curves between B3LYP and CASPT2//B3LYP are similar, the latter surface has roughness, which is reflected by the discontinuous switch from B3LYP (structure optimization) to CASPT2 (excitation energy calculation).

Figure 6.

(A) Excitation energy profile by CASPT2//B3LYP calculations with respect to C6-C7 dihedral angle. (B) Ground state potential energy surface for the C6–C7 dihedral angle calculated by B3LYP (dash) and CASPT2//B3LYP (solid).

Regarding the equilibrium angle in vacuo, the B3LYP optimization exhibited −37.1° for 6-s-cis, which was consistent with CASPT2//B3LYP, and −170.6° for the 6-s-trans, with a −5° CASPT2 correction. On the other hand, the parameter set used in the present MD simulation, which was taken from a previous study28, generated a 6-s-cis orientation with a dihedral angle of −19° in vacuo and was closer to the planar form. The dihedral angle centered at −44° obtained by the QM/MM optimization in protein (Fig. 5B) is markedly different from that of MM and QM in vacuum, which implies that the retinal-protein interaction essentially determined the orientation of the β-ionone ring, consistent with a previous study by Tajkhorshid et al.28.

Next, we consider the characteristics of the potential barrier. The structural and color tuning models have been discussed assuming a potential barrier symmetric to −90°. However, the peaks were around −100° for B3LYP and −110° for CASPT2//B3LYP, thereby deviating from 90°. Regarding the height of the potential barrier, the energy difference between 6-s-cis and the top of the barrier was 20.4 and 13.9 kJ/mol in B3LYP and CASPT2//B3LYP, respectively. Cembran et al.40 reported that in case of CASPT2//CASSCF, the energy difference between −68.2° (6-s-cis) and −90.0° of the dihedral angle is 2.28 kJ/mol, which is much smaller than that of B3LYP. The torsion parameter used in the present MD simulation, which was taken from a previous study28, generated an energy barrier of 45.0 kJ/mol between 0° and 90° of the dihedral angle, and the steric energy in the planar orientation of the β-ionone ring was replaced by the Lennard-Jones interaction. In the present MD simulation, only seven transitions between 6-s-cis and 6-s-trans were observed for SM during 48.0-ns simulation, whereas more than 50 transitions were observed for p501 and DM. Apparently, the statistical ensemble dose not sufficiently evaluate the precise free energy difference between the two conformers. To estimate the free energy difference, enhanced sampling techniques such as umbrella sampling51 and Replica Exchange Molecular Dynamics (REMD)52, may be helpful.

The CASPT2//B3LYP profile for the PSBR excitation energy (Fig. 6A) exhibits a maximum peak at −95°, which corresponds to an almost orthogonal β-ionone ring conformation, while the peak shifts toward larger angles on the ground state surface (Fig 6B). In this regime around −90°, this observation indicates that the shape of the excitation energy profile is predominantly determined by the shape of the excited-state energy profile and not by the ground-state energy profile. Regarding the double minimum on the ground-state surface, the calculated optical absorption wavelength values were 548 (6-s-cis) and 554 nm (6-s-trans). The latter is red-shifted by 6 nm (0.04 eV) relative to the former. It is worth comparing with the previous CASPT2//CASSCF study for all-trans reporting that 6-s-cis conformer has excitation energy of 0.27 eV larger than that of 6-s-trans53. In the present study, the difference in the excitation energy between two conformers is much smaller than the previous study.

Then, we cannot definitely attribute the blue shift of DM to the isomerization in C6–C7, despite the clear correlation between the dihedral angle and excitation energy, because there remains an uncertainty regarding the accuracy, which is affected by the incoherent switching from B3LYP to CASPT2//B3LYP. In fact, the single-peaked distribution (Fig. 3) of DM’s excitation energy obviously differs from the double-peaked distribution (Fig. 5) of its dihedral angle.

Although this section is devoted for the C6–C7 dihedral angle, the retinal molecule exhibits various configurations in the statistical ensemble, where excitation energy should be markedly different from one another. For more comprehensive analysis in the multidimensional conformational space of the rhodopsin molecule, we employed the principal component analysis and identified major determinants of the color tuning mechanism. (See Structural determinant for color tuning: Principal Component Analysis).

Retinal-counterion interaction

Figure 7, A–D show distance distribution between the nitrogen atom (N16) in the Schiff base and the oxygen atoms (the closer atom is Oε1 and farther atom is Oε2) of the counterion, which is, Glu-113, obtained from MD simulation (Fig. 7A) and QM/MM optimization (Fig. 7, B–D), respectively.

Figure 7.

Distance distributions between the nitrogen atom of the Schiff-base and oxygen atoms in Glu-113. (A) MD simulation. The closer (solid) and farther (dashed) oxygen atoms in Glu-113. The histograms are based on the trajectories extracted every 1.0 ps from the MD simulations for p501 (black), SM (red), DM (green). (B–D) B3LYP//CHARMM optimization. The closer (solid) and farther (dashed) oxygen atoms in Glu-113 for p501 (B), SM (C), DM (D).

On the whole, the average distance between the Schiff-base nitrogen atom and the closest oxygen atom of Glu-113 was about 0.8 Å shorter than that observed in the X-ray structure of bovine rhodopsin. As previously shown by Coto et al.41, a 0.1 Å increase in the distance results in a decrease of about 4.3×10−2 eV in the excitation energy, leading to a red shift. Thus, a 0.8 Å proximity generates a 3.4×10−1 eV increase in the excitation energy, which accounts for a major part of the 4.9×10−1 eV overestimation of the excitation energy during simulation. It is quite possible that the proximity is caused by the CHARMM force field parameter, as a previous study also reported an extremely blue-shifted absorption wavelength, in which the MM layer of the geometry was due to the CHARMM force field, and CASPT2 was used for the excitation energy calculation54. As shown in Figure 7A, the proximity of the counterion can actually be observed in MD simulation using the CHARMM force field. The rest 0.9×10−1 eV of the overestimation in the excitation energy probably caused by the missing polarization effect. In the present study, the retinal environment do not redistribute fixed point charges in response to the retinal polarization upon photoexcitation. In the case of bacteriorhodopsin, the estimates for the decrease of excitation energy due to the polarization effect are 0.09–0.19 eV with polarization force field55–56, and 0.07 and 0.05 eV with DFTB57 and PBE058, respectively.

In Figure 7A, the distance distribution of DM was similar to that of p501, whereas it was markedly different from that of SM, where the Oε1 atom moved away from the Schiff base, while the Oε2 atom approached closer to it. In Figure 7, B–D, the distributions of p501 and DM obviously differ from each other, as well as from that of SM.

The interactions between the two pairs of atoms N16–Oε1 and N16–Oε2 were evaluated at the MM level in the B3LYP//CHARMM optimization, because the counterion was not included within the QM layer during optimization. Notably, optimizations with the same parameters for the three pigments made a difference to their distance distributions, although MD simulations did not. This implied that QM treatment of the retinal generated a structural difference around the retinal, which allowed for the proximity of the counterion.

After QM/MM optimization, the average distance of the Oε1 and Oε2 atom of DM became shorter than that of p501 by 0.02 Å and 0.06 Å, respectively. It is possible that the retinal counterion interaction is responsible for the blue shift in DM spectrum, however, this small distance difference does not fully account for the excitation energy shift of DM. In addition, no difference in λmax could be observed between p501 and SM, despite the different patterns of their distance distributions. For SM, it is likely that the separation effect of the Oε1 atom and proximity effect of the Oε2 atom probably cancel out. To verify more specifically the contribution of this interaction, we performed principal component analysis as described below.

Protein environment

Before discussing the protein environment’s effects on the excitation energy, it is important to understand the protonation state problem at Glu-181 and the features of the computational method. Although Glu-181 has been inferred to be responsible for the rhodopsin photoreaction59, its protonation in the dark state remains an ongoing debate21–23. According to the MD simulation by Rohrig et al., in term of the structural stability the deprotonated model is preferable at least at MM level21. On the other hand, most previous QM calculations for excitation energy support the protonated model because the protonated model reproduces well the experimental λmax. Therefore, from the neutral point of view, the deprotonated models were employed in the present MD simulation to refine the homology model, and the optimization process adopted ONIOM-ME where the model dependency on MM layer can be suppressed, because ONIOM-ME does not take into account the polarization effects on the retinal by the protein differently from ONIOM-EE.

As mentioned, the polarization effects of the protein environment were not the major factors for the λmax shifts reproduced in the present study. However, the effects of the protein environment were incorporated through steric interactions, which affected the configurations of the retinal and its counterion. Actually, the orientations of amino acid residues making up the retinal pocket in DM are different from those of the other two pigments (Fig. 4). This implies that the orientations of these amino acid residues are susceptible to mutational information at the amino acid site 195. However, the effective structural difference around site 195 among the three pigments could not be definitively identified because considerable structural fluctuations smeared out subtle differences between the pigments. Much longer MD simulations are required to decrease the statistical error.

When comparing λmax values from the calculations with and without external charges, incorporating the protein external charges generates red shifts of 0.04 to 0.07 eV, consistent with a study by Coto et al.41, although the shift magnitudes are smaller. Therefore, except for the counterion, the protein environment decreases the total excitation energy.

Structural determinants for color tuning: Principal component analysis (PCA)

As mentioned previously, the excitation energy is sensitive to various correlated structural changes. To identify the effective structural determinants involved with the λmax shifts of DM, we used PCA. PCA is conventionally applied to mode analysis where a variance-covariance matrix is constructed with the degrees of freedom of protein, and is subsequently diagonalized. The use of PCA in the present study is somewhat unconventional in that the excitation energy is considered as an additional dimension for the protein conformational space. They no longer correspond to the actual oscillations in the protein.

Equations

For a system consisting of N atoms in three-dimensional space, the variance-covariance matrix is defined as follows, where M is the number of configurations, X (i=1, 2, ..., 3N) is the atomic coordinates of the m-th configuration, and X̄ is the statistical average of X.

In this study, a variance-covariance matrix was constructed in terms of 3N-6 internal coordinates (i.e. bond length, angle, dihedral angle) of retinal and counterion. (Table S1). Thus, the total number N of heavy atoms became 25. In addition, the excitation energy was also considered as an extra “coordinate” to examine the correlation between geometry and the absorption wavelength. Thereby the variance-covariance matrix was constructed in the 70-dimensional space spanned by the excitation energy (n=1), 24 bond lengths (n=2, 3, ..., 25), 23 angles (n=26, 27, ..., 48), and 22 dihedral angles (n=49, 50, ..., 70). We used 100 molecular configurations for each pigment, and the total number of configurations, M, became 300. Another expression for is where In this expression, contributions from each internal coordinate were normalized equivalently, except for i=1 where the are substituted by σ/w. Here we introduce a weight factor w (in this study we set w=5) so that the correlation between geometry and the absorption wavelength is illustrated well with the first few principal components. Next, the matrix was diagonalized as: where U is an orthonormal transformation matrix, whose jth column vector is the jth principal axis, and λ is the i-th eigenvalue of .

For the mth configuration, the principal components p(m) is given by, Then, a element u of vector is a coordinate corresponding to the excitation energy. We transformed the first two axes and into and via the linear combination as follows where represents the direction of steepest-descent and does that of the zero-gradient for the excitation energy.

Results

The weight factor and linear combination was used so that the fast principal mode would be most responsible for the excitation energy in Eq. (3) and (6). Figure 8 illustrates elements (black line Fig. 8) and the structural shift (red line Fig 8) of the unit vector along the first principal axis. The pattern of bond length alternation (BLA) demonstrates Schiff base proximity has the larger effect on the shift with the smaller structural change in contrast to the β-ionone ring proximity (Fig. 8A). Note that the excitation energy shift is accompanied not only the conjugated chain but also by C5-C4 a middle part of the β-ionone ring and C6-C1, a branch of the conjugate chain, as if BLA is recognized up to C5-C4 also for C6-C1, although their length are same as that of typical single bond. For the bond angles (Fig. 8B), the role the Schiff base proximity is more conspicuous, where the first four angles are essential, whereas the overall contribution to the shift is smaller than from BLA. Dihedral angles also exhibit alternation pattern (Fig. 8C). From the Schiff base to the 11-cis bond the dihedral angles around the double bonds have positive correlation with the shift, in contrast to those around single bonds. Although the large torsion at the C6-C7 angle for DM is recognized (red line in Fig. 8C), the contribution to the shift is relatively small (black line in Fig. 8C). For more detailed analysis of the effect of isomerization between 6-s-cis and 6-s-trans, 300 points in the 70-dimensional space were projected onto the plane spanned by the first two principal axes (Fig. S1). DM configurations clustered into two groups along the vertical axis. Since, the two clusters of DM in Figure 5 corresponds to the cis and trans isomers, because the vertical axis represents the distortion of the dihedral angle C6-C7 (Table 2S). Therefore the C6-C7 isomerization is excluded from a direct factor of DM blue shift. However it is possible that the isomerization may be involved indirectly in the color tuning through couplings with other mode.

Figure 8.

Correlation between the excitation energy and retinal backbone geometry: (A) bond lengths (B) bond angles, and (C) dihedral angles. The shift of excitation energy was expressed by a unit vector in the plane spanned by the first two principal components. For this multidimensional vector, the elements of internal coordinates are indicated by the labels on the horizontal axes in A, B and C. The component of the vector for each internal coordinate is shown by the black line (labels are indicated along the vertical axis on the left had side), while the corresponding change in the internal coordinate from the average value is shown by the red line (labels are indicated on the right hand side).

The effects of counterion on the shift are demonstrated by the vector components for bond lengths, N16-Oε1 (−0.044) and N16-Oε2 (0.024), bond angles, C15-N16-Oε1 (−0.029) and C15-N16-Oε2 (−0.018), and dihedral angles, C14-C15-N16-Oε1 (−0.029) and C14-C15-N16-Oε2 (0.019). Obviously, the movement of Oε1 against Schiff base is more dominative than that of Oε2. Then, the proximity of the counterion to the Schiff base proton causes blue shift in agreement with the widely accepted view of the counterion effect on the shift.

We conclude that the retinal-counterion proximity and the increase of BLA causes cooperatively the blue shift of DM.

Spectral broadening

In general, optical absorption spectral lines are broadened as a result of contributions from both homogeneous and inhomogeneous components. Although the two components cannot be distinguished precisely, the former are mainly related to the finite lifetimes of the excited states, whereas the latter results from the difference in site energy attributed to protein-chromophore interaction, namely, solvation heterogeneity responsible for color tuning15. In the case of bovine rhodopsin, the optical absorption spectrum has a full width at half-maximum (FWHM) of 4,250 cm−1 (0.53 eV) at 77 K60. Based on hole-burning spectroscopy at room temperature, the contribution of homogeneous broadening to FWHM was estimated to be 340 cm−1 (0.042 eV), while that of inhomogeneous broadening to be 1,800 cm−1 (0.22 eV)15.

Because of the difference in site energy, the statistical ensemble is indispensable for the precise estimation of λmax. However, most computational studies have not directly addressed this issue because an accurate excitation energy calculation requires considerable computational resources, even for a single point calculation.

In 2006, Hoffmann et al. performed OM2-MRCI calculations in parallel with QM/MM dynamics simulation for bacteriorhodopsin and phoborhodopsin. They reproduced the line width of 0.26 eV for bacteriorhodopsin using 1000 configurations obtained by a 100-ps MD simulation14. Later, CASPT2 calculations were performed for 10 configurations of photoactive yellow protein, and the structural origin of the static disorder was analyzed in detail61.

In the present study, we applied CASPT2 calculations for 100 different conformations obtained from 48-ns MD simulations for p501, SM, and DM and analyzed the static disorder. The FWHMs of the Gaussian fitting curves for p501, SM, and DM were 0.12, 0.15, and 0.11 eV, respectively (Fig. 3).

The inhomogeneous broadening effect is known to depend upon temperature, for instance 0.09 eV at 1.5 K, and 0.22 eV at 298 K15. In the present study, 48-ns MD simulations for structural sampling at 300 K were followed by structural optimization. The coverage of the conformational space explored within a finite MD simulation is limited, and optimization decreases the variety in a sampled structure. Therefore, the narrower FWHMs of three pigments compared to experimental results are possibly due to the limited sampling space.

Conclusions

We investigated the long-range synergetic color tuning mechanisms in an ancestral conger rhodopsin, p501, at the atomic level using homology modeling, molecular dynamics simulation, and electronic state calculations. The molecular dynamics simulations demonstrated a remarkable distortion of the β-ionone ring in the double mutant, N195A/A292S. Additional analysis illustrated that the retinal structural change was coupled to the spatial rearrangement of atoms comprising the retinal pocket, indicating that the effect of mutation far from the retinal was facilitated by another mutation near the retinal and somehow propagated via the retinal pocket to the chromophore.

In addition, we performed electronic state calculations for p501, A292S, and N195A/A292S. Using statistical ensembles of QM/MM optimized structures, the intrinsic excitation energy of retinal chromophore was evaluated. PCA was used to identify the determinants for the color tuning mechanisms of the visual pigments: bond length alternation and counterion displacement.