Insights into the Tunnel Mechanism of Cholesteryl Ester Transfer Protein through All-atom Molecular Dynamics Simulations

Cholesteryl ester transfer protein (CETP) mediates cholesteryl ester (CE) transfer from the atheroprotective high density lipoprotein (HDL) cholesterol to the atherogenic low density lipoprotein cholesterol. In the past decade, this property has driven the development of CETP inhibitors, which have been evaluated in large scale clinical trials for treating cardiovascular diseases. Despite the pharmacological interest, little is known about the fundamental mechanism of CETP in CE transfer. Recent electron microscopy (EM) experiments have suggested a tunnel mechanism, and molecular dynamics simulations have shown that the flexible N-terminal distal end of CETP penetrates into the HDL surface and takes up a CE molecule through an open pore. However, it is not known whether a CE molecule can completely transfer through an entire CETP molecule. Here, we used all-atom molecular dynamics simulations to evaluate this possibility. The results showed that a hydrophobic tunnel inside CETP is sufficient to allow a CE molecule to completely transfer through the entire CETP within a predicted transfer time and at a rate comparable with those obtained through physiological measurements. Analyses of the detailed interactions revealed several residues that might be critical for CETP function, which may provide important clues for the effective development of CETP inhibitors and treatment of cardiovascular diseases.

density lipoprotein (LDL) cholesterol (2,3). In plasma, cholesteryl ester transfer proteins (CETPs) transfer cholesteryl esters (CEs) from high density lipoprotein (HDL) to low density lipoprotein and very low density lipoprotein (VLDL) (4). Human plasma CETP is a 476-amino acid hydrophobic glycoprotein. Genetic deficiency of CETP increases HDL-C levels and decreases LDL cholesterol levels (5,6). As a result, CETP has become a promising drug target for the prevention and treatment of cardiovascular diseases. To date, four synthetic CETP inhibitors, torcetrapib, dalcetrapib, anacetrapib, and evacetrapib, have been evaluated in Phase III clinical studies in the past decade (7)(8)(9)(10). The clinical data have shown that all CETP inhibitors elevate HDL-C levels, and nearly all CETP inhibitors decrease LDL cholesterol levels (7)(8)(9)(10). Although anacetrapib is still being tested in clinical studies, the failures of torcetrapib (due to off-target effects) (7) and dalcetrapib and evacetrapib (due to their unexpected low efficacy) (8,11) highlight that the detailed mechanism underlying the CETP-mediated transfer of CE molecules between HDL and LDL/VLDL needs to be elucidated.
Three hypotheses of the CETP mechanism have been proposed: a shuttle mechanism in which CETP collects CEs from HDL and then delivers them through the aqueous environment to LDL/VLDL (12); a ternary complex mechanism in which CETP simultaneously interacts with HDL and LDL/VLDL, thus forming a ternary complex that induces CE transfer from HDL to LDL/VLDL (13); and a modified ternary complex mechanism that implicates a CETP dimer instead of a monomer (14).
The crystal structure of CETP revealed a banana-shaped molecule composed of N-and C-terminal ␤-barrel domains, a central ␤-sheet, and an ϳ60-Å-long hydrophobic central cavity (15). Two pores near the central ␤-sheet domain occupied by two phospholipid molecules have been suspected to be the gates between which the central cavity would interact with the aqueous environment or lipoproteins, thus favoring the shuttle mechanism (15).
Experimental studies on the kinetics of the plasma proteincatalyzed exchange of phosphatidylcholine and CE between plasma LDL and HDL have suggested that the exchange occurs by a sequential mechanism involving an HDL-CETP-LDL ternary complex (13). Recent electron microscopy (EM) studies have shown the presence of a ternary complex among CETP, HDL, and LDL/VLDL (16) in which the CETP N-terminal distal end interacts with HDL via a hydrophobic interaction (17). Independently, molecular dynamics (MD) simulations have revealed that the CETP distal ends are structurally flexible (18), and the N-terminal end penetrates into the HDL surface and takes up a CE molecule (19). A tunnel mechanism has been proposed because the CETP crystal structure displays a long hydrophobic central cavity containing two CE molecules (16).
Here, we used all-atom MD simulations to evaluate the tunnel mechanism by testing whether a long tunnel could be formed and is large enough to allow CE to pass through the entire CETP molecule.

Experimental Procedures
Strategy for Setting Up the Simulation System-In the tunnel mechanism, CETP interacts with HDL and LDL simultaneously and forms a ternary complex (see Fig. 1, A and B). EM experiments have suggested that the flexible N-terminal ␤-barrel domain of CETP penetrates into the phospholipid monolayer on the HDL surface via hydrophobic interactions (16,17) and that the distal end reaches the HDL CE-rich pool. However, the flexible C-terminal ␤-barrel domain penetrates into the phospholipid monolayer of the LDL surface and interacts with the inner LDL triglyceride (TG)-rich pool (16). To test whether the ternary complex could form a tunnel within CETP that allows a CE molecule to transfer through the entire CETP molecule, we constructed an MD simulation system to match the observed EM ternary complex as closely as possible. In this system, a CETP molecule was embedded into an aqueous environment with its two distal ␤-barrel domain ends penetrating into the two opposite phospholipid monolayers that envelop the CE and TG pools (Fig. 1, C and D).
Simulation System Components-To reach equilibrium for the above simulated system within a practical period of time, our strategy was to separately equilibrate each component, i.e. the CE and TG pools, CETP, the phospholipid monolayers, and their simple combinations. In brief, the CE and TG pools were constructed from a non-CHARMM force field equilibration (20,21) with a periodically expanded cross-sectional area of 96 ϫ 96 Å. The CE pool and TG pool were then submitted for equilibration for 74.0 and 9.5 ns under a CHARMM force field (22)(23)(24)(25) and physiological conditions using all-atom MD simulations. The equilibrated phospholipid monolayers were obtained from an equilibrated 1-palmitoyl-2-oleoyl-snglycero-3-phosphocholine (POPC) bilayer (26) that was further equilibrated for 20.0 ns by using a CHARMM force field (24,25). All equilibrations were conducted with the Nanoscale Molecular Dynamics version 2 (NAMD2) software package (from the University of Illinois at Champaign-Urbana) (27) at 310 K and 1 atm. During the equilibrations, the 310 K temperature was maintained through Langevin dynamics (28) with a damping coefficient of 5 ps, whereas the 1-atm pressure was maintained by using the Langevin piston Nose-Hoover method (29) with a piston period of 100 fs and a decay time of 50 fs. Periodic boundary conditions and a cutoff distance of 12 Å for van der Waals interactions were applied, and the particle mesh Ewald method (30) (grid spacing Ͻ1 Å) was used to compute the long range electrostatic interactions.
The equilibrated CETP was obtained by equilibrating the CETP crystal structure (Protein Data Bank code 2OBD) (15) embedded into a cubic box containing 66,020 TIP3P (Transferable Intermolecular Potential 3 Point) water molecules and 15 Na ϩ atoms (for neutralizing the CETP surface charges). The water box boundary was at least 25 Å away from the CETP surface. The missing hydrogen atoms within the CETP crystal structure were recovered by using the AutoPSF module of the Visual Molecular Dynamics (VMD) software package (from the University of Illinois at Champaign-Urbana) (31). The system, containing a total of 206,041 atoms, was subjected to energy minimization via a total of 20,000 steps to remove the atomic clashes. Within the first 10,000 steps, the protein backbone atoms were fixed. In the next 10,000 steps, the protein backbone atoms were constrained under a force constant of 5 kcal/ mol/Å 2 . The energy-minimized system was subsequently heated from 0 to 310 K over 31 ps to initiate the all-atom MD simulations. The system at 310 K and 1 atm reached equilibration after 5.4 ns; the constraints on the protein were removed after the first 0.4-ns simulation.
Assembling the Simulation System-The assembly process was conducted through the following three steps: (i) sandwiching the above equilibrated CE pool between two equilibrated POPC monolayers before further equilibration, (ii) sandwiching the above equilibrated TG pool with two monolayers, and (iii) allowing the two CETP ␤-barrel domain distal ends to penetrate into each of the above phospholipid sandwiches. The depths of penetration were ϳ3-4 nm according to transmission EM measurements (16) in which the molecules clashing with CETP were removed. The simulation system, consisting of a total of ϳ330,000 atoms, was then obtained after filling the gap between the two phospholipid sandwiches with equilibrated water (with 15 Na ϩ atoms for neutralization).
Equilibration of the Simulation System-The equilibration of the above system was achieved via three steps. (i) The overlapping molecules were manually shifted away from CETP, and then the system was subsequently subjected to 10,000 steps of energy minimization and a 0.2-ns MD simulation. This 0.2-ns simulation included a step of heating from 0 to 310 K within 62 ps and a simulation step at 310 K and 1 atm and a constraint of the heavy atoms (non-hydrogen atoms) of both water and CETP. (ii) The constraints on the water heavy atoms (oxygen) were gradually released via 10,000 steps of energy minimization followed by a 0.2-ns MD simulation. Within this step, an intermediate stage of constraint was used in which the water oxygen atoms were allowed to transfer only in parallel to the lipid monolayers under an initial force constant of 10 kcal/mol/Å 2 . (iii) The constraints on the CETP heavy atoms were gradually released via 15,000 steps of energy minimization and a 0.4-ns MD simulation. In this step, three intermediate stages were used to prevent deformations in the protein structure, i.e. releasing the side chain constraints, releasing non-C␣ atom constraints, and gradually releasing the constraints on the C␣ atoms (10 kcal/mol/Å 2 ). After releasing all constraints, the system was subjected to 40-ns MD simulations at 310 K and 1 atm with NAMD2 software to achieve equilibration (27). The MD simulation used to achieve this equilibrated system was repeated three times.
Equilibration Analysis of the Simulation System-The equilibration analyses were conducted using VMD (31) by monitoring the changes in the following parameters. (i) The root mean square deviation (RMSD) was calculated by the spatial changes in CETP C␣ atoms relative to their initial positions. (ii) The size of CETP was computed from the radius of gyration of the C␣ atoms. (iii) The molecular volume was measured by using a grid size of 0.25 Å. (iv) The volumes of the CE and TG pools were estimated on the basis of the distance between the attached POPC monolayers. (v) The average distance among POPC molecules was determined from the first peak of the phosphorus radial distribution function.
The system was suggested to be equilibrated after 20 ns on the basis of the convergence analyses of the CETP and lipid structures. The CETP convergence was indicated by the RMSD (a plateau of ϳ2 Å), the radius of gyration of the CETP C␣ atoms, and the CETP molecular volume. The CETP structure with the lowest C␣ RMSD (compared with the averaged structure within the last 20 ns of the simulation) was used to analyze the internal cavities and pore positions, which were calculated with the fpocket program with a minimum ␣ sphere radius of 3 Å (32). The small cavities and pores containing fewer than 15 ␣ spheres were discarded. For comparison, the cavities in the crystal structure were also identified through the same procedure. Lipid convergence was suggested on the basis of the volumes of the CE and TG pools and the average distances of the POPC molecules in each monolayer. No systematic drift was observed during the last 20 ns of the simulation.
Determination of the CE Transfer Pathway-Given that the CE transfer from HDL to LDL is on an approximately second time scale, a driving force was used to detect the transfer pathway within a practical time period. A CE located near but not directly contacting the distal end of the N-terminal ␤-barrel domain (the minimum distance between CE and CETP was more than 2.4 Å, the diameter of the hydrogen van der Waals surface) was selected as a representative molecule to probe the CE transfer pathway within the equilibrated CETP after the CEs and phospholipids it contained were removed. A step-by-step method for pathway searching was conducted as described below. When the selected CE molecule was pulled toward the center of a "pore" located in the distal end of the CETP N-terminal ␤-barrel domain under an example force of 8 kcal/mol/Å (applied onto the CE acyl chain end), the pore became larger (ϳ7 Å in diameter) and deeper (calculated by the MOLE 2.0 program (33)). The center of the deeper pore was then used as a new target to guide the driving force required to pull CE even further into CETP. After sequentially repeating the above procedures, the CE molecule was able to completely transfer through the entire CETP and exit (where the minimum distance between CE and CETP was more than 2.4 Å) from the C-terminal distal end. This process showed that CE can be transferred through a CETP molecule.
Transfer Time under Different Driving Forces-To test how different driving forces influence the CE transfer process, 18 different forces (6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22, and 23 kcal/mol/Å) were applied to the same CE molecule. Under each driving force, the MD simulations were repeated four times. Therefore, a total of 72 simulations were performed in this test. These simulations yielded a total of 72 different transfer times. The relationship between the driving force and transfer time was analyzed with MATLAB.
Estimation of the CE Transfer Time through Radiolabeling Experiments-To the best of our knowledge, there are no experimental data on the CE transfer time. However, the CE transfer time can be estimated from the experimental data of the CE transfer rates. Radioactive experiments have shown that the CE transfer rate with CETP is ϳ75.7 Ϯ 1.5 nmol/h/g (34) and that the CE transfer rate with plasma is ϳ64 Ϯ 5 g/h/ml (35), which corresponds to ϳ1.14 -1.54 CE molecules/s/CETP. This calculation is based on a plasma CETP concentration of ϳ1.75 g/ml (34), a CETP molecular mass of ϳ73 kDa, and a CE molecular mass of ϳ651 Da. Thus, the transfer time is ϳ0.65-0.88 s/CE molecule/CETP.
Analyses of the Hydrophobicity and Diameter of the CE Transfer Pathway-The hydrophobicity of the CE transfer pathway was calculated by using a weighted averaging method and the hydrophilicities of the contact residues. The contact residues were defined as those CETP residues within a distance of less than 2.4 Å from the transferred CE at each particular pathway position. The weight of a contact residue was calculated as the ratio of the solvent-accessible surface area of the residue to the total solvent-accessible surface area of all contact residues. The hydrophilicity of a contacting residue was defined by the Kyte and Doolittle (36) scale. The diameter of the CE transfer pathway was measured with MOLE 2.0 (33) with a 3.0-Å probe radius and a 1.25-Å interior threshold.
Analyses of the Interaction Energies of CE-CETP and CE-Lipid along the CE Transfer Pathway-The interaction energy between the transferred CE and CETP at each particular point along the pathway was calculated by using the pair interaction function of NAMD2. The periodic boundary conditions, cutoff distance of 12 Å for van der Waals interactions, and particle mesh Ewald method for calculating long range electrostatic interactions were used for this calculation. The interaction energy between the transferred CE and all other lipids at each particular point along the pathway was calculated by the same method.

Results and Discussion
Simulation System Setup-To simulate CE transfer from HDL to LDL via CETP, a simulated system was established on the basis of the results of transmission EM experiments (16,17). Transmission EM experiments have shown that the CETP N-terminal ␤-barrel domain penetrates into the HDL surface via a hydrophobic interaction (16,17), whereas the C-terminal ␤-barrel domain penetrates into the LDL/VLDL surface; the depths of penetration are ϳ3-4 nm (16). A representative transmission EM image in which a CETP penetrates into both HDL and LDL is shown in Fig. 1A. To simulate this CETP interaction with HDL and LDL (Fig. 1B), (i) the N-terminal ␤-barrel domain of an equilibrated crystal structure of CETP (15) was allowed to penetrate ϳ35 Å into an equilibrated POPC monolayer adhering to an equilibrated CE pool. (ii) The C-terminal ␤-barrel domain was allowed to penetrate ϳ32 Å into an opposite POPC monolayer adjacent to a TG pool (Fig. 1, C and D).
(iii) The middle portion of CETP was surrounded by equili-brated and neutralized water molecules (Fig. 1D). This system contains a total of ϳ330,000 atoms.
Conformational Changes of CETP-The above system was subjected to all-atom MD simulations for 40 ns to reach an equilibrium state under physiological conditions (Fig. 2, A and  B). After the simulation was repeated three times, a root mean square fluctuation (RMSF) analysis of the C␣ atoms showed that CETP presented higher stability when it was embedded in the lipid monolayers (ϳ4% C␣ atoms showed an RMSF greater than 1.4 Å) than when it was in solution (ϳ11% above 1.4 Å) (18). However, the distal ends of the ␤-barrel domains, helix X (Glu 465 -Ser 476 ), and the loop connecting the two ␤-barrel domains (Asp 240 -Arg 259 ) remained relatively flexible compared with the crystal structure (Fig. 2, C and D). Both distal ends exhibited distinguishable conformational changes (Fig.  3A). In the N-terminal distal end, loop ⍀6 (Gln 155 -Trp 162 ) shifted by a maximum distance of ϳ6.2 Å (on the basis of the C␣ shift) at residues Gln 155 and Gly 156 , which directly enlarged internal cavity C1, consequently connected it to the nearby distal end pore, and indirectly connected cavities C2 and C3 (Fig. 3, B and C). These changes are similar to results from a recent MD simulation (19) where the ⍀5 and ⍀6 loops separate upon penetrating into HDL and open the N-terminal distal end to take up a CE molecule. This conformational change may be a result of the hydrophobic interaction between the CETP N-terminal distal end and the hydrophobic core (17). In the C-terminal distal end, a portion of helix BЈ (Ile 405 -Met 412 ) and the flexible regions ⍀1 (Ser 286 -Thr 322 ), ⍀2 (Lys 347 -Thr 362 ), and ⍀3 (Lys 392 -Ser 404 ) bend along the ␤-barrel domain center away from the longitudinal axis of CETP with a maximum angle of ϳ6° (Fig. 3A, right panel). This bending directly generates a pore (P2) in this distal end (Fig. 3C), indirectly enlarges the central cavity in the central ␤-sheet domain, subsequently connects the central cavity to pore P2, and ultimately forms a large tunnel linked to the TG pool (Fig. 3C).
CE Transfer Strategy-To test whether a CE molecule could transfer through an entire CETP under conditions as close as possible to the real experimental conditions, an all-atom simulation method, i.e. a steered molecular dynamics simulation, was used. Because the time required to transfer a CE molecule from HDL to LDL is on an approximately second time scale (described under "Experimental Procedures"), all-atom simulations on such a large simulation system for nearly a second are quite challenging for the current advanced supercomputers. Although some coarse grain and implicit model methods can speed up the simulation by simplifying a group of atoms as one FIGURE 1. Simulation system for studying CE transfer. A, a representative negative staining EM image shows that CETP bridged HDL and LDL, forming a ternary complex (shown in schematic in B). C, a simplified simulation system (shown in schematic) was used to simulate the ternary complex to elucidate CE transfer from HDL to LDL at an atomic level. D, the N-terminal ␤-barrel domain of the CETP was inserted ϳ35 Å deep into a POPC monolayer adhered to a CE pool, whereas the C-terminal ␤-barrel domain penetrated ϳ32 Å inside an opposing POPC monolayer attached to a TG pool. The region between the two opposing lipid monolayers was filled with water molecules. The POPC headgroups and fatty tails are colored yellow and green, respectively, and the CE, TG, CETP, and water molecules are colored pink, gray, blue, and orange, respectively. The CE molecule is highlighted using van der Waals spheres. united bead or the solvent/lipid as a continuous medium and are highly successful in simulating protein dynamics, modeling, design, and prediction (37,38), this simplification may also reduce the accuracy of the simulation and distort the protein structure. For example, an implicit model for both the solvent and lipids has been shown to fail to maintain the structure of the transmembrane helix dimers of ErbB1/B2 and EphA1, potentially causing kinking, bending, or even twisting within the structure of these helices (39). Therefore, an all-atom simulation is still required to simulate the detailed process of CE trans-  fer. However, we must find a practical solution tailored to the limitations of current computational capabilities regarding simulation time.
Our intermediate solution was to introduce a driving force on the transferred CE molecule to speed up the CE transfer process on the basis of the following reasons. (i) The transfer process involves no chemical reaction or external energy. (ii) The chemical properties of HDL and LDL are similar (40,41). (iii) CE transfer from HDL to LDL is likely to be a kinetic process (13). We believe that the differences in the physical properties of HDL and LDL, such as the inner pressures and CE concentrations, play a key role in generating a driving force for CE transfer. Experimental measurements have shown that HDL and LDL have a similar surface tension (42)(43)(44). According to the Young-Laplace equation, the smaller the particle, such as HDL, the higher the internal pressure over a larger particle, such as LDL, under the same surface tension. This difference in internal pressure would generate a force to drive CE transfer.
To avoid a potential error induced by using a specific force, a series of driving forces were also applied to the transferred CE. By statistically analyzing the simulation results under different driving forces, we may gain a new understanding of the detailed transfer process and mechanism.
Probing the CE Transfer Pathway-Under a representative force of 11 kcal/mol/Å, a CE molecule was pulled to the surface pore at the CETP N-terminal distal end (prior to this simulation, two phospholipids and two CE molecules inside this CETP were removed). When CE penetrated the pore (Fig. 4A), a deeper pore was generated via the rearrangement of the nearby side chains of CETP. Through the continued pulling of CE deeper into the pore (Fig. 4B), additional conformational changes were induced in CETP (Fig. 4C), i.e. a roughly counterclockwise rotation of the N-terminal ␤-barrel domain, as proposed previously (16). This rotation connected the pore to the nearby internal cavities C2 and C3 and formed a connection to the central cavity within the central ␤-sheet domain. By sequentially pulling the CE through CETP, this molecule eventually exited from the C-terminal distal pore after a total simulation time of ϳ0.5 ns (Fig. 4A; before and after transfer, there is no physical contact between CE and CETP). This process was repeated four times (Figs. 4E and 5A), and each time, we confirmed that a CE molecule could be transferred through the entire CETP molecule (Figs. 4A and 5A and supplemental Movie 1).
Although the CE orientations during the four repeated simulations (indicated by the direction of the C10 -C19 bond shown in Fig. 5A) were flexible before entering the CETP and after exiting from the CETP, the CE steroid ring exhibited relatively constrained orientations that were parallel to each other during CE transfer within the N-and C-terminal ␤-barrel domains (Fig. 5, A and B). Notably, CE made an ϳ90°turn near the central ␤-sheet, a result consistent with the observed ϳ90°d ifference in orientations between the two CEs within the crystal structure (15).  (Fig. 5C). Statistical analysis of the relationship between the driving force and transfer time showed that the negative log of the transfer time presented a linear relationship with the log of the driving force with an R-factor of 0.96 and a slope of 2.75 (Fig. 5D).

Analysis of CE Transfer
From the above relationship between the driving force and transfer time, we predicted that the CE transfer time should be ϳ0.008 -0.03 s/CE molecule (i.e. ϳ33-125 CE molecules/s/ CETP) under an estimated driving force of ϳ0.018 -0.029 kcal/ mol/Å. The exact experimental driving force on the transferred CE has not been measured. However, we were able to estimate the driving force by using the method described below. The driving force should be calculated on the basis of the different chemical potentials between HDL and LDL. The chemical potential is defined as U ϩ PV Ϫ TS where the U is the internal energy (related to the chemical reaction and phase transition); P and V are the pressure and volume, respectively; and T and S are the temperature and entropy, respectively. Given that the chemical components within HDL and LDL are similar (neutral lipid core, i.e. CEs and TGs, surrounded by a phospholipid monolayer and amphipathic apolipoproteins) and that there is no relevant chemical reaction, temperature change, phase transition, or net CE volume change during CE transfer, the differences in the pressures and entropies would dominate the difference in chemical potential. HDL should have a higher chemical potential than LDL because of its higher internal pressure (discussed above) and lower entropy (i.e. higher CE concentration).
From the internal pressure differences and the dimensions of the CETP tunnel (ϳ6 Å measured using the MOLE 2.0 program, which is similar to that of the steroid ring of CE, i.e. ϳ6 Å in width and ϳ4 Å in thickness (15)), we calculated that the driving force is ϳ0.018 -0.029 kcal/mol/Å. The inner pressure difference was calculated by using the Young-Laplace equation (P ϭ 2␣/R where ␣ is the surface tension and R is the radius of the particle) and the following experimental results. (i) The surface tensions of HDL and LDL are both ϳ0.020 -0.033 newton/m according to previous surface balance and surface radioactivity experiments (42)(43)(44). (ii) The average diameters of HDL and LDL are ϳ100 (45) and ϳ220 Å (46), respectively.
On the basis of the above driving force, the transfer times were predicted to be ϳ0.008 -0.03 s/CE molecule (i.e. ϳ33-125 CE molecules/s/CETP). This transfer time is shorter than the value that was calculated from the radiolabeled CE experiments with the plasma CETP concentration (34,35), i.e. ϳ0.65-0.88 s/CE molecule or ϳ1.14 -1.54 CE molecules/s/ CETP molecule. Because the time spent for CETP travel, sensing, interacting, and penetrating both HDL and LDL was not included in the predicted transfer time from our simulations, a longer time should be expected for CE transfer under physiological conditions.
Interaction Energy between CE and CETP-Calculation of the free energy is an ideal method to analyze the system energy. This free energy calculation, particularly the portion attributed to entropy (47), is practically impossible in our simulated system containing ϳ34,000 molecules with multiple dimensions of freedom (such as the internal bending and rotational degrees of freedom of the transferred CE molecule) because of the limitation of computational power. Although a fast and easy estimation of free energy can be achieved by using an implicit solvent model, the model may also distort the protein structure as described previously (39). Therefore, we chose the most reliable all-atom simulation result to calculate two specific interaction energies (i.e. CE-CETP and CE-lipid) as a part of our energy analyses.
Under a representative driving force of 11 kcal/mol/Å, the all-atom simulation results were used to compute the interaction energy between the transferred CE and CETP against the CE transfer pathway (Fig. 6, A and B, blue lines). The energy distribution showed that CE exhibited a higher CE-CETP interaction energy, which quickly decreased when CE penetrated into the CETP tunnel. In the middle of the tunnel, this interaction energy was relatively flat except for a relatively high energy barrier that was observed near the central ␤-sheet where the CE steroid ring turned ϳ90°to pass through this portion of the tunnel (Fig. 5, A and B).
Interaction Energy between CE and Lipids-Similarly, we used the same all-atom simulation result to analyze the interaction energy between the transferred CE and all other lipids (Fig. 6, A and B, orange lines). The energy distribution throughout the transfer pathway showed that CE exhibited a lower energy when it was close to either side of the neutral lipid pools and exhibited a higher energy when it was located between the neutral lipid pools.
The combined energy of the CE-CETP and CE-lipid showed a much more flat distribution than did either the CE-CETP or CE-lipid interaction. Notably, the combination energies were similarly low at the positions before CE penetration or after CE exit from the CETP tunnel, thus suggesting that CE transfer required nearly no initial potential energy. When CE was transferring into the tunnel, the combination energy was dominated by the CE-CETP interactions in which the highest interaction energy was around the central tunnel region. These combination energy distribution profiles were similar to one another under different driving forces. The central energy barrier was maintained at the same position even after the driving force was decreased to 6 from 23 kcal/mol/Å (Fig. 6A), which corresponded to a ϳ50-fold increase in simulation time. The similarities in energy profiles and barrier locations reflected that our simulations may reveal intrinsic information about CE transfer through the CETP tunnel.
Residues Associated with the Energy Barrier and the Energy Wells-Although the combined energy distribution was relatively low and flat (Fig. 6, A and B, black line) compared with either the CE-lipid or CE-CETP energy, a relatively high energy barrier was observed when CE was located near the central tunnel region. To determine the residues that might contribute to this energy barrier, we conducted the following analyses based on hydrophobicity and geometrics.
The hydrophobicity analysis showed that the entire tunnel/ pathway is relatively hydrophobic (Fig. 6C), thus providing the condition required for CE transfer through CETP. The geometric analysis showed that the tunnel diameter near the central region was only ϳ6.0 Ϯ 0.6 Å. This narrow region (called the "neck" of the tunnel in the crystal structure (15)) is similar to or slightly larger than the width of a CE molecule (ϳ6 Å wide and ϳ4 Å thick (15)), thus suggesting that it may be possible to transfer an entire CE molecule through this region (Fig. 6D). We found that the residues surrounding this narrow region (within ϳ10-Å distance range), including Ile 15 , Leu 23 , Ala 202 , Ile 205 , Leu 206 , Phe 263 , Phe 265 , and Met 433 (Fig. 6D, central  panel), participated in physical contacts with the CE molecule (distance within 2.4 Å, the diameter of a hydrogen van der Waals surface), which might contribute to this local high energy barrier.
In addition to the above energy barrier, a low energy well was observed within both the N-terminal and C-terminal ␤-barrel domains (Fig. 6B). The Phe 115 , Arg 158 , and Phe 167 residues contributed the greatest energy to the N-terminal energy well (Fig.  6D, left panel), which may attract CE molecules from the CE pool into the CETP tunnel. The Phe 301 and Met 412 residues contributed the greatest energy to the C-terminal energy well (Fig. 6D, right panel), which may orient or facilitate the rotation of the CE steroid ring to allow exit from the C-terminal pore (Fig. 5A). Because the residues identified above can form physical contacts or strong interactions with the transferred CE, mutation of these residues may change their contacts or interactions with CE, thus affecting the CE transfer rate. Consistent with this hypothesis, experiments have shown that mutating the Ile 15 and Met 433 at the CETP central narrow region to the bulky residue tryptophan (i.e. I15W and M433W) and the F301D mutation at the C-terminal distal end decreases CETP transfer activities by ϳ60, ϳ40, and ϳ20%, respectively (15). For the other identified residues, there have been no reported results from mutational experiments to the best of our knowledge. Additional mutational analyses of these residues are needed to further verify the tunnel transfer mechanism. Notably, the mutations of several residues, such as L468Q (48), T138Y, and L457W (15), have shown effects on CE transfer. However, those residues were not identified in our energy and physical interaction analyses. The failure to identify these residues may be because only the mutated residues may form direct interactions with the transferred CE, or their mutation may generate a secondary effect on other nearby residues that interact with the transferred CE.
CE Transfer Mechanisms-Our previous EM studies have demonstrated the existence of a ternary complex of HDL-CETP-LDL/VLDL (16), suggesting a tunnel mechanism for CETP. The simulations indicated that it is possible for a hydrophobic CE molecule to transfer through an entire CETP via a continuous hydrophobic tunnel between the two distal ends of CETP (Fig. 6E). This tunnel might facilitate CE uptake into the tunnel without using energy. Upon entering the tunnel, hydrophobic interactions among the CE fatty tails and residues Phe 115 and Phe 167 may mediate the penetration of the CE molecule further into the tunnel. During subsequent penetration of the hydrophobic CE molecule through this hydrophobic tunnel, a narrow region might act as a rate-limiting site, leading to a prominent energy barrier for CE transfer. The internal pressure difference between HDL and LDL might provide the additional energy required for CE to pass through this energy barrier. Upon passing through this region, no further major energy barrier would prevent the exit of CE from the C-terminal opening. Residues Phe 301 and Met 412 may even assist the exit of CE by facilitating the rotation of the CE steroid ring. This observed CE flow from HDL to LDL requires no external energy input, thus suggesting that CE may transfer through the CETP tunnel by diffusion.
The prominent interaction energy barrier for CE transfer through CETP is ϳ21 kcal/mol (Fig. 6B) (estimated under a driving force of 11 kcal/mol/Å and transfer time of ϳ0.5 ns and without including entropy). To understand this barrier, we estimated the energy barrier for CE transfer into an aqueous phase through a lipid monolayer. It has been reported that the peak free energy barrier is ϳ140 kJ/mol (i.e. ϳ33.5 kcal/mol) for transfer of a cholesteryl oleate through a 1,2-dimyristoyl-snglycero-3-phosphocholine lipid monolayer into an aqueous phase (49,50) in which the interaction energy (without including the entropy) has been estimated to be ϳ28.6 kcal/mol. This energy barrier is higher than our energy barrier. Moreover, because the transfer time (under a driving force of 11 kcal/ mol/Å) is 10 9 -10 10 times faster than the transfer time mea-sured experimentally (approximately seconds) (34,35), the energy barrier under the experimental transfer time should be even much lower than ϳ21 kcal/mol, thus suggesting that it would be even easier for the CE molecule to transfer through the CETP molecule than through the lipid monolayer into an aqueous phase. Although the tunnel model was favored by our test, some other experiments have not favored this model. For example, by incubating small peptides of the CETP C terminus with HDL (49), small micelle-like structures have been generated through secondary structure disorder-to-order transitions of these small peptides. The authors of this previous study have proposed a small micelle-like structure transfer model in which those micelle-like structures mediated CE transfer from HDL to LDL/VLDL. In this proposed model, the authors have not explained how those micelle-like particles, which contain no protein, could function to sense, target, and bind to LDL/VLDL (instead of HDL, other plasma vesicles, and cells in plasma) for directional CE transfer. Moreover, these micelle-like structures have not been observed after incubation of full-length CETP with HDL (16,17). A recent, independent EM study has concluded that CE transfer does not require a ternary tunnel complex with CETP according to negative evidence, i.e. the lack of an observation showing the existence of the CETP-LDL complex or the HDL-CETP-LDL complex (51). The weakness of the logic in this study is that this failure to observe the complex has been used as evidence of its non-existence instead of being used as non-supporting evidence.
Although our simulations support the possibility of CE transfer through the CETP tunnel, a study directly observing the CE molecule transfer process or more convincing experiments are still necessary to validate this mechanism. An indirect validation study could be designed by using new CETP inhibitors to target the CETP hydrophobic tunnel, particularly the narrow central region and distal end pore regions, to introduce energy barriers for CE passage. Our insights into the CE transfer pathway through CETP may provide important clues for the design of new CETP inhibitors to efficiently treat cardiovascular diseases.
Author Contributions-G. R. initiated and designed this project, and D. L., S. Z., X. Z., M. R., and L. Z. modified this project. D. L. performed the modeling, MD simulations, and data analyses. D. L., M. R., and G. R. interpreted the data and drafted the initial manuscript, and D. L., X. Z., M. R., L. Z., G. R., and S. Z. revised the manuscript.