Molecular dynamics study of Cl− permeation through cystic fibrosis transmembrane conductance regulator (CFTR)

The recent elucidation of atomistic structures of Cl− channel CFTR provides opportunities for understanding the molecular basis of cystic fibrosis. Despite having been activated through phosphorylation and provided with ATP ligands, several near-atomistic cryo-EM structures of CFTR are in a closed state, as inferred from the lack of a continuous passage through a hydrophobic bottleneck region located in the extracellular portion of the pore. Here, we present repeated, microsecond-long molecular dynamics simulations of human CFTR solvated in a lipid bilayer and aqueous NaCl. At equilibrium, Cl− ions enter the channel through a lateral intracellular portal and bind to two distinct cationic sites inside the channel pore but do not traverse the narrow, de-wetted bottleneck. Simulations conducted in the presence of a strong hyperpolarizing electric field led to spontaneous Cl− translocation events through the bottleneck region of the channel, suggesting that the protein relaxed to a functionally open state. Conformational changes of small magnitude involving transmembrane helices 1 and 6 preceded ion permeation through diverging exit routes at the extracellular end of the pore. The pore bottleneck undergoes wetting prior to Cl− translocation, suggesting that it acts as a hydrophobic gate. Although permeating Cl− ions remain mostly hydrated, partial dehydration occurs at the binding sites and in the bottleneck. The observed Cl− pathway is largely consistent with the loci of mutations that alter channel conductance, anion binding, and ion selectivity, supporting the model of the open state of CFTR obtained in the present study. Supplementary Information The online version contains supplementary material available at 10.1007/s00018-022-04621-7.


Molecular System
The structural model of phosphorylated, ATP-bound zebrafish E1372Q CFTR (PDB: 5W81) was used [1]. Residues missing from the PDB structure include the loops connecting each TMD-NBD pair (406-436, 1182-1202), the R-region (642-843), the extracellular loop between TM7 and TM8 (887-917; also known as the extracellular loop 4), and the segment at the C-terminal end of NBD2 (1459-1485). These missing segments were not modelled in this study, so that our structural model consisted of five peptide chains. All 5 of these chains were acetylated at the N-termini and amidated at the C-termini into primary amides. Mg-ATP moieties at the NBD interface were retained in the simulations. Unlike in the human CFTR PDB structure (6MSM), cholesterol, phospholipids, and the helix at the TMD-NBD interface were absent in the zebrafish CFTR PDB structure.
The PPM server of Orientations of Proteins in Membranes (OPM) database was used to determine the starting position of the lipid bilayer with respect to the CFTR protein [2].
The CFTR protein with Mg-ATP bound was then embedded in a POPC bilayer using InflateGRO2 [3]. The Perl script of InflateGRO2 used was slightly modified to enable the use of a hexagonal simulation box. The hexagonal periodic unit cell configuration was chosen with starting dimensions: a = b = 11 nm, c = 18 nm, α = β = 90°, and γ = 120°.
This procedure resulted in a total of 262 POPC molecules added around the protein.
GROMACS 2018 was used to add SPC216 water to the system [4]. Water molecules within the lipid bilayer and the transmembrane region of the channel pore were removed using VMD [5]. GROMACS was used to add 150 mM excess NaCl to the system along with neutralizing the charge of the system with the same ions.

Simulation setup and protocol
All MD simulations were conducted using GROMACS 2018 [4]. The CHARMM36 forcefield was used for protein, lipids, ions, and ATP, together with the TIP3P water model [6][7][8][9]. Simulations were run in the NpT ensemble (T = 300 K, p = 1 atm) at 2 fs integration timesteps. Constant temperature was maintained using the Nosé-Hoover thermostat (τT = 0.5) [10,11]; constant pressure was maintained using the Parrinello-Rahman barostat (τp = 2.0) [12,13]. Semi-isotropic pressure coupling was used, with isothermal compressibility set to 4.5×10 -5 bar -1 both in the xy-plane and along the z-axis. Nonbonded interactions were calculated using Verlet neighbor lists [14,15]. Lennard-Jones interactions were cut off at 1.2 nm and a force-based switching function with a range of 1.0 nm was used. The particle-mesh Ewald (PME) method was used to compute electrostatic interactions with a real-space cut-off of 1.2 nm [16,17]. The LINCS algorithm was used to constrain covalent bonds involving H atoms [18].
The simulation system was first subjected to steepest descent energy minimization until maximum force dropped below 1000 kJ/mol/nm. Random velocities were generated at the beginning of the NpT equilibration phase, which was conducted in three 10-ns stages, successively with protein heavy atoms, protein backbone atoms, and protein Cα atoms restrained (with force constants of 1000 kJ/mol/nm 2 in x, y, and z directions).
Following the NpT-equilibration phase, the production phase of the run was conducted in NpT ensemble at the same temperature and pressure without applied restraints. New random velocities were generated at the beginning of production runs. Ten 1-μs-long simulations were produced for this system.

Analysis of durations of ion binding at sites 1 and 2 2.1 Determination of Clion binding durations
Ions were considered to be in the channel pore if they were within the boundaries of the cylinder constructed based on the coordinates of the Cα atoms of pore-lining helices, as described in the main Methods section. Only time periods during which Clions resided continuously in the channel pore (without exiting) for at least 50 ns were considered for the following analysis. Sites 1 and 2 are defined as -40 Å ≤ z ≤ -30 Å and -20 Å ≤ z ≤ -10 Å, respectively. We also defined three transition regions: -42.5 Å ≤ z ≤ -40 Å, -30 Å ≤ z ≤ -20 Å, and -10 Å ≤ z ≤ -5 Å. Whenever an ion departed from a binding site and entered a transition region before returning to the original binding site it came from, the time spent in the transition region was added to the binding event at the site.

Calculation of time constants of ion binding
Assuming ion unbinding to be a first-order exponential decay process, the probability of an ion that was bound at t=0 to remain bound to the site at time t is P(t) = , where is the time constant. It follows that is the probability of the ion remaining bound for a duration t' such that 0 < t' < t. Hence, where p(t) is the probability density function of the binding duration t. It follows that (2) Using the method described in 2.1, the durations of all binding events are collected. The binding durations were summarized in a probability density histogram, which is an estimate of the probability density function p(t) (see Fig. S2). We computed two estimates of average binding duration using two different approaches. The first estimate ( ) corresponds to the average binding durations obtained from the above analysis. The average duration of binding is an estimate of , since the expectation value of binding duration <t> is: The second estimate ( ) was obtained by exponential fit to the probability density histogram using equation (2). Note that these methods of estimation assume that the probabilities of Clbinding at sites 1 and 2 are mutually independent.