Pro219 is an electrostatic color determinant in the light-driven sodium pump KR2

Color tuning in animal and microbial rhodopsins has attracted the interest of many researchers, as the color of their common retinal chromophores is modulated by the amino acid residues forming the chromophore cavity. Critical cavity amino acid residues are often called “color switches”, as the rhodopsin color is effectively tuned through their substitution. Well-known color switches are the L/Q and A/TS switches located in the C and G helices of the microbial rhodopsin structure respectively. Recently, we reported on a third G/P switch located in the F helix of the light-driven sodium pumps of KR2 and JsNaR causing substantial spectral red-shifts in the latter with respect to the former. In order to investigate the molecular-level mechanism driving such switching function, here we present an exhaustive mutation, spectroscopic and computational investigation of the P219X mutant set of KR2. To do so, we study the changes in the absorption band of the 19 possible mutants and construct, semi-automatically, the corresponding hybrid quantum mechanics/molecular mechanics models. We found that the P219X feature a red-shifted light absorption with the only exception of P219R. The analysis of the corresponding models indicate that the G/P switch induces red-shifting variations via electrostatic interactions, while replacement-induced chromophore geometrical (steric) distortions play a minor role. However, the same analysis indicates that the P219R blue-shifted variant has a more complex origin involving both electrostatic and steric changes accompanied by protonation state and hydrogen bond networks modifications. These results make it difficult to extract simple rules or formulate theories for predicting how a switch operates without considering the atomistic details and environmental consequences of the side chain replacement.


Supplementary Note 2. Phase I: Input File Generator of a-ARM
The first phase of a-ARM (Figure 4b) allows either the automatic or semi-automatic computeraided preparation (i.e., via a five-step command-line procedure) of a 3D structure in PDB format (without hydrogens), which contains information on: the monomeric protein structure, automatically selected as Chain A but user-customizable, including the retinal proton Schiff base (rPSB) chromophore and excluding membrane lipids and non-functional ions (step 1, Section 2.2.1 of Ref. [1]); the mutant(s) automatically produced via side-chain replacement using MODELLER [2] (step 3); the protonation states for all the ionizable residues based on an algorithm that analyze pK a and partial charges using PROPKA [3], automatically assigned but user-customizable (step 4, Section 2.2.3 of Ref. [1]); the positions of Cl − /Na + external counterions needed to neutralize both IS and OS, based on an energy minimization procedure using PUTION [4], optimized automatically and not user-customizable (step 5, Section 2.2.4 of Ref. [1]); and an independent file containing the list of amino acid residues forming the cavity hosting the rPSB, determined automatically with Fpocket [5] but user-customizable (step 2, Section 2.2.2 of Ref. [1]). The resulting PDB structure (PDB ARM ) plus cavity file constitute the so-called a-ARM input, for the QM/MM model generator phase (see below), that according to the parameters selected (via command-line) in steps 1-5 can be considered as a-ARM default or a-ARM customized . While the former refers to a fully automatic input generation, which uses default parameters as suggested by the code (see above), the latter allows the command-line assisted user-customization of some of them.

-Model customization
The a-ARM customized customized approach, that is specifically used in cases where the default choices produce QM/MM models that are not suitable for the reproduction of trends in absorption properties, can be performed by adopting well-defined guided-procedure easily replicable. This protocol, documented in Refs. [1] and [6], involves three phases which only concern the ionization states of the ionizable residues: i) when the pH is > 6 we modify the ionization states by setting the pH to 5.2 in step 4 and ii) we check the protonation state of the main and secondary conterions of the rPSB, and if our analysis give them both ionized we neutralize the secondary conterion and iii) in case the model generated in step ii does not reproduce the experimental absorption maxima, then we exchange the secondary and main counterion ionization states. Supplementary Figure 3 illustrates the customization of the WT-KR2 model, described in the main manuscript, in which steps i) and ii) were employed. Table 2 Calculated pKa computed with PROPKA3.1. pKa computed at pH 5.2 and 7.0 and partial charges obtained in "step 4" of the input file generator of a-ARM. The only residue sensitive to pH change is Glu160. In the PROPKA calculations, the retinal chromophore is not included in the PQR file.

Supplementary Note 3. Phase II: QM/MM model generator of a-ARM
The second phase (Figure 4c) allows the automatic generation of ground-state (S 0 ) QM/MM models for rhodopsins ( Figure 4a) and the subsequent computation of the maximum absorption wavelength (λ max a ) via vertical excitation energy (ΔE S1-S0 ) calculations. The procedure is described as follows: -Classical molecular dynamics simulations a-ARM input is pre-processed using classical Molecular Dynamics (MD) simulations. First, the positions of crystallographic/comparative waters are optimized and the hydrogens for waters and polar residues are added by using DOWSER [7]. Then, the hydrogens for the rest of the protein and chromophore are added and their positions are optimized by a Molecular Mechanics (MM) energy minimization using GROMACS [8]. A second MM energy minimization is performed, this time on the side-chains (backbone atoms are fixed at the crystallographic/comparative positions) of the residues belonging to the chromophore cavity sub-system. The resulting structure is employed as an input to generate N=10 independent simulated annealing/MD relaxations at 298 K, each starting with a different randomly chosen seed to warrant independent initial conditions that allow to explore the possible relative conformational phase space of the cavity residue sidechains and chromophore. In the ARM MD approach, that uses GROMACS [8] and AMBER [9] force field, all side-chains of the Lys-QM and chromophore cavity (including cavity waters) subsystems are relaxed, while the backbone is fixed at the crystallographic/comparative structure. The Lys-QM subsystem is described by using a MM parametrization and partial charges computed as AMBER-like Restrained Electrostatic Potential (RESP) charges, which are specifically parametrized for each employed isomer of the chromophore (e.g., 11-cis, all-trans and 13-cis rPSB). [4] Such parameters are reported in the Supporting Information of Ref. [1]. Moreover, the default heating, equilibration, and production times for the MD (selected via benchmark calculations in Ref. [4]) are 50, 150, and 800 fs, respectively, for a total length of 1 ns. For each of the N=10 replicas, the frame closest to the average structure of the 1 ns simulation is selected as the starting geometry (i.e., guess structure) for constructing the corresponding QM/MM model.

-QM/MM calculations
Each of the 10 replicas (i.e., frame extracted from the N=10 independent MDs) is processed by a particular QM/MM approach implemented into the [Open]Molcas/TINKER [10,11]] interface [12], where the electrostatic embedding scheme used to describe the interaction between the QM and MM parts of the Lys-QM sub-system (see Figure 4a), involves an unconventional treatment called Electrostatic Potential Fitted (ESPF) [13]. In the ESPF method, the QM part of the chromophore directly interacts with the MM electrostatic potential through one-electron operators whose expectation values represent the QM charge distribution of the chromophore. Notice that S8 mutual polarization effects between the QM and MM sub-systems are not considered. Although this issue can, in principle, be solved by employing a polarizable embedding method, we have not adopted/benchmarked such technologies in this first version of our specialized QM/MM models since they are still under development in the QM/MM area (see for instance Ref. [14]). In addition, the QM/MM frontier is treated within a link atom approach whose position is restrained according to the Morokuma scheme, and it is placed across the covalently bonded lysine Cε-Cδ bond (where Cε is a QM atom). The charges of the covalently linked lysine are modified by setting the Cδ charge to zero to avoid hyperpolarization and redistribute the residual fractional charge on the most electronegative atoms of the lysine, thus ensuring a +1 integer charge of the Lys-QM layer. All the 63 Lys-QM atoms (i.e., 62 atoms + linker atom) are free to relax during the QM/MM calculation. By employing such an approach, the procedure to obtain an ARM QM/MM model consisting of N=10 replicas, can be described as follows. First, to complete the pre-processing step, a geometry optimization at the Hartree-Fock (HF) level is performed (HF/3-21G/AMBER). Then, another geometry optimization is carried out this time modeling the QM sub-system with the multiconfigurational complete active space self-consistent field (CASSCF) at the 2-roots single-state, (CASSCF(12,12)/6-31G(d)/AMBER level). This follows an energy correction at the multiconfigurational second-order perturbation theory (CASPT2) to recover the missing dynamical electron correlation associated with the CASSCF description. Thus, a 3-roots state-average CASPT2 that uses the three-root stage-average CASSCF(12,12)/6-31G(d)/AMBER as the zeroorder reference wavefunction, is computed (CASTP2(12,12)/6-31G(d)/AMBER). Ultimately, each model replica corresponds to an equilibrated gas-phase and globally uncharged monomer QM/MM model and it is associated with a calculated between S 0 →S 1 ΔE S1-S0 . The final a-ARM result is the average of the 10 ΔE S1−S0 values. A detailed explanation of the a-ARM protocol workflow is provided in Refs. [1] and [6]. Total Energies calculated at the CASSCF/AMBER//CASPT2/6-21G(d) level. First vertical excitation energy (∆ES1−S0), maximum absorption wavelength (λ a max), transition oscillator strength (fOsc), and difference between calculated and experimental data (ΔΔE S1-S0 Exp,a-ARM ). Standard deviation for the N=10 replicas (σN) is presented as subindex Experimental a-ARM (N=10) a-ARM (N=1) replica with ΔE S1-S0 a-ARM closest to the average Exp,a-ARM CASPT2 S0 CASPT2 S1 ΔE S1-S0 a-ARM λ max a,a-ARM fOsc ΔΔE S1-S0

Supplementary Note 5. Comparison of WT-KR2 models built with original ARM
Although the conclusions derived from the study of Inoue et. al. [33] can be qualitatively compared with the results presented in this work, notice that they cannot be quantitatively compared since i) the QM/MM models were constructed from a different X-ray structure (i.e., 3X3C), using the original ARM [15] that featured ii) manual input file generation (i.e., handmade and not reproducible counterion placement, different chromophore cavity, etc.) and iii) a different methodological approach for the generation of the mutant side-chain. Remarkably, we have verified that the side-chain conformation of the P219T mutant selected automatically in a-ARM  Table 4 demonstrates that although the computed magnitudes are different, the signs are the same, suggesting that both versions of the protocol produce consistent results.
As has been established in Ref. [1], the MAE is computed inside the a-ARM framework to estimate the ability of the protocol to predict photophysical properties and then be able to compare the results in trends for heterogeneous sets of rhodopsin variants. For instance, in the main text we have compared the results in vertical excitation energy obtained for the benchmark set in Ref. [1], with the results obtained in this work for the KR2 set. -

Mean absolute deviation
As has been established in Ref. [1], the MAD is computed inside the a-ARM framework as a measure of dispersion, that represents how much the absolute errors between computed and experimental values in photophysical properties in the data set (! ! , see definition above) are likely to differ from their MAE. The absolute value is used to avoid deviations with opposite signs cancelling each other. The MAE is calculated by using the following formula: Where n is the number of rhodopsin variants. See Section 6.2.1.

-Weighted average
The weighted average ($̅ ) takes into account the varying degrees of importance of the numbers in a data set. It is equal to the sum of the product of the weight (wi) times the data number (xi) divided by the sum of the weights:

Supplementary Note 7. Parallelism between computed and experimental values: weighted average
In order to evaluate the parallelism between calculated and experimental data, we have computed the trend deviation defined in the previous section as where Δ max,X WT,Exp is the difference between experimental λ max of each of the P219X mutants with respect to the experimental value of the WT, while Δ max,X WT,a-ARM is the difference between a-ARM computed λ max of each of the P219X mutants with respect to the a-ARM computed value of the WT.
According with both Δmax,X WT,Exp and Δmax,X WT,a-ARM values, reported in Supplementary Table 5 and Figure 5d, we have classified the 19 P219X KR2 variants into two clusters, namely "blueshifting" and "red-shifting" clusters. The former is composed of the only blue-shifted variant, P219R, while the latter is composed of the other 18 variants. The main purpose of such cluster definition is to threat the data of variants of a same cluster in a weighted average fashion, instead of discussing individual values. We felt that such a treatment is necessary since, as explained in the main manuscript, the <2.0 kcal mol -1 observed Δmax,X WT,Exp variations among members of the red-shifting cluster includes, in most cases, too small (we set a threshold at ≤1 kcal mol -1 ) for the Δmax,X WT,Exp trend to be safely reproduced by a a-ARM QM/MM model. In order to locate the cluster center, we selected a weighted average to give more importance to the most frequent deviations seen in our small data sample. Accordingly, the weighted average of the red-shifting cluster, for both computed and experimental data, was computed following the next procedure (notice that the results are not substantially different from the results obtained with a standard average. See Supplementary Table 5.) -Computation of the individual weights for the red-shifting cluster: In order to compute the individual weights for each of the members of the red-shifting cluster, a histogram for both Δ max,X WT,Exp and Δ max,X WT,a-ARM was generated, using the data (in kcal mol -1 ) presented in Supplementary Table 5. Then, the entire range of values was divided into a series of intervals, using a bin width of 0.2 kcal mol -1 , and the weights were calculated as the frequency of the data in the corresponding interval.
Compute the weighted average value for the red-shifted cluster. By using the data of the histograms of Supplementary Figure 5, we assigned a weight for both the Δ max,X WT,Exp and Δ max,X WT,a-ARM values for each of the members of the red-shifting cluster, as reported in Supplementary Table  5. Finally, we computed the weighted average of the red-shifting cluster, for both experimental and computed data, by using the equation shown in Supplementary Note 6, as -1.38 and -0.97 kcal mol -1 , respectively. Since the definition of such weighted average values mean the average difference in vertical excitation energy between the center of the red-shifting cluster and the wild type (WT-KR2), we add such values to the experimental and computed data of the WT-KR2. In this way, we obtained the coordinate [538.3,519.8] that we plot as black point in Figure 5b.      Table 6 Weighted average values for both experimental and computed vertical excitation energy, as well as for the different components, for the red-shifted cluster. The standard average is also reported. Values are presented in kcal mol -1 .

Supplementary Note 10. Structural comparison between P219X mutants and WT-KR2
Supplementary Figure 6 3D structures for the WT-KR2 and its P219X (X= A, C, D, E, F, G, H, I, K, L, M, N, Q, R, S, T, V, W,  Y) variants. The chromophore and the covalently linked Lys are presented as green and blue balls, respectively, whereas each X mutant is presented as orange sticks. Purple shadow illustrates the reference WT structure.

Supplementary Note 12. Limitations and pitfalls of a-ARM
In spite of the encouraging outcome of the photochemical studies based on a-ARM, additional work is necessary to generate a tool that can be systematically applied to larger arrays of rhodopsins. In this section, we describe the current methodological limitations of a-ARM in terms of both input file generation and QM/MM model building.
Regarding the input generation, three main issues have to be tackled: o Assignment of the protonation states: There are two aspects which limit the confidence in the automation of the ionizable state assignment described above. The first is that, due to the fact that the information provided by PROPKA [3] is approximated (i.e., the retinal chromophore is not included in the PQR file), the computed pK a Calc value may, in certain cases, be not sufficiently realistic. The second aspect regards the assignment of the correct tautomer of histidine. This amino acid has +1 charge when both the - o Automatic prediction of side-chain conformation for mutants: In order to achieve a successful technology for systematically predicting mutant structures, a level of accuracy of the a-ARM models superior to the one currently available is needed. To deal with that, in this work we attempt at the improvement of the mutations routine by replacing SCWRL4 (i.e., a backbone-dependent rotamer library) with MODELLER a software based on comparative modeling.
Accordingly, their simplified definition makes the ARM models more exposed to potential pitfalls with respect to more complex QM/MM models. Such possible pitfalls can be summarized as: 1) lack of a proper description of the protein environment (membrane + explicit solvent), 2) rigid protein backbone and non-cavity side-chains, 3) approximated protonation states for ionizable residues, 4) missing description of any mutual polarization effects between the QM and MM sub-systems, that can be accounted for by polarizable embedding using a polarizable force field. Since polarizable force fields are technologies still under development in the QM/MM area (see for instance Ref. [14]), we have not adopted/benchmarked them in this first version of our specialized QM/MM models.
Considering points 1-4 above, the different properties computed by ARM are expected to be affected by a systematic error. The current research of our research group is aimed, also, to overcome these points, while maintaining reasonable computational costs, or estimating the errors due to them. Nevertheless, according to the philosophy of the ARM protocol, the main focus of ARM is the ability to reproduce observed trends in vertical excitation energies, as specified in section "Validation, capabilities, and potential applications of a-ARM" in the manuscript main text and illustrated in Supplementary Figure 7.

Supplementary Figure 7 a-ARM protocol validation. (A) Computed excitation energies
ΔE S1-S0 in both kcal mol −1 (left axis) and eV (right axis) for a set of 44 rhodopsin variants. The employed protein structures were obtained from X-ray crystallography (left panel) or through comparative (homology) modeling (center panel). Two sets of variants for bovine rhodopsin (Rh) and bacteriorhodopsin (bR) are also reported (right panel). All computed data were obtained using the a-ARM default approach (blue up-turned triangles), and specific models were refined with the a-ARM customized (gold squares) approach. Experimental data, as energy difference corresponding to the wavelength of the absorption maxima, are also reported (red down-turned triangles). (B) Differences between computed and experimental excitation energies ΔΔE Exp S1-S0 in both kcal mol −1 (left axis) and eV (right axis). The trend deviation of the set, obtained after customization, is 0.7 ± 0.5 kcal mol −1 (0.03 ± 0.02 eV) and the mean absolute error (MAE) is 1.0 kcal mol −1 (0.04 eV). Reproduced with permissions from [16]. Copyright 2020 Wiley Online Library.

Supplementary Note 13. Modification of the mutation routine (step 3) of a-ARM
It is well-known that the success of in silico modeling of point mutations in proteins, relies on the selection of a robust methodology for the prediction of the side-chain conformation of the replaced amino acid. [17,18,19,20,21,22,23,24,25,26] [27,28,29] Both original [4] and updated [1,12,6] versions of the ARM protocol employ the software SCWRL4 [30] to predict the side-chain conformation of the mutated residues. Such a prediction is based on backbone- bacteriorhodopsin (bR) [32] and KR2 rhodopsins. [33] The current research work is, however, the first attempt to use the a-ARM rhodopsin model building protocol to systematically mutate a single residue with each of the remaining 19 essential amino acids. More specifically, we performed single point mutations of the residue P219, in the KR2 rhodopsin, that is located near the β-ionone ring of the rPSBAT (see Accordingly, and also consistently with the structurally "conservative" approach of the a-ARM protocol where our models are designed to retain information from the X-ray crystallographic or comparative structures, our methodology replaces only the side-chains of the mutated residues keeping the backbone atoms at fixed positions. To this aim, the optimization of the mutated side-chains is obtained using a combined approach which alternates conjugate gradient minimizations and short MD mod simulations with simulated annealing (for further details see Ref. [2]). As described by Feyfant et al. [34], this intends to minimize a scoring function including homology-derived restraints, force field energy terms (CHARMM22), and a statistical potential for non-bonded interactions. Notice that the MD mod used in this step differs from the MD employed in the QM/MM model generator phase that is described in Supplementary Note 3.
In the above procedure, the script mutate_model.py uses a default initial condition or "seed" (variable rand_seed= -49837) for the MD mod simulation. Therefore, since Modeller is deterministic, if such seed value is not modified the MD mod run will always produce the same side-chain conformation when a certain template is used as input. In order to sample more extensively the conformational space of a mutated residue and evaluate its effect on the vertical excitation energy (∆# !"#!$ % ), our customized setup produces multiple rotamers of the same mutated side-chain by providing the script with different initial seeds (i.e., initial velocities) for Supplementary Figure 9 illustrates the procedure for selecting the rotamer from the three evaluated models. Furthermore, Supplementary Figure 10 shows the performance of all the possible models generated for the P219X set, while Supplementary Figure 11 reports the 20 models selected for the analyses on color tuning presented in the manuscript.  Generate the a-ARM QM/MM models for each of the 3 rotamers and compute the vertical excitation energy ( E a ARM S0 S1 ) Compute, for the wild-type (WT) and each rotamer (rotX), the di↵erence between experimental and computed values: 1) Exp variant E a S0 S1 = E Exp S0 S1 -E a ARM S0 S1 2) rotX = ( Exp rotX E a S0 S1 -Exp W T E a S0 S1 ) Select the rotamer that features the lowest rotX value (preferring, blueshifted values) as the representative ARM QM/MM model have the same protonation states than the WT template and their side-chain is modelled with the rotamer 1 (dark orange bars). The side-chain of the customized models could be modelled with the rotamer 2 (green bars) or the rotamer 3 (light orange bars). Most of the customized models exhibit the same protonation states than the WT, with exception of P219K and P219R marked with a star.
-Limitations and pitfalls of side-chain predictor o Insufficient description of possible cavity rearrangements after mutation: The procedure for modeling the side-chain conformation comprises a short MD mod , where the produced side-chain is allowed to relax, whereas the rest of the cavity residues, waters, chromophore and protein environment remain fixed at crystallographic/comparative structure. Therefore, possible local steric/electronic rearrangements of the residues of the chromophore cavity surrounding the mutated residue are not correctly described.
Although during the QM/MM generation phase of a-ARM the geometry of this sidechain along with the side-chain of the residues in the chromophore cavity are refined via a more sophisticated Molecular Dynamics (MD) (see Section S1.2.1), in some cases this step would not be sufficient to achieve a proper description of the impact of the new side-chain on the protein environment.