Free energy landscape remodeling of the cardiac pacemaker channel explains the molecular basis of familial sinus bradycardia

The hyperpolarization-activated and cyclic nucleotide-modulated ion channel (HCN) drives the pacemaker activity in the heart, and its malfunction can result in heart disorders. One such disorder, familial sinus bradycardia, is caused by the S672R mutation in HCN, whose electrophysiological phenotypes include a negative shift in the channel activation voltage and an accelerated HCN deactivation. The outcomes of these changes are abnormally low resting heart rates. However, the molecular mechanism underlying these electrophysiological changes is currently not fully understood. Crystallographic investigations indicate that the S672R mutation causes limited changes in the structure of the HCN intracellular gating tetramer, but its effects on protein dynamics are unknown. Here, we utilize comparative S672R versus WT NMR analyses to show that the S672R mutation results in extensive perturbations of the dynamics in both apo- and holo-forms of the HCN4 isoform, reflecting how S672R remodels the free energy landscape for the modulation of HCN4 by cAMP, i.e. the primary cyclic nucleotide modulator of HCN channels. We show that the S672R mutation results in a constitutive shift of the dynamic auto-inhibitory equilibrium toward inactive states of HCN4 and broadens the free-energy well of the apo-form, enhancing the millisecond to microsecond dynamics of the holo-form at sites critical for gating cAMP binding. These S672R-induced variations in dynamics provide a molecular basis for the electrophysiological phenotypes of this mutation and demonstrate that the pathogenic effects of the S672R mutation can be rationalized primarily in terms of modulations of protein dynamics.

ing from the apo-CBD, the C-linker forms a stable IR tetramer, and the tonic HCN auto-inhibition is removed (9).
The coupling between cAMP binding, CBD activation, and channel gating is of central importance to the physiological function of HCN (9,14). Mutations that perturb or disrupt this coupling cause pathological conditions for the heart and pacemaker current (11,(15)(16)(17)(18)(19). For example, alteration of the highly conserved Ser-672 residue in the HCN4 CBD by the heterozygous S672R mutation results in bradycardia, i.e. abnormally low resting heart rates that are reduced by ϳ30% (17,18). Electrophysiological profiling using integral HCN4 channels revealed that the S672R HCN4 variant is activated at voltages that are ϳ10% more negative than wild type in either the absence or the presence of endogenous cAMP (18). The negative S672R versus WT voltage shift was independently confirmed using chimeric constructs in which the HCN2 trans-membrane channel was combined with the HCN4 IR (20) , indicating that S672R represents a loss-of-function mutation with enhanced resistance to channel opening, i.e. S672R causes a reduced inward current flow in the diastolic phase, explaining the slower than normal heart rate observed in patients affected by the mutation (17,18). Furthermore, another S672R-included change consistently revealed by electrophysiology is a reduction on the deactivation time constant, i.e. the S672R mutation causes a significant acceleration of HCN channel deactivation (18,20).
As a first step toward elucidating the molecular mechanism for the observed electrophysiological phenotypes of the bradycardia S672R mutation, the structure of the S672R HCN4 IR mutant was solved in the presence of cAMP (Fig. 1C) (20). However, the crystal structure of the cAMP-bound S672R HCN4 IR tetramer did not reveal any major change relative to wild type (Fig. 1C, r.m.s.d. 0.41 Å), with the exception of the BBR (Fig. 1C), which exhibited increased local disorder for four loop residues, as indicated by missing electron density and side chain re-packing of selected surrounding residues ( Fig. 1D) (5,8,20). No appreciable variations in subunit assembly as a result of the mutation were detected (5,8,20). Although these initial structural results provide an essential framework for further investigations of the S672R variant, several questions about the molecular mechanism underlying the S672R bradycardia-inducing mutation still remain open.
It is still not clear how the S672R mutation perturbs the apo-CBD, and in general the CBD dynamics, which has been recognized as a key determinant of the auto-inhibitory CBD function in HCN gating and in other cAMP-dependent systems (9,13,(21)(22)(23)(24)(25)(26)(27)(28)(29). One of the aspects of CBD dynamics that is most relevant for auto-inhibition is the equilibrium between inactive and active conformations (Fig. 1, B and E) (30 -38). The coupling of this two-state dynamic conformational equilibrium with the cAMP-binding equilibrium results in a four-state thermodynamic cycle (Fig. 1E), which provides the simplest allosteric model to rationalize the cAMP-dependent modulation of HCN channels, i.e. a shift to less negative activation voltages, as well as to rationalize the effect of disease-related mutations. Considering that the S672R mutant leads to a shift to more negative activation voltages, i.e. an effect opposite to cAMP, while still preserving similar K1 ⁄ 2 values for integral HCN4 channels (17,18), we hypothesize that the S672R mutation compromises the cAMP-dependent activation either by giving rise to a partial reversal of the two-state auto-inhibitory equilibrium toward the inactive conformation (Fig. 1, E and F) and/or by partially stabilizing additional states, distinct from the fully active or fully inactive states, in which the activation function of the CBD is compromised. To test these hypotheses, here we comparatively analyze the S672R and WT HCN4 CBDs in both apo-and holo-forms using NMR spectroscopy, which reports on dynamics at residue resolution and over multiple time scales (30,39). Furthermore, we examined the S672R versus WT relaxation and chemical shift changes, which are exquisitely sensitive atomic reporters of subtle but functionally relevant allosteric perturbations that are often elusive to other structural techniques (27,30).

S672R mutation leads to pervasive perturbations of the HCN4 CBD in both Apo-and Holo-forms
As a first step toward understanding how the bradycardialinked S672R mutation affects the apo-and holo-HCN4 CBD, we assigned the 1 H- 15 N HSQC spectra of the S672R HCN4(563-724) construct (supplemental Fig. S2) in the absence of cAMP ( Fig. 2A) and in the presence of saturating amounts of cAMP (Fig. 2B). The comparison with the corresponding spectra of WT HCN4(563-724) (9) reveals that the S672R mutation results in pervasive changes for both apo-and holo-forms (Fig. 2, A-C), well beyond that previously anticipated based on X-ray structure comparisons (Fig. 1C) (20). In both apo-and holo-forms, the most pronounced chemical shift changes are observed not only at the site of the mutation but extend to loci located within a 15-Å radius from Ser-672 (Fig. 2, C and D), including the ␤2-3 loop, the PBC, the BBR, and, in the case of the holo sample, the C-terminal lid as well (Fig. 2C). In addition, several residues throughout the ␤-subdomain and the C-terminal helix in both apo-and holo-forms exhibits notable ppm changes (Ͼ 0.1 ppm; Fig. 2, C, E, and F).
Some of the largest S672R versus WT ppm variations arise either from losses of key interactions mediated by the Ser-672 hydroxyl, as in the case of Glu-625 and adjacent residues in the ␤2-3 loop (Fig. 1D), or from local structural changes, as in the case of the BBR that exhibits prominent S672R versus WT r.m.s.d. values based on the available crystal structures (Fig.  2D). However, most residues outside of the BBR display only relatively subtle S672R versus WT structural changes, as confirmed by marginal r.m.s.d. values (Fig. 2D) (8,20), and by the lack of significant S672R versus WT differences in the secondary structure profiles computed based on the secondary chemical shifts (supplemental Fig. S3). These observations suggest that the S672R versus WT ppm variations of Fig. 2C may also reflect mutation-dependent modulations in dynamics. In fact, chemical shifts of CBDs are known to report on the position of the dynamic auto-inhibitory equilibrium of CBDs (Fig. 1E), i.e. they encode for the relatively inactive versus active state populations. Specifically, NMR peak positions reflect populationweighted averages of the respective ppm values for the pure inactive and active CBD states, which exchange rapidly on the chemical shift NMR time scale (9,40). A simple but effective means to extract state populations from average chemical shifts is the projection analysis (9,21,28,29,30,41), in which changes in HSQC cross-peak positions arising from the mutation are evaluated relative to a reference vector defined by the WT apo and WT holo samples. For example, Fig. 3, A-C, shows how the residue-specific S672R holo versus WT holo vectors are analyzed in terms of the angles defined with the reference vectors ()and of the normalized projections onto the reference vectors (X). The cos() values (Fig. 3, A and B) report on the direction of the mutation-induced shifts in the inhibitory equilibrium (i.e. cos() Ͻ0 for inactivation) and on deviations from the linearity expected based on pure two-state equilibria (i.e. ͉cos()͉ Ͻ1), whereas the X values (Fig. 3C) quantify the degree of mutationinduced fractional inhibition or activation (9, 21, 28 -30, 41).

S672R mutation shifts the dynamic inhibitory equilibrium of the HCN4 CBD toward an inactive but partially De-correlated state
The chemical shift projection analysis (CHESPA) of holo-S672R shows that negative cos() values are observed for the A, topology of tetrameric HCN4 channels. The TMD is connected to the CBD through an ␣-helical domain known as the C-linker, which also serves as the main point of contact between protomers within the gating tetramer. Ser-672 and cAMP are shown as gray and black spheres, respectively. B, overlay of the apo (red; PDB code 2MNG) and cAMP-bound (gray; PDB code 3OTF) CBDs of WT HCN4. Select regions (i.e. N3A, BBR, and Lid) that undergo cAMP-dependent conformational changes are indicated with dashed boxes. cAMP and Ser-672 are depicted as in A. C, structural overlay of holo-WT (black, PDB code 3OTF) and holo-S672R (blue, PDB code 4HBN) HCN4 CBDs. D, changes in key contacts mediated by Ser-672 (black ribbon, yellow sticks) and S672R (blue ribbon, cyan sticks). Polar contacts for Ser-672 and S672R are depicted with green and red dashed lines, respectively. Selected residues that form interactions with S672(R) or that experience side-chain reorientations are shown as sticks. cAMP is shown as white spheres. E, four-state thermodynamic cycle to model allostery for the HCN CBD, in which the auto-inhibitory equilibrium is coupled to the cAMP-binding equilibrium through the active versus inactive state selectivity of cAMP. A hypothetical mechanism for the bradycardia-induced S672R mutant posits that the autoinhibitory CBD equilibrium is shifted by the mutation further toward the inactive state (red lines). F, simplified free energy diagram illustrating that the hypothesized mutation-induced shift of the auto-inhibitory equilibrium toward the inactive state may arise either from a stabilization of the inactive state and/or from a destabilization of the active state. This scheme assumes that the state-specific association constants are not significantly affected by the mutation. majority of the residues (Fig. 3A), with a cos() distribution highly skewed toward Ϫ1 (Fig. 3B), reflecting a clear prevalence of mutation-induced shifts in cross-peak positions from the WT holo toward the WT apo, i.e. the auto-inhibitory state. However, the extents of these shifts, i.e. the magnitudes of the X values, are highly residue-dependent (Fig. 3B). Residue-specific variations in X values are expected for sites proximal to either the mutation and/or cAMP, such as the BBR and the ␤2-3 loop (Fig. 3A), because these loci are affected by nearest neighbor effects (NNEs). Nevertheless, in our case residue-specific variations in X values are also observed for sites more removed from the mutation and cAMP, such as the N-terminal helices and the C-terminal region (Fig. 3C). For example, the X values of CTL residues 714 -724 range from Ϫ0.2 to ϳϪ1 with an average X value of Ϫ0.4 (Fig. 3, C and E), whereas in the N-terminal N3A significantly lower X values close to Ϫ0.1 are mea-sured (Fig. 3, C and F). The variations in residue-specific X values represent a clear deviation from the pattern of nearly residue-independent X values expected for these regions based on the two-state model (Fig. 1E), indicating that the response of the holo-CBD to the S672R mutation cannot be fully recapitulated by a simple two-state auto-inhibitory equilibrium (Fig.  1E). Deviations from the two-state model are also confirmed by the appearance of ͉cos()͉ values Ͻ1 for several residues, including those not affected by NNE (Fig. 3, A and D) (27), pointing to the presence of non-linear patterns (Fig. 3G). Overall, the CHESPA analysis of holo-S672R reveals that S672R causes a partial reversal of the inhibitory equilibrium toward a less activated state (Fig. 3B), but such reversion is not as correlated as expected based on a purely two-state model because it is significantly more pronounced in the C-terminal lid than in other regions largely insensitive to NNEs (Fig. 3, C, E, and F). Overlay of NH-HSQC spectra for WT (black) and S672R (blue) HCN4(563-724) in the apo (A) and cAMP-bound (B) samples. C, CCS differences between the WT and S672R variant for the apo (red) and cAMP-bound (blue) samples. The black circles along the bottom of the plot depict residues within 15 Å of the mutation site. The protein's secondary structure is depicted along the top as white (␣-helices) and gray (␤-strands) rectangles, and the BBR and several other motifs are highlighted with gray background. D, residue-specific local WT versus S672R r.m.s.d. values (vertical bars) computed in MOLMOL (56) by aligning three adjacent residues from the WT and S672R cAMP bound X-ray structures (PDB codes 3OTF (8) and 4HBN (20), respectively), as was described previously (27,41). The solid and dashed black horizontal lines depict the average and average ϩ S.D. r.m.s.d. values, and the horizontal red line represents a 15-Å distance from S672R. The blue and red lines represent the backbone amide nitrogen distances from residue S672(R) for the cAMP bound structure (4HBN) and the apo NMR ensemble (2MNG), respectively. The mutation site (S672R) is highlighted with a red asterisk in both C and D. E and F, residues with CCS changes greater than the dashed line in C were mapped onto the crystal structure as spheres and surface for both the apo (E) and the cAMP-bound (F) samples.
To investigate whether the mutation-induced de-correlation is unique to the holo sample or applies also to the apo-form, we extended the CHESPA analysis to the apo-S672R sample (Fig. 3, I-L). Although the S672R holo-CHESPA analysis clearly shows a shift toward inactivation, this shift is not expected to affect the auto-inhibitory S672R apo-CHESPA. The S672R-induced auto-inhibitory shift arises from a partial stabilization of the inactive state and/or from partial destabilization of the active . The red box highlights a subset of lid residues with markedly negative X values, pointing to a significant shift toward the inactive state. Selected residues reported in E-H are labeled. D, 3D map of CHESPA outliers, defined as residues with either ͉cos͉ Ͻ0.9 or ͉X͉ Ͼ1.1. Such outliers reveal deviations from the two-state equilibrium, such as those arising from nearest neighbor effects or population of a third state. E, representative NH-HSQC cross-peaks for a residue from the C-terminal lid region with a holo-S672R X value close to Ϫ0.4 (C, red box). F, similar to E, but for a representative residue with a holo-S672R X value close to Ϫ0.1. G, representative residue (Val-642) from the S672R holo spectra that exhibited a non-linear chemical shift change relative to the WT apo-and holo-states. H, representative residues from the S672R apo spectra compared with WT apo-and holo-states that displayed non-linear chemical shift changes. This non-linearity is also evident from the CHESPA analysis of the apo-S672R HCN4 CBD, shown in I and J, revealing a nearly uniform cos() distribution (I) between the Ϫ1 and 1 extremes. Selected residues reported in E-H are labeled in I. K, X values for apo-S672R residues. L, similar to D but for the apo-S672R CHESPA outliers. state (Fig. 1F). However, in either case, the WT apo inhibitory equilibria of our construct are known to be already almost quantitatively skewed toward the inactive state ( Fig. 1E) (9,40). Hence, the Ser-672-induced auto-inhibitory shift is expected to be negligible under our experimental conditions, suggesting that the S672R apo-CHESPA analysis should reflect primarily de-correlative effects of the mutation, i.e. deviations from the two-state model. This prediction is confirmed by the cos() values measured for S672R apo (Fig. 3I), exhibiting a broad scattered distribution spanning the full Ϫ1 to 1 range (Fig. 3J), which is in marked contrast to the inactivation-skewed pattern observed for S672R holo (Fig. 3B). Furthermore, the abundance of residues with ͉cos()͉ Ͻ1 (Fig. 3, H and J) and the scatter detected in the residue-specific X values (Fig. 3K) confirm that the perturbations caused by S672R cannot be modeled exclusively based on a two-state model (Fig. 1E). These deviations from the twostate transition point to mutation-induced de-correlations and possibly to enhanced dynamics. The latter were gauged through S672R versus WT comparative 15

Comparative S672R versus WT analysis of internal dynamics for the Apo-CBD reveals that the S672R mutation results in a significant overall enhancement in the amplitude of the picosecond to nanosecond internal motions with local increases in millisecond to microsecond dynamics near cAMP phosphate recognition sites
The overall S 2 distributions observed in the absence of cAMP for the WT and S672R CBD reveal that the mutation results in a marked net shift toward lower S 2 values (Fig. 4A). This result points to increased amplitudes of picosecond to nanosecond motions in S672R versus WT apo-CBD. Furthermore, considering that order parameters for picosecond to nanosecond motions are a proxy for conformational entropy (30), the S672R-induced shift toward reduced S 2 values suggests a widening of the well in the free-energy landscape minimum representing the apo-CBD (Fig. 5A). The S672R versus WT dynamics enhancements of the apo-CBD are not limited to the picosecond to nanosecond time scale, but they include also microsecond to millisecond motions (Fig. 4K). In particular, the ␤2-␤3 loop proximal to the mutated site ( Fig. 1D) displays some of the greatest residue-specific enhancements of dynamics, as observed for both low and high frequency spectral density functions (supplemental Fig. S4; Fig. 4, C and K). Specifically, residue Glu-625, which forms a hydrogen bond with Ser-672 ( Fig.  1D), exhibits a remarkable ϳ3-fold increase in J(0) (supplemen-tal Fig. S4A; Fig. 4C), consistent with a large increase in microsecond to millisecond dynamics. In addition, several residues adjacent to Glu-625 display enhanced picosecond to nanosecond dynamics ( Fig. 4K; supplemental Fig. S4C). These data suggest that by disrupting the interaction between Ser-672 and Glu-625, the S672R mutation destabilizes the entire ␤2-␤3 loop. Because the ␤2-␤3 loop lies adjacent to the PBC, we checked whether this destabilization could extend to this cAMP-binding motif. Although the dynamic characterization of the apo-PBC is sparse due to line broadening, several residues at the edges of the PBC do indeed exhibit enhanced microsecond to millisecond dynamics ( Fig. 4G; supplemental Fig.  S4A), confirming that the ␤2-␤3 enhancement of dynamics does at least partially extend to the PBC, which is a primary cAMP recognition site. So next we checked whether the S672R causes significant changes in dynamics at the other main cAMP contact sites, i.e. the BBR and the lid.
As to the BBR, whereas the 645-650 region in the WT apo-CBD displays higher J( N ϩ H ) values than adjacent residues, pointing to flexibility in the picosecond to nanosecond time scale, no further significant enhancements of dynamics could be detected for S672R apo under our experimental conditions (supplemental Fig. S4C; Fig. 4E). Although BBR residues Lys-645 and Lys-648 did exhibit slight decreases in J(0) ( Fig. 4E; supplemental Fig. S4A), Lys-648 displayed quenched dynamics at the N frequency (supplemental Fig. S4B), and for either residue no significant changes at the J( Nϩ H ) spectral density were detected. In contrast, the apo C-helix exhibited quenching of microsecond to millisecond dynamics as seen by decreased J(0) values for Thr-705 and Ala-707 ( Fig. 4L; supplemental Figs. S4A and S6A). Meanwhile, our probes for the BBR and lid do not support a major net S672R versus WT enhancement in picosecond to nanosecond dynamics in these regions (Fig. 4, E and I). However, it should be considered that these regions are already flexible in the picosecond to nanosecond time scale in the WT apo-CBD, as shown by a significantly increased J( N ϩ H ) value relative to those expected for overall tumbling (supplemental Fig. S4C). Beyond the cAMP-binding sites, significant dynamic changes were found at the interface between ␣and ␤-subdomains and even in parts of the N-terminal C-linker (Fig. 4, K and L; supplemental Fig. S4). The net effect of these pervasive variations in dynamics in the apo-CBD is that the S672R mutation results in an enhancement in the amplitude of fast (picosecond to nanosecond) motions (Fig. 4A) and in the associated entropy of the conformational ensemble sampled by the apo-CBD (Fig. 5A). To evaluate the effect of the S672R mutation on the holo-CBD dynamics, the comparative S672R versus WT 15 N-relaxation analyses were extended to the data acquired in the presence of excess cAMP (Fig. 4, B, D-J, M, and N; supplemental Figs. S5 and S6).

In the Holo-CBD, the S672R mutation results in localized enhancements in the dynamics of key cAMP recognition and gating regions (PBC and Lid)
Unlike the apo-CBD (Fig. 4A), in the case of the holo-CBD no net S672R versus WT shift to larger amplitude picosecond to nanosecond motions is observed in the overall distribution of S 2 order parameters (Fig. 4B). However, several localized S672R versus WT enhancements of dynamics were detected also for the holo-CBD (Fig. 4, H, J, M, and N; supplemental Figs. S5 and S6). The most dramatic increase in dynamics was found in two loci within the cAMP-binding pocket, i.e. PBC and lid, but significant enhancements of dynamics were also identified in selected sites within the N3A and C-linker (Fig. 4, H, J and M; supplemental Fig. S5). The PBC, which is nearly completely assigned in the cAMP-bound state and better sampled than in the apo-state, experiences large enhancements of microsecond to millisecond dynamics, as indicated by several significant increases in the J(0) (Fig. 4H; supplemental Fig. S5A). Unlike the apo-CBD, the mutation induced increase in dynamics for the PBC did not extent to the adjacent ␤2-␤3 loop (Fig. 4D). Similarly to the apo-CBD, for the BBR, i.e. the other cAMP contact region in the ␤-subdomain, only a few significant changes in dynamics were recorded, and these exhibited a mixed pattern, including both enhancements and quenches ( Fig. 4F; supplemental Fig. S5).
Besides the PBC, the other region of the holo-CBD experiencing major dynamic enhancements is the C-helix and the adjacent C-terminal lid segment ( Interestingly, for the S672R holo sample the picosecond to nanosecond dynamics in the lid region, as probed by the average J( Nϩ H ) spectral densities for the C-helix and CTL, lie in-between those of the WT apo and WT holo samples (supplemental Tables S2 and S3). This result is in agreement with the fractional activations observed by CHESPA for this region, which also shows a partial reversion toward the apo-inactive state (Fig. 3C). Together, the 15 N-relaxation data and the CHESPA comparative analyses for the holo-CBD consistently point to a partial mutation-induced disengagement of the CTL, which is part of a key gating element for cAMP. Furthermore, the mutation also enhances the dynamics of key cAMP recognition elements, such as the PBC (Fig. 4, G and H), suggesting that S672R may lead also to altered on-off cAMP exchange kinetics. This hypothesis is also supported by the overall S672R versus WT net shift toward lower order parameters, i.e. higher conformational entropy, observed for the apo-CBD (Fig. 4A).
The wider free energy well of the mutant is in agreement with reduced free energy barriers for the apo-holo exchange (Fig.  5A), corroborating the hypothesized S672R versus WT acceleration of cAMP-binding kinetics. To further test this hypothesis, we analyzed the cAMP-binding dynamics utilizing proteinbased and ligand-based NMR experiments.

S672R mutation accelerates the on/off rates of the HCN4 CBD-cAMP complex
The cAMP-binding kinetics were probed by comparative analyses of chemical exchange regimes in HSQC spectra (Fig. 5, B-F) and ROESY buildup curves (Fig. 5, G and H). The HSQC analysis involved monitoring residue-specific transitions from slow to intermediate and from intermediate to fast for apo-holo exchange regimes in HSQC spectra acquired at sub-stoichiometric ligand/protein ratios, as illustrated in Fig. 5B. At a given protein/ligand stoichiometric ratio, the effective apo-holo exchange constant, k ex , is defined by Equation 1, A particular value of k ex results in different exchange regimes for different residues, depending on the respective differences in Hz between the apo and holo resonance frequencies (⌬ A-H ). If k ex Ͻ Ͻ ⌬ A-H for a selected residue, then that amino acid will appear in the slow exchange regime, which is characterized by two peaks appearing at the positions of the apo-and holo-states with intensities proportional to the populations of ligandbound and -unbound protein (Fig. 5B, blue peaks). However, under the same experimental conditions, and therefore with the same k ex , for other residues that experience lower ⌬ A-H values, it is possible that k ex Ϸ ⌬ A-H , i.e. the residues fall within the intermediate exchange regime in which peaks broaden significantly, often beyond detection (Fig. 5B, black peak). For amino acids with lower ⌬ A-H so that k ex Ͼ Ͼ ⌬ A-H , the fast exchange regime applies, and only a single peak is observed whose chemical shift is a population weighted average of the apo and holo chemical shifts (Fig. 5B, red peak). The difference in residue-specific exchange regimes provides an opportunity for a simple estimation of k ex by assessing the range of ⌬ A-H values that result in intermediate chemical exchange. For this purpose, the ⌬ A-H values were sorted in decreasing order, and the chemical exchange regime was determined for each residue to define the boundaries between slow, intermediate, and fast exchange (Fig. 5, D and F). The k ex was then assessed based on the range of ⌬ A-H values that resulted in intermediate chemical exchange (Fig. 5, D and F). Overall, Fig. 5, D and F, illustrate how the ranked ⌬ A-H frequency differences define a "built-in" scale for gauging k ex . This approach is not as quantitative as  Table S1. C-J, two-dimensional reduced spectral density plots for the apo-WT (black) and S672R (red) HCN4 CBD. The following four regions were probed: the ␤2-␤3 loop (aa 621-636), BBR (aa 639 -652), extended PBC (aa 659 -676), and Lid region (aa 699 -724), for both apo (top row) and cAMP-bound (bottom row) samples. The Lid is further dissected into the Lid and CTL in supplemental Fig. S6. The solid black line represents the expected spectral densities for an isotropic rigid rotor, and the black and red ellipsoids are centered on the average spectral density values with axes defined by Ϯ2 S.D. for the WT and mutant, respectively. K, map of apo residues with S672R versus WT enhanced picosecond to nanosecond (red spheres and surface) and microsecond to millisecond (orange spheres and surface) dynamics. L, map of apo residues with S672R versus WT quenched picosecond to nanosecond (blue sphere and surface) and microsecond to millisecond (cyan sphere and surface) dynamics. M and N, as in K and L except for the cAMP-bound sample.
other methods (e.g. N z -exchange spectroscopy (26,42)), but it is applicable to both HCN4 WT and mutant, unlike N z -exchange experiments, as explained in the supplemental Materials and Methods.
The HSQC-based k ex assessment was performed for both the WT and S672R HCN4 CBD in the presence of sub-stoichiometric amounts of cAMP (Fig. 5, C-F). Evidence that k ex increases in going from S672R to WT is already provided by residues that The diagonal peaks indicated in the spectra correspond to the H8 proton of cAMP in its bound (H8b) and free (H8f) states. The cross-peak for the chemical exchange from bound to free H8 is also indicated as H8b 3 f. Black contours indicate positive peaks; red contours indicate negative peaks. The H8b 3 f preserves the same sign as the diagonal, as expected for chemical exchange cross-peaks. H, plot of H8b 3 f versus H8b intensity ratios versus ROESY mixing times. The initial slopes of the buildups based on the first order approximation are reported in the figure. For the S672R mutant, the bound ligand diagonal peak was detectable only up to a ROESY mixing time of 30 ms due to rapid decay. In contrast, the WT bound H8 signal was detectable up to a ROESY mixing time of 50 ms. The solid black line depicts the slope of the two WT points shown in the plot, and the dashed line depicts the slope that includes two additional points at mixing times of 40 and 50 ms. experience similar apo versus holo frequency differences but fall in the slow exchange regime in WT and in the fast exchange regime in S672R. This is the case, for example, for Ala-597 as shown in Fig. 5, C and E. In WT HCN4, two separate peaks are detected for Ala-597 (Fig. 5C), consistent with slow apo-holo exchange, although in S672R a single peak is observed (Fig. 5E), consistent with fast apo-holo exchange. Considering that the ⌬ A-H values of WT and S672R are comparable (Fig. 5, C and E), the slow versus fast exchange regimes observed for Ala-597 in WT versus S672R, suggests that the mutation results in an increased k ex value. The two-peak versus one-peak patterns shown in Fig. 5, C and E, are the basis for determining whether a given residue experiences slow or fast apo-holo exchange (see supplemental Materials and Methods and supplemental Fig. S8 for further details) so that the plots of the sorted ⌬ A-H values are partitioned into the three different kinetic regimes (Fig. 5, D and F). The comparison of Fig. 5, D and F, reveals marked WT versus S672R differences in the positions of the apo-holo exchange regimes within the scales defined by the sorted ⌬ A-H values (Fig. 5, D and F). For WT, the slow to fast transition occurs in the 30 -40 Hz range (Fig. 5D), and for S672R, the intermediate regimes fall in the 150 -260 Hz window (Fig. 5F), confirming that the mutation leads to a significantly enhanced k ex value.
The S672R-induced k ex enhancement revealed by the HSQC analyses (Fig. 5, B-F) reflects a primarily kinetic effect, as no significant S672R versus WT differences were observed for the thermodynamic affinities of the HCN4 CBD, with K d values in the 2-10 M range for both mutant and WT (Fig. 6). Specifically, based on supplemental Equation S3, the assessed k off values for WT and S672R are expected to fall in the 23-31 and 113-202 s Ϫ1 ranges, respectively, pointing to a major S672Rinduced acceleration in the cAMP on/off dynamics. This result was independently confirmed through ligand-based NMR experiments, i.e. ROESY spectra (Fig. 5, G and H), which complement the protein-based NMR data, i.e. the semi-quantitative HSQC analysis outlined above. The ROESY data were acquired in the presence of excess ligand ,and slow exchange was observed between the free and bound ligand peaks, as confirmed by ROESY cross-peaks with the same sign as the diagonal (Fig. 5G). The initial slopes of the mixing time buildup measured for the bound-to-free exchange cross-peak fall in the 18 -32 and 110 -134 s Ϫ1 ranges for WT and S672R, respectively (Fig. 5H). These values are consistent with the k off ranges independently estimated based on the HSQC data and confirm that the S672R mutation leads to an at least ϳ3-fold enhancement in the on-off exchange rates of cAMP.

Bradycardia-related S672R mutation remodels the free energy landscape and the dynamic profiles critical for the cAMPdependent modulation of HCN4 channels
The comparative NMR analyses presented here (Figs. 2-6; supplemental Figs. S3-S8) have revealed how the S672R bradycardia-inducing mutation remodels the free energy landscape (FEL) for the cAMP-dependent regulation of the HCN4 channels, as illustrated in Fig. 7A. The FEL remodeling caused by S672R (Fig. 7A) is recapitulated by two main perturbations. First, the mutation results in a shift of the auto-inhibitory CBD equilibrium toward inactivation, as shown by the CHESPA data (Fig. 3), consistent with the stabilization of inactive conformations and/or the destabilization of the active state (Fig. 7A, vertical red arrows). Second, the S672R variant leads to a widening of the apo-inactive free energy basin (Fig. 7A, horizontal red  arrows), as indicated by a redistribution of the apo-CBD order parameters for fast (picosecond to nanosecond) motions toward lower values (Fig. 4A). The broadening of the apo-inactive conformational sub-set is also confirmed by the CHESPA data (Fig. 3), which indicate that the mutation-induced shift to  HCN4(563-724) with cAMP monitored by STD NMR. Data were fitted with a rectangular hyperbola for a 1:1 binding stoichiometry providing a K d in the 2-10 M range. The fractional saturation was measured by normalizing the STD amplification factor (9). The corresponding data for WT human HCN4(563-724) (9) is shown (black dots) for the convenience of comparison. The low micromolar affinity reported here for the S672R variant is in agreement with the range of K d values of 0.8 and 11.2 M measured by fluorescence anisotropy and isothermal titration calorimetry by Xu et al. (20), respectively. In addition, considering that the apo-WT state is already primarily sampling the inactive state, based on the scheme in Fig. 1F, a mutation-induced shift toward inactivation is not expected to result in major variations in the observed affinity for cAMP. Figure 7. Proposed mechanism of S672R-induced bradycardia based on free energy landscape re-modeling. A, schematic free energy landscape for the allosteric cAMP-dependent modulation of the HCN4 CBD and its perturbation by the S672R mutation (red arrows). The horizontal red arrows denote the S672R versus WT enhanced conformational entropy of the inactive state as indicated by the S672R versus WT decreased S 2 parameters of the apo samples (Fig. 4A). The vertical red arrows illustrate that the S672R mutation results in a CBD shift toward inactivation either by destabilizing the active state and/or by stabilizing the inactive state. B, summary of critical S672R versus WT changes in dynamics. Only the ␤-core and the C-terminal ␣-subdomain (B and C helices and the adjacent CTL) are shown for simplicity. Red wavy lines indicate S672R versus WT enhanced dynamics in the HCN4 CBD. According to the proposed model, the enhanced dynamics at sites lining the cAMP-binding pocket as well as the disengagement of the CTL in the cAMPbound sample lead to an acceleration of the kinetics of cAMP binding and unbinding. The red dot indicates the site of the S672R mutation. The red outline and red shading denote areas with enhanced S672R versus WT dynamics, leading to an increased k off rate for cAMP (red arrow).
inactive states cannot be fully captured by a simple two-state exchange model, i.e. rather than a single cooperative transition, S672R causes a partial de-correlation, which accounts for the observed residue dependence of the fractional activation values (X) and for the ͉cos()͉ Ͻ1 value observed at multiple sites (Fig.  3). Considering that the apo-CBD samples primarily inactive conformations, the broadening of the free energy well caused by S672R (Fig. 7A) is expected to lower the free energy barrier for apo-holo exchange (Fig. 5A). This prediction was confirmed by the observation that S672R accelerates the cAMP-binding kinetics by a factor of at least 3-fold, as measured by both protein-and ligand-based NMR methods (Fig. 5, C-H).
The mutation-induced enhancement of dynamics is not limited only to the fast (picosecond to nanosecond) time scale, but it extends also to the microsecond to millisecond dynamic window (Fig. 4, C-N). The S672R variant exhibits enhanced millisecond to microsecond dynamics at sites critical for cAMP recognition and binding, such as the PBC (Fig. 7B). The PBC enhancement of millisecond to microsecond dynamics was detected in both apo-and holo-forms and extends to the adjacent ␤2-3 loop in the apo-CBD, whereas for the holo-CBD it affects the C-terminal lid, referred to as CTL (Fig. 7B). The CTL is part of the lid region critical for the cAMP-dependent HCN channel activation (14) and for gating cAMP binding in WT (Fig. 1, B and C). However, in the holo-S672R, the CTL becomes partially disengaged, as consistently indicated by both the negative CHESPA fractional activation values (Fig. 3, B, C, and E) and the concurrent mutation-induced enhancement of both fast and slow dynamics ( Fig. 4; supplemental Fig. S6D and Table  S3). The enhanced lid dynamics in S672R versus WT is also in full agreement with the faster binding kinetics observed for the mutated variant. Overall, the picture emerging from the NMR investigation presented here suggests that the S672R-induced perturbations extend well beyond what was previously anticipated based on X-ray crystallography of the cAMP-bound intracellular region (20). The crystal structure of S672R and WT holo-HCN4 IR tetramers differs primarily at the BBR, whereas NMR chemical shift changes show that for both apoand holo-forms S672R affects multiple sites within an ϳ15-Å radius, spanning not only the BBR but also several other loci in the ␤-subdomain and lid regions (Fig. 2, D-F).

S672R-induced FEL remodeling provides an initial molecular explanation for the electrophysiological changes underlying familial bradycardia
The S672R-induced remodeling of the free energy landscape for the cAMP-dependent modulation of HCN (Fig. 7) provides a molecular framework to rationalize the two main electrophysiological signatures of this bradycardia-related mutation, i.e. the negative shift in the activation voltage and the faster deactivation. The former phenotype may arise as a result of a reduction in cAMP affinity (20) and/or of a constitutive S672Rinduced variation in the position of the auto-inhibitory CBD equilibrium toward inactivation (17,18). Our data support the notion that the second mechanism significantly contributes to the negative S672R-induced shift in the activation voltage observed by electrophysiology. Considering that activation voltages of full-length HCN4 have been shown to correlate with the position of the auto-inhibitory CBD equilibrium as probed by CHESPA (9), the partial inactivation observed for the S672R variant (vertical red arrows in Fig. 7A) is consistent with the negative S672R versus WT change of the activation voltage observed by electrophysiology for the full-length channel (18,20). Furthermore, the position of the auto-inhibitory CBD equilibrium in reconstituted apo full-length channels may be affected by the presence of dimers known to be populated in the absence of cAMP (43,44). If dimerization causes a partial shift of the WT apo-CBD toward activation, the S672R mutation may result in an appreciably increased sampling of inactive states also in the absence of cAMP, explaining the constitutive negative shifts in activation voltages observed for apo-channels (18).
The other putative mechanism, i.e. a reduction of cAMP affinity, does not appear to be sufficient alone to explain the observed negative S672R-induced shift in the activation voltage without any constitutive mutation-induced partial inactivation. This conclusion is supported by several independent observations. First, S672R is known to cause a negative shift in activation voltages even in the absence of cAMP (18). Second, we did not observe any major S672R versus WT increases in the cAMP K d values measured under our experimental conditions by STD for the monomeric HCN4 CBD, as both K d values fall in the 2-10 M range (Fig. 6). Third, no significant S672R-induced enhancements in the EC 50 values were reported for the fulllength HCN4 channels either (18). Fourth, although S672R versus WT K1 ⁄ 2 increase was reported for either HCN2 or hybrid HCN2/4 channels (20), a significant mutation-induced reduction in integral HCN4 affinity would decrease the sensitivity to physiological levels of cAMP and result in chronotropic incompetence, which, however, was not observed in patients affected by the S672R mutation (17,18). These observations consistently suggest that the negative S672R-induced shift in the activation voltage cannot be fully accounted for by a reduction of cAMP affinity alone without any constitutive S672Rinduced variation in the position of the auto-inhibitory CBD equilibrium.
The other electrophysiological signature of the bradycardiarelated mutation, i.e. the faster deactivation, is explained at least in part by the S672R-induced acceleration of the cAMP off kinetics, as seen by our combined analyses of protein-and ligand-based NMR experiments (Fig. 5). This result is in agreement with the S594R versus WT k off enhancement reported for the unbinding of a fluorescent cAMP analog from full-length mHCN2 channels and is rationalized by both the S672R-induced widening of the apo free energy well (Fig. 5A) and the enhancement of dynamics observed in the holo-CBD for the PBC and the lid, two critical elements of the cAMP-binding sites (Fig. 7). Unlike the PBC and the lid, the remaining cAMP recognition region, i.e. the BBR, does not appear to be subject to significant mutation-induced enhancements of dynamics (Fig. 4).
In conclusion, we have shown how the bradycardia-related S672R mutation remodels the FEL for the cAMP-dependent modulation of HCN4 channel gating. The S672R-induced FEL remodeling is recapitulated by two main effects. First, S672R results in the constitutive stabilization of inactive versus active states, contributing to the negative S672R versus WT shift in the activation voltage observed by electrophysiology for HCN4 channels. Second, S672R leads to enhanced dynamics at sites pivotal for cAMP binding, such as the PBC and the lid, which in turn leads to an acceleration of the cAMP-unbinding kinetics, explaining at least in part the faster deactivation reported for S672R versus WT HCN4 channels. Hence, our results provide a mechanistic explanation of the main electrophysiological determinants of the familial bradycardia diagnosed in patients carrying the S672R mutation. In addition, the NMR approaches outlined here will be useful to map how other disease-related mutations remodel functionally relevant free energy landscapes in other systems as well (45)(46)(47)(48).

HCN4 sample preparation
The expression and purification of HCN4(563-724), referred to here as HCN4 CBD for simplicity, was adapted from a previously established protocol (9). A rare codon optimized SUMO-HCN4 fusion construct was subcloned into a pET302NT-His vector (Invitrogen) and expressed in Escherichia coli BL-21(DE3) cells in M9 minimal media isotopically enriched with [ 15 N]ammonium chloride (Sigma). Expression was induced with 0.5 mM isopropyl ␤-D-1-thiogalactopyranoside (Thermo Fisher Scientific) at an A 600 of 0.6 and incubated for a further 16 -18 h at 20°C. Cells were lysed with a cell disruptor, and the cell lysate was centrifuged at 20,000 ϫ g for 1 h. The cell lysate was then filtered with 0.8/0.2-m syringe filters (VWR Scientific) and purified by Ni-NTA affinity chromatography following the same protocol as in Akimoto et al. (9). Following proteolytic cleavage of the His 6 -SUMO tag with His 6 -tagged TEV protease, the sample was purified again with Ni-NTA affinity chromatography to remove both the SUMO tag and the TEV protease. The cleaved protein was then diluted 3-fold with Sp Wash buffer (10 mM MES, 10 mM KCl, and 10 mM 2-mercaptoethanol, pH 6.0) and further purified by cation exchange chromatography with a 5-ml HiTrap Sp column (GE Healthcare). Following injection, the protein was eluted using a 100% gradient over 100 min with Sp elution buffer (10 mM MES, 1 M KCl, and 10 mM 2-mercaptoethanol, pH 6.0). The purified protein was then dialyzed overnight in NMR buffer (20 mM MES, 75 mM KCl, 2 mM EDTA, 2 mM EGTA, 1 mM DTT, and 0.02% sodium azide, pH 6.5) and stored at 4°C. HCN4(563-724) S672R was created by site-directed mutagenesis using a variation of the QuikChange protocol. The reaction used the KOD Hot Start master mix (Novagen) with primers synthesized in-house. The reaction protocol followed the manufacturer's specifications for the KOD Hot Start kit. Following the reaction, parental DNA was degraded with DpnI for 1 h at 37°C, and the reaction product was transformed in E. coli top10 competent cells. Plasmids were later isolated from single colonies grown in liquid LB media using a GeneJET Plasmid Miniprep kit (Thermo Fisher Scientific) and sequenced to confirm the mutation. HCN4(563-724) S672R was expressed and purified similarly to the WT with the exception that the cation exchange gradient was changed to 100% Sp elution buffer over 25 min.

NMR spectroscopy
NMR experiments were acquired on a Bruker Avance 700 MHz spectrometer equipped with a 5-mm TCI cryo-probe, unless otherwise specified. All experiments were acquired at 300 K in the NMR buffer defined above with 5% D 2 O unless otherwise specified. Spectra were processed either with NMRpipe (49) or directly in Bruker's Topspin and analyzed in NMRFAM-SPARKY (50). Chemical shift assignments of S672R were obtained through spectral comparison with WT spectra (9), if no ambiguities were present, and through triple resonance experiments (i.e. HNCO, HNCA, and HN(CO)CACB) otherwise. The Talosϩ server was used to predict the protein's secondary structures using HN, C␣, C␤, and CЈ chemical shifts and to cross-check assignments based on acceptable backbone torsional angles (51).

Chemical shift analyses
Compounded chemical shift differences (⌬CCS) were computed as described previously (9,27,29,40,41), using Equation 2, where the ␦ symbols represent the proton (H) and nitrogen (N) chemical shifts (in ppm) for states 1 and 2. The CHESPA protocol has been described elsewhere (9,23,29). Briefly, the WT apo and WT holo cross-peaks for a given residue define a reference vector A, whereas the S672R (either apo or holo) and its corresponding WT state define a perturbation vector B. The cos and fractional activation (X) values are then computed for residues with mutant versus WT ⌬CCS Ͼ0.05 ppm, using Equations 3 and 4.

NMR spin relaxation measurements and reduced spectral density mapping
The longitudinal and transverse relaxation rates, as well as the { 1 H-15 N} steady state NOEs were measured for 150 M WT and S672R HCN4(563-724) samples in their apo-form and cAMP-bound form (2 mM cAMP). Pseudo-3D versions of the inversion recovery and CPMG pulse sequences modified to include a water flip-back pulse and sensitivity enhancement were used to measure T 1 and T 2 rates, respectively (52,53). Both experiments included 1024 complex points with a spectral width of 14.28 ppm in the t 2 dimension and 256 complex points with a spectral width of 31.82 ppm in the t 1 dimension. The T 1 experiment was acquired with 16 scans, 128 dummy scans, a recycle delay of 2.5 s, and the following relaxation delays: 100 ms (ϫ2); 200, 300, and 400 ms (ϫ2); and 500, 600, 800, and 1000 ms (ϫ2) (where ϫ2 indicates duplicate spectra). The T 2 experiment was acquired with 32 scans, 128 dummy scans, a recycle delay of 1.8 s, and the following CPMG relaxation delays: 8 (9), and for the S672R mutant the HN NOE were measured similarly to WT. In brief, the HN-NOE experiment used a 10-s recycle delay with a 5-s proton saturation period. NOE experiments were acquired in 10 -12 sets of saturated and unsaturated spectra, which were co-added and processed similarly to the HSQC spectra noted above. The HN-NOEs were then computed from the intensity ratios of saturated/non-saturated peaks.
Error bars were determined from replicate measurements. T 1 and T 2 experiments were processed in Topspin using cosine squared window functions and residual water deconvolution. Spectra were then transferred to Sparky and fitted to an exponential decay to extract T 1 and T 2 rates. Reduced spectral densities were computed as described previously (24,25). Significant changes in dynamics were identified by computing the average spectral density values for all residues along with their standard deviations (S.D.) and by excluding from subsequent averaging calculations any residues outside the average Ϯ 1 S.D. This computation of progressively trimmed averages and standard deviations was completed three times, after which any further iteration resulted only in marginal changes, i.e. ϳ1%. At this point, the average spectral density was considered to arise primarily from overall molecular tumbling as opposed to internal motions.
The assessment of spectral densities reflecting overall tumbling served as a reference to identify significant S672R versus WT variations in internal dynamics based on the following criteria. (i) The spectral density of a residue for either WT or S672R has to be outside the average Ϯ 5 S.D. (ii) The S672R versus WT spectral density change for a given residue must exceed five times the S.D. obtained from the overall tumbling assessment or, if greater, the sum of the respective experimental errors. Residues that met such criteria were deemed subject to significant S672R versus WT changes in internal dynamics and mapped onto the crystal structure of S672R HCN4. Residues undergoing S672R versus WT enhancement (suppression) in the picosecond to nanosecond backbone dynamics were identified through increases (decreases) in J( Hϩ N ), whereas residues with S672R versus WT enhancement (suppression) in millisecond to microsecond backbone dynamics were identified through increases (decreases) in J(0), typically with values above those computed for overall tumbling, in the absence of J( Hϩ N ) changes, as well as through 850 versus 700 MHz R 2 comparisons. The T 2 relaxation rates were also measured on a Bruker Avance 850 MHz spectrometer equipped with a TXI probe to further examine residues with significant contributions from R ex arising from dynamics in the millisecond to microsecond time scale. The pulse sequence employed for this analysis was a pseudo-3D version of the CPMG experiment with sensitivity enhancement and a water flip-back pulse. Spectra were acquired with 2048 complex points and a spectral width of 14 ppm in the direct dimension and 256 complex points and a spectral width of 36 ppm in the indirect dimension. Six CPMG relaxation delays were acquired in duplicate: 16 Fig. S7. The difference in R 2 rates at 850 and 700 MHz was then com-puted, and significant changes in microsecond to millisecond dynamics were identified using an approach similar to the one described above for the reduced spectral density analysis. The S 2 order parameters were calculated from J(0) and J( N ϩ H ) as described previously (54,55). The S 2 distributions were modeled through non-linear fitting with both one and two Gaussian curves. Further details are available in supplemental Table S1.

Measurements of binding thermodynamics and kinetics
Binding affinities of cAMP for the WT and S672R HCN4 CBD were determined by measuring binding isotherms through the STD amplification factor, as explained previously (9). For the assessment of the binding kinetics through proteinbased experiments, N z experiments could not be implemented due to weak holo peaks at cAMP concentrations sufficiently low to ensure slow exchange (supplemental Materials and Methods). Hence, a simpler and less quantitative approach was taken based on the analysis of gradient-and sensitivityenhanced { 15 N, 1 H} HSQC spectra, as explained under "Results." These HSQC datasets were acquired with 1024 complex points and a spectral width of 14.1 ppm in the directly detected 1 H dimension and 128 complex points and a spectral width of 31.8 ppm in the indirectly detected 15 N dimension. Spectra were processed employing a 60°shifted sine squared bell window function, zero filling, and linear prediction in the indirect dimension. Titrations were acquired with 100 M total protein concentrations and total cAMP concentrations of 25, 50, 100, 200, 400, 600, 1000, and 1500 M. The HSQC-based assessment of the cAMP-binding kinetics was complemented by acquiring 2D-ROESY spectra, recorded with mixing times in the 15-50 ms range in the presence of excess ligand. Spectra were acquired with spectral widths of 12 ppm in both dimensions and 2048 and 512 complex t 2 and t 1 points, respectively. The strength of the continuous wave spin-lock strength was set at 2.5 kHz. The cAMP off kinetics was probed through the initial slope of the plot of the ratio of the bound-to-free (b 3 f) ROE cross-peak intensity to the bound ligand diagonal peak intensity versus the ROESY mixing times.