A Multi-Scale–Multi-Stable Model for the Rhodopsin Photocycle

We report a multi-scale simulation study of the photocycle of the rhodopsins. The quasi-atomistic representation (“united atoms” UA) of retinal is combined with a minimalist coarse grained (CG, one-bead-per amino acid) representation of the protein, in a hybrid UA/CG approach, which is the homolog of QM/MM, but at lower resolution. An accurate multi-stable parameterization of the model allows simulating each state and transition among them, and the combination of different scale representation allows addressing the entire photocycle. We test the model on bacterial rhodopsin, for which more experimental data are available, and then also report results for mammalian rhodopsins. In particular, the analysis of simulations reveals the spontaneous appearance of meta-stable states in quantitative agreement with experimental data.

and Table S1 report the known phenomenology about the BR and MR photocycles. For each state and transition, the absorption wavelengths and transition times are reported. For the first transition the two different times refers to the chromophore isomerization and subsequent relaxation respectively. For BR also the free energies are reported in the table (referred to the BR rest state, the first in the table), measured in different pH conditions. The plots report correlation between these quantities. Figure S1a shows a clear correlation between the barrier height (measured as difference beteween the activated state and the starting one, green and red dots = experimental data in different conditions, see caption) and the transition time, as expected. This correlation is reproduced in simulations (black dots).  Table S1 with their references. Times are in log scale.
(a) (b) Conversely, the correlation between the transition wavelengths and the transition times ( Figure S1b) does not show a clear trend. The optical properties are, in fact, more related to the state of the retinal (cis-trans), of the hydrogen bond network surrounding it and the protonation state of the Schiff base [1][2][3][4][5]. Table S1. Adsorption wave-length, the transition times, and relative energies of each state of the photo-cycle of BR (upper table) and of MR (lower table. Data from [6][7][8][9][10]

S2. Reference States Choice
The reference structures for the photo cycle states were chosen as follows. For each optical state, we searched in the PDB [11,12] repository the available structures. Among these, we chose the best compromise between resolution and completeness. When possible, we also excluded the mutants and choose crystallization conditions as far as possible natural in terms of temperature and illumination. We also needed the same number of residues for each of the cycle photo-state. Thus when necessary the structures were completed using data from other structures of the same state. The structures used and some of their features (PDB code, reference, publication year) are reported in Table S2, and represented in Figure S2. Table S2. List of data used to build the reference structures for the photocycle states. The PDB codes are reported (in case a second structure is used to complete the first, the list of the AA used is also reported), the resolution, and the publication year of the structure.

S3. Force Field Details and Parameters
The explicit form of each of the FF term is the following: The functional forms used for the us depend on the specific bead type and class of interaction, as specified in Tables S3 and S4. The functional forms are chosen among the following: Constraint = holonomic constraint (maintained by the SHAKE algorithm [13,14] The optimized parameters and functional forms of the multi-scale model for single states are reported in Table S4. Table S4. FF parameters of the multi-scale model for BR. The distances (r0 rcut) are in Å, the angles (θ, ϕ) in deg, the costants kθ, Kteth and A in Kcal/mol and for the 12-6 potential a and b are in Å 12 and Å 6 respectively. The value of the rcut for the separation between local and non local part of the non bonded interactions is uniformly set at 8.5 Å.

FF Term Analytical Form Parameters
Opsin constraint

S4. Interface Parameters Optimization Procedure
Interface interactions between the active site and the rest of the protein need a de novo parameterization. These interactions involve Ul, Uθ, UΦ of the retinal-lys side chain linkage, a number of non-bonded Uloc interactions between the retinal and residues facing it. The bonding potentials are parameterized combining the partially biased model and the UA model: the structural parameters (equilibrium distances, angles and dihedrals) are taken from the reference structures, while elastic constants are obtained by the UA parameters.
The Uloc term is represented as a sum of Morse potentials (see previous section). The parameters of this term are fitted on energy profiles evaluated by atomistic simulations using the OPLS FF, as a function of the independent variables, by constraining it to sampled values and optimizing the systems with respect to all the other coordinates. As an example, Figure S2 illustrates the procedure to evaluate the Morse term representing the hydrogen bond between the water and an aspartic acid and between two water molecules. In this case the independent variable is the distance between the beads representing oxygen and the aspartic acid (Cα of asp). As it can be seen from the Figure S3, the Morse potential fits very well the atomistic energy profile. The obtained value for the width parameter α is little less than 2 Å −1 , larger than the other local interactions (see Table S2) as an effect of the deeper and shorter range hydrogen bonding interactions, with respect to other local ones. The value obtained for ε is further optimized. In fact, as an effect of the resolution reduction and of the topology of the local interaction, in our model each hydrogen bond is described by a sum of several Morse interaction. This implies a scale factor to be applied to ε, which was estimated to 0.005 by comparing the multi-scale model with atomistic simulations in the interface region. Figure S3. Fitting procedure to parameterize the hydrogen bonding CG-UA or UA-UA interactions (a,c): atomistic representations of the elements involved in the hydrogen bond (water and aspartic acid in (a) and two waters in (c)). (b) energy profile evaluated with the OPLS FF, as a function of the distance between beads (red dots) and its fit with a morse potential (black line).

S7
The hydrogen bonds formed by water molecule W402, namely with Asp85, Asp212 and Lys216 where treated explicitly and individually, thus in this case no rescaling is needed. In these bonds, water acts as a donor or as acceptor [5]. We added other three potentials to the Uloc term, Uh-bond, described by a Morse potential with ε and α depending on the type of hydrogen bond created according to the procedure described above. In particular if the water molecule acts as a donor, the binding is higher, see Table S2. The final values for the parameters are reported in Table S2.

S5. Multi-Stable Models Building
As explained in the main text, except for the first transition, the multi-stable models are obtained by combining in couple subsequent models for the single states. The model is built in order to reproduce the relative energies of the states taken from experiment, δ. The bi-stable models are composed of the sum of a number of double well potentials, thus the amount of δ must be distributed by these. Table S5 describes how this distribution is performed. First, the number of term treated with the double well form in place of the single well one is limited to those whose equilibrium values difference overcomes a threshold (fixed at 10% for distance depending interactions and to 5% for angles and dihedrals), while a single well is used otherwise. This limits the use of double well potentials to the terms which substantially contributes to the difference of the structure. This criterion is applied to bonding terms and to local non-bonded. It is to be noted that if a given couple i-j is present only in one of the two states, then it is represented with a double well Morse where the second well has its minimum corresponding to the value of the generic non local-non bonded potentials. Then the value of δ is distributed among those double well interactions which are mainly involved in the transition, which are reported in Table S5, namely those included in the active site. The physical meaning of this procedure is that the energy differences between the different states are mainly due to conformational diversity arising from the active site.   Figure S4 reports the energy variation of vs. time during the single state dynamics simulations. The rest state simulation is 2 μs long and the other states are 50 ns long. As noted in the main text, the N state simulation shows a change in conformation of the turn between helix E and F, an "opening",

S8
which results in a change in energy in the dihedral term of 15 Kcal/moL. The duration of the event is 33 ns. All other states do not show noticeable conformational changes during dynamics.

S7. Details of the State Transitions
For the analysis of the different state transitions it is useful to define a reasonable "reaction path" for each of them. For the BR-K and L-M transitions, this can be recognized to be the dihedral angle around the 13th bond in the retinal chain, namely that describing the bond isomerization. For each transition, we determined a series of configurations at fixed reaction coordinates letting all the others S9 to relax. We remark that these configurations approximately interpolate between the starting and final states, and are used only for rough evaluation of the energetics during the transition and test in a simple way the effectiveness of the double well approach. The energy of the intermediate configuration was evaluated in each case with three FFs: that of the starting state, that of the final state and the double well FF connecting the two states. Sample results of energy profiles are reported in Figure S5.
The total potential energies (red) and temperatures (black) along the transition simulations are reported in Figure S6. As noted in the main text, small or null barrier transitions (e.g. K-L or BR*-K) occur fast or spontaneously. Other transitions are induced by forcing the system along the above determined paths (forcing the coordinate to increase smoothly), up to the top or the barrier. The energy profiles are then re-evaluated along the simulations and reported in the main text. During all transitions, the temperature is maintained constant at 300K with the Berendsen thermostat (time constant 1 fs).