On the other hand, HDX-MS experiments produce richer information. A typical experiment starts by equilibrating a protein in H2O at room temperature under physiological conditions (pH 7–8). Then, the protein is diluted with excess D2O for the HDX to occur. At various time points, a small quantity of solution is sampled. The HDX reaction is quenched in the sample by adding acid to lower pH to 2.5, and by cooling it to 0◦C. Proteins in the sample are then digested using acidic proteases (such as pepsin) that are active under quenching conditions. This proteolytic digestion generates numerous peptides, which are portions of the protein typically 6–20 amino acids in length. The sample is then introduced into a chromatography system, to separate the peptides and automatically send them for MS analysis. This analysis allows identifying the peptides generated by the proteolytic digestion and quantifying their deuterium uptake. As the digestion and MS analysis are repeated at various time points, HDX-MS experimental results are usually reported as a set of deuterium-uptake kinetic curves for various peptides (Huang and Chen, 2014).

A crucial technical aspect of HDX-MS experiments is known as back-exchange. This is the process by which the deuterium atoms incorporated by the peptides exchange back to hydrogens. This happens when the sample is prepared for MS analysis because all the required steps (quenching, enzymatic digestion, desalting, chromatographic separation) are performed in H2O solution. On the one hand, back-exchange is beneficial because it enables fast-exchanging side-chain positions to revert to hydrogens, which greatly facilitates the MS identification of peptides by limiting mass changes to amide groups (Wei et al., 2013). On the other hand, back-exchange can become detrimental if slower-exchanging amide groups start reverting to hydrogens, which means losing the information generated by the experiment (Mayne, 2016). To mitigate this problem, all experimental steps have to be performed rapidly, at low temperature.

Unfortunately, back-exchange of amide groups cannot be totally avoided, which affects several aspects of HDX-MS experiments. First, depending on the kind of performed analysis, the measurements produced by the mass spectrometer might have to be corrected for back-exchange (Engen et al., 2011). Second, because terminal positions of a polypeptide chain are more susceptible to back-exchange than other positions, the analysis of peptide-level deuterium-uptake curves has to account for it. More specifically, if the HDX experienced by a given peptide is considered as the average HDX undergone by its amino acids (as done in Section 3.1), the first two amino acids in the chain have to be ignored (Konermann et al., 2011; Huang and Chen, 2014). Indeed, after digestion, the first amino acid of the peptide becomes an amine-terminus, therefore losing its deuterium; as a result, the second amino acid usually undergoes back-exchange as well (Mayne, 2016).

### 2.2. Hydrogen Exchange Estimated from Protein Structure

Numerous theoretical models have been suggested to formalize a relationship between local and/or global structural properties of a protein and the level of hydrogen exchange it undergoes locally. However, none of these models has yet been largely accepted by the scientific community. Several of them have also shown limitations in a recent evaluation study (Skinner et al., 2012b). In this section, we mention the ideas that prevailed in the early days of the research on hydrogen exchange mechanisms, and introduce various models proposed during the past 10 years.

Early attempts to connect hydrogen-exchange mechanisms with protein structure, in the 1970s, were based on accessibility or penetration models. A common view was that solvent-accessible hydrogens located at the protein's surface would exchange rapidly, and that buried hydrogens would exchange more slowly. In other words, protection from exchange was thought to be positively correlated with atom burial or, equivalently, negatively correlated with solvent penetration in the protein matrix. However, it is now well recognized that atom burial is not the primary factor in characterizing hydrogen exchange (Konermann et al., 2011; Skinner et al., 2012b). Indeed, hydrogen-bonded amide groups at the surface can exchange as slowly as deeplyburied amide groups. A variant of this early model of hydrogen exchange based on solvent penetration became popular in the 1980s: hydrogen exchange was thought to be positively correlated with solvent accessibility surface area (SASA). Although this correlation is in general relatively weak (Skinner et al., 2012b; Radou et al., 2014), it has been observed for surface loops of non-globular proteins (Truhlar et al., 2006). This model has been used in qualitative studies of hydrogen exchange (Petruk et al., 2013), sometimes including rigidity properties for increased accuracy (Sljoka and Wilson, 2013).

To explain the fact that even solvent-exposed hydrogens can exchange very slowly, several protein properties have been investigated. For example, there have been some attempts to show that hydrogen exchange is modulated by electrostatic effects on the relative acidity of amides (Anderson et al., 2008; Avbelj and Baldwin, 2009; Hernández et al., 2009; LeMaster et al., 2009). Although this appears to be true in specific cases, in general, no correlation can be expected between protection from hydrogen exchange and changes in relative acidity of amides evaluated via electrostatic calculations (Skinner et al., 2012b). On the other hand, participation in hydrogen bonds is usually recognized as a strong determinant of protection from hydrogen exchange (Skinner et al., 2012b). However, approaches that consider only hydrogen bonding to explain protection from exchange, such as those described in Ma and Nussinov (2011) and Park et al. (2015), are not expected to generalize well. Therefore, some attempts have been made to combine several factors, such as N–H coupling constants and residue fluctuation (Brand et al., 2007).

The most successful approaches to date have been those that combine packing density with various properties related to protein dynamics. On the one hand, some approaches, such as the COREX family of tools, have attempted to link hydrogen exchange to large segmental unfolding reactions (Hilser et al., 2006; Wrabl et al., 2011; Liu et al., 2012). However, a drawback of COREX is that it heavily relies on SASA for doing so. On the other hand, other approaches have attempted to link hydrogen exchange to local interactions (Wu et al., 2009; Gogonea et al., 2010; Craig et al., 2011). Among them, the approach we have adopted in our work relies on the combined evaluation of hydrogen bonding and packing density (Vendruscolo et al., 2003; Best and Vendruscolo, 2006; Gsponer et al., 2006; Kieseritzky et al., 2006; Radou et al., 2014). It is based on a phenomenological equation approximating hydrogen-exchange protection, which is detailed in Section 3.1. Of note, there has been an attempt to predict the coefficients of this phenomenological equation from a protein's amino acid sequence (Tartaglia et al., 2007). Other methods have similarly focused on estimating structural parameters related to hydrogen exchange, directly from protein sequence (Dovidchenko et al., 2009; Lobanov et al., 2013).

### 3. MATERIALS AND METHODS

### 3.1. Phenomenological Approximation of Hydrogen Exchange

As mentioned in Section 2.1, the levels of hydrogen exchange observed in different parts of a protein are known to be partly influenced by its local structure. Several theoretical models have been proposed to formalize a relationship between a protein's conformation and the corresponding hydrogen exchange (see Section 2.2). However, none of them benefits from a consensus of the scientific community, and several of them have shown limitations (Skinner et al., 2012b). Among these models, we chose the one that seemed the most promising, based on its performance in a recent comparative study (Skinner et al., 2012b) and on the number of publications in which it features (Vendruscolo et al., 2003; Best and Vendruscolo, 2006; Gsponer et al., 2006; Kieseritzky et al., 2006; Tartaglia et al., 2007; Radou et al., 2014).

The model we use to estimate hydrogen exchange from a protein's conformation relies on the definition of a phenomenological expression to approximate the protection factors (cf. Section 2.1) of the protein's residues (Vendruscolo et al., 2003). In this theoretical model, it is assumed that protection from hydrogen exchange results from the presence of hydrogen bonds involving amide groups and from the packing density of atoms around these amide groups. More precisely, the protection factor of residue i in conformation C, Pi(C), is derived from the phenomenological expression

$$
\ln P\_i(\mathbf{C}) = \beta^{\text{h}} N\_i^{\text{h}}(\mathbf{C}) + \beta^{\text{c}} N\_i^{\text{c}}(\mathbf{C}) \ , \tag{1}
$$

where N h i (C) is the number of hydrogen bonds formed by the amide hydrogen of residue i, and N c i (C) is the number of so-called "atom contacts" (which is used to quantify packing density) involving residue i. Parameters β h and β c were estimated by fitting experimental hydrogen-exchange data from seven proteins, which lead to: β h = 2 and β c = 0.35 (Best and Vendruscolo, 2006).

Instead of being estimated from a single conformation, hydrogen exchange can also be estimated from a conformational ensemble. In that case, protection factors are computed as ensemble averages. Given a set of conformations, S, the protection factor of residue i with respect to S is derived from

$$\ln P\_i(\mathcal{S}) = \frac{1}{|\mathcal{S}|} \sum\_{\mathcal{C} \in \mathcal{S}} \ln P\_i(\mathcal{C}) \ . \tag{2}$$

The way hydrogen bonds and atom contacts are accounted for has changed over the years, following the evolution of the theoretical model (Vendruscolo et al., 2003; Best and Vendruscolo, 2006). Additionally, not all the details of the methodology have been published. Building on this model, we define hydrogen bonds and atom contacts in the following way:

The residues' protection factors derived from Equation (1) can be directly compared to protection factors obtained from an HDX-NMR experiment. On the other hand, HDX-MS experiments produce deuterium-uptake curves of peptides extracted from a protein. Therefore, a similar kind of data has to be derived from the protein's structure to allow for a comparison with HDX-MS data. For that, we consider that the deuterium uptake of a residue follows pseudo-first-order kinetics (Brier and Engen, 2008; Konermann et al., 2011; Huang and Chen, 2014). Knowing that Pi = k int i / k obs i , the fraction of deuterium incorporated by residue i at time t can be expressed as

$$d\_i(t) = 1 - \exp(-k\_i^{\text{obs}} \ t) = 1 - \exp(-(k\_i^{\text{int}}/P\_i) \ t) \ . \tag{3}$$

As k int i is known (Bai et al., 1993; Connelly et al., 1993), di(t) can be derived from the protein's conformation by calculating Pi . The deuterium uptake of a peptide can be considered as an average over the residues it contains. Therefore, the fraction of deuterium incorporated by peptide j at time t is

$$D\_{\dot{\jmath}}(t) = \frac{1}{n\_{\dot{\jmath}}} \sum\_{i=1}^{n\_{\dot{\jmath}}} d\_{\dot{\imath}}(t) \ , \tag{4}$$

where nj is the number of residues containing an exchangeable amide hydrogen in peptide j (Radou et al., 2014). Note that, in addition to the N-terminal amino acid and to prolines, we systematically exclude from the average the second amino acid (even if it contains an amide group) because of back-exchange (see Section 2.1) (Konermann et al., 2011; Huang and Chen, 2014). Using Equation (4), one can obtain deuterium-uptake curves for various peptides, from any protein conformation.

### 3.2. Goodness-of-Fit between Structurally-Derived and Experimental HDX Data

Using the HDX prediction model presented in Section 3.1, one can derive HDX data from a protein's conformation and compare it to the experimental HDX data. Then, assessing the goodnessof-fit between structurally-derived and experimentally-observed HDX data can be done as follows:

uptake at time t, and D obs j (t) is the experimentally-observed deuterium uptake at time t. This histogram can also be aggregated into an average difference over all peptides (as done in Section 4.3).

### 3.3. Conformation Providing the Best Fit to Experimental HDX Data

The question that remains is: which conformation should the HDX data be derived from to obtain a good fit to the experimentally-observed HDX data? Several studies have shown that conformations reported in the PDB (and more specifically crystal structures) do not provide good estimates for experimental HDX data (Radou et al., 2014; Devaurs et al., 2016). This can be explained by the very nature of HDX data: as it reflects the inherent flexibility of a molecule, in theory, it cannot be accurately predicted from a single conformation. Therefore, it was suggested that hydrogen exchange should be estimated from an ensemble of conformations extracted from an MD simulation, to account for the variability of a protein's structure (Best and Vendruscolo, 2006). Our previous study shows that this methodology also has limitations: better estimates of the experimental HDX data can sometimes be obtained from a single conformation extracted from a conformational ensemble produced by an MD simulation than from the whole ensemble (Devaurs et al., 2016). This shows that, in the context of the structural analysis of experimental HDX data, it is relevant to try and fit this data using a single conformation.

In this work, using computational methods that can sample protein conformations, we aim to obtain a single conformation that can help analyze the experimental HDX data collected for a protein. As PDB conformations produced by X-ray crystallography do not generally provide good estimates for experimental HDX data, they are usually not the best choice for a structural analysis of this HDX data. In spite of this, in our experiments, we systematically evaluate the goodnessof-fit achieved when comparing experimentally-observed HDX data against HDX data derived from a PDB conformation. This provides a baseline against which other methods can be compared. The two methods we evaluate in this study are MD simulations and coarse-grained conformational sampling.

#### 3.3.1. MD Simulations

In this study, all MD simulations were performed with the GROMACS v4.6.5 package (Pronk et al., 2013) using the GROMOS96 (53a6) force field. A cubic box was defined with at least 9 Å of liquid layer around the protein (the exact dimensions were different for each protein), using SPC water model and periodic boundary conditions. An appropriate number of sodium (Na+) and chloride (Cl−) counter-ions were added to neutralize the system, with final concentration of 0.15 mol/L. The algorithms v-rescale (τt = 0.1 ps) and parrinello-rhaman (τp = 2 ps) were used for temperature and pressure coupling, respectively. Cutoff values of 1.2 nm were used both for van der Waals and Coulomb interactions, with Fast Particle-Mesh Ewald (PME) electrostatics. For all MD simulations, the production stage was preceded by (i) three steps of Energy Minimization (alternating steepest-descent and conjugate gradient) and (ii) eight steps of Equilibration. The Equilibration stage started with position restraints for all heavy atoms (5,000 kJ−1mol−1nm−1 ) and a temperature of 310 K, for a period of 300 ps, to allow for the formation of solvation layers. The temperature was then reduced to 280 K and the position restraints were gradually reduced. This process was followed by a gradual increase in temperature (up to 300 K). Together, these Equilibration steps represent the first 500 ps of each simulation. During the production stage, the system was held at constant temperature (300 K) without restraint. The MD simulations were run on various high-performance computers, using between 32 and 144 threads, depending on the size of the protein; the production stage lasted between 150 and 300 ns (additional protein-specific information is provided in Section 3.4). Then, we estimated HDX data as an average over the ensemble of conformations produced by a simulation. We also derived HDX data from every single conformation extracted from such a conformational ensemble.

#### 3.3.2. Structured Intuitive Move Selector (SIMS)

In this paper, we propose a new methodology to obtain a better fit to experimental HDX data, using conformations produced by a coarse-grained conformational sampling approach. For that, we use a computational framework, called Structured Intuitive Move Selector (SIMS), that was developed to explore a protein's conformational space (Gipson et al., 2013). This framework integrates methods known as sampling-based motion-planning algorithms, initially proposed in the field of robotics to randomly explore high-dimensional spaces (Hsu et al., 1999; ¸Sucan and Kavraki, 2010). Using these methods, exploring a protein's conformational space consists of incrementally building a graph whose nodes are conformations and whose edges represent potential transitions between them (Al-Bluwi et al., 2012; Gipson et al., 2012). SIMS follows a "coarse-grained" approach, similarly to MD-like methods using coarse-grained force fields (Davtyan et al., 2012), Monte-Carlo-based simulations (Sim et al., 2012; Boomsma et al., 2013), methods using elastic network models (López-Blanco and Chacón, 2016), or other roboticsinspired conformational sampling methods (Devaurs et al., 2013, 2015).

In SIMS, the exploration starts from a known conformation of the protein (usually, a crystal structure available in the PDB) and aims at producing new conformations by perturbing existing ones. Conformational sampling involves perturbations of the protein's structure, referred to as protein moves. These moves are common perturbation strategies, such as loop sampling, rigid-body motion (i.e., fix one loop's end and move the other end), random perturbation of backbone dihedral angles, and overall energy minimization. To implement these moves and calculate molecular energy, SIMS relies on the Rosetta modeling software (Das and Baker, 2008; Kaufmann et al., 2010). Additionally, SIMS involves an energy threshold to filter the conformations it generates. Note that, by varying this threshold, SIMS can be made more permissive than a typical MD simulation, with respect to the energy of the conformations it generates.

SIMS involves an internal-coordinate representation of proteins in which bond lengths and bond angles are assumed to be constant. Additionally, taking into account the planarity of peptide bonds, the associated torsion angles are restricted to their trans conformation (i.e., ω = 180◦ ). In SIMS, a protein's conformation is represented by a vector of backbone (ϕ, ψ) dihedral angles. Side chains are not explicitly modeled in a conformation, but they are automatically optimized by Rosetta when a move is performed. As a result, a protein composed of N + 1 residues is modeled with 2N degrees of freedom. Such a coarse-grained model has long been shown to provide a good approximation of a protein's behavior (Levitt, 1976).

In SIMS, proteins are decomposed into fragments on which moves are applied. Fragments are specific sets of residues that can be defined automatically, based on secondary structure, or that can be chosen by a domain expert. A fragment can be a protein domain, a single secondary structure element (or several of them), a single residue, or several (non-necessarily contiguous) residues. Using these fragments, one can favor the sampling of specific regions of the protein during conformational exploration. Indeed, based on how flexible some regions are expected to be, fragments are assigned probabilities to be perturbed during conformational sampling (Gipson et al., 2013). These probabilities can reflect available expert knowledge of the protein, reported experimental data (such as B factors) or predicted data resulting from a computational analysis (Fox et al., 2011). In some of the experiments presented here, we use discrepancies between experimentally-observed and structurallyderived HDX data to define these probabilities and therefore guide conformational exploration.

Our experimental methodology can be summarized as follows: First, we use SIMS to perform a conformational exploration starting from the crystal structure of a protein, without using any sampling bias. From the ensemble of conformations generated by SIMS, we determine which conformation provides the best estimates of the HDX data, using Equations (1)–(4). If a good fit is obtained, no additional run of SIMS is performed. If the goodness-of-fit is too low, we run SIMS again, using the largest discrepancies between experimentally-observed and structurallyderived HDX data as a sampling bias: the protein regions where these discrepancies are the largest are assigned higher probabilities to be sampled. We repeat this process a given number of times, or stop when a conformation providing a good fit to the HDX data is obtained.

A single run of SIMS lasted 24 h and was performed on four threads of a 3.6 GHz Intel i7-4790 quad-core CPU. For small proteins, we ran SIMS only once, but for the largest one, we ran SIMS five consecutive times (see Section 4 for more details). For comparison, if the aforementioned MD simulations were run on the same computer, 24 h of computation would yield only 5– 15 ns of simulation, therefore requiring days to weeks for a whole simulation, depending on protein size.

### 3.4. Studied Proteins and Experimental HDX Data

First, we use two small proteins (CI2 and Im7) to illustrate the concepts involved in our methodology. As they have been extensively studied, they represent useful benchmarks. Then, we analyze two medium-size proteins (SN and C3d) that represent more challenging targets for our methodology. Note that we consider two kinds of HDX data: HDX-NMR for CI2, Im7 and SN; and HDX-MS for C3d.

#### 3.4.1. Chymotrypsin Inhibitor 2 (CI2)

We consider a truncated form of chymotrypsin inhibitor 2 (CI2) composed of 64 amino acids (PDB 1TM1), where residue 1 corresponds to residue 20 of the full protein. The main secondary structure elements of CI2 are the following: residues Val13 to Asp23 form an α-helix; residues Gln28 to Pro33 and residues Arg46 to Val51 form two β-strands. As a simple system for folding studies, CI2 has been the subject of several HDX-NMR experiments (Itzhaki et al., 1997; Neira et al., 1997). Protection factors for more than half of CI2's residues have been reported. However, as done in other studies (Best and Vendruscolo, 2006), we only use the protection factors associated with local hydrogen-exchange mechanisms characteristic of CI2's native state. Therefore, we only consider the following 14 residues (whose protection factors are given in parentheses): Leu8 (8.1), Val9 (9.9), Val13 (7.2), Ala16 (7.1), Lys17 (6.6), Lys18 (8.2), Gln22 (9.5), Ala27 (6.7), Gln28 (8.2), Asp52 (8.5), Asn56 (8.4), Ala58 (9), Gln59 (10.5), Val63 (7.4). Note that these protection factors are given as ln P, based on published exchange rates (Itzhaki et al., 1997). Three trajectories of CI2 were obtained by running MD simulations with a 150 ns production stage.

#### 3.4.2. Bacterial Immunity Protein Im7

The bacterial immunity protein Im7 is a single-domain α protein composed of 86 residues (PDB 1AYI). Im7's native state comprises four α-helices: residues Glu12 to Lys24 (I), residues Asp32 to Thr45 (II), residues Thr51 to Tyr56 (III), and residues Glu66 to Asn79 (IV). Helices I and II form an N-terminal helical hairpin, and helix IV is located along the open end of this hairpin. Im7 has been shown to fold through an on-pathway intermediate whose structure is significantly different from that of its native state (Gorski et al., 2004). In this non-native state, helices I, II and IV are conserved, but helix III is not formed. A computational analysis of this intermediate state has shown that helices I, II, and IV are not organized as they are in the native state (Gsponer et al., 2006). This analysis was based on protection factors of 26 residues, obtained via an HDX-NMR experiment aimed at characterizing Im7's folding intermediate (Gorski et al., 2004). Here, we consider the same residues (whose protection factors are given in parentheses): Asp9 (8), Tyr10 (9.2), Thr11 (10.2), Val16 (11.5), Gln17 (11.6), Leu18 (11.2), Glu21 (8.5), Glu23 (6), Leu37 (6.1), Leu38 (7.5), Phe41 (6), Val42 (10.3), Leu53 (3.5), Ile54 (3.6), Tyr55 (4.4), Tyr56 (5.6), Gly67 (8.1), Val69 (8.8), Ile72 (9), Lys73 (9.5), Glu74 (9), Trp75 (8.8), Arg76 (9.9), Ala77 (9.8), Ala78 (8), Lys85 (6.8). These protection factors are given as ln P, based on published exchange rates (Gorski et al., 2004). Three trajectories were obtained by running MD simulations with a 200 ns production stage.

#### 3.4.3. Staphylococcal Nuclease (SN)

Micrococcal nuclease, or Staphylococcal nuclease (SN), is a mixed α/β protein composed of 149 amino acids organized in two domains (PDB 1SNO). The first domain (residues 1– 98) belongs to the oligonucleotide/oligosaccharide-binding-fold (or OB-fold) superfamily. It consists of a five-stranded β-barrel with Greek key topology, capped by an α-helix (residues Gly55 to Glu67) located between the third and fourth strands. The five β-strands are: residues Lys9 to Ala17, residues Thr22 to Tyr27, residues Gln30 to Leu36, residues Ile72 to Phe76, residues Gly88 to Ala94. The second domain (residues 99–149) contains two α-helices: residues Val99 to Arg105, and residues Glu122 to Lys134. SN also contains two minor β-strands. HDX-NMR experiments have been performed on a double mutant of SN with similar structure but increased stability, to characterize its native state (Skinner et al., 2012b). This allowed measuring hydrogenexchange rates for most residues and deriving corresponding protection factors. Here, we use 100 of these protection factors: residues of the N and C terminals that are missing from the crystal structure (PDB 1SNO) are not considered. Note that protection factors were reported as log10 P, instead of ln P (Skinner et al., 2012b).

#### 3.4.4. Complement Protein C3d

C3d is a fragment of the complement component C3 (Nagar et al., 1998; Hammel et al., 2007). It is a single-domain α protein composed of 297 residues (PDB 2GOX), where residue 1 corresponds to residue 991 of the full C3 molecule. C3d contains twelve α-helices and five 310-helices that are organized into an α-α barrel where most consecutive helices alternate between the inside and the outside. Based on previous notations (Nagar et al., 1998), the core of the barrel consists of the following six parallel α-helices: α1 (residues Glu22 to Thr41), α3 (residues Thr86 to Leu102), α5 (residues Lys149 to Ala164), α8 (residues Ser196 to Met209), α10 (residues Gln236 to Leu253), and α12 (residues Ser278 to Asp295). It is surrounded by another set of six parallel helices (running anti-parallel to those of the core) comprising one 310-helix, T1 (residues Ala7 to Leu13), and five α-helices: α2 (residues Leu49 to Arg70), α4 (residues Ser107 to Lys121), α7 (residues Ser174 to Asn189), α9 (residues Pro215 to Thr223), and α11 (residues Phe256 to Gln269). In previous work, we performed an HDX-MS experiment and several MD simulations on C3d, to characterize its native state (Devaurs et al., 2016). In this paper, as in our previous study, we use the deuterium-uptake data obtained for 81 peptides extracted from C3d.

### 4. RESULTS

We now report the results we have obtained for the four proteins introduced in Section 3.4. First, a comparative analysis of CI2 and Im7 sheds light on the issues encountered when trying to fit experimental HDX data with the different methods presented in Section 3.3. Then, we examine two medium-size proteins: SN and C3d.

### 4.1. Chymotrypsin Inhibitor 2 vs. Bacterial Immunity Protein Im7

Our first set of results aims at highlighting differences and similarities that exist between two possible scenarios: (i) the case where the HDX data and the crystal structure describe the same state of the protein, and (ii) the case where they describe two different states. As mentioned in Section 3.4.1, the HDX-NMR data (i.e., protection factors of residues) obtained for CI2 is characteristic of its native state, whose structure has been described (PDB 1TM1). On the contrary, the HDX-NMR data gathered for Im7 is known to characterize a non-native folding intermediate (see Section 3.4.2) that is structurally different from Im7's native state (PDB 1AYI).

The comparison between these two scenarios is illustrated in **Figure 1**. The native conformations of CI2 and Im7, as reported in the PDB, are depicted in green using the ribbon model. The bar charts show that the HDX data derived from the PDB conformations (i.e., the crystal structures reported in the PDB) using Equation (1) does not match well the experimentallyobserved HDX data: the average difference between structurallyderived and experimentally-observed protection factors (see Section 3.2) is close to 3. Although this is not surprising in the case of Im7 (because the HDX data and the crystal structure describe different states), it is important to note that the HDX estimates are equally bad in the case of CI2 (even though the HDX data and the crystal structure describe the same state).

For both CI2 and Im7, we performed three MD simulations. We observe that, as suggested in Radou et al. (2014), deriving HDX data from the ensemble of conformations extracted from each MD simulation leads to a better fit to experimentallyobserved HDX data than if the PDB conformation is used. However, it also appears that the best fit is usually obtained with a single conformation selected within these MD ensembles, which confirms our previous results on C3d (Devaurs et al., 2016). The bar charts in **Figure 1** show the differences between the experimentally-observed and structurally-derived protection factors, when deriving these protection factors from the MD conformation providing the best fit. It is clear that using this MD conformation yields a better fit to the experimental data than using the PDB conformation, but not drastically. In the case of Im7, the limited improvement was expected: our MD simulations were meant to sample the native state; they were not long enough to observe a transition to the folding intermediate. Even in the case of CI2, we will show that better results can be obtained.

We used SIMS to sample the conformational space of CI2 and Im7, starting the exploration from their PDB conformations, without any bias. From the sets of conformations generated by SIMS, we extracted the conformation yielding the best fit between structurally-derived and experimentally-observed HDX data. The bar charts in **Figure 1** show differences between experimentally-observed and structurally-derived protection factors, when deriving them from the SIMS conformation with the best fit. This conformation yields a significantly better fit to experimental HDX data than the PDB and MD conformations. The SIMS conformations for CI2 and Im7 are depicted in red using the ribbon model in **Figure 1**. In the case of CI2, the SIMS conformation is very similar to the PDB conformation: differences occur mostly in side-chain positions and not in backbone structure. This was expected because the HDX data and the crystal structure describe the same state. This also highlights the strong impact that even small structural variations can have when estimating protection factors with Equation (1). In the case of Im7, the SIMS conformation providing the best fit to the experimental HDX data is significantly different from the PDB conformation, which confirms that the HDX data and the crystal structure describe different states.

### 4.2. Staphylococcal Nuclease

A recent evaluation study of various models for deriving hydrogen exchange from a protein's structure involved HDX-NMR data gathered for SN (see Section 3.4.3) (Skinner et al., 2012b). The study concludes that, at least for SN, none of the evaluated models can produce HDX data that fits well the experimentally-observed HDX data. The best results are achieved by the model based on Equation (1), with a correlation coefficient R 2 = 0.51 between the structurally-derived and experimentally-observed protection factors. That study follows the methodology in Best and Vendruscolo (2006), estimating the protection factors of SN's residues using an ensemble of conformations extracted from an MD simulation. However, that study does not consider estimating HDX data from the PDB conformation alone, using Equation (1). Interestingly, using this PDB conformation, we obtained a correlation coefficient R 2 = 0.69 between the structurally-derived and experimentallyobserved protection factors (see Section 3.2). Note that better results can be achieved with our novel methodology.

We used SIMS to explore the conformational space of SN, starting from its PDB conformation, without introducing any bias. From the ensemble of conformations generated by SIMS, we extracted the conformation providing estimates of protection factors that best fit the experimental HDX data. This yields a correlation coefficient R 2 = 0.78 between the structurallyderived and experimentally-observed protection factors, as shown in **Figure 2**. Importantly, the SIMS conformation is very similar to the PDB conformation: only small structural differences are observed at the backbone level (see **Figure 2**). This confirms that the HDX data and the crystal structure both describe SN's native state.

### 4.3. Complement Protein C3d

In previous work (Devaurs et al., 2016), we performed an HDX-MS experiment on C3d and obtained deuterium-uptake curves for 86 peptides. As in that previous study, we restrict the current analysis to the 81 peptides whose data is the most reliable. This HDX data is expected to describe the native state of C3d when present alone in solution. However, once again, deuterium-uptake curves derived from the PDB conformation of C3d (PDB 2GOX) using Equations (1)–(4) do not fit well the experimental data (see **Figure 3**). The average difference between the experimentally-observed and structurallyderived deuterium-uptake curves across all peptides is 1.23 (see Section 3.2). Discrepancies are especially significant in the region of C3d comprising residues Met191 to Ala242. As shown in Devaurs et al. (2016), this is not due to structural differences between the native state of C3d and the conformation observed during the HDX-MS experiment, but rather to the limitations of predicting HDX data using crystal structures.

We carried out three MD simulations to sample the variability of C3d's native state (Devaurs et al., 2016). Using the ensemble

conformations produced by SIMS that provide the best estimates of the experimental HDX data are colored in red. The bar chart at the bottom shows the average difference (across residues) between experimentally-observed and structurally-derived HDX data (i.e., protection factors), when deriving this data using conformations in the PDB (green), conformations produced by MD which best fit the HDX data (blue), or conformations generated by SIMS which best fit the HDX data (red). The bar charts in the middle show these differences for all residues separately.

of conformations extracted from each simulation allows deriving deuterium-uptake curves of peptides that fit the experimental data better than when using the PDB conformation. However, an important conclusion of our previous study is that: using a single conformation extracted from these MD ensembles produces even better results (Devaurs et al., 2016). The conformation providing the estimates of deuterium-uptake curves that best fit the experimental HDX data is referred to as the MD conformation. It yields a decrease in the average difference (0.89) between structurally-derived and experimentally-observed HDX data. Despite the improvement in goodness-of-fit, large discrepancies remain (see **Figure 3**), especially in the region [Met191-Ala242] of C3d (Devaurs et al., 2016).

To sample C3d's conformational space more extensively, we carried out the following iterative process with SIMS: using the PDB conformation as input, we ran SIMS once without

introducing any bias; then, we ran SIMS four times, using the discrepancies between structurally-derived and experimentallyobserved HDX data as a sampling bias. This bias is introduced in the following way: at the end of each run, we select the conformation generated by SIMS providing estimates of deuterium-uptake curves that best fit the experimental HDX data, and we determine the regions of C3d where discrepancies are the largest; then, in the following run, these regions are assigned higher probabilities to be sampled (cf. Section 3.3.2). This SIMS-based iterative process generated a conformation providing estimates of deuterium-uptake curves that fit well the experimental HDX data (see **Figure 3**). Using this SIMS conformation, the average difference between the experimentally-observed and structurally-derived deuteriumuptake curves across all peptides decreases to 0.6. Importantly, the SIMS conformation is very similar to the PDB conformation: all the helices forming the α-α barrel are conserved; only two short helices have unfolded. The α-α barrel of the SIMS conformation (radius of gyration: 19 Å) is only slightly wider than the α-α barrel of the PDB conformation (radius of gyration: 18 Å). This confirms that the HDX data and the crystal structure both describe C3d's native state. It also confirms that the native state of C3d is relatively stable, with little flexibility recorded by the HDX-MS experiment. Finally, note that using the SIMS conformation can provide an improved structural analysis of C3d, by refining the HDX data from the peptide level to the residue level (Devaurs et al., 2016).

### 5. DISCUSSION

We chose to include in our methodology the phenomenological approximation of hydrogen-exchange protection, as expressed by Equation (1), because it seemed to be the most promising HDX prediction model. Indeed, it performed best at predicting experimental HDX data, when compared to several other models (Skinner et al., 2012b). Even though the goodness-offit achieved with this model was not impressive (Skinner et al., 2012b), our work demonstrates how it can be used successfully.

Our results clearly show that using the conformation of a protein as reported in the PDB does not provide good estimates of experimental HDX data. This confirms what was observed in previous similar studies (Radou et al., 2014; Devaurs et al., 2016). This was also indirectly acknowledged when this HDX prediction model was first proposed (Vendruscolo et al., 2003). In an attempt to consider structural dynamics, it was suggested that HDX data should be derived as an average over an ensemble of conformations produced by an MD simulation (Best and Vendruscolo, 2006).

We have indeed observed that computing an average over an MD ensemble provides better estimates of experimental HDX data than using a single PDB conformation. However, as shown in our previous work (Devaurs et al., 2016), this study confirms that using a single conformation carefully extracted from the MD ensemble usually provides even better estimates. In other words, the MD conformation that provides the best estimates within the MD ensemble performs generally better than the whole ensemble. Note that we do not claim that this MD conformation constitutes a better representation of a protein's state than a PDB conformation or an MD ensemble. In theory, the best estimates for experimental HDX data would be obtained by computing an average over an ensemble of conformations best representing a protein's state and its inherent flexibility. However, in the same way as estimates derived from two similar conformations can significantly differ, estimates derived from two similar ensembles

can be very different. In practice, it is thus more convenient to generate many conformations and select the one providing the best estimates than to find the best conformational ensemble.

The fact that numerous conformations have to be generated in order to obtain good estimates of experimental HDX data, and that a PDB conformation is not enough, is also linked to weaknesses of the HDX prediction model based on Equation (1). The first limitation of this model is its lack of robustness: it is very sensitive to small variations in the protein structure. As well illustrated by the case of CI2, two conformations that are very similar at the backbone level and present differences only in their side-chain conformations can produce very different HDX estimates. The second limitation of this model is that it only partially reflects the mechanisms underlying hydrogen exchange. For example, it does not consider any dynamic aspect of proteins. Therefore, it could be interesting to develop a more accurate model by accounting for additional structural and dynamic properties of proteins (Skinner et al., 2012a). Since such a model has not been proposed yet, we believe it is best to compensate for the weaknesses of the current model by performing conformational sampling.

### 6. CONCLUSION

When performing a structural analysis of HDX data collected for a protein, a premise to an accurate analysis is to use a conformation that matches this data. Several studies, including ours, show that crystal structures reported in the PDB are not a good choice because they often provide bad estimates of experimental HDX data. Because HDX data reflects the inherent flexibility of a protein, a conformational ensemble should ideally provide better estimates than a single conformation. However, our work has shown that this is not always the case with a conformational ensemble produced by an MD simulation. Therefore, it is perfectly justified to try and fit experimental HDX data using a single conformation. In this paper, we have shown that this can be done using a coarse-grained conformational sampling tool to explore a protein's conformational space. The specific tool we used, called SIMS, yields a conformational ensemble from which one can extract a conformation providing a good fit to the experimental HDX data. Note that we do not claim that a conformation produced by SIMS is a better representation of a protein's state than its crystal structure. Besides the improved accuracy, another advantage of using SIMS is its efficiency: a conformation providing a good fit to experimental HDX data can be obtained at a fraction of the computational cost of running a traditional MD simulation. Finally, we believe that other conformational sampling methods could produce similar results, in terms of accuracy and efficiency. The achievement of our study mostly consists of revealing the technicalities that must be addressed for such methods to be successful.

Our methodology relies on the use of an HDX prediction model defining how to derive HDX data from a protein's structure. This model is based on a phenomenological approximation of the protection factors of a protein's residues. Despite its limitations, this model enables our methodology to successfully produce a conformation fitting the experimental HDX data. Another interesting benefit of this model is that, besides the validation of experimental HDX data, in the case of HDX-MS experiments, it offers the possibility to refine the HDX data from the peptide level to the residue level (Radou et al., 2014; Devaurs et al., 2016). This has the potential to enhance applications of the HDX-MS technique (Pirrone et al., 2015).

As part of our future work, we intend to apply our methodology to larger proteins, to evaluate its scalability. Since coarse-grained conformational sampling scales better than MD, we expect our methodology to be even more beneficial with large proteins. Additionally, we plan to investigate several useful applications of this work. First, as demonstrated with Im7, our method can be used to obtain a structural model of a non-native state of a protein when only its native state is described in the PDB and only HDX data is available for this non-native state. Second, although we applied our method only to cases where the experimental HDX data was expected to characterize a single protein conformation because a single conformer was assumed to be present in solution, it could be applied to more complex cases, where several conformers are involved. Indeed, if structurallyderived HDX data better fits experimentally-observed HDX data when deriving it from a small set of structurally-different conformations (i.e., two or three, or a handful of conformations) than when deriving it from a single conformation, we can suspect that several protein conformers are present together in solution.

### AUTHOR CONTRIBUTIONS

DR, JL, LK, MM suggested the initial idea behind this work. DD conceived the study, including the methodological choices and the selection of proteins, and implemented the scripts required for the analysis. MP performed the HDX-MS experiment. DA performed the MD simulations. DD carried out the analysis of the results and wrote the paper. DA, DR, JL, LK, MM, MP reviewed and commented on the paper.

### FUNDING

This work was supported in part by the National Science Foundation under Grant CCF 1423304, as well as the National Institutes of Health under Grants R21CA209941, AI068730, and AI030040. Computational simulations were run on equipment supported in part by the Data Analysis and Visualization Cyberinfrastructure funded by NSF under Grant OCI 0959097, on equipment supported by the Cyberinfrastructure for Computational Research funded by NSF under Grant CNS 0821727, as well as on equipment supported in part by the Big-Data Private-Cloud Research Cyberinfrastructure MRI-award funded by NSF under grant CNS 1338099 and by Rice University.

### ACKNOWLEDGMENTS

The authors are thankful to Emanuele Paci and Gaël Radou for providing information on their implementation of the HDX prediction model based on the phenomenological approximation of protection factors. All protein structures are depicted by images produced with the PyMOL Molecular Graphics System, Version 1.8 Schrödinger, LLC.

#### Devaurs et al. Conformational Sampling for HDX Data Fitting

### REFERENCES

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Devaurs, Antunes, Papanastasiou, Moll, Ricklin, Lambris and Kavraki. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Hybrid Methods in Iron-Sulfur Cluster Biogenesis

Filippo Prischi 1 and Annalisa Pastore2, 3 \*

1 School of Biological Sciences, University of Essex, Colchester, UK, 2 Maurice Wohl Institute, King's College London, London, UK, 3 Molecular Medicine Department, University of Pavia, Pavia, Italy

Hybrid methods, which combine and integrate several biochemical and biophysical techniques, have rapidly caught up in the last twenty years to provide a way to obtain a fuller description of proteins and molecular complexes with sizes and complexity otherwise not easily affordable. Here, we review the use of a robust hybrid methodology based on a mixture of NMR, SAXS, site directed mutagenesis and molecular docking which we have developed to determine the structure of weakly interacting molecular complexes. We applied this technique to gain insights into the structure of complexes formed amongst proteins involved in the molecular machine, which produces the essential iron-sulfur cluster prosthetic groups. Our results were validated both by X-ray structures and by other groups who adopted the same approach. We discuss the advantages and the limitations of our methodology and propose new avenues, which could improve it.

#### Edited by:

Slavica Jonic, UPMC Sorbonne Universités, France

#### Reviewed by:

Michael Sattler, Technische Universität München, Germany Ivaylo Ivanov, Georgia State University, USA

> \*Correspondence: Annalisa Pastore<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biophysics, a section of the journal Frontiers in Molecular Biosciences

Received: 14 November 2016 Accepted: 23 February 2017 Published: 13 March 2017

#### Citation:

Prischi F and Pastore A (2017) Hybrid Methods in Iron-Sulfur Cluster Biogenesis. Front. Mol. Biosci. 4:12. doi: 10.3389/fmolb.2017.00012 Keywords: frataxin, NMR, molecular complexes, small angle X-ray scattering, structural biology, iron-sulfur cluster machinery, hybrid methods

### INTRODUCTION

Biophysical approaches that make the combined and integrated use of different methodologies are named "hybrid techniques."  A particularly useful application of hybrid techniques is the use of a combination of high and low resolution techniques which first target the local structure of a molecule (a domain or a complex component) and then reconstruct the full picture of the assembly (the so-called cut-and-paste approach) (Grishaev et al., 2005, 2008; Parsons et al., 2008; Deshmukh et al., 2013). Hybrid methods have, for instance, been successfully introduced to gain structural insights of complexes with different sizes which would be unaffordable if approached by only one technique (Wüthrich, 2001). One of the very first examples of hybrid methods was our study based on a combination of small angle scattering (SAXS) and nuclear magnetic resonance (NMR) to approach the arrangement of the domains of titin, a giant muscle modular protein containing more than 300 copies of two all-β sequence motifs, the fibronectin type 3 and the immunoglobulinlike modules (Improta et al., 1998). More recently, Michael Sattler and co-workers extended the approach to the study of RNA-protein interactions (Gabel et al., 2006; Madl et al., 2011; Huang et al., 2014). In 2010, we implemented a robust methodology which brings together NMR, SAXS, in site directed mutagenesis, ITC and other techniques to study weak complexes. This methodology has proven particularly effective for proteins of the iron-sulfur (FeS) clusters biogenesis machine, a highly conserved and essential metabolic pathway (Zheng et al., 1998). These proteins share important features which make particularly useful the application of hybrid methods to their study:

all the components of this cellular machine (i) tend to form transient interactions, making co-purification of the complexes difficult to impossible; (ii) have different binding affinities from each other; (iii) compete for the same binding sites; (iv) have different likelihood to crystallize, which often results in proteins forming crystals alone and not as part of the complex. In addition, many of the complexes are, although relatively large for NMR studies, too small for cryo-electron microscopy studies (Nogales and Scheres, 2015). As a result, high-resolution structures of most of these protein complexes are still not available. Here, we review our approach, discuss its successes and clarify the limitations. We also suggest ways to circumvent specific problems.

### THE PARADIGMATIC EXAMPLE OF THE IRON-SULFUR CLUSTER BIOGENESIS COMPONENTS

Present ubiquitously in nearly all life forms, FeS clusters are protein inorganic prosthetic groups involved in a multitude of biological functions, such as electron transfer, gene expression regulation, thiolation, photosynthesis, nitrogen fixation, metal trafficking, substrate binding, DNA repair/replication and RNA modification (Johnson et al., 2005; Mettert and Kiley, 2015). FeS clusters are formed from iron ions and inorganic sulfide. Due to the toxic nature of these elements, formation of intracellular FeS clusters does not occur spontaneously, but all organisms have evolved protein machineries for the production of clusters. The FeS cluster assembly (ISC) system is a highly conserved factory found both in prokaryotes and eukaryotes and capable of providing FeS clusters to a wide range of apo-proteins. In particular, the eukaryote ISC machine is found in the matrix space of mitochondria and is distinct from the system that produces the clusters in the cytosol (Lill, 2009; Rouault, 2015). In E. coli, which is most studied as a model system because of its lower complexity, the ISC machine is composed of eight genes clustered in an operon, iscRSUA-hscBA-fdx-iscX (Takahashi and Nakamura, 1999; **Figure 1**). The operon is controlled by IscR, a transcriptional repressor (Schwartz et al., 2001) followed, in the order, by genes coding for a cysteine desulfurase (IscS) (Schwartz et al., 2000), a scaffold protein upon which clusters are built (IscU) (Agar et al., 2000), an A-type carrier with unclear function (IscA) (Krebs et al., 2001; Ollagnier-de-Choudens et al., 2001), a co-chaperone/chaperone system that is thought to facilitate cluster transfer from IscU to the final acceptor (hscA and hscB) (Chandramouli and Johnson, 2006), an electron donor ferredoxin (Fdx) (Yan et al., 2013b) and a protein with unknown function (IscX or YfhJ) (Pastore et al., 2006). This system constitutes the so called core assembly machine. The formation of FeS clusters by the core machine starts with the production of S 0 from L-cysteine by IscS, followed by reduction of S0 to S2− by Fdx (Yan et al., 2013b) and ends with the incorporation of Fe2+ or Fe3+ and formation of a [2Fe-2S] cluster on IscU (Agar et al., 2000). It is still unclear how the iron is delivered to the system.

Among these proteins, the crucial ones are IscS (NFS1 in human) and IscU (ISCU in human). IscS is a pyridoxal 5 ′ -phosphate (PLP)-dependent desulfurase. PLP is not only necessary for the catalysis of L-Cys to L-Ala, but is also important for structural stabilization. CD spectra of recombinant IscS without PLP revealed in fact that the protein is completely unfolded, albeit proteolytically stable and not prone to aggregation (Prischi et al., 2010b). IscU, a 10 kDa protein, is predominantly in a monomeric state in solution and binds to IscS to accept the sulfur, which will form the clusters. The function of IscS and IscU are regulated by the protein frataxin (FNX). This is an essential protein highly conserved both in prokaryotes (where takes the name of CyaY) and eukaryotes where it is present in mitochondria. FXN was first identified for its connection to Friedreich's ataxia (Campuzano et al., 1996), a progressive neurodegenerative disease caused by an expansion of a GAA trinucleotide repeat within the first intron of the FXN gene, which results in reduced levels of FXN (Campuzano et al., 1996, 1997). Studies on the yeast frataxin homolog (YFH1) helped to understand that reduced levels of FXN causes loss of function of FeS cluster containing enzymes, increased amount of free radicals and iron deposits in mitochondria (Babcock et al., 1997; Foury and Cazzalini, 1997; Koutnikova et al., 1997; Rötig et al., 1997). Proteins from the FXN family bind weakly ferrous ions (Kd 4 µM) and ferric ions (Bou-Abdallah et al., 2004). These features are strongly conserved: human FXN is able to bind Fe2+ and Fe3+ in a similar way (Yoon and Cowan, 2003). Ability to weakly bind iron could be in agreement with the hypothesis that the protein functions as an iron chaperone, but the way FXR binds iron is unusual. The FXN fold, which is composed of two α-helices packed against an anti-parallel β-sheet (Cho et al., 2000), does not share any similarity with any other known iron binding proteins, like ferritins, ferredoxins or hemoglobins (Harrison and Arosio, 1996). It is also unusual that iron coordination occurs solely through carboxylate residues and no conserved histidine, cysteine, or tyrosine - residues usually found in iron binding motifs - are present in frataxins (Nair et al., 2004). Finally, cation binding is highly unspecific since, in addition to iron, frataxins bind to diamagnetic Ca2+, Zn2+, Lu3+, and paramagnetic ions Mn2+, Co2+, Gd3+, Eu3+, and Yb3+ (Nair et al., 2004). Twenty years have passed since these initial studies which have made clear that FXN is connected to FeS cluster formation, but the exact function of FXN remains elusive. Different theories have been proposed: (i) FXN is the iron chaperone that delivers Fe2+ or Fe3+ to IscU (Yoon and Cowan, 2003); (ii) FXN acts as a scavenger that is able to sequester mitochondrial iron through formation of high-molecular-weight aggregates and to maintain it in a bioavailable form (Adamec et al., 2000). We have proposed a third hypothesis which is currently the most accredited: (iii) FXN acts as an iron sensor that regulates the amount of FeS cluster formed to match the concentration of the available acceptors (Adinolfi et al., 2009). Our model, which proposes a completely new function of FXN, is based on studies that rely on the demonstration that FXN binds to the IscS/IscU complex in an iron dependent manner (Prischi et al., 2010a). To gain information on this ternary complex, we adopted a hybrid approach, which relied on NMR, SAXS, site directed mutagenesis, molecular docking and molecular dynamics simulations.

### THE DEVELOPMENT OF A HYBRID METHOD

The rationale of our hybrid method develops through the following logic steps (**Figure 2**):

### Step 1: Identifying the ISC Interactome

The network of interactions between IscS, IscU, and CyaY was probed by NMR spectroscopy. We exploited the well-known concept that the spectrum of a molecule is very sensitive to the chemical environment. Protein-protein interaction cause changes in the chemical environment of the reporter nucleus. This means that titration of a protein with another molecule results in shifts in the position of some or all resonances in the spectrum or Chemical Shift Perturbation (CSP) (Roberts, 1993; Zuiderweg, 2002), which can then be used to map the regions involved in the interaction. Typically, we titrated a 15N labeled component of the isc operon with another un-labeled protein. However, the resonance line widths depend inversely on the tumbling time (Bloembergen et al., 1948) and, thus, the larger the complex, the broader are the line widths up to spectral disappearance. Many of the ISC components have sizes well within the limits of NMR observation except for IscS, which is an obligate dimer of 90 kDa, and the chaperone HscA. This meant that we could alternatively titrate the low size proteins (i.e., adding unlabeled protein A into 15N labeled protein B and, viceversa, unlabeled protein B into 15N labeled protein A) and map the interacting site on both proteins. The case was quite different when adding the 90 kDa IscS to a smaller component. In this case the result would strongly depend on the regime of exchange of the complex.

The most common NMR experiment used to measure CSP is the two-dimensional 15N heteronuclear single-quantum coherence NMR ([1H,15N]-HSQC NMR), a method that allows the detection of correlations between 15N nucleus and 1H nucleus which are covalently bound. Titration of IscS into 15N labeled IscU caused complete disappearance of IscU signal from the [1H,15N]-HSQC NMR spectra at a 1:0.7 IscU:IscS molar ratio (Prischi et al., 2010b) without previous CSP. Absence of detectable CSP for these titrations and disappearance of the IscU signal indicates binding but also suggests that the process is under an intermediate-to-slow exchange regime in the NMR time range (**Figures 3A–C**). The exchange regime is the rate kex at which a nucleus switches from one conformation to another (in this case a "free state" to a "bound state"). The NMR linewidths depend on the populations of each state, the relative values of the exchange rate kex and the chemical shift difference ∆ν. In the slow exchange regime (kex << |∆ν|), signals from both states are observed at their distinct chemical shifts, intensities and linewidths; if the regime is fast (kex >> |∆ν|), a single peaks will be observed at the chemical shift between free and bound conformations weighed according to the populations; if it is in an intermediate regime (kex ≈ |∆ν|), a single peak is observed between the two states but due to the presence at the same time of the free state and the bound state, the resulting resonance is broadened (Kleckner and Foster, 2011). In our case IscU alone corresponds to the "free state", while the IscU-IscS complex is the "bound state". Since the spectra are completely unperturbed until we reach a 1:1 ratio IscU:IscS, we can deduce that this is not in a fast exchange regime in the NMR time range and we can deduce that the process is an intermediate-to-slow exchange regime. We would expect then to see peaks for both the free state and the bound one. However, the high molecular weight of the complex causes that the bound state is outside the limit of NMR observation and we do not observe it. This did not allow us to map the interaction surface of IscU on IscS, a problem often observed in the NMR studies of complexes. We also did not observe CSP when titrating directly IscU and CyaY in the absence of IscS but in this case the spectra of the two proteins, individually labeled in turn, where completely unaffected. This meant no

direct interaction in contrast with studies carried out on the human and yeast proteins, where a direct interaction between the scaffold protein and FXN was observed (Yoon and Cowan, 2003; Correia et al., 2009; Leidgens et al., 2010), suggesting a different behavior of the bacterial proteins. The difference could perhaps be ascribed to the lack of the N-terminal extension which, in eukaryotes, is part of the mitochondrial signal and absent in prokaryotes. Finally, when titrated with IscS, the spectrum of CyaY remains visible and shows clear CSP, which allowed us to map the interaction on a specific surface (**Figures 3D–F**). CyaY interacts with IscS using a negatively charged surface area localized on α1, β1 and α1β1 and β1β2 loops (Adinolfi et al., 2009). Interestingly, this negatively charged surface is the same involved in iron binding (Yoon and Cowan, 2003; Bou-Abdallah et al., 2004; Nair et al., 2004). We then tested for competition between CyaY and IscU binding on IscS by [1H,15N]-HSQC NMR spectra titrating 15N labeled IscU with up to an equimolar amount of unlabeled IscS with unlabeled CyaY. Presence of competition should cause dissociation of 15N labeled IscU from IscS, resulting in the reappearance or increase of the NMR signal, in a way proportional to the amount of competitor added. We did not observe competition with IscU (Adinolfi et al., 2009; Prischi et al., 2010a). We could thus conclude that both CyaY and IscU bind to IscS, but not each other and obtain the surface of interaction on CyaY from NMR only.

CSP data did rule out the presence of a direct interaction or competition between CyaY and IscU, but this did not automatically exclude binding between the two when in the presence of IscS. We titrated 2H, 15N double-labeled CyaY with unlabeled IscU and IscS up to a 1:1:1 molar ratio. 2H labeling reduces spin-spin relaxation, a parameter inversely proportional to the linewidths of the resonance in the spectrum. This results in narrower linewidths and thus provide higher resolution (Gardner and Kay, 1998). We observed a new set of spectral perturbations, which we attributed to a direct contact of the residues involved with IscU. Once mapped onto CyaY structure, these residues clustered on the anti-parallel β-sheet surface of CyaY (**Figure 4A**). In particular, the conserved Trp61 in CyaY was found to be involved in the interaction with IscU. These data are in agreement with studies on human FXN, where it was shown that the exposed side chain of Trp155 (equivalent to CyaY Trp61) is indispensable for FXN-ISU binding (Correia et al., 2009; Leidgens et al., 2010).

To map the surface of interaction on IscS we used site directed mutagenesis. We designed mutations of IscS targeting solvent exposed residues. We aimed to abolish interaction with ISC components, while keeping IscS stable and functional. We titrated 15N labeled IscU and CyaY with five different IscS mutants, i.e., IscS\_R220E/R223E/R225E, IscS\_I314E/M315E, IscS\_K101E/K105E, IscS\_E334S/R340S and IscS\_R39E/W45E (Prischi et al., 2010a). IscS\_R220E/R223E/R225E triple-mutant, in which a positively charged patch formed mainly by arginines close to the dimer interface was inverted in charge, does not bind CyaY (**Figure 4B**; Prischi et al., 2010a). This strongly supported the assumption that binding between these proteins is driven by electrostatic interactions and is in agreement with our competition studies. Differently, in IscS\_I314E\_M315E we inserted two charged residues into an uncharged/hydrophobic patch. This mutant has a reduced affinity for IscU, which, due to a change in the exchange regime, caused chemical shift perturbation of the 15N-labeled IscU HSQC spectrum. This not only allowed us to identify IscS interacting surface, but also to identify the IscU residues involved in IscS binding (**Figure 4A**; Prischi et al., 2010a). All other IscS mutants behaved like the wild type in titration experiments and provided us with solid controls (Prischi et al., 2010a).

#### Step 2: Restrained Docking Simulation

To gain a visual impression and understand the relative orientation of proteins in the complex, we generated models of the central ISC machine using NMR restrained molecular docking simulations. We used the docking software HADDOCK (Dominguez et al., 2003). HADDOCK can incorporate NMR or other distance restraints and implement them as "ambiguous interaction restraints" (AIRs). The software forces the protein interfaces to come together without imposing a particular

FIGURE 3 | Comparison among close-ups from HSQC spectra of 15N uniformly labeled IscU. (A) From left to right, IscU alone; (B) IscU in the presence of unlabeled IscS at a molar ratio of 1:0.5; (C) The same as in (B) but at a molar ration of 1:1. No detectable CSP of IscU peaks and concomitant disappearance of the IscU spectra when saturated with IscS suggest that the process is under an intermediate-to-slow exchange regime in the NMR time range. For comparison, The spectra were all recorded at 25◦C and 600 MHz and visualized using the same signal to noise ratio. (D) CyaY alone; (E) CyaY in the presence of unlabeled IscS at a molar ration of 1:0.5; (F) The same as in (E) but at a 1:1 molar ratio. The spectra were all recorded at 25◦C and 600 MHz and visualized using the same signal to noise ratio.

orientation. Using the AIRs that we determined experimentally during Step 1, we obtained different families of complexes, which differed by details but all reported IscU bound on the opposite tips of the IscS dimer, roughly close to the N-terminus. CyaY was instead consistently located near the cavity that contains the active site of IscS, spatially close but not overlapping with IscU. While these results could have been sufficient for having a first rough model of the IscS-IscU-CyaY complex, we felt that further validation was needed to confirm independently the relative positions of the three proteins.

### Step 3: Validation through SAXS Data

The generated models were then experimentally verified and re-scored using SAXS data. SAXS is a solution technique that allows to study the shape, conformation and assembly state of proteins and, more in general, macromolecular complexes (Mertens and Svergun, 2010). Despite being a low-resolution technique, SAXS is well suited for the study of flexible systems and intrinsically disordered proteins (Wang et al., 2011), which are major limitation in X-ray crystallography and cryo-EM. It also allows the study of proteins in solution in nearly-native conditions. SAXS experiments are not time consuming. Recent hardware improvements allow high-throughput studies (Round et al., 2008). A monochromatic X-ray beam is scattered by the protein sample in solution. At low (below 0.1 Å−1 ) and medium momentum transfer (s) (between 0.1 Å−1 and 0.25/0.3 Å −1 ) scattering angle, we obtain different information about the system (**Figure 5A**). Above 0.3 Å−1 the noise masks the signal and above 0.5 Å−1 data are collected at wide angle. This technique is not called SAXS anymore, but WAXS (Wide Angle X-ray Scattering) (Graewert and Svergun, 2013). At low s it is possible to extrapolate the radius of gyration Rg, which provide information about the size of the protein (Grant et al., 2015; Kikhney and Svergun, 2015). In order to obtain a reliable Rg, it is important to

have a precise measurement of the sample concentration before SAXS measurements. Medium s provides information about the shape of the protein. More precisely it is possible to obtain the Dmax, which provides measurement of the maximum dimension of the protein (Svergun, 1992).

From these measurements it is possible to build a lowresolution envelop, which represents the shape of the protein studied. The shape is reliable only when the system under study is monodisperse. In fact, it is always possible to obtain envelops from SAXS data, but poly-disperse samples do not

FIGURE 5 | SAXS profiles of the complexes. (A) The X-ray scattering patterns from IscS (1), IscU (2), CyaY (3) binary complexes IscS/IscU (4) and IscS/CyaY (5) and ternary complex IscS/CyaY/IscU (6). Plots display the logarithm of the scattering intensity as a function of momentum transfer (s). At the bottom it is highlighted with a blue box the low s region of the SAXS curve, in red the medium and in green the high. The experimental data are displayed as dots with error bars, the scattering from typical ab initio models computed by DAMMMMIF as full lines and the calculated curves from the high-resolution (for proteins alone) and rigid body models (for complexes) computed by CRYSOSOL/SASREF as dashed lines. The successive curves are displayed down by one logarithmic unit for clarity (figure adapted from Prischi et al., 2010a). (B) Table summarizing SAXS data. Rg is the radius of gyration; Dmax is the maximum size of the particle; MMSAXS is the molecular mass calculated from SAXS data; MMexp experimental molecular mass of the solute and χab and χrb values for the fit curves from ab initio models and from high resolution models (for proteins alone) and rigid body modeling (for complexes) using CRYSOL/SASREF, respectively.

generate envelops that represent the real shape of the protein. For example, in the study of PERK N-terminal domain, the protein was in a dynamic equilibrium between a dimer and a tetramer (Carrara et al., 2015). In this case it is not possible to obtain protein shape information, but SAXS is still informative. It had to be assumed that the resulting shape was a weighted average of dimer and tetramer shapes. The factor of weight had to be obtained from independent techniques, such as analytical ultracentrifuge (AUC), from which the relative populations of PERK dimer and tetramer in solution were estimated. The SAXS curve was then back-calculated using PERK dimer and tetramer crystal structures, weighted according to their relative abundance in solution, and fitted on experimental data. Agreement between experimental and back-calculated SAXS curves provided a confirmation that the PERK oligomeric structures were not a crystallographic artifact, but representative of the oligomeric state of PERK in solution (Carrara et al., 2015).

Luckily, IscS, IscU, and CyaY are all mono-disperse in solution. It was thus possible to obtain reliable information about their shapes from SAXS only. We collected SAXS data for each of the individual components, as well as for the binary (IscS-IscU, IscS-CyaY) and tertiary (IscS-IscU-CyaY) complexes (Prischi et al., 2010a). As previously mentioned, due to the dynamic nature of the Fe-S cluster machinery, the binding affinities of IscU and CyaY for IscS are relatively low: KdIscU-IscS = 1.3 ± 0.2 µM and KdCyaY-IscS = 18.5 ± 2.4 µM (Prischi et al., 2010a). We were able to isolate IscU-IscS complex using Size-Exclusion Chromatography (SEC) (Prischi et al., 2010b), but not CyaY-IscS and IscS-IscU-CyaY. In these cases we directly mixed proteins in solution prior data collection. Knowledge of relative Kd allowed us to estimate the optimal proteins ratios in order to maximize formation of the (Prischi et al., 2010a). It is worth mentioning that, despite not being available when we collected our data, a new methodology, particularly useful when collecting SAXS data on protein complexes, is Size-Exclusion Chromatography in line with SAXS (SEC–SAXS) (Mathew et al., 2004). SEC-SAXS is useful for separating pure systems that are under monomer-oligomer equilibrium or to further purify the sample before SAXS data are collected (particularly indicated for low stability protein which tend to form soluble aggregates). SEC–SAXS is available as a continuous-flow sample delivery option at BioCAT (Advanced Photon Source, U.S.A.) (Mathew et al., 2004), SWING (Soleil, France) (David and Perez, 2009), the SAXS beam line at the Australian Synchrotron, BM29 (ESRF, France), BL23A1 (NSRRC, Taiwan), B21 (Diamond, U.K.) and P12 (DESY, Hamburg) (Blanchet et al., 2015). SEC-SAXS has however limitations and it shouldn't be used as a purification step (Jeffries et al., 2016).

Ab initio envelops were generated using the software DAMMIF (Franke and Svergun, 2009). The high-resolution PDB structures 1P3W (Cupp-Vickery et al., 2003) for IscS and 1SOY (Nair et al., 2004) for CyaY were fitted into the SAXS envelops by rigid body modeling. Two different structures were used for IscU: one solved by NMR (PDB ID 1Q48) (Ramelot et al., 2004) and one by X-ray crystallography (PDB ID 2Z7E) (Shimomura et al., 2008). The two structures have a similar overall secondary structure content, but while in the NMR structure the first 25 residues are in a random coil conformation, the crystal structure is more compact and the N-Terminus forms a α-helix which makes contacts with α3 and the α5α6 loop. The quality of fitting of a 3D structure on a SAXS envelop can be visually ascertained, and can be more accurately estimated using the χ 2 (Svergun, 1999). χ 2 tells us how well the back-calculated scattering intensity from a 3D structure fits the experimental SAXS data. Fitting of the isolated CyaY and IscS resulted excellent, with a χ 2 of respectively 1.01 and 1.09 and an estimated molar mass of respectively 12 kDa ± 4 kDa (expected 12.231 kDa) and 85 kDa ± 10 kDa (expected 90.180 kDa) (**Figures 5A,B**). Of the two structures, 1Q48 fitted better the SAXS data collected for the isolated IscU in agreement with the dynamic nature of isolated IscU in solution (Kim et al., 2009; Prischi et al., 2010b), with a χ 2 of 1.03 and an estimated molar mass of 13 kDa ± 4 kDa (expected 13.849 kDa) (**Figures 5A,B**). It is also strongly recommended to check the residuals of the difference between experimental and back calculated SAXS data (i.e., whether these are random and not systematic).

### Step 4: Experimental Validation of the Models

The software DAMMIF (Franke and Svergun, 2009) was used to generate envelops from SAXS data. DAMMIF, a fast version of DAMMIN (Svergun, 1999), carries out an ab initio shape determination by simulated annealing using a single phase Dummy Atom Model (DAM). The DAM is represented by a tightly packed group of beads, which mimic, but do not resemble, real atoms. Each bead has a known scattering pattern and the software puts beads together so that the accumulated scattering resembles the experimental data. The software used to generates back-calculated curves and fit them on experimental data is the CRYSOL software (Svergun et al., 1995). CRYSOL requires a 3D structure/object as an input and then, taking into account the contribution for each atom, it evaluates the scattering intensity.

Despite having a 10–20 Å resolution (2π/smax), the SAXS envelops of the binary and tertiary complexes resulted evidently different from those of the single components. We first tested whether SAXS data were sufficient to generate meaningful binary and tertiary complexes using the SASREF software (Petoukhov and Svergun, 2005). SASREF tries to build the quaternary structure of a complex using the structures of the subunits and the solution scattering data. It is particularly useful because it can work with multiple data set(s), which allows working with SAXS data from sub-complexes and creating contrast series. SASREF build the complex structure without steric clashes using a simulated annealing protocol, which minimize differences between the experimental scattering data and the back-calculated SAXS curve of the model being built.

We inputted in SASREF the SAXS data and the highresolution structures of the single components but the complexes obtained with this approach did not generate reliable models, since they were not in agreement with our binding studies. Instead, we docked our HADDOCK structures into SAXS envelops. HADDOCK models were used to generate back-calculated SAXS curves, which were fitted on experimental data. Based on χ 2 , we selected the best fitting model, which was an experimentally verified model of FeS machinery complexes. Selecting the "best fitting model" could be problematic if two similar HADDOCK models have small orientation differences, which are clearly not distinguishable at SAXS resolution. In this context, HADDOCK is particularly well suited, because it first generates a maximum of 1,000 structures and then groups them according to their relative RMSD. By aligning these generated models using the interface residues of the first molecule, the RMSD (more correctly called interface-ligand RMSD) is calculated for the interface residues (less than 10 Å distance from the first molecule) of the second molecule (Dominguez et al., 2003). In our case, all structures HADDOCK grouped within the same group had RMSD < 7.5 Å. For each group, we used the structures with overall lower energy (evaluated by HADDOCK). Analysis of the binary complexes confirmed that CyaY sits near the IscS dimer interface and the active site, while IscU is located on the periphery of the IscS dimer and is aligned with the long axis of IscS. Accordingly, the IscS-IscU (Rg = 35 Å, Dmax = 121 Å) envelop is more elongated than the IscS alone (Rg = 31 Å, Dmax = 109 Å), while the IscS-CyaY envelop is more globular (Prischi et al., 2010a; Yan et al., 2013b; **Figure 4B**).

### Step 5: Comparison with X-Ray Crystal Protein Complexes

A limitation of this procedure is the absence of a tool to predict/model major structural changes upon formation of a complex. HADDOCK can simulate small conformational changes during the molecular dynamics refinement, but the final model strongly depends on the initial 3D structures provided: HADDOCK assumes a key-in-the-lock model. If a protein has significantly different structures in the free and bound states, HADDOCK (like any other protein docking software) will fail to generate a reliable model. For the IscS-IscU complex, we found that the model generated from HADDOCK did not fit the SAXS envelop as well as the single components did. However, a crystal structure of the IscS-IscU complex (PDB ID 3LVL) (Shi et al., 2010) became available while we were carrying out our studies. This structure is in perfect agreement with our NMR and mutant binding data and fits the SAXS envelop better than the HADDOCK model. This is due to IscU going through a structural rearrangement upon binding, with a formation of a α-helix in the N-terminus, similar to the one seen in 2Z7E (Shimomura et al., 2008). IscU has an optimal orientation for FeS cluster formation, with the surface containing three conserved cysteines pointing toward Cys328 in IscS loop (Shi et al., 2010). The distance between IscS active site and IscU is around 12 Å, suggesting the presence of major conformational changes happening during FeS cluster formation.

We then used 3LVL (Shi et al., 2010) for modeling the tertiary complex, IscS-IscU-CyaY. Interestingly, the model confirmed that, despite not being able to interact between each other directly, CyaY and IscU can interact once bound to IscS. This structure helped us to explain an inhibitory effect of CyaY on FeS cluster formation: enzymatic studies had showed that the tertiary complex is "less dynamic" than the binary ones with CyaY creating an additional anchoring point between IscS and IscU (Prischi et al., 2010a). This is in agreement with the observation that CyaY binding increases the affinity of IscU for IscS thus reducing the dissociation rates of the complex (the koff for the disassembly of the IscS/IscU complex is 0.8 s−1 in the absence of CyaY, vs. 0.006 s−1 in the presence of CyaY).

### A STEP FORWARDS: MOLECULAR DYNAMIC SIMULATIONS

From our studies it emerged that the dynamic nature of the ISC proteins is a key factor in their functions. To feature this dynamical behavior, we thus complemented our previous data with extensive (400 ns) molecular dynamic (MD) simulations (di Maio et al., 2017) of the IscS-IscU complex, both in the presence and in the absence of CyaY. We showed that the binary IscS-IscU complex is stably folded in line with our SAXS evidence (Prischi et al., 2010a), but IscU adopts a likely functionally relevant pivotal motion around the interface with IscS. This means that, despite being firmly attached to IscS, IscU maintains some degree of flexibility upon complex formation, which can be connected to their low binding affinity and the need of IscU to deliver FeS cluster to protein acceptors. At the same time, the pivotal motions observed in the MD simulations suggest that IscS-IscU interface is "fluid," with IscU side chains at the interface being trapped in several local minima. This was confirmed by NMR experiments (di Maio et al., 2017).

During the trajectory, the IscS catalytic loop containing Cys328 moves spontaneously and shifts from a mostly 310-helical structure to a β-turn/310-helix equilibrium, bringing Cys328 from 12 Å to 9 Å from the FeS cluster binding site on IscU (di Maio et al., 2017). This is in agreement with the IscS-IscU Xray structure of the A. fulgidus (PDB ID 4EB5) (Marinoni et al., 2012), which brilliantly captured the delivery stage of FeS cluster from IscS to IscU. In 4EB5, the FeS cluster is bound to the Cys of the IscS catalytic loop and is about 6 Å away from the IscU FeS cluster binding site (Marinoni et al., 2012).

In agreement with our previous studies (Prischi et al., 2010a), the simulations showed that the tertiary complex IscS-IscU-CyaY is more stable and that CyaY reduces the structural fluctuations of the IscS-IscU complex (di Maio et al., 2017). The most striking feature of the complex is the absence of motions of the IscS catalytic loops (one for each protomer) over the same timescale, due to CyaY steric hindrance and a salt bridge between CyaY Arg53 and IscS Glu334. This model brings us back to the beginning of this review, where we described the possible roles of FXN. The model we generated recapitulated in an elegant way our enzymatic data and provides a mechanistic explanation of how CyaY slows down FeS cluster formation (Adinolfi et al., 2009; Prischi et al., 2010a).

### EXTENSION OF THE METHODOLOGY TO OTHER ISC COMPLEXES

The approach described here has now been adopted by others (Kim et al., 2014, 2015) also to elucidate other ISC complexes,

FIGURE 6 | Ribbon representation of Fdx, YfhJ, and IscS interaction. IscS monomers are colored in pale cyan and light green with side chains of residues mutated indicated explicitly: R112E/R116E (orange), R220E/R223E/R225E (red). PLP in IscS active site is shown in black. Side chains of holo-Fdx residues exhibiting CSP (Ile54, Val55, Gln68, Glu69, Asp70, Asp71, Met72, Leu73, Asp74, Lys75, Ala76, Trp77, Gly78, Leu79, Glu80, Glu82) are explicitly shown in red. Fdx is loaded with a [2Fe-2S] cluster. YfhJ is colored in light blue and side chains of residues exhibiting CSP (Leu3, Lys4, Glu10, Ile11, Glu13, Ala14, Asp17, Leu58, Trp61, Leu62, Asp63, Glu64) are explicitly shown in blue.

#### TABLE 1 | Hybrid method breakdown.

Flowchart of the pros and cons of the different techniques part of the hybrid method adopted for the study of Iron-sulfur cluster machinery.

increasing the robustness of the methodology. A study of the complex between IscS and YfhJ was published (Kim et al., 2014). YfhJ behaves similarly to CyaY as it is able to bind both Fe (II) and Fe (III) using an electrostatic negative surface, which is the same area involved in IscS binding (**Figure 6**; Pastore et al., 2006). YfhJ also competesfor the same site of CyaY on IscS in agreement with previous mutation studies (Shi et al., 2010).

We have ourselves recently applied this hybrid method to model the IscS complex with Fdx, a FeS cluster dependent protein which is known to provide electrons to cellular reactions. Fdx is not-functional and devoid of tertiary structure in the absence of the cluster (Yan et al., 2013a). As for IscU, the spectrum of labeled Fdx disappears completely upon addition of unlabeled IscS. To circumvent the problem, we titrated 2H, 15N doublelabeled holo-Fdx with IscS using [2H,15N]-SOFAST HMQC NMR experiments (Yan et al., 2013b). This experiment requires a shorter acquisition time compared to HSQC and is thus more suitable for unstable samples. We could then identify the residues of Fdx involved in IscS binding cluster, which reside near a uniform acidic patch on the α2-α3 loop (**Figure 6**; Yan et al., 2013b).

NMR competition studies revealed that Fdx and CyaY compete for the same site of IscS (Yan et al., 2013b). A Fdx-IscS SAXS verified model showed that Fdx sits in a position similar to that of CyaY near the active site. This was utterly validated by creating a new IscS mutant (IscS\_R112E/R116E), which interferes with Fdx binding (Yan et al., 2013b). Assuming that the two proteins exploit their functions in different times during the cluster biogenesis, competition could represent a fascinating regulation mechanism. Superimposition of the Fdx-IscS and IscS-IscU models reveals that the Fdx C-terminus (which contains two key residues for electron transfer reactions, Tyr101 and His105) points toward the interface between IscS and IscU. This nicely explains how, after production of S0 from Lcysteine by IscS, Fdx could reduce S0 to S2− (Yan et al., 2013b, 2015).

To add surprise to surprise, we have more recently shown that also the co-chaperone HscB binds to IscS in the same binding pocket, a result further validated by cross-linking experiments (Puglisi et al., 2016). This implies a picture in which IscS acts as a central platform on which several of the other bacterial ISC proteins bind and typically form 1:1 complexes (Pastore et al., 2006; Adinolfi et al., 2009; Prischi et al., 2010a; Yan et al., 2013b). It remains for us to understand why and how several different components of the same pathway compete for the same site. We suggested that this is a regulatory significance, which could operate through allosteric responses and involve the binding sites on each of the protomers present in the IscS dimer.

#### CONCLUSIONS

In conclusions, we have in this review gone through a methodology (**Table 1**), which has allowed us to gain information on a 110 kDa complex with hybrid techniques. The method can in principle be applied also to larger complexes. The most successful cases are anyway those which involve an appreciable charge of shape of the complex, leading to a clear difference of the SAXS envelop between the isolated components and the complex. Limitations are currently dictated by the number of restraints available and by their distance tolerance: restraints which can allow a tolerance of more than 11 Å, as it is the case for cross-linking studies, are informative but only if several distances are available. It would also be useful to develop HADDOCK and other software to deal with the specific problems

#### REFERENCES

of hybrid methods. Some attempts along this line have already been made but more effort would be welcome in the future. It appears particularly useful, in a future perspective, to flank NMR and SAXS studies to other techniques, such as fluorescence, isothermal calorimetry, AUC and cross-linking to obtain more complete and complementary information. As a word of caution though, very good care should anyway be paid to the validation of the results. False positives can be easily obtained if assuming the presence of the wrong species in solution. It remains nevertheless clear that hybrid methods have open a new perspective to the size and complexity of the complexes which can be studied by Structural Biology and, more importantly to the possibility of tackle not only stable and rigid assemblies but also weakly interacting dynamical machines.

### AUTHOR CONTRIBUTIONS

All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

The work was supported by MRC (U117584256).

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Prischi and Pastore. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Chlorophyll-Derivative Modulation of Rhodopsin Signaling Properties through Evolutionarily Conserved Interaction Pathways

#### Kristina N. Woods 1 \*, Jürgen Pfeffer 2 and Judith Klein-Seetharaman3

1 Lehrstuhl für BioMolekulare Optik, Ludwig-Maximilians-Universität, München, Germany, 2 Bavarian School of Public Policy, Technical University of Munich, München, Germany, 3 Warwick Medical School, University of Warwick, Coventry, United Kingdom

Retinal is the light-absorbing chromophore that is responsible for the activation of visual pigments and light-driven ion pumps. Evolutionary changes in the intermolecular interactions of the retinal with specific amino acids allow for adaptation of the spectral characteristics, referred to as spectral tuning. However, it has been proposed that a specific species of dragon fish has bypassed the adaptive evolutionary process of spectral tuning and replaced it with a single evolutionary event: photosensitization of rhodopsin by chlorophyll derivatives. Here, by using a combination of experimental measurements and computational modeling to probe retinal-receptor interactions in rhodopsin, we show how the binding of the chlorophyll derivative, chlorin-e6 (Ce6) in the intracellular domain (ICD) of the receptor allosterically excites G-protein coupled receptor class A (GPCR-A) conserved long-range correlated fluctuations that connect distant parts of the receptor. These long-range correlated motions are associated with regulating the dynamics and intermolecular interactions of specific amino acids in the retinal ligand-binding pocket that have been associated with shifts in the absorbance peak maximum (λmax) and hence, spectral sensitivity of the visual system. Moreover, the binding of Ce6 affects the overall global properties of the receptor. Specifically, we find that Ce6-induced dynamics alter the thermal stability of rhodopsin by adjusting hydrogen-bonding interactions near the receptor active-site that consequently also influences the intrinsic conformational equilibrium of the receptor. Due to the conservation of the ICD residues amongst different receptors in this class and the fact that all GPCR-A receptors share a common mechanism of activation, it is possible that the allosteric associations excited in rhodopsin with Ce6 binding are a common feature in all class A GPCRs.

Keywords: Terahertz spectroscopy, protein allostery, small-molecule allosteric agonist, Chlorine Compounds, protein dynamics

### INTRODUCTION

The dimmest habitats on earth appear at night and in the depths of the ocean (Warrant, 2004). The greatest challenge for vision in these habitats is capture of photons, and the way these photons are post-processed. The primary photoreceptor in eyes, rhodopsin, is a major target for adaptation to different light conditions. Rhodopsin as the prototypical member of the GPCR-A family adopts an

#### Edited by:

Isabelle Callebaut, Centre National de la Recherche Scientifique (CNRS), France

#### Reviewed by:

Simone Brogi, University of Siena, Italy Steffen Wolf, Albert Ludwigs University of Freiburg, Germany

> \*Correspondence: Kristina N. Woods<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biophysics, a section of the journal Frontiers in Molecular Biosciences

Received: 20 February 2017 Accepted: 28 November 2017 Published: 12 December 2017

#### Citation:

Woods KN, Pfeffer J and Klein-Seetharaman J (2017) Chlorophyll-Derivative Modulation of Rhodopsin Signaling Properties through Evolutionarily Conserved Interaction Pathways. Front. Mol. Biosci. 4:85. doi: 10.3389/fmolb.2017.00085

**59**

overall organization of seven transmembrane (TM) helices that form a bundle. Within the GPCR family there is a large sub-group of opsins, representing opsin sequences in different photoreceptor cell types and organisms. All opsins covalently bind retinal, a vitamin A derivative, at the interface between the transmembrane (TM) and extracellular domains (ECDs). Visual signal transduction is initiated by photon-induced isomerization of 11-cis retinal to all-trans retinal. This event is sensed by the TM domain which undergoes a conformational change that results in the activated state, Metarhodopsin II (Meta-II) via Meta-I. Meta-II—unlike dark state rhodopsin—binds and activates the G protein, transducin Gt , ultimately leading to receptor hyperpolarization. Signal desensitization is initiated by phosphorylation of the C-terminus of rhodopsin by rhodopsin kinase, followed by binding of arrestin, preventing further binding of Gt to rhodopsin. In vitro, Meta-II decays to opsin and free retinal, the half-life of which depends on the lipid environment (Farrens and Khorana, 1995). In addition, a storage form of rhodopsin, Meta-III, emerges in parallel with Meta-II, from the Meta-I state. The Meta-III state of rhodopsin—unlike Meta-II—has a protonated retinal Schiff base and decays into opsin and free retinal on significantly longer time scales (Heck et al., 2003b,a; Vogel et al., 2003, 2004; Stehle et al., 2014).

One mechanism of adaptation to dim habitats is spectral tuning. Spectral tuning refers to adjustment of the absorbance maxima in the spectra of the photoreceptors. Spectral tuning can be hard-wired by genetic variation of the rhodopsin sequences changing the interactions between the retinal and the protein ("opsin shift") (Nathans, 1990). This allows adaptation to the wavelengths that are maximally transmitted under different environmental conditions, such as sun-light in the shade, sun-light penetrating water, moon-light or bioluminescence generated by deep-sea fish (Douglas et al., 1998; Fishkin et al., 2004). Genetic adaptation is also the mechanism by which cone and rod opsins absorb at different wavelengths (Douglas et al., 1999). However, spectral tuning can also be achieved by interaction of the rhodopsins with small molecules. In xanthorhodopsin, a bacteriorhodopsin-like proton pump in the halophilic eubacterium Salinibacter ruber, a carotenoid is bound in addition to retinal (Balashov et al., 2005). Light energy absorbed by the carotenoid is transferred to the retinal with a quantum efficiency of ∼40% (Balashov et al., 2005) and light is funneled to the retinal similar to the photosynthetic light harvesting complex. In the deep-sea fish Malacosteus niger (dragon fish), porphyrin/chlorophyll derivatives, including chlorin e6 (Ce6), may act as photosensitizers rendering its rhodopsin sensitive to wavelengths that other deep-sea fish respond to as a result of genetic adaptations of their retinal binding pockets or retinal replacement (Douglas et al., 1998, 1999; Isayama et al., 2006). An investigation of different chlorophyll derivatives in rod outer segments on bleaching rates of the chromophore by red light highlighted Ce6 as the molecule with the strongest photosensitizing effect on bovine rhodopsin (Washington et al., 2004). The hypothesis of a binding pocket for Ce6 and energy transfer to retinal resembling that of light harvesting in photosynthesis was proposed (Washington et al., 2004; Isayama et al., 2006). The photosensitizing effects of chlorophyll-derivatives observed in salamander and deep-sea fish can be reproduced qualitatively using bovine rhodopsin, in vitro (Washington et al., 2004). Furthermore, electroretinogram recordings in mice injected with Ce6 suggested that the sensitivity of rhodopsin can be broadened to blue and red light in the presence of Ce6 (Washington et al., 2007). Indeed, the photosensitizing effect of Ce6 and other porphyrin compounds in vision has been reported as a side-effect during photodynamic therapy (Kimura, 1987).

To understand by what Ce6 binding enhances bleaching rates of rhodopsin in red light, we used a novel approach which we recently developed to probe dynamics and allostery in rhodopsin (Woods et al., 2016): through a combination of molecular dynamics simulations, evolutionary sequence analysis and Terahertz (THz) spectroscopy we show that Ce6 binding excites evolutionarily conserved communication pathways in rhodopsin that establish connections between the ligand-binding site and the rest of the receptor.

### RESULTS

#### Experimental Detection of the Conformational Ensemble Dynamics and Intermolecular Changes in Rhodopsin When Bound with the Allosteric Modulator Ce6

#### Global Fluctuations of Dark-State Rhodopsin Bound with Ce6

The <100 cm−1 region of the infrared spectrum is sensitive to global, internal fluctuations that describe the intrinsic dynamics of the receptor. These globally, correlated thermal fluctuations provide a mechanism for sampling the ensemble of conformations that comprise the free energy landscape of possible receptor conformations (Frauenfelder et al., 1991). Hence, the modes detected in the experimental THz spectrum in this region provide direct information about the sampling of conformational substates in rhodopsin. An inspection of the low frequency modes of dark-state rhodopsin in the <100 cm−1 spectral region in **Figure 1A** in the presence and absence of Ce6 reveals that the region of the spectrum above 50 cm−1 is dramatically altered when contrasting the two states of the receptors. For instance, in the Ce6-bound receptor in **Figure 1B** there is a prominent peak at 80 cm−1 and a general increase in the absorption peak intensity at a frequency above 50 cm−1 when contrasted with the dark-state receptor that is not bound with Ce6 (**Figure 1** and Woods et al., 2016). The differences in the spectra of the two receptors suggest that binding of Ce6 alters the conformational ensemble dynamics in rhodopsin.

In our previous work on the inactive-receptor in the unbound state (Woods et al., 2016), we have found there is an equilibrium of both inactive and active-state protein conformational fluctuations in the dark-state protein. This conformational heterogeneity in the inactive receptor indicates that rhodopsin samples a diverse set of functional structures even before any activation event has taken place. In particular, the heterogeneity of global structural fluctuations detected experimentally is

spectral region.

intimately tied with multiple pre-existing allosteric associations in the inactive receptor. Peaks at approximately 75 and 55 cm−1 in the experimental spectrum of the unbound receptor were identified as rhodopsin global fluctuations in an inactive-type conformation. Specifically, the peak at 55 cm−1 was found to be associated with a retinal (polyene chain) torsional fluctuation that is coupled with a protein, global backbone torsion, whereas the peak at 75 cm−1 is associated with a retinal torsional oscillation that is coupled with collective out-of-plane protein side-chain fluctuations (Woods et al., 2016). In addition to the inactivestate conformational fluctuations, we also uncovered a weaker band at approximately 40 cm−1 in the experimental spectrum that was later found to be associated with transient interactions of receptor amino acid side-chains that supported a more activelike rhodopsin conformation. These retinal-induced transient interactions of the more active-like receptor conformation were hypothesized to serve as a necessary precursor in the mechanism that eventually leads to the active-state receptor. Therefore, the enhancement of the ∼75 cm−1 mode in the Ce6-bound receptor

spectral region and in the (D) 100–20 cm−1

suggests that Ce6 stabilizes a specific conformational state of rhodopsin—the ground state.

#### Ce6-Induced Intra- and Inter-Protein Interaction Networks in Dark-State Rhodopsin

In the 100–250 cm−1 region of the experimental spectrum we detect motions in rhodopsin that reveal more localized intermolecular interactions (Woods, 2014b) such as interhelical contacts as well as helical interactions with the solvent. In general, peaks in the experimental spectrum in the 100–160 cm−1 are related to protein intra- and intermolecular interactions. For example in **Figure 1A**, the peaks at approximately 150 and 140 cm−1 in the dark-state spectrum of the unbound receptor were found to be associated with interhelical and solvent-induced H-bonding interactions in our earlier study (Woods et al., 2016), while peaks ≥170 cm−1 are predominately associated with solvent-solvent interactions of water molecules in the receptor hydration shell. A comparison of dark-state rhodopsin in the unbound and unbound state in **Figure 1A** suggests that

line) in the (C) 240–100 cm−1

binding of Ce6 alters the inter- and intra-protein interactions in rhodopsin but has a much stronger effect on an extensive network of solvent H-bonds that stabilize the ground state of the receptor (Woods, 2014a).

### MD simulation of Ce6-Induced Long-Range Correlated Fluctuations and the Effect on Dark-State Rhodopsin Global Motions

Allostery in proteins enables the activity of one site in a protein to modulate function at another spatially distinct region (Hawkins and McLeish, 2006; Fenwick et al., 2011; Motlagh et al., 2014). Recent experimental and computation investigations on a number of proteins and enzymes have demonstrated that allosteric signal transmission is mediated by protein local structural fluctuations (Whitten et al., 2005; Daily and Gray, 2007; Pandini et al., 2013). To assess if Ce6 allosterically affects retinal-protein interactions in rhodopsin, we carried out MD simulations of dark-state bound rhodopsin in the presence and absence of Ce6. We find that Ce6 binds only weakly on the cytoplasmic surface (Supplementary text, Ce6–ligand binding affinity and Figures S2, S3) of the receptor. Experimental fluorescence measurements of Ce6 binding to rhodopsin also reveal micromolar binding affinity (unpublished results). Despite the weak binding, Ce6 has a strong effect on the long-range correlated fluctuations of the bound state when contrasted with the unbound receptor. This can be deduced from our analysis of the MD simulation induced localized structural fluctuations (LSFs) and the consequent collective dynamics in the receptor that arise from the addition of Ce6 to rhodopsin. For instance, **Figure 2** reveals the development of long-range correlated fluctuations in dark-state rhodopsin with Ce6 that significantly involve contributions from the intracellular loops close to where Ce6 binding takes place. Particularly, Ce6-induced long-range correlated fluctuations involving intracellular loop 2 (CL2) and extracellular loop 2 (EL2) are coupled with a structurally conserved collection of aromatic and polar residues

at the extracellular end of H4 (Pro170, Pro171, Tyr175, Ser176, Arg177, and Tyr178) leading from the ligand-binding pocket to the EC domain. This group of GPCR-A wide conserved longdistance fluctuations modulates the H-bonding environment of residues and solvent molecules lining the retinal ligandbinding pocket (**Figure 2B**). Specifically, side-chain fluctuations of residues such as Glu122, Trp126, and Phe261 in Ce6 rhodopsin have increased (H-bonding) interactions with the retinal β-ionone ring due to the CL2–EL2 correlation when contrasted with unbound receptor (Figure S1b). The Ce6 induced correlated fluctuations reorient Glu122 so that the there is a rearrangement of the H-bonding network surrounding the ring such that it stabilizes the retinal in a defined conformation (**Figure 2** and Figure S4). In our previous work (Woods et al., 2016) we have conjectured that the ability of the retinal to assume multiple conformational states is intertwined with the formation of multiple, distinct signaling pathways in both the inactive and active receptor. In line with what was observed in the experimental detection of the global modes of the dark-state receptor in **Figure 1**, we find from the MD analysis of the Ce6 induced structural fluctuations in rhodopsin, that binding of the allosteric modulator assists in stabilizing the receptor in a specific conformation, namely the ground state because of the modified electrostatic interactions in the vicinity of the retinal ring that accompany Ce6-binding. This conclusion is further supported by a mapping of the low-frequency torsional dynamics of the retinal from the dark-state MD simulations of rhodopsin in both the unbound and bound states. The absence of the 40 cm−1 mode in the MD retinal torsional spectrum of Ce6-bound rhodopsin in **Figure 3A** clearly shows that Ce6 stabilizes only a subset of the receptor conformations. Both the 80 and 65 cm−1 represent the effect of retinal torsional dynamics from interactions with the inactive-like conformation of rhodopsin, whereas the 40 cm−1 mode reflects the influence on the retinal from protein interactions in a more active-like conformation (Woods et al., 2016).

Ce6-induced long-distance induced fluctuations involving CL1 and CL3 are instrumental in altering the interhelical packing within the receptor core. Correlated fluctuations involving these two loops modulate the dynamics of TM1and TM2 (**Figure 2**) and subsequently destabilize a conserved cluster of water molecules near Gly90 in helix 2 that maintain the shape of the ligand-binding pocket. In particular, the Ce6-induced distortion of the binding pocket weakens interactions in the ligand-binding region that stabilize residues forming the receptor hydrophobic core such as Gly114, Ala117, Thr118, and Gly120– Glu122 and concurrently, also disrupts a conserved network (Angel et al., 2009) of water molecules near Gly90 in helix 2 that stabilizes the Glu113 salt bridge of the protonated Schiff base (PSB). These Ce6-induced changes within the receptor ligandbinding pocket allow Glu113 to move further away from the PSB. This perturbation of the PSB H-bonding network increases the amplitude of the dynamic structural fluctuations of the counterion and affects the thermal stability of the entire receptor. The reduced electrostatic interactions from residues surrounding the PSB, resulting from the movement of Glu113 away from the PSB, decreases the energy gap between the ground and excited

FIGURE 3 | The torsional spectrum of the retinal from rhodopsin MD simulations in (A) the dark-state (blue line) and in the dark-state when bound to Ce6 (green line) and (B) in Meta-II (cyan line) and Meta-II when bound to Ce6 (purple line). The peaks at 15 and 25 cm−1 in (B) are related to torsional fluctuations of the C-20 methyl group near the terminus of the polyene chain and a collective chain-twisting oscillation, respectively in Meta-II whereas the peak close to 60 cm−1 in (B) is associated with a chain torsion coupled with a retinal ring bending motion of a more inactive-type rhodopsin structure. In all cases, the torsion of the retinal is defined by the angle created by the C5-, C9-, and C13- methyl groups.

state of the receptor (Rajamani et al., 2011). The reduction of the energy difference between the two states also reduces the barrier for activation; hence, it reduces receptor thermal stability (Figure S4). Numerous previous studies on rhodopsin (Lin et al., 1998; Rajamani et al., 2011; Imamoto and Shichida, 2014) have also shown that the reduction in the energy difference between the ground and excited state reduces the thermal stability of the receptor and is also directly correlated with the red-shift of the λmax of the receptor. Although, it is important to point out that previous thermal denaturing studies (Balem et al., 2009) on inactive-state rhodopsin found an increase in the helical content of the receptor when bound to Ce6. In other words, the helical regions of the receptor are stabilized in the dark-state receptor when bound to Ce6 (when contrasted with the unbound receptor). These results are in direct contrast to both the MD simulation studies and (THz) experimental studies carried out in this investigation. Although, the disparity in interpretations could be a consequence of what the distinctive experimental methods measure. The melting data is primarily sensitive to changes in the helical content of the receptor, whereas the THz data (≤100 cm−1 ) is sensitive to the global changes in the receptor as a whole. The initial findings from this investigation suggest that the principal changes in thermal stability arise from modifications around the ligand-binding site in the Ce6-bound receptor and these alterations influence the global stability of the receptor.

The deformation of the retinal ligand-binding pocket due to the Ce6-induced correlated fluctuations of CL1 and CL3 is further stabilized by the rearrangement of conserved network of water molecules in EL2. An analysis of the CL1-CL3 Ce6 induced correlated dynamics in **Figure 2** also reveals alterations in solvent-protein interactions that include the water molecules shared between the side-chains of Glu181 and Ser186 that directly connect the dynamics of the EL2 with the retinal bindingpocket (Figure S5). The CL1-CL3 coordinated dynamics shift Glu181 further away from Ser186 by means of the introduction of an additional water molecule into the retinal binding pocket that also links the backbone atoms of Glu181 and Ser186 (while retaining the water coordinated side-chain linkage of the two residues). This results in the overall rearrangement of EL2 through the correlated dynamics of helix 3 via the conserved Cys110-Cys187 disulfide bond. The correlated movement shifts Ile189 (as well as the entire β4 loop of EL2) in the direction of the receptor N-terminus. The end result of the altered Ce6 induced dynamics of dark-state rhodopsin is a more open ligandbinding pocket that is created by the larger separation between the retinal and extracellular region of the receptor. The increased distance between Ile189 and the retinal disrupts the H-bonding network of the PSB but in doing so, would also likely lead to a decrease in receptor thermal stability (associated with higher levels of dark-noise) as well as increase the rate of hydrolysis of the PSB. In fact, earlier experimental studies (Imamoto and Shichida, 2014; Yanagawa et al., 2015) focusing on the thermal activation rate of rhodopsin have revealed a direct correlation between the thermal stability of the dark-state of the receptor and that of the lifetime of the active-state intermediate Meta-II. In this regard, both Glu122 and Ile189 were identified as two residues that play a crucial role in suppressing thermal fluctuations in the retinal-binding pocket in rhodopsin, which ultimately imparts the low dark noise characteristic of the receptor.

It is also interesting to return to the experimental detection of the global modes of the dark-state receptor (**Figure 1A**) and to consider the observed changes in the conformational ensemble dynamics that were detected with Ce6-binding. Previous experimental studies probing the conformational stability of Ce6 in (dark-state) rhodopsin have indicated that Ce6 stabilizes the helices of the receptor. This earlier work is in line with the experimentally detected population shift in the conformational ensemble dynamics of the Ce6-bound receptor (when contrasted with the unbound receptor) that we have detected in this investigation on rhodopsin (**Figure 1A**). But interestingly, is in direct contrast to the identified changes in the ligandbinding pocket dynamics of the bound receptor that we have uncovered from MD simulation. In the latter case, the Ce6 induced dynamics appears to create a closer potential energy surface between the ground and excited state (Ala-Laurila et al., 2004; Hofmann and Palczewski, 2015) of the retinal in Ce6 rhodopsin. This reduction in the energy difference between the two states would, in principal, create a thermally unstable ground state that consequently would shift the equilibrium toward the excited-state of the receptor. The disparity in interpretations that arise when examining the local dynamics of the retinal-binding pocket with approaches that map the global characteristics of the receptor suggests that induced long-range structural fluctuations may also play an important role in the Ce6 mechanism of spectral tuning in rhodopsin.

### Experimental Detection of the Induced Changes in Global Dynamics and the Allosterically Coupled Motions Associated with Activation in Meta-II Rhodopsin When Bound by Ce6

Global Dynamics and Thermal Stability in Meta-II-Ce6 The global motions of Meta-II rhodopsin in the presence of Ce6 are dramatically altered as compared to dark rhodopsin with and without Ce6 in **Figure 1D**. In contrast to the dark-state unbound spectrum in **Figure 1B**, there are no clearly resolved vibrational bands in the <100 cm−1 region of the experimental spectrum. Based on previous work on rhodopsin (Woods et al., 2016), this likely indicates that the major protein modes excited after isomerization are red-shifted to very low frequencies. Despite the lack of spectral features, one can still deduce information about the stability of the receptor states based on the general shape of the spectrum. For instance, it is apparent that the Ce6-bound Meta-II state in **Figure 1C** is far more oscillatory above 70 cm−1 when contrasted with the unbound Meta-II state. The oscillatory nature of the Meta-II-Ce6 spectrum indicates that there is instability in the internal modes of the receptor. Moreover, the increased instability detected in the global modes of the Ce6-bound receptor also implies a general decrease in the thermal stability of the Ce6-bound activated-state of the receptor.

#### Experimental Detection of Alterations in the Allosteric Interaction Network Associated with Activation in Meta-II-Ce6

An inspection of the higher frequency spectra of the Meta-II states of rhodopsin in **Figure 1C** strongly suggests that there are major modifications in the interhelical interactions in Meta-II-Ce6 when compared with the unbound active-state receptor. For example, the absence of the ∼140 and 150 cm−1 modes in the Ce6 receptor bound spectrum suggests that both the interhelical packing and the solvent-protein interactions in the bound receptor are dramatically altered in Meta-II-Ce6. Both modes (at 140 and 150 cm−1 ) are strongly influenced by retinal torsional dynamics associated with both the β-ionone ring and the polyene chain, indicating that the coupling between the retinal ligand and its immediate protein environment is considerably changed in the Ce6-activated state. In our earlier study (Woods et al., 2016) of Meta-II rhodopsin we determined that the vibrational modes at approximately 110, 120, and 130 cm−1 are spectral markers of the water-mediated pathway of activation in Meta-II that connects the dynamics in the ligandbinding pocket with the G-protein binding site in the ICD. The vibrational modes represent anharmonic, solvent-mediated fluctuations of protein backbone and side-chain atoms that form allosteric sites in Meta-II. The water-mediated allosteric sites create a coherent, defined pathway of signal propagation in Meta-II that connect the dynamics of the retinal taking place in the ligand-binding pocket with intermolecular interactions taking place in the receptor intracellular domain. The fact that both spectra in **Figure 1C** have prominent vibrational bands at about 110, 120, and 130 cm−1 implies that activation in the Ce6 bound receptor still takes place, but based on the other prominent differences in the spectrum of the distinct Meta-II states, it is likely that the activation mechanism in Meta-II-Ce6 is somehow altered. One clue of how the pathway of activation in Meta-II-Ce6 may differ from Meta-II can be discerned from the solvent—solvent interaction region of the experimental spectrum (>170 cm−1 ) in **Figure 1C**. A comparison of the spectra makes it readily apparent that solvent interaction network in Meta-II is disrupted with Ce6 binding. The dramatic drop in intensity in Meta-II-Ce6 in the 200–240 cm−1 spectral region suggests that the addition of the allosteric modulator (Ce6) disrupts the solvent H-bonding network that both stabilizes and aids in the efficiency of the active Meta-II intermediate binding with the G-protein. Although interestingly, in the same solvent-interaction region, there is a new peak in the Meta-II-Ce6 spectrum centered at 190 cm−1 that possibly hints at a new set of solvent interactions that support an alternative pathway of activation in the Ce6 bound receptor. In previous THz studies on globular proteins (Woods, 2010, 2014a) we have identified an equivalent peak at approximately 190 cm−1 in the protein-solvent coupling region of the experimental spectrum that describes interaction dynamics directly tied with the formation of long-range communication channels in the protein. Particularly, in the previous instances we found that the anharmonic dynamics of the solvent molecules within the hydration shell coupled with protein motions to promote long-range coherence pathways in the protein threedimensional structure.

#### MD simulation of Meta-II bound with Ce6

Localized Structural Fluctuations in Meta-II-Ce6 and the Disruption of the Active-State Activation Pathway A network representation of the Meta-II LSFs from the MD simulations of Ce6-bound and -unbound rhodopsin structures are shown in **Figures 4A,B** and a mapping of the LSF components onto the two Meta-II structures is shown in **Figures 4C,D**. One of the distinguishing features that we detect in the Meta-II-Ce6 LSF when contrasted with unbound Meta-II is a disruption in the structural overlap of interactions that creates an activation pathway from the ligand-binding region to the cytoplasmic surface. Specifically, in the Meta-II-Ce6 LSF in **Figure 4B** we find that there is a disconnect between the network of interactions linking rearrangements taking place in the CWxP motif (in the retinal-binding pocket) with a conserved intracellular pathway of intermolecular associations consisting of: a conserved network of internal water molecules (Nygaard et al., 2010), conserved residues in TM1 and TM2, and residues comprising the NPxxY motif in helix 7. This conserved pathway dynamically connects residues from the protein interior (in contact with the retinal) with residues in the cytoplasmic surface that are crucial for GT binding. In **Figure 4** and Figure S6, a comparison of the unbound and bound receptor structures from the MD simulation reveals that key water molecules involved with the activation pathway have been displaced in the Meta-II-Ce6 structure. The water-mediated Hbonding connection between Ala260 and Asn302 is disrupted in Meta-II-Ce6 and consequently, the connection between Trp265 in helix 6 and Asn302 in helix 7 is weakened. The modification in the conserved water-mediated network in Meta-II-Ce6 effects the rearrangement of H-bonding associations at the intracellular side of the receptor that form the G-protein binding site in the active state. Specifically, the rupture of the Ala260–Asn302 watermediated H-bond linkage diminishes the correlation between the rotation "toggle" switch of the CWxP motif comprised of conserved residues Tyr268, Trp265, and Phe261 on helix 6 and the structural changes that take within the intracellular core of the receptor involving Met257, Tyr306, and Tyr223 that lead the "breaking" of the ionic-lock between Arg135 on helix 3 and Glu247 on helix 6. The ionic-lock connects the intracellular side of TM3 and TM6 in the active state and breaking of the lock is a crucial step in forming the G-protein binding site. In the MD simulations carried out in this study, we find evidence that the ionic-lock is not fully stabilized in the open form (broken) in the active-state of the Ce6-bound receptor (Figure S7).

#### Ce6-Induced Correlated Fluctuations and the Creation of an Altered Activation Pathway in Meta-II-Ce6

The reason for the disruption in the activation pathway becomes more apparent when analyzing the induced interactions and correlated motions that accompany Ce6 binding in Meta-II-Ce6. In **Figure 5**, we find that Ce6-induced correlated fluctuations introduce new LSFs in the vicinity of the receptor ligand-binding pocket. An illustration is provided in a comparison of the correlated dynamics in both receptors involving helix 4 near the ligand-binding pocket in **Figures 4C,D**. In the unbound, activestate receptor there is a strong correlation between residues in helix 4 with residues residing in the N-terminus region of the receptor. The helix 4–N-terminus correlated motion is associated with stabilizing the activation pathway that connects the dynamics taking place within the retinal pocket with the dynamics occurring in the cytoplasmic side of the receptor. For instance, in a previous computational study on squid rhodopsin it was revealed that Ala167 (in helix4) in addition to Ala304 (helix 7) and internal water molecules in the intracellular region of helix 7 are instrumental in creating a maximally connected H-bonding pathway that links the active-state retinal binding pocket with the cytoplasmic region (Bondar et al., 2011).

The absence of the (helix 4 - N-terminus) long-range connection in Meta-II-Ce6 suggests that the Ce6 bound receptor may support an alternate pathway for activation. For instance, in **Figure 4A** we observe that the binding of Ce6 excites new collective fluctuations in Meta-II that involve residues in CL3

cartoon representation of the 3-D structure of (C) Meta-II and (D) Meta-II-Ce6.

near the Ce6-binding site that are coupled with collective fluctuations of residues in both EL2 and EL3. The Ce6-induced dynamical fluctuations have a prominent effect on the shape of the ligand-binding pocket and consequently, the allosteric interactions that determine the signal communication pathway from the ligand-binding site to the rest of the protein. Specifically, a comparison of the LSFs of the unbound and bound active-state receptor in **Figures 4A,B** reveals that Ce6-induced fluctuations disrupt long-range correlations between EL1 and residues that constitute the hydrophobic core of the receptor, particularly Ala117, Thr118, and Gly120–Glu122. In the unbound receptor, this long-range correlation is associated with maintaining the shape of ligand binding pocket (**Figures 4A,C**). The disruption EL1–receptor core correlated dynamics in Meta-II-Ce6 promotes a deformation of the ligand binding region that consequently allows residues surrounding β-ionone ring, such as Glu122 and Pro215 on helix 5, to shift closer to the retinal as well as to the extracellular side of helix 6. This shift simultaneously alters a stable network of polar interactions within the core of the receptor that connects the ligand-binding site (involving helices 3, 5, and 6) with the rest of the protein. Explicitly, the shift of helix 5 ligand-binding residues promotes a shift in the hydrophobic interactions involving Gly121 and Leu125 on helix 3 and Phe261 on helix 6. The adjustment of the hydrophobic packing interaction allows the intracellular side of the helix (helix 3) to move closer to helix 6. Further, the helix 3 shift is connected with a correlated fluctuation involving the motion of residues Phe261, Trp265, and Tyr268 on helix 6 with that of EL2, such that the movement of these residues promotes an overall tilt of the IC side of H6 and as a consequence also reduces the distance between TM3 and TM6 (**Figure 6**). The outcome is interhelical packing changes in the core of the molecule that shift TM2 and TM8 upward in the direction of the C-terminus of the receptor and tilts TM7 inward toward the hydrophobic core. Consequently, the altered ligand-binding cavity shape also supports a 7TM structure with a weakened but still associated