Potential of mean force calculations on an L-type calcium channel model

Gabriela Barreiro1, Cristiano Ruch Werneck Guimarães2 and Ricardo Bicca de Alencastro

Physical Organic Chemistry Group, Departamento de Química Orgânica, Instituto de Química, UFRJ, Bloco A, Sala 609, Cidade Universitária, Ilha do Fundão, CT, RJ 21949-900, Rio de Janeiro, Brazil

1 To whom correspondence should be addressed. Present address: Department of Chemistry, Wesleyan University, Middletown,CT 06459, USA. E-mail: gbarreiro{at}wesleyan.edu


    Abstract
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
To understand the mechanisms of Na+/Li+ permeation at submicromolar Ca2+ concentrations, Na+/Li+ blocking at higher Ca2+ concentrations (10-6–10-4 M) and Ca2+ permeation at millimolar Ca2+ concentrations, we used our recently described L-type calcium channel model. For this purpose, we obtained potential of mean force (pmf) curves for the position change of one Na+ and one Ca2+ ion inside the channel and for the position change of a second Ca2+ ion when the EEEE locus is coordinated to Ca2+. The pmf curves suggest that (i) at submicromolar Ca2+ concentrations, because of the low velocity of Ca2+ entry in the channel, monovalent ion flux occurs; (ii) at Ca2+ concentrations between 10-6 and 10-4 M, thermodynamic equilibrium between the channel and Ca2+ is achieved; as the coordination of Ca2+ with the locus is more favorable than the coordination of Na+, the monovalent ion flux is blocked; and (iii) to put a second Ca2+ ion inside the channel at an appropriate rate, the Ca2+ concentration should reach millimolar levels. Nevertheless, the entry of a second Ca2+ is thermodynamically unfavorable, indicating that the competition of two Ca2+ ions for the locus leads to Ca2+ permeation.

Keywords: Ca2+ permeation/FDTI method/L-type Ca2+ channel/molecular model/potential of mean force calculations


    Introduction
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
The family of voltage-gated calcium channels (VGCCs) is responsible for the major influx pathway for Ca2+ in different types of cells. Ca2+ ions play important roles in a wide variety of processes such as muscle contraction, cell proliferation, gene expression, neurotransmitter release and other secretion processes and cell death (Catterall, 1988Go; Imoto, 1993Go; Varadi et al., 1995Go, 1999Go).

One low (T type) and five high VGCC types (L, N, P, Q, R) have been studied through pharmacological and electrophysiological studies. Although, up to now, only the pore-forming subunits of three members of T-type calcium channels have been cloned (Perez-Reyes, 1998Go; Perez-Reyes et al., 1998Go; Lee et al., 1999Go; Lacinova et al., 2000Go; Lory et al., 2000Go; McRory et al., 2001Go), the L-type subfamily has been extensively characterized by biochemical approaches. These studies revealed that the L-type calcium channel complex is a heteropentamer consisting of {alpha}1, ß, {alpha}2/{delta} and {gamma} subunits (Figure 1Go) (Varadi et al., 1995Go, 1999Go; Randall and Benham, 1999Go). Ca2+ ions reach the cytoplasm by passing selectively through the pores of the {alpha}1-subunit (Catterall, 1988Go; Imoto, 1993Go; Varadi et al., 1995Go, 1999Go; Randall and Benham, 1999Go), formed by four repeating motifs (MI–MIV), each comprising six hydrophobic segments (S1–S6). A highly conserved segment connecting the S5 and S6 transmembrane domains in each motif, termed the P loop or ‘SS1–SS2’ region, is responsible for calcium selectivity in the pore region (Figure 1Go) (Catterall, 1988Go; Varadi et al., 1999Go). More specifically, four closely aligned glutamate residues, the EEEE locus, localized in the SS2 segment of each P loop segment, molecularly express the calcium selectivity of VGCCs. These negatively charged side chains form the Ca2+ binding site (Yang et al., 1993Go; Varadi et al., 1995Go, 1999Go).



View larger version (46K):
[in this window]
[in a new window]
 
Fig. 1. Molecular structure of the VGCC complex [Adapted from Varadi et al.(1999)Go] with zoom of the proposed transmembrane arrangement of the VGCC {alpha}1-subunit. The dark and light gray cylinders are the SS1 and SS2 segments, respectively [Adapted from Catterall (1988)Go].

 
Two different theories describe the Ca2+ binding site. The first, called the multiple-binding site model, suggests that the EEEE locus is divided into two high-affinity Ca2+ binding sites (Almers and McCleskey, 1984Go; Hess and Tsien, 1984Go; McCleskey, 1999Go). The second, on the other hand, suggests the existence of a single binding site of high affinity (Yang et al., 1993Go; McCleskey, 1999Go). Site-directed mutagenesis studies, in which the E residues were mutated to D, Q, A or K, demonstrated that E residues in the SS2 regions collectively contribute to the high-affinity binding of Ca2+ (Kim et al., 1993Go; Tang et al., 1993Go; Klöckner et al., 1996Go; Varadi et al., 1999Go). Moreover, Yang et al. showed that the neutralization of any of the E residues to a Q shifted the blocking of the monovalent current to a higher Ca2+ concentration (Yang et al., 1993Go). These results, therefore, suggest that the EEEE locus in L-type VGCCs must form a single high-affinity binding site.

Since the three-dimensional structure of VGCCs is not available, recent studies using different theoretical techniques focused on understanding the highly selective mechanism of calcium channels (Doughty et al., 1995Go; Guy and Durrell, 1995; Boda et al., 2000Go, 2001Go, 2002Go; Nonner et al., 2000Go). Some of them (Corry et al. 2000Go; Lipkind and Fozzard, 2001Go; Zhorov et al., 2001Go) used the KcsA bacterial potassium channel crystal structure as a template to build a model of the L-type calcium channel. Although some structural and experimental evidence favor the idea that the KcsA channel is an evolutionary predecessor of sodium and calcium channels (Lipkind and Fozzard, 2000Go, 2001Go), very important differences between potassium and calcium channels suggest that the selectivity filters of both channels are structurally dissimilar. Among these, the following may be mentioned.

  1. The EEEE locus, localized in the ß-turn region of the SS2 segment of each P loop form the Ca2+ selectivity filter (Yang et al., 1993Go; Varadi et al., 1995Go, 1999Go). In K+ channels, on the other hand, the selectivity filter is lined by the carbonyl oxygen atoms of the TVGYG signature sequence of each repeat (Choe et al., 1999Go; Snyders, 1999Go). Their side chains face outward into the protein (Lipkind and Fozzard, 2000Go). For example, the side-chain of the Y residue of each repeat points away from the pore and interacts with other conserved aromatic residues of the pore helices (Snyders, 1999Go).
  2. Yang and co-workers (Yang et al., 1993Go) and Ellinor and co-workers (Ellinor et al., 1995Go), in mutagenesis studies, verified the distinct contributions of each of the E residues to high-affinity divalent cation binding. They concluded that the pore would not be arranged in a symmetrical ring of four strictly equivalent E residues and proposed that the non-equivalence of these amino acids might arise from asymmetry in the positions of their {alpha}-carbon atoms. To explain this non-equivalence, Varadi and co-workers suggested that the E residues that form the ion selectivity filter should be arranged in two close but non-equivalent planes, occupying trans positions (Varadi et al., 1999Go), unlike those of K+ channels or nicotinic acetylcholine receptor channels, where symmetrically arranged rings of amino acid residues are capable of imparting ion selectivity (Varadi et al., 1995Go).
  3. As discussed by Varadi and co-workers (Varadi et al., 1995Go, 1999Go), the SS1 region for Ca2+ channels is composed of amino acids that are predicted to be in an {alpha}-helical arrangement. In recent work (Barreiro et al., 2002Go), we evaluated the dynamic stability of the proposed {alpha}-helical arrangement for the SS1 segments using molecular dynamics (MD) simulations and verified that during the time span of the simulations the {alpha}-helix structures were very stable. Thus, the concept of a ß-hairpin structure for the pore-lining region, as in K+ channels, is probably incorrect.
  4. Following the framework of K+ channels, the amino acids that would form the turn between SS1 and SS2 in Ca2+ channels, valine and isoleucine, are the least likely amino acids to form turns (Varadi et al., 1995Go, 1999Go).
  5. Another significant argument is the strict consensus sequences that may form a bend structure in Ca2+ channels are located in the SS2 region (MEGW, GEDW, FEGW and ATGE in motifs I, II, III and IV, respectively) (Varadi et al., 1995Go, 1999Go).
  6. The outer mouth of the Ca2+ channel vestibule is assumed to be much wider than that of the KcsA channel (Lipkind and Fozzard, 2000Go, 2001Go).
  7. The lengths of the SS1 segments in Ca2+ channels are much greater than those in K+ channels (Schetz and Anderson, 1993Go).
  8. There is a poor sequence identity of transmembrane segments between L-type VGCCs and the KcsA channel (Zhorov, 2001Go).

Because of this, one expects the predicted three-dimensional geometry for the P loop of Ca2+ channels to be different from that of K+ channels. Accordingly, the L-type VGCC pore model recently derived by us (Figure 2Go) (Barreiro et al., 2002Go) was built based purely on experimental information regarding Ca2+ channels and theoretical results obtained from docking and MD simulations. The agreement between the results of our simulations and the experimental observations regarding Ca2+ channels (Yang et al., 1993Go; Ellinor et al., 1995Go; Varadi et al., 1995Go, 1999Go; McCleskey, 1999Go) was very good. More specifically, the results of the MD simulations in our L-type VGCC pore model with one and two Ca2+ ions at the binding site are in accordance with mutation studies, which suggest that the EEEE locus in the L-type calcium channel must form a single high-affinity binding site. These results suggest that the Ca2+ permeation through the channel would indeed be derived from competition between two ions for the single binding site. Furthermore, the experimentally observed blocking of the Na+ flux at micromolar Ca2+ concentrations, probably due to the occupancy of the single high-affinity binding site by one Ca2+, was also reproduced by our model.



View larger version (129K):
[in this window]
[in a new window]
 
Fig. 2. The VGCC pore model with one Ca2+ ion at the binding site.

 
The main interest of this work was to verify if our molecular model of the L-type VGCC pore could contribute to the understanding of the experimental results obtained by Yang and co-workers (Yang et al., 1993Go), discussed in detail by McCleskey (McCleskey, 1999Go). In other words, we intended to test the ability of our model to explain the mechanism of Na+/Li+ permeation at submicromolar extracellular Ca2+ concentrations, Na+/Li+ current blocking at higher Ca2+ concentrations ([Ca2+] {approx} 10-6–10-4 M) and Ca2+ permeation at millimolar Ca2+ concentrations. To accomplish this, we used our L-type VGCC pore model to obtain potential of mean force (pmf) curves, (i) for the position change of a Na+ ion inside the channel, (ii) for the position change of a Ca2+ ion inside the channel and (iii) for the position change of a second Ca2+ ion when the EEEE locus is already coordinated to one Ca2+ ion.


    Methods
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
Potential of mean force calculations

To obtain the pmf curves for the position change of an ion, Na+ or Ca2+ was placed in nine different positions within the channel model. Only its x coordinate (7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 36.5, 37.5 and 38.5 Å) was modified, while the y and z coordinates (Figure 2Go) were kept the same. On the other hand, to obtain the pmf curve for the position change of a second Ca2+ ion with the EEEE locus already coordinated to one Ca2+ ion, the second ion was placed in seven different positions within the channel (7.5, 12.5, 17.5, 22.5, 27.5, 32.5 and 34.5 Å) while the first ion was kept at the x coordinate value of 37.5 Å (the EEEE locus region). In other words, the second ion was placed 30, 25, 20, 15, 10, 5 and 3 Å away from the Ca2+ ion already coordinated to the EEEE locus.

To obtain the free energy variation ({Delta}A) for the position change of an ion inside the channel, we built the thermodynamic cycle exhibited in the center of Figure 3Go. In this case, the free energy variation associated with the annihilation of an ion in a given position inside the channel ({Delta}A4), subtracted from the free energy variation associated with the annihilation of an ion in a different position inside the channel ({Delta}A2), gives the free energy variation associated with the displacement of an ion between these two positions ({Delta}A1) (Figure 3Go). The quantity {Delta}A3 is zero because the ion is annihilated in both the initial and final states. The advantage of this procedure when compared to standard pmf calculations (Beveridge and DiCapua, 1989Go) is that larger changes in the ion position can be obtained with a reduced number of simulations. The non-physical processes of annihilation of an ion were simulated using the finite difference thermodynamic integration (FDTI) method (Mezei, 1987Go).



View larger version (24K):
[in this window]
[in a new window]
 
Fig. 3. Thermodynamic cycles used to obtain the potential of mean force curves. {Delta}A1 and {Delta}A1(real), free energy variations associated with the displacement of an ion between two positions in the model and real systems; {Delta}A2 and {Delta}A2(real), free energy variations associated with the annihilation of an ion in a given position inside the model and real systems; {Delta}A3, in the model system and {Delta}A3(real), in the real system, are zero because the ion is annihilated in both the initial and final states; {Delta}A4 and {Delta}A4(real), free energy variations associated with the annihilation of an ion in a different position inside the model and real systems; {Delta}Atransf1 and {Delta}Atransf2, free energy variations associated with the transfer of the model system with one ion in different positions to the real environment; {Delta}Atransf3 and {Delta}Atransf4, free energy variations associated with the transfer of the model system with one annihilated ion in different positions to the real environment.

 
The finite difference thermodynamic integration (FDTI) method

As we have defined the thermodynamic states throughout this paper by constant number of particles (N), volume (V) and temperature (T), we will be discussing free energy changes in terms of Helmholtz free energy (A). As discussed above, we simulated non-physical processes of annihilation of an ion to obtain potential of mean force curves. Since only potential energy terms were changed between the initial and the final states, there is no kinetic contribution to the free energy variation. Thus, the perturbation method (PM) equation (Beveridge and DiCapua, 1989Go; Jorgensen, 1989Go; Kollman, 1993Go) may be written as


(1)

where U({lambda}i) and U({lambda}i+1) are the potential energy functions for the states {lambda}i and {lambda}i+ 1, respectively and < >{lambda}i refers to an average over the ensemble of configurations generated using U({lambda}i). Another equivalent expression for {Delta}A can be employed in free energy difference calculations. It is referred to as the thermodynamic integration (TI) method (Beveridge and DiCapua, 1989Go; Kollman, 1993Go), where {Delta}A is obtained by the equation


(2)

where is the average value of the potential energy function derivative with respect to the coupling parameter at fixed values of {lambda}.

The FDTI algorithm (Mezei, 1987Go) combines the formalisms of PM and TI. FDTI computes numerically the derivatives of the free energy with respect to the coupling parameter at fixed values of {lambda}, followed by numerical integration using a quadrature scheme (Equation 3Go). In other words, differently from PM, the number of {lambda} to be used in the free energy calculation is not dependent on the degree of overlap between the initial and final phase spaces.


(3)

where n is the number of quadrature points, {Delta}{lambda}i is a parameter that depends on the numerical integration scheme and {delta}{lambda} is the increment used to compute the numerical derivatives.

The long-range electrostatic interaction problem

Since our L-type VGCC pore model is a finite system, the free energy variation associated with the annihilation of an ion inside the channel will always lack contributions from the medium outside the finite system (Åqvist, 1990Go, 1992Go; Hansson and Åqvist, 1995Go). However, it is reasonable to assume that these contributions are similar for different positions of one ion inside the channel. Thus, the contributions from the medium outside the finite system to {Delta}A2 and {Delta}A4 would tend to cancel out, and, consequently, would not be so important to {Delta}A1, the free energy variation associated with the displacement of an ion inside the model. This can be better visualized in Figure 3Go, which shows that the free energy variation associated with the annihilation of an ion in a given position inside the real system [{Delta}A4(real)], subtracted from the free energy variation associated with the annihilation of an ion in a different position inside the real system [{Delta}A2(real)], gives the real free energy variation associated with the displacement of an ion between these two positions [{Delta}A1(real)]. The quantity {Delta}A3(real) is zero because the ion is annihilated in both the initial and final states. By using the thermodynamic cycle shown in Figure 3Go, {Delta}A1(real) may be rewritten as


(4)

where {Delta}Atransf1 and {Delta}Atransf2 are the free energy variations associated with the transfer of the model system with one ion in different positions to the real environment and {Delta}Atransf3 and {Delta}Atransf4 are the free energy variations associated with the transfer of the model system with one annihilated ion in different positions to the real environment. In Equation 4Go, we can safely consider that {Delta}Atransf3 is equal to {Delta}Atransf4 because the ion is already annihilated in this case. If the same is true for {Delta}Atransf1 and {Delta}Atransf2 or, at least, if {Delta}Atransf1 {approx} {Delta}Atransf2, then {Delta}A1(real) {approx} {Delta}A2 - {Delta}A4. In any event, this last equation only holds if all atoms inside the model interact with each other.

Simulation conditions

Before the pmf calculations, all systems generated by placing one ion (Ca2+ or Na+) in different positions inside the channel were energy minimized. Following this procedure, the maximum derivative in each case was less than 0.01 kcal/mol.Å. Then, each system was submitted to an MD simulation of 60 ps (30 ps for the equilibration stage and 30 ps for the data collection) in order to relax its structure to different positions of the ion. The coordinates of the ion were kept fixed in both the energy minimization and MD simulation steps. For the model with two Ca2+, the ion already coordinated to the EEEE locus was harmonically restrained to its position, using a high force constant during this whole preparation process. The last configuration of each MD simulation was selected as an initial configuration to be used in the pmf calculations.

Each alchemical annihilation of one ion in different positions was calculated by the FDTI method. To avoid numerical instability resulting from the presence of a finite charge on ions bearing very small Lennard–Jones parameters, the electrostatic decoupling approach (Bash et al., 1987Go) was used. The electrostatic and van der Waals contributions were scaled down linearly using the {lambda} parameter in separate MD runs of 120 ps each, giving a total of 240 ps of simulation to annihilate completely an ion. First, only the charge was annihilated. Then, the ‘neutral’ ion was transformed into a dummy atom. The coordinates of the ‘disappearing’ ion were kept fixed during the annihilation process. For the model with two Ca2+, the ion already coordinated to the EEEE locus was harmonically restrained to its position using a high force constant during the annihilation process, while the ‘disappearing’ Ca2+ was kept fixed. The number of quadrature points, n, used in each MD run was six. The values of {Delta}{lambda}i ({Delta}{lambda}1 = 0.0857, {Delta}{lambda}2 = 0.1804, {Delta}{lambda}3 = 0.2340, {Delta}{lambda}4 = 0.2340, {Delta}{lambda}5 = 0.1804, {Delta}{lambda}6 = 0.0857) and the quadrature points ({lambda}1 = 0.03377, {lambda}2 = 0.16940, {lambda}3 = 0.38069, {lambda}4 = 0.61931, {lambda}5 = 0.83060, {lambda}6 = 0.96623) (see Equation 3Go) were automatically calculated by the Gaussian–Legendre quadrature method (Press et al., 1986Go). The increment {delta}{lambda} used in Equation 3Go was 0.005.

The equations of motion were integrated every 1.0 fs using the Verlet Leapfrog algorithm (Hockney, 1970Go). In all simulations, the ensemble average of the Boltzmann factor, e{Delta}U/kT, was evaluated through an MD run of 10 ps after 10 ps of equilibration at each {lambda}i, giving a total of 120 ps of simulation for each transformation and 240 ps for a complete annihilation. For further analysis, the trajectory was sampled every 0.5 ps. The convergence at each {lambda}i was verified by plotting the Boltzmann factor vs time. The FDTI method computes {Delta}Ai/{delta}{lambda} at each {lambda}i sampling, both forward [{Delta}A({lambda}i -> {lambda}i + {delta}{lambda})/{delta}{lambda}] and backward. [{Delta}A({lambda}i -> {lambda}i - {delta}{lambda})/{delta}{lambda}]. {Delta}Ai/{delta}{lambda} at {lambda}i is the first quantity plus the negative of the second one, divided by two. Similar values for [{Delta}A({lambda}i -> {lambda}i + {delta}{lambda})/{delta}{lambda})]and [–{Delta}A({lambda}i -> {lambda}i - {delta}{lambda})/{delta}{lambda}] are an additional measure of convergence.

All MD runs were performed in an NVT ensemble (T = 310 K). During MD runs, the temperature was maintained at 310 K via the velocity-scaling algorithm (Woodcock, 1971) at the equilibration stage and via a weak coupling to an external temperature bath with a time constant of 0.1 ps (Berendsen et al., 1984Go) at the data collection stage. Acetyl and N-methylamine blocking groups were used to cap the truncation points in the model. To avoid anomalous behavior in the channel structure, owing to the absence of the discarded amino acids, the methyl carbon atoms of the acetyl and N-methylamine blocking groups were kept fixed. In addition, to avoid evaporation of water, a 5 Å layer of water molecules surrounding the whole system was kept fixed during all simulations (Figure 2Go).

Since there is no charge variation in the transformation of a ‘neutral’ ion into a dummy atom, a spherical residue-based non-bonded interaction cutoff of 10 Å was applied to save computational time. To ‘turn off’ the interactions smoothly, we employed a quintic spline function from 8.5 to 10 Å. In this manner, discontinuities in the potential energy surface were minimized. Whenever an atom moved more than half the length of the buffer region (between 10 and 11 Å), the neighbor list was updated. This ensured that no atoms outside the buffer region were able to move close enough to interact. On the other hand, the free energy variation associated with the annihilation of charges would always lack contributions from the medium outside the finite system. As mentioned above, in this case, all atoms inside the model had to interact with each other. Thus, we applied the double-cutoff methodology (Auffinger and Beveridge, 1995Go). The first cutoff had a 10 Å radius, while the radius of the second cutoff was long enough to include interactions of all atoms with every other atom.

All calculations were performed using molecular modeling software from MSI (USA), running on a Silicon Graphics O2 R10000 workstation. The InsightII program (version 97.0) was employed as a graphical interface for the construction and visualization of molecular structures. The all-atom CVFF force field (Dauber-Osguthorpe et al., 1988Go), within the Discover program (version 2.9.7), was employed in energy minimization and MD simulations.


    Results and discussion
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
Pmf curves for the position change of one Ca2+ ion and one Na+ ion inside the unoccupied L-type VGCC pore model

The pmf curves for the change in position of one Na+ ion and one Ca2+ ion in the unoccupied L-type VGCC pore model are shown in Figure 4aGo and b. As expected, the pmf curves point out the EEEE locus as the region of coordination, more precisely the free energy minimum was located in x coordinate values of 37.5 and 36.5 Å for Na+ and Ca2+, respectively. The pmf curves also show that the interaction of one Ca2+ ion with the binding site is more favorable than the interaction of one Na+ by ~140 kcal/mol, with respect to the entrance of the channel (x coordinate value of 7.5 Å). This result shows that our model reproduces the high selectivity of the channel for Ca2+ and suggests that this high selectivity derives from strong electrostatic interactions between one ion with charge 2+ and the EEEE locus with charge 4–.



View larger version (14K):
[in this window]
[in a new window]
 
Fig. 4. Potential of mean force curves (a) for the position change of an Na+ ion inside the channel, (b) for the position change of a Ca2+ ion inside the channel and (c) for the position change of a second Ca2+ ion when the EEEE locus is already coordinated to one Ca2+ ion.

 
Interestingly, Figure 4a and bGo show a free energy barrier inside the channel (located at an x coordinate value of 17.5 Å) for both Na+ and Ca2+. This probably derives from the existence of hydrophobic amino acids in this region (Figure 5Go) and suggests that the entry of ions into the channel is kinetically controlled. With respect to the entry into the channel, the free energy barriers for Na+ and Ca2+ are 12.96 and 34.46 kcal/mol, respectively. Therefore, while the entry of one Na+ into the channel is much more favorable kinetically than the entry of one Ca2+, the entry of one Ca2+ is much more favorable thermodynamically than the entry of one Na+.



View larger version (66K):
[in this window]
[in a new window]
 
Fig. 5. Hydrophobic amino acids responsible for the free energy barrier located at an x coordinate value of 17.5 Å: (a) for Na+; (b) for Ca2+. The other amino acids and water molecules are not shown for better visualization.

 
Pmf curves for the position change of one Ca2+ ion inside the L-type VGCC pore model when the EEEE locus is already coordinated to one Ca2+ ion

Figure 4cGo shows the pmf curve for the change in position of a second Ca2+ ion in the L-type VGCC pore model when the EEEE locus is already coordinated to one Ca2+. The pmf curve shows that the entry of a second Ca2+ ion into the channel is very unfavorable thermodynamically. More specifically, the free energy minimum located near the EEEE locus is ~50 kcal/mol less stable than the free energy minimum near the entrance of the channel. This less stable free energy minimum is located at an x coordinate value of 32.5 Å or 5 Å away from the Ca2+ ion coordinated to the EEEE locus. Although the EEEE locus bears a negative charge of 4–, this last result shows that it is able to accommodate only one Ca2+ ion at a time.

Figure 4cGo shows a free energy barrier for the entry of the second Ca2+ ion in the channel, also located at an x coordinate value of 17.5 Å. As discussed above, this is probably due to the presence of hydrophobic amino acids in this region. This barrier is, however, much higher in free energy (114.62 kcal/mol with respect to the entrance of the channel) than the free energy barrier for the entry of one Ca2+ or one Na+ ion when the EEEE locus is unoccupied. This last result suggests that the entry of a second Ca2+ ion into the channel is considerably slower than the entry of an ion (Ca2+ or Na+) when the channel is unoccupied.

In Figure 4Go, it could be argued that the values of the free energy barriers and free energy minima in the pmf curves are overestimated. This is probably derived from non-equivalence between {Delta}Atransf1 and {Delta}Atransf2 (Figure 3Go). In other words, the contributions to {Delta}A2 and {Delta}A4 from the medium outside the finite system would indeed vary with the position of the ion inside the channel and, consequently, would be important to {Delta}A1. However, in spite of this overestimation, as may be noted below, the pmf curves give valuable interpretations to the experimental studies performed by Yang and co-workers.

A possible interpretation for the mechanism of ion permeation in L-type VGCCs

As mentioned above, the main goal of this work was to verify if our molecular model of the L-type VGCC pore could contribute to the understanding of the experimental studies performed by Yang and co-workers (Yang et al., 1993Go), discussed in detail by McCleskey (McCleskey, 1999Go). In their experimental studies, Yang and co-workers verified that at submicromolar extracellular Ca2+ concentrations ([Ca2+] < 10-6 M), Ca2+ channels were highly permeable to monovalent cations, such as Na+ and Li+. As described previously, Figure 4a and bGo show the existence of a free energy barrier inside the channel, located at an x coordinate value of 17.5 Å. Whereas the free energy barrier for Na+ is 12.96 kcal/mol, for Ca2+ it is 34.46 kcal/mol. The smaller value for Na+ suggests that its entry into the unoccupied channel is considerably faster than the entry of Ca2+. In this manner, if the kinetics of entry of Ca2+ into the channel depend on its extracellular concentration, its velocity of entry at extracellular Ca2+ concentrations <10-6 M should be considerably slower than the velocity of entry of Na+. Therefore, Ca2+ ions would be almost absent from the high-affinity binding site, allowing the flux of monovalent ions.

Yang and co-workers also verified that when the extracellular concentration of Ca2+ increases (10-6 M < [Ca2+] < 10-4 M), the Na+ flux is no longer observed. Again, assuming that the kinetics of Ca2+ entry depend on its extracellular concentration, an increase in its concentration to 10-6–10-4 M would raise the velocity of entry of Ca2+. In this manner, Ca2+ would be able to enter the channel at a rate comparable to the rate of entry of Na+. Since the kinetics would no longer be an issue, thermodynamic equilibrium between the channel and extracellular Ca2+ could be achieved. As shown by the pmf curves in Figure 4a and bGo, the coordination of one Ca2+ ion with the EEEE locus is more favorable than the coordination of one Na+ by ~140 kcal/mol. This result shows that when Ca2+ enters the channel it occupies the EEEE locus and blocks the monovalent ions flux through the channel. We arrived at the same conclusion through MD studies in our L-type VGCC pore model in the presence of the ion pair Na+/Ca2+ at the EEEE locus (Barreiro et al., 2002Go).

As discussed by Yang et al. (Yang et al., 1993Go) and McCleskey (McCleskey, 1999Go), when the extracellular Ca2+ concentration reaches millimolar levels ([Ca2+] > 10-3 M), which corresponds to physiological concentrations, it is possible to observe Ca2+ flux through the channel. Figure 4cGo shows the existence of a free energy barrier inside the channel for the entry of a second Ca2+ ion, located also at an x coordinate value of 17.5 Å. However, this barrier is much higher (114.62 kcal/mol) than the free energy barrier for the entry of Ca2+ (34.46 kcal/mol) and Na+ (12.96 kcal/mol) when the EEEE locus is unoccupied. This shows that the entry of a second Ca2+ into the channel is considerably slower than the entry of an ion when the channel is empty. Assuming again that the kinetics of Ca2+ entry depend on its extracellular concentration, these results seem to indicate that in order to put a second Ca2+ ion inside the channel at a reasonable rate, the Ca2+ concentration outside the cell should reach millimolar levels ([Ca2+] > 10-3 M). Nevertheless, Figure 4cGo shows that the entry of a second Ca2+ is thermodynamically unfavorable by ~50 kcal/mol, indicating that the EEEE locus is able to accommodate only one Ca2+ ion at a time. Thus, when two Ca2+ ions occupy simultaneously a single high-affinity binding site, both compete with each other for the E residues of the locus, decreasing the affinity of the pore for Ca2+. However, to have Ca2+ permeation through the channel, the energy barrier for the Ca2+ ion closest to the cytoplasm to leave the EEEE locus towards the cell interior should be lower than the energy barrier for the second Ca2+ ion to leave the EEEE locus towards the outside of the cell. We demonstrated that this is exactly the case through MD simulations in our L-type VGCC pore model with two Ca2+ ions at the EEEE locus (Barreiro et al., 2002Go). We showed that the coordination with the EEEE locus of the Ca2+ ion closest to the cytoplasm was almost lost at around 65 ps of simulation. Conversely, we observed that during the MD simulation the full coordination of the second Ca2+ ion with the single high-affinity binding site started to be re-established.

Conclusions

In this work, we used our recently derived L-type VGCC pore model to try to understand the mechanisms of Na+/Li+ permeation at submicromolar extracellular Ca2+ concentrations, Na+/Li+ current blocking at higher Ca2+ concentrations ([Ca2+] {approx} 10-6–10-4 M) and Ca2+ permeation at millimolar Ca2+ concentrations. Our interpretation of the pmf curves suggests that at submicromolar extracellular Ca2+ concentrations, because of the low velocity of Ca2+ entry into the channel, Ca2+ ions would be almost absent from the EEEE locus, leading to a flux of monovalent ions. At extracellular Ca2+ concentrations between 10-6 and 10-4 M, on the other hand, Ca2+ is able to enter the channel at a more competitive rate compared with the entry of Na+ and thermodynamic equilibrium between the channel and extracellular Ca2+ can be established. Consequently, the monovalent ion flux is blocked because the coordination of one Ca2+ ion with the EEEE locus is much more favorable than the coordination of one Na+. Finally, at millimolar extracellular Ca2+ concentrations, a second Ca2+ ion is able to enter the channel. However, as the entry of a second Ca2+ is thermodynamically unfavorable, when two Ca2+ ions occupy simultaneously the single high-affinity binding site, both compete with each other for the E residues of the locus, leading to Ca2+ permeation. In conclusion, it may be noted that the pmf curves obtained with our model give valuable interpretations to the experimental results obtained by Yang and co-workers.


    Notes
 
2 Present address: Department of Chemistry, Yale University, New Haven, CT 06520, USA Back


    Acknowledgments
 
This research received partial financial support from the Brazilian agencies CNPq and FAPERJ (grants E26/151494/99 and E26/151898/2000). G.B. acknowledges CAPES, Brazil, and C.R.W.G. acknowledges, FAPERJ, Brazil, for scholarships.


    References
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
Almers,W. and McCleskey,E.W. (1984) J. Physiol., 153, 585–608.

Åqvist,J. (1990) J. Phys. Chem., 94, 8021–8024.[ISI]

Åqvist,J. (1992) J. Mol. Struct. (THEOCHEM), 256, 135–152.

Auffinger,P. and Beveridge,D.L. (1995) Chem. Phys. Lett., 234, 413–415.[CrossRef][ISI]

Barreiro,G., Guimarães,C.R.W. and Bicca de Alencastro,R. (2002) Protein Eng., 15, 109–122.[Abstract/Free Full Text]

Bash,P.A., Singh,U.C., Langridge,R. and Kollman,P.A. (1987) Science, 236, 564–568.[ISI][Medline]

Berendsen,H.C.J., Postma,J.P.M., van Gunsteren,W.F., DiNola,A. and Haak,J.R. (1984) J. Chem. Phys., 81, 3684–3690.[CrossRef][ISI]

Beveridge,D.L. and DiCapua,F.M. (1989) Annu. Rev. Biophys. Biophys. Biochem., 18, 431–492.[CrossRef][ISI][Medline]

Boda,D., Busath,D.D., Henderson,D. and Sokotowski,S. (2000) J. Phys. Chem. B, 104, 8903–8910.[CrossRef][ISI]

Boda,D., Henderson,D. and Busath,D.D. (2001) J. Phys. Chem. B, 105, 11574–11577.[CrossRef][ISI]

Boda,D., Busath,D.D. and Henderson,D. (2002) Appl. Surf. Sci., 7834, 1–3.[CrossRef]

Catterall,W.A. (1988) Science, 242, 50–61.[ISI][Medline]

Choe,S., Kreusch,A. and Pfaffinger,P.J. (1999) Trends Biochem. Sci., 24, 345–349.[CrossRef][ISI][Medline]

Corry,B., Allen,T.W., Kuyucak,S. and Chung,S.-H. (2000) Biochim. Biophys. Acta, 1509, 1–6.[ISI][Medline]

Dauber-Osguthorpe,P., Roberts,V.A., Osguthorpe,D.J., Wolff,J., Genest,M. and Hagler,A.T. (1988) Proteins: Struct. Funct. Genet., 4, 31–47.[ISI][Medline]

Doughty,S.W., Blaney,F.E. and Richards,W.G. (1995) J. Mol. Graphics, 13, 342–348.[CrossRef][ISI][Medline]

Ellinor,P.T., Yang,J., Sather,W.A., Zhang,J.-F. and Tsien,R.W. (1995) Neuron, 15, 1121–1132.[ISI][Medline]

Guy,H.R. and Durell,S.R. (1995) J. Gen. Physiol. Ser., 48, 1–16.

Hansson,T. and Åqvist,J. (1995) Protein Eng., 8, 1137–1144.[Abstract]

Hess,P. and Tsien,R.W. (1984) Nature, 309, 453–456.[ISI][Medline]

Hockney,R.W. (1970) Methods Comput. Phys., 9, 136–211.

Imoto,K. (1993) FEBS Lett., 325, 100–103.[CrossRef][ISI][Medline]

Jorgensen,W.L. (1989) Acc. Chem. Res., 22, 184–189.[ISI]

Kim,M-S., Morii,T., Sun,L.-X., Imoto,K. and Mori,Y. (1993) FEBS Lett., 318, 145–148.[CrossRef][ISI][Medline]

Klöckner,U., Mikala,G., Schwartz,A. and Varadi,G. (1996) J. Biol. Chem., 271, 22293–22296.[Abstract/Free Full Text]

Kollman,P.A. (1993) Chem. Rev., 93, 2395–2417.[ISI]

Lacinova,L., Klugbauer,N. and Hofmann,F. (2000) Gen. Physiol. Biophys., 19, 121–136.[ISI][Medline]

Lee,J.H., Daud,A.N., Cribbs,L.L., Lacerda,A.E., Pereverzev,A., Klockner,U., Schneider,T. and Perez-Reyes,E. (1999) J. Neurosci., 19, 1912–1921.[Abstract/Free Full Text]

Lipkind,G.M. and Fozzard,H.A. (2000) Biochemistry, 39, 8161–8170.[CrossRef][ISI][Medline]

Lipkind,G.M. and Fozzard,H.A. (2001) Biochemistry, 40, 6786–6794.[CrossRef][ISI][Medline]

Lory,P., Monteil,A., Chemin,J., Leuranguer,V., Bourinet,E. and Nargeot,J. (2000) Therapie, 55, 249–254.[ISI][Medline]

McCleskey,E.W. (1999) J. Gen. Physiol., 113, 765–772.[Free Full Text]

McRory,J.E., Santi,C.M., Hamming,K.S.C., Mezeyova,J., Sutton,K.G., Baillie,D.L., Stea,A. and Snutch,T.P. (2001) J. Biol. Chem., 276, 3999–4011.[Abstract/Free Full Text]

Mezei,M. (1987) J. Chem. Phys., 86, 7084–7088.[CrossRef][ISI]

Nonner,W., Catacuzzeno,L. and Eisenberg,B. (2000) Biophys. J., 79, 1976–1992.[Abstract/Free Full Text]

Perez-Reyes,E. (1998) J. Bioenerg. Biomembr., 30, 313–318.[CrossRef][ISI][Medline]

Perez-Reyes,E., Cribbs,L.L., Daud,A., Lacerda,A.E., Barclay,J., Williamson,M.P., Fox,M., Rees,M. and Lee,J.H. (1998) Nature, 391, 896–900.[CrossRef][ISI][Medline]

Press,W.H., Flanney,B.P., Teukolsky,S.A. and Vetterling,W.T. (1986) In Numerical Recipes. The Art of Scientific Computing. Cambridge University Press, Cambridge.

Randall,A. and Benham,C.D. (1999) Mol. Cell. Neurosci., 14, 255–272.[CrossRef][ISI][Medline]

Schetz,J.A. and Anderson,P.A.V. (1993) Biol. Bull., 185, 462–466.[Abstract/Free Full Text]

Snyders,D.J. (1999) Cardiovasc. Res., 42, 377–390.[CrossRef][ISI][Medline]

Tang,S., Mikala,G., Bahinski,A., Yatani,A., Varadi,G. and Schwartz,A. (1993) J. Biol. Chem., 268, 13026–13029.[Abstract/Free Full Text]

Varadi,G., Mori,Y., Mikala,G. and Schwartz,A. (1995) Trends Pharmacol. Sci., 16, 43–49.[CrossRef][ISI][Medline]

Varadi,G., Strobeck,M., Koch,S., Caglioti,L., Zucchi,C. and Palyi,G. (1999) Crit. Rev. Biochem. Mol. Biol., 34, 181–214.[Abstract/Free Full Text]

Woodock,L.V. (1971) Chem. Phys. Lett., 10, 257–261.[CrossRef][ISI]

Yang,J., Ellinor,P.T., Sather,W.A., Zhang,J.-F. and Tsien,R.W. (1993) Nature, 366, 158–161.[CrossRef][ISI][Medline]

Zhorov,B.S., Folkman,E.V. and Ananthanarayanan,V.S. (2001) Arch. Biochem. Biophys., 393, 22–41.[CrossRef][ISI][Medline]

Received June 4, 2002; revised October 22, 2002; accepted January 17, 2003.





This Article
Abstract
FREE Full Text (PDF)
Alert me when this article is cited
Alert me if a correction is posted
Services
Email this article to a friend
Similar articles in this journal
Similar articles in ISI Web of Science
Similar articles in PubMed
Alert me to new issues of the journal
Add to My Personal Archive
Download to citation manager
Search for citing articles in:
ISI Web of Science (2)
Request Permissions
Google Scholar
Articles by Barreiro, G.
Articles by de Alencastro, R. B.
PubMed
PubMed Citation
Articles by Barreiro, G.
Articles by de Alencastro, R. B.