Unravelling the Turn‐On Fluorescence Mechanism of a Fluorescein‐Based Probe in GABAA Receptors

Abstract GABAA (γ‐aminobutyric acid type A) receptors are ligand‐gated ion channels mediating fast inhibitory transmission in the mammalian brain. Here we report the molecular and electronic mechanism governing the turn‐on emission of a fluorescein‐based imaging probe able to target the human GABAA receptor. Multiscale calculations evidence a drastic conformational change of the probe from folded in solution to extended upon binding to the receptor. Intramolecular ππ‐stacking interactions present in the folded probe are responsible for quenching fluorescence in solution. In contrast, unfolding within the GABAA receptor changes the nature of the bright excited state triggering emission. Remarkably, this turn‐on effect only manifests for the dianionic prototropic form of the imaging probe, which is found to be the strongest binder to the GABAA receptor. This study is expected to assist the design of new photoactivatable screening tools for allosteric modulators of the GABAA receptor.


Table of contents
Section S1 Selection of the density functional 2 Tables S1, S2 Section S2 Ground state of the fluorescent probe 1 5 Figure S1 Section S3 Emission of the fluorescent probe 1 6 Figure S2 Section S4 Setting up of the systems 7 Table S3 Section S5 MD simulations 8 Figures S3, S4, S5, S6 Binding energy analysis  11  Tables S4, S5 Emission. Emission energies and associated orbital characters are investigated with the best performing functionals for absorption, M06 and B3LYP, [1][2][3] as well as CAM-B3LYP, [4] benchmarked against ADC (2). We first optimized the minimum energy geometry on the S1 PES of 2a and 2b

Section S6
with TD-DFT/def2-TZVP and LR-PCM. As for absorption, the emission energies were calculated with the default LR-PCM implicit solvation as well as with the corrected approach cLR-PCM for B3LYP and M06.
As it can be seen in Table S2, the TD-DFT calculations for 2a is in line with the results of ADC(2), with transition character (π Xan ← π OG * ) and energies by +0.44 eV (similar to the shift seen with the absorption calculations). The errors for the non-emissive 2b are less pronounced: ADC(2) gives a value of 760 nm (1.63 eV), again in reasonable agreement with B3LYP (841 nm, 1.47 eV). We stress that for 2b we found a surprisingly low energy S1 state with ADC(2) (1.63 eV). Out of Table S1. Calculated S1 absorption energies (Eab in eV and nm), oscillator strengths (fosc, arb. u.), and state orbital characters of 2a and 2b using TD-DFT and ADC (2) in implicit water solvation.  (2) emission energy.

Method
Choice of the level of theory. According to the results obtained above, we consider that B3LYP is an adequate level of theory to perform single point TD-DFT calculations on the OG488 fluorophore, also within the GABAAR. We used the default LR-PCM model, as we do not find appreciable changes when using the corrected LR approach. The results of the static calculations on 2a and 2b are presented in Table 1 of the main manuscript.

Section S2 GROUND STATE OF THE FLUORESCENT PROBE 1
The Cartesian coordinates (see repository: https://phaidra.univie.ac.at/o:1383031) of the fluorescent probe gabazine (Gzn) conjugated with Oregon Green (OG488) [14] (1) were generated as discussed in the Introduction of the main manuscript, with two of its more probable protonation states based on a pH value of 7.4 (see Figure S1 and Scheme 1) and taking the pharmacophoric points for the binding of GABA into account. The relevant protonatable moieties are the carboxylate of Gzn, the iminium cation of Gzn, the carboxylate of OG488, and the xanthene of OG488. The latter either has a hydroxy group in the monoanionic form (1b) or is fully deprotonated in the dianionic form (1a) (Scheme 1). The ground state geometry for both forms of 1 was optimized using DFT with the B3LYP [15] functional, the def2-SVP [6] basis set, the Grimme D3 dispersion correction [7] and PCM water solvation. The calculations were carried out in the Gaussian 16, v. A.01 suite. [5]

Section S3 EMISSION OF THE FLUORESCENT PROBE 1
The geometries of the first excited state S1 of 1a and 1b were optimized using with TD-B3LYP-D3/def2-SVP/LR-PCM(water), which is the protocol described in Section S1, apart from the basis set size. The choice of the B3LYP functional is justified from the comparison with other DFTfunctionals against ADC(2) as method of reference on 2a and 2b. We note that TD-DFT has been extensively used in the literature to describe the electronic structure of fluorescein-based [37,38] and organic photoprobes. [39,40] The results for the S1 to S0 emission for the whole fluorescent probes 1a and 1b can be seen in Table 1 in the main text. To further validate the DFT results also for 1, we also calculated the emission of 1b with ADC(2)/def2-SVP using the B3LYP-optimized geometry with implicit water solvation (COSMO), as implemented in TURBOMOLE v7.0. [12] We obtain an emission energy of 1.38 eV with an oscillator strength (fosc ) of 0.010, which is comparable to that obtained with B3LYP-D3 (  Figure S2. Natural transition orbitals (NTOs) involved in the S1 → S0 emission of 1b (the folded conformation found in water). The calculations are done on the optimized S1 geometry at the ADC(2)/def2-SVP/COSMO(water)//TD-B3LYP-D3/def2-SVP(PCM(water) level of theory.

Section S4 SETTING UP OF THE SYSTEMS
The Cartesian coordinates of the receptor were obtained from the cryogenic electron microscopy (cryo-EM) solved structure of the human full-length heteromeric GABAAR(α1/β3/γ2) in complex with the anxiolytic alprazolam (Xanax), the transmitter GABA, and the megabody Mb38 (PDB id. 6HUO). [16] This structure, among other reasons mentioned in the "Results and Discussion" section, was used in this work, because the structure of the same receptor in the presence of another antagonist (bicuculline, see Table S3, PDB id. 6HUK), like Gzn, presented a distorted orientation of α1Arg67 by one of the aromatics in bicuculline. Thus, to preserve the former residue as one of the pharmacophoric points for the binding of Gzn at the beginning of our simulation, we selected the structure in the presence of GABA. Alprazolam, GABA, and Mb38 were removed for our simulations. Since the secondary structure of the intracellular domain (ICD) of all subunits is unknown (α1, amino acids 313-391; β3, amino acids 308-423; γ2, amino acids 323-400) the peptide linker of sequence Ser-Gln-Pro-Ala-Arg-Ala-Ala used experimentally by Miller et al. [17] was introduced as replacement for each of the subunits ICD. [18,19] The protonation states of the titratable residues on the receptor were estimated at pH 7.4 using the H++ server. [20][21][22][23]  in the upper membrane) after reorientation with the OPM server. [24] The system was solvated with TIP3P water [25] molecules with a box of size 130x125x195 Å using the server CHARMM-GUI. [26][27][28][29][30] In addition to 21 neutralizing Clions, a concentration of 0.15 mM of brine (210 Na + and Clions) was added to reproduce the experimental conditions of Sakamoto et al. [14] Before adding the probe, the unbound system was relaxed by means of standard short MD simulations in AMBER18. [31] The placing of the probe inside the GABA binding sites was done with the hypothesis of preserving the two pharmacophoric points of the binding of GABA: the salt-bridges between GABA and residues α1Arg67 and β3Glu155.

Section S5 MD SIMULATIONS
Force field parameters for the fluorescent probe were computed on the DFT(B3LYP-D3/def2-SVP) optimized structure of 1a and 1b using the program Antechamber implemented in AmberTools19. [31] AMBER atom types and RESP charges were used. The receptor was described via the AMBER ff14SB force field [33] and the lipids via the Lipid17 force field. [34] TIP3P water [25] parameters and the ion parameters of Li/Merz were used. [35] Periodic boundary conditions were employed with a cutoff of 10 Å for non-bonded interactions and the Particle Mesh Ewald method for long-range electrostatics. All the MD simulations were run using the AMBER18 software package. [31] The ground state optimized probe forms were neutralized with one or two Na + ions for the 1b and 1a form, respectively, and solvated in a box of TIP3P water [25] with a minimum distance of 40 Å from the solute to the boundaries of the box using tLEaP. [31] Then, the system was minimized in three stages: protons, then solvent molecules and counterions, and finally, the whole system.
Harmonic restraints of 200 kcal mol -1 Å -2 were imposed to the rest of the atoms in each stage.
After that, the system was thermally equilibrated for 20 ps from 100 K to 300 K using the Berendsen thermostat and restraining the solute with 20 kcal mol -1 Å -2 . Finally, and before production, the system was equilibrated in 6 steps from NVT to NPT ensemble for 38 ps. The system was further simulated for 100 ns with the Langevin thermostat with a time step of 0.2 fs.
The simulation of the receptor in the membrane was executed following a similar protocol as for the probe in water. The protocol only deviates in the heating and equilibration process: The system was heated 200 ps from 0 K to 300 K using the Langevin thermostat with harmonic restraint imposed on all the atoms of receptor, probes, and lipids (10 kcal mol -1 Å -2 ); the equilibration was executed in 10 steps for 5 ns in total. Finally, the system was further simulated for 100 ns.
Several parameters of the MD simulation, apart from the ones pictured in the main work, were analyzed. Figure S3 shows the distances d1-d3 along the 100 ns simulation. The choice of the initial receptor structure (PDB id. 6HUK vs. 6HUO) was validated in Figure S4. Figure S5 pictures snapshots of 1a and 1b at the GABA binding site B of the GABAAR(α1/β3/γ2) and Figure S6 displays the root-means square deviation along the 100 ns simulation.  6HUK [16] (magenta), PDB id. 6HUO [16] (purple), and one snapshot from the MD simulation of 1a (green) (B) Plot of the root-means square deviation (RMSD, Å) of the backbone carbon atoms of the receptor along the MD simulation of 1a of 100 ns against PDB id. 6HUK (magenta line) and PDB id. 6HUO (purple line). The RMSD of 6HUK against 6HUO is shown as a red dot.

Section S6 BINDING ENERGY ANALYSIS
The binding free energies of 1a and 1b in binding sites A and B were calculated using the MM/ISMSA program. [36] It estimates the solvent-corrected binding energies (Table S4) as well as their per-residue decomposition (Table S5). The values reported in the results are averages. The last 15 ns of the MD simulations were considered for the calculation of the binding energies in accordance with their root-means square deviation (RMSD) plots ( Figure S6).

Section S7 QM/MM MD SIMULATIONS
As QM region, the OG488 part of 1a (41 atoms) was treated using the third order semiempirical density-functional tight binding method DFTB3 [41][42][43] implemented in AMBER18. [44] The environment (solvent, protein, and lipid membrane) was treated as point charges in the calculation of the energy of the QM part. First, taking the 100 ns MD simulations of 1a bound to GABAAR and in a box of water molecules, one snapshot per system was randomly selected from the major conformational clusters (58% in site A in the membrane context and 58% in aqueous solution).
The cluster analysis was performed using the k-means algorithm implemented in cpptraj for 3 clusters. Second, QM/MM MD simulations were run for 20 ps in the ground state (electronic state S0) for both systems (see Scheme S1). SHAKE and Ewald were deactivated in the QM region and a cut-off of 15 Å (water)/30 Å (GABAAR) was applied from the center of the QM region for the electrostatic term. Third, from the former trajectories on the ground state, 11 snapshots from the last 10 ps were taken and used as input geometries for a 500 fs Born-Oppenheimer MD simulation on the S1 excited state (11 simulations for each system). The QM region, the OG488 part of 1a, was treated with TD-DFT (TD-BP86-D3/def2-SVP) [45] using the interface AMBER18:ORCA. [46,47] The SlowConv algorithm and the deactivation of the SOSCF algorithm was required to get SCF convergences. 8 roots were included in the calculation and the energy of the potential energy surface (PES) of the S1 state was used for the propagation. We also tested the use of B3LYP as DFT-functional for the QM/MM MD propagation on the S1 PES. See Table S6 for the comparison. Although we obtain higher energies with B3LYP compared to BP86, we obtained similar results for oscillator strength (fosc), charge transfer numbers (CT), participation ratio (PR) numbers for the initial/final electron delocalization, and transition character (not shown).
We conclude that the geometry of the molecule is well enough described by BP86. The emission spectrum for 1a was calculated using 5 snapshots from the last 250 fs of each of the 11 S1 excited state Born-Oppenheimer MD simulations in aqueous solution and the protein context (55 snapshots each, see Scheme S1, see repository for xyz-structures: https://phaidra.univie.ac.at/o:1383031). Each snapshot was then used for a TD-DFT single point calculation using Gaussian 16, v. A.01 [5] (B3LYP-D3/def2-SVP) and an electrostatic cut-off of 55 Å (see Section S1 for the discussion about the DFT-functional selection). The QM region is composed of the whole probe in aqueous solution and of only the fluorophore in the protein context. The BP86 functional is sufficient for an efficient S1 propagation, but not for the calculation of the emission spectrum, as it delivers drastic red-shifts. The emission energy is calculated as the vertical excitation energy from ground to first excited state in the optimized S1 geometry. The oscillator strength is used as an approximation for the emission intensity. [48][49][50] The spectra were convoluted by Lorentzian functions weighted by the oscillator strength with a full width at half maximum of 0.10 eV using TheoDORE [51] (see Figure 5A). To one more time validate the use of the B3LYP functional, a second spectrum of 1a was convoluted based on M06-D3/def2-SVP calculations (30 snapshots). The comparison between B3LYP and M06 is shown in Figure S7.
The spectrum with the M06 functional is shifted to higher energies compared to B3LYP, but the characteristics of the spectrum stay consistent. Scheme S1. Schematic workflow of how the 55 snapshots for the emission spectrum of the probe at the receptor were obtained using classical MD, QM/MM MD and QM@MM methods.