Conformational free-energy landscapes of a Na+/Ca2+ exchanger explain its alternating-access mechanism and functional specificity

Secondary-active transporters catalyze the movement of myriad substances across all cellular membranes, typically against opposing concentration gradients, and without consuming any ATP. To do so, these proteins employ an intriguing structural mechanism evolved to be activated only upon recognition or release of the transported species. We examine this self-regulated mechanism using a homolog of the cardiac Na+/Ca2+ exchanger as a model system. Using advanced computer simulations, we map out the complete functional cycle of this transporter, including unknown conformations that we validate against existing experimental data. Calculated free-energy landscapes reveal why this transporter functions as an antiporter rather than a symporter, why it specifically exchanges Na+ and Ca2+, and why the stoichiometry of this exchange is exactly 3:1. We also rationalize why the protein does not exchange H+ for either Ca2+ or Na+, despite being able to bind H+ and its high similarity with H+/Ca2+ exchangers. Interestingly, the nature of this transporter is not explained by its primary structural states, known as inward- and outward-open conformations; instead, the defining factor is the feasibility of conformational intermediates between those states, wherein access pathways leading to the substrate binding sites become simultaneously occluded from both sides of the membrane. This analysis offers a physically-coherent, broadly transferable route to understand the emergence of function from structure among secondary-active membrane transporters.


Supporting Figures
).The selections include Cα atoms (green spheres), backbone N and O atoms (blue and red spheres), and sidechain carbon atoms (green sticks).Note these selections are approximately symmetric, i.e. they include equivalent residues in both topological repeats.(C,D).The residues that exhibit changes in water accessibility and lipid contacts are shown as grey and brown volumes, respectively, with water and lipid molecules in their vicinity shown with density isosurfaces (blue and yellow mesh, respectively).Residues that exhibit changes in protein-protein H-bonding and salt bridges are shown as sticks.The number of lipid-tail and H-bond contacts between two groups of atoms in a given simulation snapshot was calculated using distance cut-off values of 5.2 and 3.2 Å, respectively.Hydrogen atoms were excluded from this analysis.

Figure S7.
Changing interaction patterns in the alternating-access mechanism of NCX_Mj.(A) Extracellular and intracellular views of representative OF and IF conformations, with 3 Na + ions bound (green spheres), highlighting water molecules within 10 Å of the ion-binding sites (red spheres).The approximate distances between TM helices involved in the opening and closing of intra-and extracellular pathways are indicated in Å. (B) Same views as panel (A), showing phenylalanine sidechains (spheres) that occupy or become released from hydrophobic pockets during the alternating access transition (F182 and F202 on the extracellular side, and F23 and F39 on the intracellular side).Protein residues that are in contact with those phenylalanine sidechains when buried in a hydrophobic pocket are also highlighted, as are lipid molecules in contact with either F23 or F182 when disengaged (yellow).(C) Close-up of the backbone of N-terminal portions of TM2 (TM2a-TM2b) and TM7 (TM7a-TM7b) in the OF and IF state, highlighting selected intramolecular H-bonds (gray lines).

Figure S8.
Global changes in hydrophobic, H-bond and electrostatic interactions alternating-access mechanism of NCX_Mj.Reported quantities are calculated for ensembles of either OF or IF structures, with 3 Na + ions bound, extracted from the free-energy basins revealed in the landscapes shown in Fig. 2. Data is also shown for an ensemble of intermediate configurations.(A) Number of contacts between water and protein apolar groups in each state, calculated as in Fig. S6, for  != 4.2 Å .The residues selected are those that change significantly their accessibility to water during the alternating-access transition, namely: 24,31,38,39,41,42,43,44,46,50,53,56,57,60,61,66,68,69,72,76,138,177,179,183,195,198,199,200,201,202,203,205,209,212,216,220,225,228,231,234,235, and 260.(B) Number of contacts between the apolar portion of the previously mentioned protein residues and the apolar portion of all other protein residues, for each state.The number of contacts between two groups of atoms indexed by  and , for a given simulation snapshot at time t, was calculated as . "# ()/ !0] %& , where r '( is the distance between any pair of atoms and r ! is a constant that depends on the type of interaction; the values used for water-protein and protein-protein contacts were 4.2 and 5.2 Å, respectively.(C) Persistence of H-bonding interactions between residues on the extracellular side of the protein, for each state; only interactions that change significantly during the transition are examined.(D) Same as (C), for residues on the intracellular side.(E) Analysis of the contribution of electrostatic interactions to the free-energy difference between the OF state and either the IF or the intermediate state, evaluated after scaling down the charges of selected atom groups (via a multiplication factor).To quantify the contribution of H-bonds and salt bridges, we scaled down the atomic charges in the polar groups of residues D21, E28, Y60, E180, K187, D197, K198, R223 and E257 (specifically, to the atoms in each of the 'charge groups' bearing a either a neutral or finite integer charge e.g. the -CH2-COO -group of glutamate or aspartate).To assess the contribution of polar residues changing accessibility during the transition, we scaled down the charges of the polar groups of T44, T57 and T203.The contribution of the binding site was evaluated by scaling down the charges of the polar groups of residues coordinating the bound Na + and water molecules, as well the ions and water molecules.To estimate the change in these free-energy differences, we resampled the BE-META trajectories (with 3 Na + bound, for the replicas with bias on ξ !, ξ "# , ξ "! or ξ "$ ; see Table S1), introducing a reweighing factor for each configuration given by  " =  {− ( " ) −  " )/ B , where the term ( " ) −  " ) is the difference between the potential-energy functions with modified and conventional charges, respectively.To quantify this deformation, we evaluated the mean Z-coordinate of the ester layer across the XY plane relative to the mean Z-coordinate of both layers, in the bulk region (at a distance larger that 45 Å from the ion-binding sites center, where the membrane is approximately flat).(B) Hydrophobic thickness of the lipid bilayer in the OF, intermediate and IF states.The thickness was evaluated as the difference between the mean Z-coordinates of the outer and inner ester layers.Contours are shown in intervals of 2 Å.All maps were calculated after first aligning the protein onto a reference frame, using the portion of the structure that is invariant when OF and IF states are compared (Fig. 2).To prevent tilting of the membrane during this process, this alignment was two-dimensional (i.e.only X and Y components of the atomic coordinates are considered).After structural alignment we calculated for each state a threedimensional occupancy map for the ester oxygen atoms of all POPC lipids in either the outer or the inner leaflet, using the Volmap tool of VMD {Humphrey, 1996 #124}.The mean Z-coordinate was calculated as the weighted average of the Z-coordinate for each value of X and Y, with a weight provided by the occupancy at that position.Values of X and Y featuring less than 1% of the overall occupancy are excluded (and correspond to the region occupied by the protein in the center of each map).

Construction and simulation of repeat-swap model of NCX_Mj
To construct the repeat-swap model of NCX_Mj shown in Fig. 1, we first generated a reference sequence by swapping the residues of the two repeats, with the specific sequence alignment provided in Fig. S1, as described elsewhere (3).We then generated 250 models with MODELLER v9.15 (4) using as template the OF X-ray structure of NCX_Mj bound to 3 Na + ions (5).To preserve the C2 pseudo-symmetry of the binding site region with respect to an axis parallel to the membrane plane, we restrained the distances between the 3 Na + ions and their coordinating oxygens as well as the oxygen-oxygen distances.The lowest-energy model was selected for further refinement.First, sidechain rotamers were optimized using SCWRL4 (6).Then, the model was energy-minimized according to the CHARMM27/CMAP forcefield (7,8) and embedded in a solvated POPC bilayer using GRIFFIN (9).An MD simulation of 800 ns was then carried out to equilibrate the model at constant temperature (298 K), pressure (1 atm) and membrane surface area (∼69 Å 2 per lipid).This simulation was carried out with NAMD 2.12(10) and the same force field.This initial equilibration involved a series of stages in which first position restraints and then RMSD restraints are gradually removed on protein sidechain and backbone atoms.Thereafter we performed additional 3 µs MD simulation without any restraint (see Fig. S4 for results).

Improved refinement of crystallographic data for NCX_Mj bound to Ca 2+ and Sr 2+
A statistical survey of the Ca 2+ -binding sites in the PDB (2) suggests that when Ca 2+ is bound to two carboxylate groups in a bidentate mode, the latter ion has a total coordination of eight.To fulfill this coordination geometry, in our previous simulation study of NCX_Mj, based on the OF X-ray structure (2), we placed Ca 2+ in the SCa site and introduced two water molecules coordinating Ca 2+ at near the Smid site.This ion coordination mode has also been used in this study.By contrast, in the original interpretation of the crystallographic data, Ca 2+ was assumed to partially occupy SCa and Smid sites, assisted by one water molecule (2) (Fig. S10).To validate or refute our previous interpretation, we carried out a new refinement of NCX_Mj bound to Ca 2+ based on the electron-density map available in the Protein Data Bank (entry 5HXR).Model building was completed with COOT (11), and structure refinement was carried out with PHENIX (12).The results are shown in Fig. S10.It is clear that the geometry we propose for the Ca 2+ state is just as compatible with the electron-density data as the original model; it is however more plausible, statistically speaking, and is also supported by recent ATR-FTIR experiments (13), which show that Ca 2+ does not reside near residue D240 at the Smid site.To further support the proposed coordination geometry, we also carried out a new refinement of crystallographic data obtained for Sr 2+ (2), which is also transported by NCX_Mj (14).This analysis clearly reveals that in the published structural model (entry 5HXS) the coordination of the Sr 2+ ion is incomplete since additional coordinating atoms ought to be present near the Smid site (Fig. S10).Indeed, a new refinement adding two coordinating water molecules in the defective region significantly improves the agreement with the crystallographic data (Fig. S10).The conclusion from this analysis is that both Ca 2+ and Sr 2+ bind only at the SCa site and in a nearly identical mode, which includes two coordinating water molecules at the Smid site.

Restraining potentials used in the bias-exchange Metadynamics calculations
In the bias-exchange Metadynamics (BE-META) simulations of NCX_Mj with 2 or 3 Na + ions bound, we applied restraining potentials on a coordination variable to preclude Na + dissociation (even though spontaneous dissociation was not observed in unbiased simulations).The coordination variable was defined as: where  # is the distance between ions and coordinating oxygen atoms and  $ = 0.24 nm.For this variable we set lower values of  $ = 4.3 for Na + ions bound to Sint and Sext sites (coordinated to 5 oxygens) and  $ = 4.75 for the ion bound to site SCa (coordinated to 6 oxygens; Fig. 1C), using a potential of the form  !( !() −  $ ) ' , applied when  !() <  $ and for which  != 1000 kJ/mol.The water molecule in the binding site (Fig. 2C) was confined with two half harmonic potentials,  ( (() −  $ ) ) , on distances, (), between the water oxygen and the Cγ atoms of N81 and D240 respectively.This potential was set with  ( = 10000 kJ/(mol nm 2 ) and applied when the above mentioned distances are larger than  $ = 0.46 nm.
A similar scheme was used in the BE-METAD simulations with Ca 2+ bound, but only for the SCa site.This potential was applied on the coordination variable defined by Eq.S1 in which  # is the distance between Ca 2+ and a coordinating oxygen atom, and  $ = 0.3 nm.For this variable we set a lower value of  $ = 7.4 using a potential of the form  !( !() −  $ ) ) , applied when  !() <  $ and for which  != 400 kJ/mol.The two water molecules at the Smid site (Fig. 2C) were also confined with a half harmonic potential,  ( (() −  $ ) ' , on the distances () between the Cγ atoms of either N81 or D240 residues with the oxygen atom of the adjacent water molecule.For this potential  $ = 0.41 nm and  ( = 10,000,000 kJ/(mol nm 4 ).To prevent spurious water penetration and facilitate convergence, we also applied a boundary potential on the coordination variable defined by Eq.S1 in which  # is the distance between Ca 2+ and all oxygen atoms of water molecules (excluding the those coordinating Ca 2+ ), and  $ = 0.3 nm.Distances,  # , beyond 0.5 nm were excluded.For this variable we set an upper value of  $ = 0.5 using a potential of the form  !( !() −  $ ) ) , applied when  !() >  $ and for which  != 400 kJ/mol.To prevent unfolding of binding site portions of TM2 and TM7 we also applied harmonic distance restraints (  ( (() −  $ ) ) , with  $ = 0.41 nm and  ( = 10000 kJ mol -1 nm -4 ) between: 1) carbonyl oxygen of A206 and side chain oxygen of S210; 2) carbonyl oxygen of A206 and amide nitrogen of S210; 3) carbonyl oxygen of V205 and amide nitrogen of T209; 4) homologous distance pairs based on inverted symmetry (A47-S51, M46-T50).

Maximum-entropy reweighting of the simulated ensembles based on HDX-MS data
In its original formulation, the HDXer (15) methods assumes that the input HDX-MS datasets are corrected for deuterium back-exchange.While this is an ideal situation, not all datasets can be produced with this additional layer of information.In particular the HDX-MS data for NCX_Mj do not include this correction (personal communication from Prof. Khananshvili).Therefore, it was necessary to develop an additional strategy that allows for reweighting ensembles on data that have not been corrected for back-exchange levels.
This new scheme assumes that two HDX-MS datasets are available, for which uncorrected deuterium fractions (with values ranging from 0 to 1) of peptide j measured at time t are denoted as  *,, -./& and  *,, -./) respectively.These datasets refer to experiments carried out under two different biochemical/physiological conditions (e.g., the presence or absence of ligands) or for two different mutants of the same protein.All HDX experiments are assumed to be performed with the same protocol so that the back-exchange contribution to a given peptide can be reasonably assumed to be the same across the two datasets.In case those datasets correspond to measurements carried out on two protein variants, only identical peptides must be considered for the analysis (i.e., peptides containing sequence variations should be excluded).
A key approximation in this approach is that back-exchange is accounted for by normalization of the deuterated fraction of a given peptide by a factor,  * , which retains approximately the same value for different timepoints (16,17) and, as previously mentioned, also across the two different HDX datasets.Provided a fully deuterated control were available, the factor  * could be experimentally determined using the relation: , where  &$$ ,  $ are the centroid masses of the fully-deuterated and non-deuterated peptide respectively,  % is the D/H fraction in the buffer and  * is the associated number of exchangeable amide hydrogens.However, in this case we assume that the experimental data are not corrected for back-exchange.Hence, the term  * is unknown and therefore is set as a free parameter during the reweighting procedure.The implication is that, rather than the absolute experimental deuteration levels, the target of the ensemble reweighting is primarily determined by the variations between the deuterium fractions of the two datasets.Hence, the sensitivity of the reweighting procedure will depend on the magnitude of these variations (in addition to the number of data points).The proposed scheme involves two concurrent ensemble reweighting calculations, each targeting one of the HDX-MS datasets, and in which the factors  * are determined at each reweighting step by best fitting to both datasets.The overall likelihood optimization function is given by the sum of the two individual contributions for each dataset (1 and 2): where  1// #4&,) is the apparent work required to reweight the ensemble and  -33 is the error distribution of the experimental data, including both HDX datasets.
To compute the work requires a model of the exchange for a given frame in the ensemble.As in our previous work, we use the exchange model of Best and Vendruscolo (1): in which  5 657# is the predicted protection factor of residue i from reweighted ensemble k,  H,5 and  C,5 are the number of H-bonds and heavy-atom contacts of the amide group of residue i, respectively, and 〈… 〉 # denotes a mean value over the reweighted ensemble k.Then, based on the maximum entropy principle, the apparent work can be expressed as (15): where  5 8 (Χ # ) and  5 !(Χ # ) are the number of H-bonds and heavy atom contacts, respectively, for the amide of residue i, at a molecular configuration, Χ # .We assume that the error distribution,  -33 , is Gaussian: As in our previous work (15), the terms  # 5 8 and  # 5 ! are iteratively determined in order to minimize the likelihood function of Eq.S2: where  denotes the iteration step,  is a parameter that sets the experimental uncertainty and  is a suitably chosen update rate.

Optional optimization of the model parameters
In the previous section we assumed that the model parameters,  8 and  !, have a fixed value, usually set to that derived by Best and Vendruscolo, which best fit a set of HDX measurements ( 8 = 2.0,  != 0.35) (1).To reduce model bias, we set  8 and  ! as optimizable parameters during reweighting, which are sampled using a Monte Carlo algorithm and averaged out during reweighting.The parameter sampling stems from choosing an error probability distribution in which  8 and  ! are integrated out: where ( 8  ! ) is the prior error distribution of the model parameters.Applying a prior distribution can be useful for sparse and noisy experimental data, as it restricts the sampling to a reasonable range of values.Here, ( 8  ! ) was set to mimic the trend of the square residuals as a function of  8 and  ! (Fig. S5B) that was reported by Best and Vendruscolo (1).Specifically, those square residuals describe the discrepancy between predicted and experimental protection factors for a set of proteins.Based on Eq.S9, the update rule for  # 5 !/8 becomes: where 〈… 〉 = !,= " denotes an average performed over  !,  8 parameter space, with probability The states with 2 Na + ions bound to NCX_Mj were generated from the states with 3 Na + ions bound through a slow alchemical transformation that annihilates the ion bound at either the Sext or the Sint site, over a 40-ns MD trajectory.A Cl -ion in the bulk solution was concurrently annihilated to maintain neutrality.The apo states (i.e., no ions) were constructed from the equilibrated states with 2 Na + ions bound by gradually annihilating both ions, again over a 40-ns MD trajectory, together with two Cl -ions in the bulk solution.The states with 2 H + bound were generated from the states with 3 Na + bound by gradually annihilating the Na + ions from the binding sites, and a Cl -ion from the bulk, while creating protonated forms of E54 and E213, over a 100-ns MD simulation.Once initial configurations had been obtained, trajectories ranging from 1 to 2 µs were generated for all occupancy states before initiating the enhanced-sampling simulations.a) Gaussian functions were added to the biasing potential acting on each variable at 5 ps intervals and were initially 2 kJ/mol in height.The height of these Gaussians functions was progressively scaled down as prescribed by the well-tempered Metadynamics method (19).Note the specified 'temperature bias' is not the simulation temperature, which was in all cases 298 K. Following a previous study (2), we used this scheme for an efficient exploration of the multidimensional collective-variable space: replicas with larger values of the temperature-bias and Gaussian widths are able to access high free-energy regions, albeit at low spatial resolution.Replicas with smaller values primarily sample low free-energy regions but at high spatial resolution.
Figure S1.Alignment of the amino-acid sequences of the two inverted topological repeats in NCX_Mj.The sequence identity is 32%.The transmembrane regions in each repeat are indicated and color-coded as in Fig. 1.Residues directly involved in ion-binding are marked with asterisks.The repeat-swap model shown in Fig. 1 was constructed using this alignment.

Figure S2 .
Figure S2.Atom selections used to define the path-collective variables employed in the bias-exchange Metadynamics simulations of the alternating-access transition (see Methods).The specific reaction coordinate associated with each selection is indicated (see also TableS1).The selections include Cα atoms (green spheres), backbone N and O atoms (blue and red spheres), and sidechain carbon atoms (green sticks).Note these selections are approximately symmetric, i.e. they include equivalent residues in both topological repeats.

Figure S3 .
Figure S3.Molecular simulation system and error estimates.(A) All-atom simulation system used to calculate the free-energy landscapes shown in Fig. 2 and Fig. 5.The transporter (orange/blue) is embedded in a hydrated POPC lipid bilayer; electrolyte ions are omitted for clarity.(B) Estimated error of the free-energy maps shown in Fig. 2, based on block analysis.(C) Estimated error of the free-energy maps shown in Fig. 5, also from block analysis.

Figure S4 .
Figure S4.Conventional MD simulations of NCX_Mj bound to 3 Na + ions.(A) Trajectories initiated in the inward-facing (IF) occluded state identified in the free-energy landscape shown in Fig. 2. The plot quantifies the average root-mean-square (RMS) difference between the Ca trace at a given simulation time and the Ca traces in an ensemble of structures extracted from each of the two main free-energy basins in that landscape (within 2 kcal/mol of minimum), i.e. the original IF state, and the outward-facing (OF) occluded state.Unstructured loop regions were omitted in this analysis.(B) Same as (A), for trajectories initiated instead in the OF occluded state.(C) Same as (A), for a trajectory initiated in the repeat-swap model shown in Fig. 1.

Figure S5 .
Figure S5.Evaluation of the predicted conformational ensembles for OF and IF states, in reference to HDX-MS data for two forms of NCX_Mj with different propensities for the OF and IF states.See also Fig. 3C.(A)The plot quantifies the degree to which a conformational ensemble that initially consists of equally weighted populations OF and IF structures (extracted from simulation) must be re-reweighted to produce optimal agreement with the experimental data obtained for each protein construct (see Additional Methods).The degree of reweighting is quantified by an 'apparent work', and the agreement with experiment by the mean-squared deviation between calculated and measured deuteration fractions for a collection of protein fragments.(B) Prior probability distribution of the model parameters used to predict hydrogen-deuterium exchange protection factors based on the forward model of Best and Vendruscolo (1).

Figure S6 .
Figure S6.Changing interaction patterns in the alternating-access mechanism of NCX_Mj.Reported quantities are calculated for ensembles of either OF or IF structures, with 3 Na + ions bound, extracted from the free-energy basins revealed in the landscapes shown in Fig. 2. (A, B) The charts enumerate the set of residues in the extracellular or intracellular side of the protein, respectively, that show a significant change in hydrophobic lipid-tail contacts or protein-protein H-bonding and salt-bridges when OF and IF stats are compared.(C, D) The charts enumerate the set of residues in the extracellular or intracellular side of the protein, respectively, that show a change in hydration number of 2-fold or more; the hydration number was defined using a cutoff distance of 4.2 Å. (C) Representative structures of the OF and IF ensembles illustrating the changes in interaction patterns quantified in panels (A,B) and (C,D).The residues that exhibit changes in water accessibility and lipid contacts are shown as grey and brown volumes, respectively, with water and lipid molecules in their vicinity shown with density isosurfaces (blue and yellow mesh, respectively).Residues that exhibit changes in protein-protein H-bonding and salt bridges are shown as sticks.The number of lipid-tail and H-bond contacts between two groups of atoms in a given simulation snapshot was calculated using distance cut-off values of 5.2 and 3.2 Å, respectively.Hydrogen atoms were excluded from this analysis.

Figure S9 .
Figure S9.Changes in membrane shape and thickness during the alternating access transition.(A) Deformation of each of the lipid bilayer leaflets relative to a flat shape, for the OF, intermediate and IF states.To quantify this deformation, we evaluated the mean Z-coordinate of the ester layer across the XY plane relative to the mean Z-coordinate of both layers, in the bulk region (at a distance larger that 45 Å from the ion-binding sites center, where the membrane is approximately flat).(B) Hydrophobic thickness of the lipid bilayer in the OF, intermediate and IF states.The thickness was evaluated as the difference between the mean Z-coordinates of the outer and inner ester layers.Contours are shown in intervals of 2 Å.All maps were calculated after first aligning the protein onto a reference frame, using the portion of the structure that is invariant when OF and IF states are compared (Fig.2).To prevent tilting of the membrane during this process, this alignment was two-dimensional (i.e.only X and Y components of the atomic coordinates are considered).After structural alignment we calculated for each state a threedimensional occupancy map for the ester oxygen atoms of all POPC lipids in either the outer or the inner leaflet, using the Volmap tool of VMD {Humphrey, 1996 #124}.The mean Z-coordinate was calculated as the weighted average of the Z-coordinate for each value of X and Y, with a weight provided by the occupancy at that position.Values of X and Y featuring less than 1% of the overall occupancy are excluded (and correspond to the region occupied by the protein in the center of each map).

Figure S10 .
Figure S10.Improved molecular refinement of the Ca 2+ and Sr 2+ coordination geometries in outwardfacing NCX_Mj based on previously published crystallographic data(2).The original models (PDB entries 5HXR and 5HXS), shown in the left, are contrasted with our proposed interpretation, shown on the right.2Fo−Fc maps contoured at 1s is shown as a grey mesh and Fo−Fc omit maps at 3s is shown as in magenta.See Supplementary Methods for further details.

Table S2 .
Gaussian widths () and temperature bias (ΔT) used in each of the replicas in the bias-exchange well-tempered Metadynamics simulations of opening and occlusion of NCX_Mj a .