Energetics of ion competition in the DEKA selectivity filter of neuronal sodium channels

The energetics of ionic selectivity in the neuronal sodium channels is studied. A simple model constructed for the selectivity filter of the channel is used. The selectivity filter of this channel type contains aspartate (D), glutamate (E), lysine (K), and alanine (A) residues (the DEKA locus). We use Grand Canonical Monte Carlo simulations to compute equilibrium binding selectivity in the selectivity filter and to obtain various terms of the excess chemical potential from a particle insertion procedure based on Widom's method. We show that K$^{+}$ ions in competition with Na$^{+}$ are efficiently excluded from the selectivity filter due to entropic hard sphere exclusion. The dielectric constant of protein has no effect on this selectivity. Ca$^{2+}$ ions, on the other hand, are excluded from the filter due to a free energetic penalty which is enhanced by the low dielectric constant of protein.


Introduction
Sodium (Na) channels can be categorized on the basis of their function, the cell in which they are found, structure of the protein (both secondary and tertiary), and the structure of the selectivity filter (SF). The SF is a narrow region of the permeation pathway, where the channel discriminates between different ions. The selectivity properties of different channels primarily depend on the type of amino acid motifs present in their SF.
The two most widely studied classes of Na channels are the neuronal (this is the one studied here) and bacterial Na channels. The SF of neuronal Na channels has a DEKA locus made of aspartate (D), glutamate (E), lysine (K), and alanine (A) residues. On the basis of their homology with L-type calcium (Ca) channels [1], these amino acids seem to face the permeation pathway. The accurate structure of the DEKA Na channels is still unknown, so theoretical studies are restricted to using models based on homologies of the known structures or on reduced models based on minimal structural information available. This is the approach used in this work, while the minimal structural information is that the SF has the DEKA locus.
The bacterial Na channels, on the other hand, have X-ray structures measured recently [2][3][4][5]. These channels include NavMs [2,5], NavAb [3,4], and NaChBac [6]. The structure of a Ca 2+ -selective mutant of NavAb is also available [7]. These channels have a lot of aspartates and glutamates in their SF. Therefore, at a first glance, they look like a Ca channel. Hydration plays an important role in the selectivity mechanisms of these channels [8], but this is not the subject of the present study.
Simulation studies for Na channels have been based on models of different resolutions. All-atom, explicit-water models are usually used when X-ray structures are available. They are generally studied using molecular dynamics (MD) simulations [4,5,9]. In the case of the DEKA channel, Lipkind and Fozzard [10] performed MD simulations to explore the Na + vs. K + selectivity for various mutants of DEKA based on extreme homology modelling.
Boda et al. [11][12][13] and Vora et al. [14] used reduced models of Na channels in the implicit solvent framework. In these models, only the SF amino acids were represented in an explicit way, while other parts of the channel protein were reduced into a dielectric body. This is the modelling level that we use in this work. An intermediate approach is that of Finnerty et al. [15], who proposed a localization method, where SF aminoacid terminal groups are localized into certain positions inspired by structural information.
The advantage of reduced models is that they permit the design of simulation setups in time and length scales that mimic experimental setups and are capable of studying a wide range of concentrations and voltage. Also, they make it possible to focus on the essential features of the system (SF structure, pore geometry, bath concentrations, voltage, etc.) and to take the effect of the remaining degrees of freedom into account in an averaged, but physically well-based manner (dielectric response as well as external constraints such as the walls of the channel and the membrane).
Simulations can also be distinguished on the basis of the fact whether they were performed in or out of equilibrium. Equilibrium simulations can study the selective binding of various ions in the SF. Monte Carlo (MC) simulations, especially in the grand canonical (GC) ensemble (Grand Canonical Monte Carlo, GCMC) are ideal tools for this purpose [11][12][13]15]. MD simulations can also be used to study selective binding [10,16,17]. Some properties of transport, however, can be extrapolated even from equilibrium simulations on the basis of the integrated Nernst-Planck equation as suggested by Gillespie et al. [18,19]. Simulating transport requires a dynamical simulation method. These can be MD simulations [4,5,9], Brownian Dynamics simulations [14,[20][21][22], and Dynamical Monte Carlo (DMC) simulations [13,23,24]. Transport can also be studied with theoretical methods such as the Energy Variational approach of Eisenberg et al. [25][26][27][28].
To extend equilibrium binding-selectivity simulations to non-equilibrium situations of steady-state ionic transport is of crucial importance because experimental data are available for ionic currents from electrophysiological measurements [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44]. The relation of the fluxes carried by the competing ions (flux ratio) defines dynamical selectivity. How binding selectivity is related to dynamical selectivity is, however, a non-trivial issue as shown by Rutkai et al. [45]. In particular, the flux is determined not only by the occupancy of a given ionic species in the channel, but also by its mobility.
The main conclusions of those papers [11,12] were that K + ions are excluded from the SF by steric repulsion, while Ca 2+ ions are excluded by an electrostatic penalty. The new aspect of this study is that we provide an energetic analysis for the phenomena described in our 2007 paper [12]. The energetic analysis is performed by separating the free energy (more exactly, the chemical potential) into various terms corresponding to various interactions such as volume exclusion, ion-ion, ion-dielectrics, selfenergy, etc., interactions. This approach was introduced by Gillespie [59,60] in his density functional studies for the Ryanodine Receptor Ca channel and was extended to three-dimensional models including inhomogeneous dielectrics using a GCMC methodology [61] on the basis of Widom's particle insertion method [62,63].
In our previous work, we analyzed the energetics of the selectivity of the L-type Ca channel [61]. In that paper, the dielectric constant of the protein ( pr ) was allowed to be different from that of the baths ( w ). It was shown that the low dielectric protein surrounding the pore focusing the electric field, and thus enhancing the electrostatics, is necessary to reproduce the micromolar selectivity observed for the L-type Ca channels [64,65]. We also extended that work for the case of a dielectric constant different inside the channel ( ch ) from that of the bath [66]. This model is a simple representation of solvation. Our results showed that solvation plays a minor role in the selectivity mechanism of the L-type Ca channel. The explanation is that the solvation penalty for Ca 2+ is balanced by stronger interactions of Ca 2+ with the SF charges. Our simulations extending this work to the DEKA locus are in progress and will be published in a subsequent paper.
In this work, however, we restrict ourselves to the case, where the dielectric constants inside and outside the channel are the same ( ch = w ). This is the model that was considered in our work [12] in 2007. The SF of the Ca channel is highly charged (EEEE locus, four glutamate residues providing −4e charge).
The DEKA filter, on the other hand, is weakly charged (−1e altogether). Therefore, it does not favor divalent ions (Ca 2+ ). Additionally, the bulky terminal group of the lysine is present, which, according to our hypothesis, is there to exclude large ions such as K + . This paper examines how these mechanisms work and their energetic basis.

Model
In our model, most of the atomic structure of the Na channel is reduced to a coarse-grained geometry (figure 1). The channel protein is represented as a continuum solid with dielectric coefficient pr . The three dimensional body of the protein is obtained by rotating the thick line in figure 1 about the r = 0 axis. The protein thus forms an aqueous pore that connects the two baths. Water in the baths and pore is described as an implicit solvent that is a continuum dielectric having uniform dielectric coefficient w = 80. The central cylindrical part of the pore (with radius R = 3.5 − 4.5 Å and length 10 Å) forms the selectivity filter that includes the only atoms of the protein that are treated explicitly. These atoms are four half-charged 'oxygen ions' O 1/2− [ figure 1 (B), red spheres] representing the charged terminal groups of the D and E residues, while a positive 'ammonium ion' NH + 4 [ figure 1 (B), blue sphere] represents the terminal group of the K residue. The alanine is ignored. The structural oxygen ions are confined to the selectivity filter (their centers are in the region r R − R i , |z| 5 Å − R i , where R i is the ionic radius), but they can move freely inside the filter.
The ions are modelled as charged hard spheres with crystal radii (see caption of figure 1). The computation of the intermolecular energy terms due to screened Coulomb potentials and interactions with polarization charges induced on the dielectric boundaries [the boundary of the protein and the electrolyte; thick line in figure 1 (A)] are described in our previous works [55,61,67]. Ions are restricted to the aqueous space of the model and cannot overlap with hard walls in the system. Figure 1 (A) shows only the small central region of the simulation cell. The entire simulation cell is a cylinder having typical dimensions of radius 40 Å and length 180 Å. The channel is embedded in a membrane region that excludes ions by hard walls as described before [55].

Method of energetic analysis
In an equilibrium GCMC simulation, the acceptance of ion insertion/deletions of ions is governed by the configurational chemical potential of the respective ionic species i defined as where k is Boltzmann's constant, T is the temperature, c i (r) is the concentration profile, µ EX i (r) is the excess chemical potential profile, c i (B) is the bulk concentration, and µ EX i (B) is the bulk excess chemical potential. Although kT c i (r) and µ EX i (r) can be different in different regions (they are position dependent), their sum is constant due to equilibrium. The bulk excess chemical potentials µ EX i (B) corresponding to the prescribed bulk concentrations c i (B) are calculated using the Adaptive GCMC method [68]. By rewriting 13601-3 It can be identified with the binding free energy of an ion moved from a bath (B) to position r of the channel [61]. If we write up equation (3.2) for Na + and K + and take the difference, we can derive that which Na + is favored over K + at location r (binding selectivity is positive if location r is selective for Na + over K + ). The corresponding term on the right-hand side containing the bulk concentrations is called 'number advantage' [59] because it expresses the advantage that an ionic species gets from outnumbering the other species in the bulk. The channel can become selective for a given ionic species for two reasons: either from the number advantage or the energetic advantage expressed by ∆∆µ EX (r).
The energetic advantage, however, contains terms due to different interactions present in the system as described in appendix A. The EX term can be divided in various ways. Here, we use the division used in our latest work [66]: (3.5) or briefly EX = HS + II + ID + SELF, where HS means hard sphere exclusion, II means interaction with the ions, ID means interactions with the dielectrics (polarization charges induced by other ions), and SELF means interactions with the polarization charges induced by the ion itself. (In the division of our earlier work [61], we used the DIEL term that contained the SELF term, namely, DIEL = ID + SELF.) We can also use the division EX = HS+MF+SC+SELF, where MF means the interaction with the mean (average) electric field of all the existing charges in the system (ionic and induced). SC expresses correlations beyond the mean field level (SC refers to 'screening') [59]. The SELF term is a one-particle term (mean-field in nature) and corresponds to the average electrostatic interaction energy of the inserted ion with its selfinduced charge. It is not included in the ID or the MF term. The SELF term corresponds to the dielectric boundary force or energy of reference [69].
The computation of all these terms can be found in our original paper [61] and in appendix A. Briefly, the total EX chemical potential can unambiguously be obtained by inserting charged hard spheres (representing the ions) in the Widom particle insertion method. Different terms of EX are computed by inserting particles interacting only through short-ranged (HS) or more direct (II) interactions and obtaining the rest as residuals. For example, it is reasonable to compute the HS term by inserting uncharged hard spheres with the same radius as the respective ion in the Widom procedure. All the remaining terms (II, ID, SELF) are electrostatic in nature and obtained by deducting the HS term from the EX term. (The separation of HS and electrostatic terms and their effect on selectivity can already be found in the work of Nonner et al. [46] in the context of the mean spherical approximation.) Similar procedures are used to separate the II and ID, as well as the MF and SC terms, as described in appendix A.
The r-dependence of various terms is computed by ion insertions into grid cells shown in figure 1.
Note that the concentration profile can be computed in two different ways. First, sampling the number of ions in a volume element, computing the average ion number and dividing by the volume of the element. This is advantageous when the concentration and/or the volume element is large so there is a large enough sample of ions. The concentration, on the other hand, can also be computed from equation (3.1) by computing the EX term from the Widom method and deducting it from the chemical potential. This approach is useful where the concentration is low. This method was used in our simulations for the DEKA channel.
Our grid is two-dimensional because we have rotational symmetry. Our profiles, therefore, are expressed in terms of the (z, r ) cylindrical coordinates. In this work, however, we show the results that are averaged over the r -coordinate where R min (z) = R(z) − R larger ion (z) is the cross-section that is accessible to the center of the larger of the competing ions, [R(z) denotes the radius of the simulation domain at z].

Results and discussion
We start our discussion with competition of ions of the same charge. Specifically, we study the selectivity of Na + over various monovalent ions. In the classical mole fraction experiment, the mole fraction of one ion (Na + , for example) is changed while keeping the total cation concentration constant (when divalent ion is present, the total ionic strength is kept constant in some studies). These results are seen in figure 5 of reference [12]. In this work, the concentration of the two competing cations in the baths is the same (50 mM), so the number advantage is zero. Figure 2 shows the various terms of the ∆µ EX i (z)-profiles for Na + and K + for protein dielectric constant pr = 10 and filter radius R = 3.5 Å. The value pr = 10 is the value fixed in our studies for the L-type Ca channel [18,19,55,56,58,61,70,71]. The value R = 3.5 Å was used in our DMC study for the DEKA Na channel to reproduce experimental data [13]. Therefore, where the EX term (or any component) is negative, it energetically favors the ionic species, so it increases the concentration of that ionic species. As also seen in figure 6 of our previous paper [12], there are peaks at the entrances of the SF and the vestibules (|z| ∼ 5 Å). In the center of the SF, on the other hand, the concentrations are low. This region forms a depletion zone for both ions, where ions have difficulty to enter. The question, therefore, is which ion is excluded less from this region. The answer is that there are more Na + than K + in the SF (the EX term is lower for Na + ), so the SF is Na + -selective.

13601-5
All the electrostatic terms (II, ID, MF, SC) are negative except the SELF term. The SELF term is repulsive because the ions are in the w = 80 region, so the sign of the induced charge on the pr | w boundary is the same as the sign of the inserted ion itself. This practically corresponds to the dielectric penalty an ion must pay when it passes the low dielectric membrane region as described in classical works [72][73][74][75]. The SELF term is slightly larger for Na + because the smaller Na + can get closer to the channel wall and can induce a larger polarization charge. The other term that is positive is the HS term describing the volume exclusion. This is the term that is very different in the case of Na + and K + ; it is larger in the case of K + . Since the size of K + ions (we talk about the dehydrated (Pauling) radius) is larger, it is more difficult to insert such an ion in the SF. Therefore, K + has a larger entropic penalty than Na + does. This difference is especially apparent in the center of the SF, where the NH + 4 (the structural ion representing the large terminal group of the lysine) profile has a peak (see figure 6 of reference [12]). Without the HS term (ions of finite size) we could not get a Na + -selective filter (against K + ) in this model.
The MF term is negative, because the SF is negatively charged. There is no space for the cations to fully neutralize the SF charge. The SC term is similar to the MF term in order of magnitude indicating that mean field theories are not sufficient to study ionic systems in crowded confined spaces such as the SF of ion channels.
The dominant term that drives Na + vs. K + selectivity is the HS term. In figures 3 and 4, therefore, only the differences of the EX and HS terms are shown for various cases. In this special case, where the number advantage is zero, the EX difference is equal to the binding affinity [see equation (3.3)], while the HS term is the dominant term of ∆∆µ EX (r). Since the differences are obtained by deducting the K + terms from the Na + terms, positive values favor Na + . Figure 3 shows the profiles for various pore radii for pr = 10 (top panels) and pr = 80 (bottom panels).
Narrower channels favor Na + even more, as expected, because it is even more difficult to find space for the large K + ions in the small SF compared to Na + . Putting it in another way, Na + vs. K + selectivity is better for narrow channels, where stronger competition is forced by the confinement and lack of space,

13601-6
Energetics of ion competition in Na channel  so the smaller size of Na + has the advantage. The binding affinity curves (left-hand panels) and the HS advantages (right-hand panels) behave similarly with small differences due to other energetic terms (see figure 2).
Another conclusion of the figure is that Na + vs. K + selectivity does not depend on the dielectric constant of the protein; the curves for pr = 10 (top panels) and pr = 80 (bottom panels) behave practically the same. Figure 4 shows the same curves but now for a fixed pore radius (R = 3.5 Å) and different monovalent cations (Li + , K + , Rb + , Cs + ) competing with Na + . The main conclusion is similar to those drawn in figure 3; the crowded SF favors the smaller ion. The pore is selective for Li + against Na + , while it is selective for Na + against the larger ions.
The protein dielectric constant does not have an effect on these profiles. Of course, the value of pr has a large effect on the individual ionic profiles and the occupancies (see figure 8 of reference [12]), but not on the relative ones that we study here.
In the second half of this section, we analyze the competition of Na + against Ca 2+ . The other usual way to study the behavior of the channel having varying electrolyte composition is to keep the concentration of one species fixed (Na + , for example) and to add another species (Ca 2+ , for example) gradually. This added salt experiment was done by Almers and McCleskey in their experiment for the L-type Ca channel [64,65]. We performed this kind of experiment in our previous simulations for the DEKA locus and its DEEA mutant, see figure 2 of reference [12].
Those simulations qualitatively reproduced the experiment of Heinemann et al. [38]. Heinemann et al. found that mutating the DEKA locus into a DEEA locus, the selectivity behavior of the channel is reminiscent to Ca channels rather than Na channels. In experiment, the current drops to half (IC50) at Ca 2+ concentration 10 −4 M, while in our simulations, the number of Na + ions drops to half at the same concentration. The explanation is that the DEEA mutation has −3e charge producing a Ca channel, but with weaker selectivity than in the case of the −4e charge (EEEE locus). The DEKA locus, on the other hand, shows Na + over Ca 2+ selectivity. This selectivity is stronger for smaller pr (see figure 10 (A) of reference [12]). The dielectric constant of the protein, therefore, has a strong effect in the case of monovalent vs. divalent competition.
In figure 5, we show the results only for two chosen Ca 2+ concentrations, 10 mM (top panels) and 40 mM (bottom panels) -both are well above the physiological values (∼1-2 mM).
The background Na + concentration is 50 mM. The Na + and Ca 2+ concentration profiles are shown for pr = 80 (left-hand panels) and 10 (right-hand panels). There are more Na + than Ca 2+ ions in the filter in the case of pr = 10 for both concentrations. A single Na + ion efficiently counterbalances the filter charge. Ca 2+ ions, on the other hand, overcharge the filter, which is electrostatically unfavorable. To counterbalance this overcharge, a Cl − would be needed, but there is no space left for it in the filter.
In the case pr = 80, on the other hand, there are more Ca 2+ ions at [Ca 2+ ] = 40 mM. The explanation is that Ca 2+ is still double charged so the SF attracts it more strongly. The overcharged filter is balanced by Cl − ions from outside the filter. In this case, this is possible because the Coulomb forces are more long-ranged and more screened than in the case of pr = 10, where the low-dielectric protein focuses the electric field. This means that the low dielectric protein is needed to exclude Ca 2+ . The energetics of this phenomenon is analyzed in figures 6 and 7. The difference in Na + vs. Ca 2+ selectivity is more clearly seen by plotting the binding selectivity curves. When this is positive, the pore is Na + -selective, while it is Ca 2+ -selective in the opposite case. The number advantages are also indicated with dashed horizontal lines. As bath Ca 2+ concentration is increased, this line and the binding selectivity curve with it are shifted downwards. The shape of the binding selectivity curves does not change much with the bath Ca 2+ concentration. We can conclude, therefore, that Na + vs. Ca 2+ selectivity does not depend on the bath Ca 2+ concentration. This is because the DEKA locus is a singly occupied SF; only one cation occupies the SF at one time (or none).
This was not true for the L-type Ca channel. That channel could be multiply occupied, so selectivity behavior was a function of Ca 2+ concentration due to correlations of cations in the filter. Furthermore, the SF of the EEEE locus became more charge neutral as Ca 2+ concentration was increased. Because of that, the MF terms decreased (see figure 7 of Boda et al. [61]). That effect is absent here; the probability that a channel becomes charge neutral does not depend on ionic concentrations, but it rather depends on entropic effects (available space in the channel given by filter radius and ion sizes).  The top panels show the pr = 10 data. The left-hand panel shows the II and ID terms (EX = HS + II + ID + SELF), while the right-hand panel shows the MF and SC terms (EX = HS + MF + SC + SELF). The EX and SELF terms are shown both in the left-hand and right-hand sides. The HS term (not shown) is close to zero because the ions are of similar sizes. The EX term is also close to zero in this case, but this is the effect of the balance of the different free energy advantage terms. The SELF term is very positive, so it favors Na + . This term is about four times larger for Ca 2+ than for Na + so it plays the role of solvation penalty in this model. Without the SELF term, we could not get a Na + selective filter (against Ca 2+ ) in this model. Both the II and ID terms (as well as the MF and SC terms, see right-hand panel) favor Ca 2+ because Ca 2+ is attracted twice as strongly by the SF charges (ionic and induced) as Na + .
The bottom panels show the pr = 80 data. Here, the ID and SELF terms are absent, because there is no dielectric boundary present. The ID term favors Ca 2+ , while the SELF term favors Na + . Since the SELF term is larger in absolute value, these two terms together (ID + SELF) still favor Na + , so the channel becomes less Na + selective in their absence.
The SC term is small for pr = 80, which means that Na + vs. Ca 2+ selectivity is chiefly a mean-field effect in this case; the O 1/2− ions attract Ca 2+ twice as strongly as they attract Na + . In the case of pr = 10, on the other hand, SC is quite large indicating a SF of higher density and correlations beyond the meanfield level (mainly, with induced charges).
Summarizing, the EX term is negative for pr = 80, so it is rather a Ca channel. The EX term is close to zero for pr = 10, which means that neither ions are favored energetically. Binding selectivity is driven by the number advantage, which results in a Na + selective channel (against Ca 2+ ) at physiological Ca 2+ concentrations (1-2 mM).

Conclusions
We analyzed the energetics of ion selectivity in the SF of the DEKA Na channels. The reduced model studied before [12] was capable of reproducing the basic characteristics of this channel. We showed that K + ions are excluded from the SF due to entropic hard sphere exclusion. The dielectric constant of the protein has no effect on this selectivity. In general, this filter favors smaller ions over larger ones.
Ca 2+ ions, on the other hand, are excluded from the filter due to a free-energetic penalty which is enhanced by the low dielectric constant of the protein. The DEKA locus works as a Na channel in the Na + vs. Ca 2+ competition by not favoring Ca 2+ . The dominant term is the number advantage in the bulk solutions. In physiological situations this mechanism suffices.
We showed that the dominant term of the energetic penalty is the SELF term, which is a dielectric penalty -the interaction of the ion with the polarization charges induced by itself. This dielectric penalty is a simple, implicit representation of solvation penalty in the framework of this model, where ch = w .
Simulations, where a different dielectric constant inside the channel is used, take solvation into account explicitly.

A. Widom particle insertion method to compute the components of the excess chemical potential
The excess chemical potential profile can be computed using Widom's particle insertion method [62,63]. We divide the simulation cell into small volume elements as described in reference [61] and insert "ghost" particles into uniformly generated positions in these volume elements. We compute the interaction energy U (r) of the "ghost" ion inserted at position r with the whole system and use it in the operation W [U (r)] = −kT ln e −U (r)/kT , (A.1) where the brackets denote GC ensemble average. If the interaction energy U (r) contains all the terms (no matter how it is divided into terms), operator W provides the full excess chemical potential µ EX i (r) = W U HS i (r) +U II i (r) +U ID i (r) +U SELF i (r) .

(A.2)
A diverging term U WALL i (r) corresponding to overlap with protein and membrane walls is omitted in this equation, because we evaluate the excess chemical potential only at the allowed positions.
The II term of the energy is obtained as U II i (r) = j i z i z j e 2 ψ II i j (r, r j ), where ψ II i j (r i , r j ) = 1 8π 0 w |r i − r j |

13601-11
describes the interaction of a unit charge at r i with the polarization charge, h j (r j , s), induced by another unit charge at r j (or vice versa). Vector s runs over the dielectric boundary B. The polarization charge is determined using our Induced Charge Computation method [55,67]. We define the terms in the excess chemical potential that correspond to different interactions as suggested by Gillespie [59]. The definition of these terms is not unique. In our previous work [61], we suggested a possible and physically well-based procedure. The HS term in the excess chemical potential is computed by inserting uncharged hard spheres into the system with the same size as the corresponding ion, but without the charge: µ HS i (r) = W U HS i (r) . (A.5) The II + ID + SELF part is the difference EX − HS. If we insert charged hard spheres into the system, but ignore their interactions with the polarization charges, we can compute an excess chemical potential term describing the ion-ion interactions including the HS interactions: W U HS i (r) +U II i (r) . The II term (that corresponds solely to the interaction with the ionic charges) is obtained by subtracting the HS term: The ID term (that corresponds to the interactions with polarization charges induced by other ions) is what remains: The SELF term is a one-particle term that corresponds to the i = j term of the ID energy in equation (A.4).
The MF term is simply the interaction with the mean electric field computed by sampling with a unit point charge as described in reference [61]. The SC term, again, is what remains: SC = EX − HS − MF − SELF.