a-ARM: Automatic Rhodopsin Modeling with Chromophore Cavity Generation, Ionization State Selection, and External Counterion Placement

Abstract

The Automatic Rhodopsin Modeling (ARM) protocol has recently been proposed as a tool for the fast and parallel generation of basic hybrid quantum mechanics/molecular mechanics (QM/MM) models of wild type and mutant rhodopsins. However, in its present version, input preparation requires a few hours long user's manipulation of the template protein structure, which also impairs the reproducibility of the generated models. This limitation, which makes model building semiautomatic rather than fully automatic, comprises four tasks: definition of the retinal chromophore cavity, assignment of protonation states of the ionizable residues, neutralization of the protein with external counterions, and finally congruous generation of single or multiple mutations. In this work, we show that the automation of the original ARM protocol can be extended to a level suitable for performing the above tasks without user's manipulation and with an input preparation time of minutes. The new protocol, called a-ARM, delivers fully reproducible (i.e., user independent) rhodopsin QM/MM models as well as an improved model quality. More specifically, we show that the trend in vertical excitation energies observed for a set of 25 wild type and 14 mutant rhodopsins is predicted by the new protocol better than when using the original. Such an agreement is reflected by an estimated (relative to the probed set) trend deviation of 0.7 ± 0.5 kcal mol^-1 (0.03 ± 0.02 eV) and mean absolute error of 1.0 kcal mol^-1 (0.04 eV).

Figure 1.

General scheme of a QM/MM ARM and a-ARM model, composed by (1) main chain (cyan cartoon), (2) chromophore rPSB (green ball-and-sticks), (3) Lys side chain covalently linked to the chromophore (blue ball-and-sticks), (4) main counterion MC (cyan tubes), (5) protonated residues GLH and ASH (violet tubes), (6) external Cl^– (green balls) counterions, (7) water molecules (tubes), and the (8) residues of the chromophore cavity subsystem (red frames). Parts 1 and 6 form the environment subsystem. Parts 2 and 3 form the Lys-QM subsystem, which includes the H-link atom located along the only bond connecting blue and green atoms. Parts 4 and 8 form the cavity subsystem. Water molecules (Part 7) may be part of the environment or cavity subsystems. The external OS and IS charged residues are shown in frame representation. This figure, and all other protein structures presented in this work, were produced using PyMOL, version 1.2.

Figure 2.

a-ARM workflow. After the selection of the protein chain, a-ARM generates the ARM input files with complete information on the chromophore cavity, protonation states, and counterion placement (see Figure 1) corresponding to points B–D of section 1. The input is then used for the execution of the original ARM, obtaining as output 10 a-ARM models along with the calculated average vertical excitation energy (ΔE_S1–S0). The parallelograms represent input or output data, the continuous line squares refer to processes or actions, and the dashed lines mean software executions. The [A] mark symbolizes fully automation, whereas the [M] mark represent manual decision. Finally, the [M/A] mark indicates situation that may be either manual or automated (see text). Notice that the software execution labeled “QM/MM calculation” is the same as in the original ARM (see ref 52). In a-ARM the production of the PDB^ARM and cavity input files takes only a few minutes.

Figure 3.

External counterion placement. Schematic representation of the procedure for the definition of the number and type of external counterions needed to neutralize the IS (A) and OS (B) surfaces of bovine rhodopsin. We also illustrate the grid generated by the PUTION code to calculate the coordinates of the Cl^– counterions in the IS (C) and the Na⁺ in the OS (D). The negatively and positively charged residues are illustrated as red and blue sticks, respectively, and the Na⁺ and Cl^– counterions as blue and green spheres, respectively.

Figure 4.

Automatic generation of mutants by using the SCWRL4 software. The code does not require any interaction with the user during execution.

Figure 5.

(A) Vertical excitation energies (ΔE_S1–S0) computed with a-ARM_default (up triangles) and a-ARM_customized (squares), along with reported ARM (circles) and experimental data (down triangles). S₀ and S₁ energy calculations were performed at the CASPT2(12,12)//CASSCF(12,12)/AMBER level of theory using the 6–31G(d) basis set. The calculated ΔE_S1–S0 values are the average of 10 replicas (see Table S3 in the Supporting Information). (B) Differences between calculated and experimental ΔE_S1–S0 ($\mathrm{\Delta }\mathrm{\Delta }{E}_{\mathrm{S}1-\mathrm{S}0}^{\text{Exp}}$). Values presented in kcal mol^–1 (left vertical axis) and eV (right vertical axis).

Figure 6.

a-ARM models. Conformational (the occupancy factor of the rotamers Asp-116 and Gln-157 are presented in parentheses) and ionization state variability for KR2 [PDB ID 3X3C] (A), BPR [PDB ID 4JQ6] (B), RCone (C) [PDB ID template 1U19], bR_AT [PDB ID 6G7H] with standard (D) and modified cavity (E). MC and SC are presented as cyan and violet tubes, respectively.