|Home | About | Journals | Submit | Contact Us | Français|
The product of transesterification of phospholipid acyl chains and unesterified cholesterol (UC) by the enzyme lecithin: cholesterol acyltransferase (LCAT) is cholesteryl ester (CE). Activation of LCAT by apolipoprotein (apo) A-I on nascent (discoidal) high density lipoproteins (HDL) is essential for formation of mature (spheroidal) HDL during the antiatherogenic process of reverse cholesterol transport. Here we report all-atom and coarse grained (CG) molecular dynamics (MD) simulations of HDL particles that have major implications for mechanisms of LCAT activation. Both the all-atom and CG simulations provide support for a model in which the helix 5/5 domains of apoA-I create an amphipathic “presentation tunnel” that exposes methyl ends of acyl chains at the bilayer center to solvent. Further, CG simulations show that UC also becomes inserted with high efficiency into the amphipathic presentation tunnel with its hydroxyl moiety (UC-OH) exposed to solvent; these results are consistent with trajectory analyses of the all-atom simulations showing that UC is being concentrated in the vicinity of the presentation tunnel. Finally, consistent with known product inhibition of CE-rich HDL by CE, CG simulations of CE-rich spheroidal HDL indicate partial blockage of the amphipathic presentation tunnel by CE. These results lead us to propose the following working hypothesis: After attachment of LCAT to discoidal HDL, the helix 5/5 domains in apoA-I form amphipathic presentation tunnels for migration of hydrophobic acyl chains and amphipathic UC from the bilayer to the phospholipase A2-like and esterification active sites of LCAT, respectively. This hypothesis is currently being tested by site-directed mutagenesis.
Of the possible mechanisms suggested to explain the atheroprotective role of apolipoprotein (apo) A-I, a process called reverse cholesterol transport is most completely understood at the molecular level. Activation of the plasma enzyme lecithin:cholesterol acyltransferase (LCAT) (1) by apoA-I is necessary for esterification of the unesterified cholesterol (UC) molecules of high density lipoproteins (HDL) to cholesteryl ester (CE) and leads to the conversion of discoidal phospholipid-rich HDL to spheroidal, CE-rich (circulating) HDL, a central step of reverse cholesterol transport (2). A fuller understanding of reverse cholesterol transport demands knowledge of the structure and dynamics of the various HDL particles and intermediates in their assembly.
In 1999 we proposed an atomic resolution antiparallel double belt helical model for discoidal HDL (3). The general features of the model, including the LL5/5 registry in which the left interfacial edge of helix 5 associates pairwise to form an antiparallel interhelical domain (3), have been confirmed by several laboratories using physical chemical methods (4–11). In the model, apoA-I monomers form a curved amphipathic α helical ring with 11/3 residues per turn—termed an α 11/3 helix and since observed in other lipid-associating proteins, such as α–synuclein (12, 13)—so that the hydrophobic surface faced inward toward the lipid disc. We have suggested that the two antiparallel helical rings are held together by six interhelical solvent-shielded salt bridges (3). Although initial derivation of this model depended critically upon certain features of the lipid-free x-ray structure (14), the model can be derived a priori from profound constraints imposed on the conformation and orientation of lipid-associated proteins by lipid (3).
In another publication, we successfully modeled the conformation of apoA-I on phospholipid bilayer discs as an amphipathic helical double belt extending the full length of the lipid-associating domain (residues 44–243) of apoA-I with N- and C-termini in direct contact (15). Conversion to smaller discs would require a conformational change in the apoA-I belt. To determine the nature of this change, we performed molecular dynamics (MD) simulations on a series of progressively smaller discoidal HDL particles produced by incremental removal of palmitoyloleoylphosphatidylcholine (POPC) (16). These results provided atomic resolution models for two of the smaller particles produced by in synthetico reconstitution of PL-rich HDL particles from POPC and apoA-I (16). Most importantly, the apoA-I amphipathic α helical double belt, held together by interhelical solvent-shielded salt bridges, independent of four different conditions of particle shrinkage, twisted into a saddle-shape to conform to the edge of the bilayer by closely approximating the X-ray structure of the close dimer pairs of the tetrameric lipid-free apoA-I (14).
However, these MD simulations have limitations: Multiple long simulations are required to increase the confidence that equilibrium has been achieved and that energy barriers have been overcome. Because of the short (less than 10 ns) simulations, we could not rule out the possibility that the saddle-shaped structures that we saw represented kinetically trapped collapsed intermediates.
We recently used temperature jump MD simulations to 500 K to scan the per-residue helix stability of apoA-I in phospholipid-rich HDL (17). We found that the conformations of the overlapping N- and C-terminal domains of apoA-I in the particles were unusually mobile, exposing hydrocarbon regions of the phospholipid to solvent; a lack of buried interhelical salt bridges in the terminal domains correlated with increased mobility. On the basis of these results, we proposed that the terminal domains of apoA-I are phospholipid concentration-sensitive molecular triggers for fusion/remodeling of HDL particles.
One approach for overcoming kinetic trapping that incorporates temperature jump MD simulations is MD simulated annealing (MDSA) (18–21). This procedure, when used to refine x-ray crystallography and NMR structures, after placing constraints on the initial protein structure in implicit solvent, applies temperature jumps (to 500–1000 K or so), followed by slow cooling to physiologic temperatures. Since we wanted to overcome energy barriers during the simulations, we applied MDSA without constraints to our particles in explicit solvent as an approach to bypassing kinetically trapped intermediates that might exist in our previous simulations. We hypothesized that the lipid in the HDL assemblies would place its own restraints on apoA-I: i) the low dielectric environment of a lipid bilayer greatly increases long range electrostatic effects; ii) the two dimensional lipid bilayer has significantly fewer degrees of freedom compared to globular proteins, simplifying the protein folding problem; iii) nonspecific protein-lipid hydrophobic interactions dominate, as opposed to much more specific protein-protein interactions. As there is no standard protocol for unconstrained MDSA, we developed optimal conditions for unconstrained MDSA (degree and length of temperature jump, and rate of cooling) through trial and error. This approach provides exploration of conformational space but is not a robust way to explore longer simulations times.
A second approach for overcoming kinetic trapping is coarse grained MD (CGMD), a form of simulation that has become an effective tool for simulating lipid dynamics over relatively long time frames (µsec to msec) (22). The reduction of the number of degrees of freedom, achieved by mapping four or even more atoms into one interaction center (CG bead) and by disregarding all hydrogen atoms, and the combined use of short-range interactions and smooth potentials to produce large integration steps, make the CG model computationally very fast.
More recently, CGMD has been applied to proteins and lipoproteins (23, 24). The type of the backbone particle depends on its secondary structure. Since different dihedral and angle parameters are used to distinguish α helixes, β strands, or random coils, it is difficult, but not impossible, to study realistic folding-unfolding events at this stage. Our approach to this weakness in the Martini force field for proteins is to use MDSA to explore conformational space, applying CGMD to each member of the resulting ensemble (24). This allows exploration of both conformational space and time.
Human LCAT transacylates the sn-2 acyl chains of POPC to UC with approximately 90% fidelity (25, 26), forming CE and lysolecithin; phosphatidylcholine (PC) molecules with longer sn-2 chains relative to sn-1 have lower fidelity (25, 26). In the first step of this reaction, LCAT has a phospholipase A2-like activity, hydrolyzing an acyl chain of PC to fatty acid (27). For soluble, shorter chain PC, the phospholipase activity occurs spontaneously; apoA-I is not required (28). For insoluble, longer chain PC, the phospholipase A2-like activity requires apoA-I (28). One function, therefore, of apoA-I is activation of LCAT by presentation of acyl chains of long chain PC to the active hydrolytic site—the rate limiting step during trans-esterification of CE by the enzyme—for binding and catalysis (29).
We previously showed that the spacing between the individual helixes of the helix 5/5 domain on discoidal HDL surfaces is dynamic (30) and suggested that the opening and closing of this domain may modulate the access of lipid substrate by LCAT (15). Martin, et. al. (10) have suggested a similar structure for the helix 5/5 domain that they term the “looped belt” model in which residues 133–146 comprise a flexible loop segment and they propose that this domain could serve to modulate LCAT activity for different-sized discoidal HDL particles. Using the combined MDSA/CGMD strategy employed in our previous simulations of spheroidal HDL (24), we report here that this combined approach provides complimentary structural and dynamic support for a detailed molecular model for LCAT activation by apolipoprotein A-I. We propose that, after attachment of LCAT to discoidal HDL, the helix 5/5 domains in apolipoprotein A-I form amphipathic (presentation) tunnels for migration of hydrophobic acyl chains and amphipathic unesterified cholesterol from the discoidal HDL bilayer to the phospholipase A2-like and esterification active sites of LCAT, respectively.
One particle containing 160 POPC, 24 unesterified cholesterols, and 2 apoA-I molecules (160:24:2) was created by generating an all-atom particle containing 2 full length apoA-I molecules in a 100%α11/3 conformation wrapped in an antiparallel double belt around a POPC bilayer of 252 lipid molecules (15). The first 32 residues of the N-terminus of apoA-I (helixes G0, G1 and G2) were rotated to orient their hydrophobic faces towards the surface of the POPC acyl chains. The resulting particle was subjected to random removal of 92 POPC molecules, followed by random insertion of 24 UC molecules into the space vacated by the deleted POPC to generate a particle with a POPC:UC:apoA-I molar ratio of 160:24:2. This initial particle was simulated at 310 K and 1 bar for 10 ns, and then subjected four times to the following MDSA protocol (all at 1 bar): heated from 310 K to 500 K in 20 ps, simulated at 500 K for 10 ns, cooled from 500 K to 420 K in 2 ns, cooled from 420 K to 400 K in 5 ns, cooled from 400 K to 310 K in 3 ns, and finally simulated at 310 K for 10 ns, giving a total duration of 30 ns. For simulations at a fixed temperature, velocity reassignments occurred every 1 ns to prevent the “flying ice cube” effect (31).
All-atom simulations were performed using NAMD (32) as described in Jones et al. (16). Each system was ionized and charge neutralized with NaCl to 0.15 M with the “Add Ions” plug-in of VMD (33). The TIP3P water model was used (34). The CHARMM 22 (35, 36) and 27 (37, 38) force fields were used for protein and lipid molecules, respectively.
The initial configuration for a single cholesteryl oleate (CO) molecule was prepared from coordinate files for cholesterol and POPC as described previously (24). Initial CE-rich particles were created by the surface loading method that involves random insertion of UC and CE molecules into spaces vacated by deleted POPC followed by MDSA. A CE-rich spheroidal particle was created by generating an all-atom particle containing 2 full length apoA-I molecules in a 100%α11/3 conformation wrapped in an antiparallel double belt around a POPC bilayer of 252 lipid molecules (3, 15). The resulting PL-rich particle was subjected to random removal of 195 POPC molecules, followed by random insertion of 16 CE and 6 UC molecules into the space vacated by the deleted POPC to generate a particle with a POPC:CE:UC: apoA-I molar ratio of 57:16:6:2. This initial particle was simulated at 310 K and 1 bar for 10 ns, and then subjected four times to the MDSA protocol as described above to produce an ensemble of four particles.
All-atom final structures of MDSA simulations were coarse grained using the CG Tools plugin of the Visual Molecular Dynamics (VMD) software version 1.8.6 (31). CGMD simulations were performed using the MARTINI CG force field for lipid (37) and protein (39). Further details about CG parameters can also be found in our previous work (22).
After mapping the final structures of all-atom MDSA simulated particles with CG beads, each coarse grained discoidal 160:24:2 and spheroidal 57:16:6:2 particle was subjected to conjugate gradient energy minimization. Then, each energy minimized particle was solvated in a cubic periodic water cell extending at least 14 Å beyond the lipid headgroups and the protein molecules. This assured at least the presence of 10 hydration layers of water between two periodic images of each simulated system. The solvated system was then subjected to an additional conjugate gradient energy minimization in order to reduce steric contacts between the water molecules and the lipoprotein complex. In order to preserve overall charge neutrality sodium ions were placed randomly using the genion command of the GROMACS software version 4.0.3 (40). Then, each final energy minimized model system was subjected to a few steps of equilibration without ions. Finally, after 10 ns of equilibration each structure was subjected to long production runs as follows:
CG discoidal (160:24:2) and spheroidal (57:16:6:2) particles were subjected to CGMD simulations of 20 µs and 18 µs at 310 K and 1 atm, respectively. CGMD simulations were performed using GROMACS version 4.0.3 (40). The force field parameters were as described in Catte et al. (15). The majority of protein residues had α-helical conformation and the remaining residues were in a β-turn conformation. All protein residues of the helix 5/5 domain were kept in the α-helical conformation in both model discoidal and spheroidal HDL particles. The final analysis was performed on the last 40% of each CGMD simulation at 310 K. In all CGMD simulations we report the effective time, which corresponds to the time of simulation multiplied by 4 (41).
The RMSDs of protein α-carbon atoms were calculated for the 160:24:2 particles over the entire duration of the four 30 ns MDSA simulations. The RMSDs were calculated with respect to the 160:24:2 particles after they had been initially simulated at 310 K for 10 ns.
The distances of POPC terminal methyls of the sn-2 and sn-1 chains to the nearest residue G129 of the apoA-I molecules were measured and averaged over the entire durations of the four 30 ns MDSA simulations. A simple VMD script was written to measure these distances.
The number of salt bridges formed between any two residues E125, K133, E136, and K140 was counted over the entire duration of the MDSA-1, MDSA-2, and MDSA-4 simulations. The “saltbr” VMD command was used to count these salt bridges.
The distances of POPC terminal methyl carbon atoms and phosphorous atoms from the nearest G129 residue were measured for selected POPC molecules over the entire duration of certain 30 ns MDSA simulations. A simple VMD script was written to measure these distances.
The distances of POPC P1-oxygen atoms to nitrogen atoms on basic protein residues were measured for selected POPC molecules over the entire duration of certain 30 ns MDSA simulations. Since the “saltbr” command of VMD only works for salt bridges between protein residues, a simple VMD script was written to measure these distances.
We first described this function, similar but not identical to a radial distribution function, in Jones, et al (17). For each trajectory frame in the last 20% of the MDSA simulations of the 160:24:2 particles, distances were measured from the nearest protein atom to each POPC or UC (using COM or a representative atom such as phosphorus), and the distribution of these distances were averaged and plotted. A simple VMD script was written to measure and average these distances. For circles, analysis shows that the distribution should be linear in distance from the annulus.
The SASA of all POPC and UC atoms within 7 Å of residue G129 was measured over the course of all MDSA simulations. This was performed by using VMD’s “measure sasa” command.
To emphasize positive versus negative potentials, the APBS plugin in VMD was set to a color scale data range of – 1.0, 1.0 with the structure coloring method set to Volume. To emphasize the regions that are most positive, the APBS plugin in VMD was set to a color scale data range of – 1.0, 20.0 with the structure coloring method set to Volume.
The starting structure equilibrated at 310K was subjected four times to the 30 ns MDSA protocol shown in Fig. 1A. Changes in average RMSD of the protein of the four resulting particles were plotted over the full 30 ns protocol of the simulations and the results shown in Fig. 1B–D. Note that the RMSD of each of the four protein models decreases slightly over the last 6 ns of the protocol, supporting structural convergence.
The four final structures created by the MDSA protocol, viewed from the helix 5/5 domain of the particles, are shown in Fig. 2A. Note that in all four particles there is central separation of the opposing individual helix 5 portions of the double belt motif (green), creating a gap exposing underlying POPC acyl chains (black and white). Fig. 2B shows that in three of the four structures the acyl chain exposure within the helix 5 gap includes one (MDSA-1) or two (MDSA-2 and MDSA-4) terminal methyl groups. In MDSA-3, for reasons to be discussed, the mid portion, rather than the terminal methyl region, of a single acyl chain is exposed to solvent. The average RMSD for helix 5 reaches a maximum and begins to decline before the end of the 500 K temperature jump phase of the MDSA protocol (Fig. S1). This is good evidence that the helix 5/5 structure is not adversely affected by the T-jump step of the MDSA protocol.
In fact, the central domain, the helix 5/5 domain along with the helix 4/6 overlaps, represents the most stable part of the double belt structure, even at 500 K, while the N- and C-terminal domains represent the most mobile and generally least stable regions of the double belt (17). We suggest that central double belt stability is due to two factors: solvent inaccessible salt bridges between the helix 2/3 and helix 7/8 domains K (17) and a complex weave of inter- and intra-helical salt bridges in the helix5/5 domain (Segrest, unpublished data). The complex dynamics of helix 5/5 salt bridges is illustrated by the stochastic formation/breaking of salt bridge shown in Fig. S2.
The structural elements of the helix 5/5 gap can be seen in Fig. 2B. As noted in a previous publication (15), the orientation of the antiparallel double helical belt places the pairwise K133 residues opposite one other in the solvent-shielded salt bridge 2-position of the α11/3 helical wheel (3, 42), so that, in extended conformations, their charged moieties would point directly at one another on the L-interfacial edge of the α 11/3 helix. As seen in Fig. 2B, in the three structures in which terminal methyl groups are exposed to solvent (MDSA-1, MDSA-2 and MDSA-4), the two opposing K133 have moved away from one another to form intrahelical salt bridges with E136 and interhelical salt bridges with E125. Another basic residue, K140, also forms salt bridges with acidic residues E125 and E136 three out of six times. This arrangement of salt bridges is responsible for the helix 5/5 gap that exposes POPC acyl chains. The long sides of the helix 5/5 gap are formed by the acyl chain moieties of the pair of K133 amphipathic basic residues; the short sides of the helix 5/5 gap are formed by the two low bulk G129 residues (Fig. 2B, magenta).
In contrast, in the structure in which no terminal methyl groups are exposed to solvent, MDSA-3, a portion of the upper helix has broken (not shown), the two opposing K133 are inverted to the lipid side of the helix 5/5 domain (Fig. 2C), one interacting with solvent on the R-interfacial edge of the α 11/3 helix and the other forming salt bridges with E125. The sum of these conformational differences results in distortion and narrowing of the helix 5 gap to prevent access to terminal methyl groups; thus MDSA-3 represents an important negative control as will be noted later.
Fig. S2 of the Supporting Information is a plot of a running count of the number of salt bridges, out of a possible eight, created by the two clusters of K133, E125, E136 and K140 during the last 10 ns of the simulations of MDSA-1, MDSA-2 and MDSA-4. These results show that the formation of salt bridges by these four residues at the physiologic temperature of 310 K is surprisingly dynamic; the number of salt bridges varies from 3–8 over 10 ns in a seemingly stochastic fashion.
The structure of the lipid side of the helix 5/5 gap is shown in Fig. 2C. The lipid side is covered by 5 Leu residues per helix (gold). Interestingly, however, the top and bottom of the lipid face of the gap is lined by A130 (white), and by the previously mentioned G129 on the solvent side, providing room for insertion of the terminal methyl moieties. Due to the distortion and narrowing of the helix 5 gap in MDSA-3 noted above, the orientation of the G129 and A130 residues relative to the gap has been changed significantly in this structure.
We used total SASA of all POPC or UC atoms located within 7 Å of G129 as a dynamic metric of the opening or closing of the helix 5/5 gap in the four MDSA simulations and the results are shown in Fig. 3. This approach suggests similar gap metrics for MDSA-2 and MDSA-4; in particular, both display a relatively constant gap area of approximately 25 Å2 over the last 10 ns. The MDSA-1 simulation, although similar to MDSA-2 and MDSA-4, ends with a larger gap of approximately 50 Å2 as the result of a break in one helix 5 between residues 140–143 that occurs at approximately 22 ns. The SASA metric confirms that the gap for MDSA-3 is essentially closed from early in the MDSA protocol; the gap area is generally less than 10 Å2 from 16 ns on.
Having shown that the helix 5/5 gap is consistently present over the last 10 ns of simulation for MDSA-1, MDSA-2 and MDSA-4, we used distance from one or both G129 as a yardstick for the degree in which acyl chain terminal methyls have inserted themselves into the gap. Using this yardstick, we examined the average count per particle of all POPC acyl chain terminal methyls within 4–10 Å of one or both G129 as a function of the time of MDSA. Fig. 4 represents plots of the dynamics of the last 10 ns of those POPC molecules that come within 6 Å of G129 more than 5% of the time in MDSA-1, MDSA-2 and MDSA-4. The terminal methyls of three of the nine POPC remain inserted into the helix 5/5 gap over essentially the entirety of the 10 ns period (one in MDSA-1 and two in MDSA-4) at an average distance of approximately 5 Å. A movie of the insertion and retention of the two terminal methyls of MDSA-4 (thick black and gray curves in Fig. 2C) is shown in Movie S3 in the Supporting Information. In MDSA-2, two terminal methyls swap places midway in the 10 ns period and a third inserts during the last 2 ns. Using a 5 Å mean for full insertion, no more than two terminal methyls are inserted at any given time.
The relative immobility of the three terminal methyls inserted into the helix 5/5 gap is illustrated by comparison with other portions of the same POPC molecules; the contralateral terminal methyls of the eight inserted POPC shown in Fig. 4 are noticeably more mobile than those inserted into the helix 5/5 gap (data not shown).
Inspection of the CGMD simulations of structures created by CG mapping of three of the four members of the 160:24:2 all-atom MDSA ensemble shows a terminal methyl pseudoatom visibly inserted into the helix 5/5 gap in approximately 25% of all trajectory frames inspected at 400 ns intervals; Fig. 5 shows structural details of one such example. Fig. 6A and B quantitatively compare the dynamics of terminal methyl insertion in the all-atom and CG simulations, respectively. Although insertion occurs in both types of simulations, it reaches 100% for a single methyl group for the all-atom simulations using an 8 Å cutoff, while at the approximately comparable distance of 6 Å for the CG simulations (zero insertion occurs at 2 and 4 Å for the all-atom and CG simulations, respectively) single methyl pseudoatom insertion reaches only 32%.
Inspection of the CGMD simulations suggests at least a partial explanation for the discrepancy; UC, in addition to acyl chains, also appears in the helix 5/5 gap. At least one UC-OH is inserted and exposed to solvent in the center of the helix 5/5 gap between the K133 pair in 64% of all trajectory frames inspected at 400 ns intervals. Fig. 7 shows structural details of one such example.
It appears from Fig. 7 that insertion of UC also occurs elsewhere in the double belt model. Inspection of the CGMD simulations shows insertion between one or the other of the two helix 7/3 pair sites (specifically between the R177/D89 solvent inaccessible salt bridges) in approximately 28% of all trajectory frames inspected at 400 ns intervals; however insertion is rare between either one of the two helix 6/4 pair sites (between the R155/E111 solvent inaccessible salt bridge pairs), occurring in less than 1% of all trajectory frames inspected at 400 ns intervals. Fig. 8A–C quantitatively compare the dynamics of UC insertion between the single helix 5/5 and each one of the pair of helix7/3 and helix 6/4 domains using distance measurements of the UC-OH groups from the charged pseudoatoms of K133, R177 and H155 respectively. Using a 6 Å cutoff, UC insertion reaches 53%, 29% and 2% for the single helix 5/5 and each one of the pair of helix 7/3 and helix 6/4 domains, respectively.
To examine the issue of the proximity of UC to the protein during MDSA simulations, we measured the annular distribution of UC and PL from the protein belt in the four 160:24:2 PL-rich particles produced by MDSA and the results are shown in Fig. 9A–C. Annular distribution functions (ADF) of distances of POPC (Fig. 9A) and UC centers of mass (COM) (Fig. 9B) or POPC phosphorus atoms (Fig. 9C) from the nearest protein atom were measured over the last 20% of the MDSA simulations of the 160:24:2 particles. A comparison of Fig. 9A with Fig. 9B shows that both POPC and UC have a pronounced annular distribution shell centered approximately 5 Å from the protein. The annular shell is noticeably sharper for UC (Fig. 9B) than for POPC (Fig. 9A), partly reflecting the more condensed nature of UC versus POPC, but surprising, nevertheless.
A plot of the ADF of POPC phosphorus atoms from the nearest protein atom in Fig. 9C shows a sharp annular distribution shell centered approximately 3 Å from the protein. This, together with analysis of POPC-protein salt bridge formation (data not shown), indicates that the annular shell of POPC-COM (Fig. 9A) is largely the result of the formation of salt bridges between the phosphorus atoms of POPC and basic protein residues on the R- interfacial edge of the apoA-I amphipathic helixes (3, 43). The fact that UC forms a sharp annular shell at least comparable to POPC was surprising since electrostatic interactions of the polar but uncharged UC molecules with the protein would be expected to be weaker than interactions between the salt bridges formed by POPC with the protein.
To understand further the mechanisms for formation of a sharp UC annular shell, we measured the distribution of UC density along the apoA-I sequence. Examination of the dynamics of UC in the all-atom simulations of the MDSA-1, MDSA-2 and MDSA-4 particles showed that no UC approached closer than 11 Å from the helix 5/5 gaps (data not shown). However, a histogram plot (Fig. 9D) of the number of UC-OH oxygens that approach to within 10 Å of each amino acid residue in the apoA-I sequence for these three simulations indicated that, while not inserting into the helix 5/5 gap in the time frame of the simulations, UC does cluster in the vicinity of the single helix 5/5 and each helix 7/3 pair. In fact, the UC density in the helix 5 region (central open arrowhead) is approximately twice that of each individual helix 7/3 pair (peripheral open arrowheads), similar to the ratio of insertion of UC into the helix 5/5 gap versus each individual helix 7/3 domain (Fig. 8). Further, similar to the CG simulations (Fig. 8), there is no clustering in the vicinity of the individual helix 4/6 domains.
The distorted helix 5/5 domain in the MDSA-3 simulation serves as a useful negative control. A histogram plot (Fig. 9E) of the number of UC-OH oxygens that approach to within 10 Å of the apoA-I sequence for the MDSA-3 simulation indicates that, while UC clusters in the general vicinity of the helix 7/3 pairs, it does not cluster in the vicinity of the helix 5 domain. These results suggest that the conformation of the helix 5/5 gap may affect the clustering of UC in its vicinity.
It is well known that LCAT is susceptible to product inhibition; CE-rich HDL particles bind LCAT effectively, and contain apoA-I and ample PC and UC substrates, but exhibit relatively low reactivity with LCAT since the CE products cannot be removed from the active site (29). We therefore examined the presence or absence of UC and CE in the helix 5/5 gap during both all-atom and CGMD simulations of CE-rich particles containing full length apoA-I with a molar ratio of 57:16:6:2 (POPC:CE:UC:apoA-I). Examination of the final structures of the four all-atom MDSA simulations of the 57:16:6:2 particles shows one CE (Fig. 10A), one UC (Fig. 10B), and two terminal methyls (Fig. 10C and D), present in the four helix 5/5 gaps.
We proceeded to create CG models by CG mapping of three of the final all-atom MDSA structures (MDSA-1, MDSA-2 and MDSA-4) and subjected these CG structures to 18 µs CGMD simulations. Inspection of these simulations showed that CE blocks the helix 5/5 gap in approximately 27% of all trajectory frames inspected at 400 ns intervals (see examples of insertion of methyl beads, UC and CE into the helix 5/5 gap in Fig. 11). Fig. 12 quantitatively compares the dynamics of CE insertion between the helix 5 pair sites using distance measurements of CE from the pseudoatom of the nearest G129; CE blocks the helix 5/5 gap in 29% of the frames using a 6 Å cutoff distance. We conclude that the results of Fig. 10–Fig. 12 showing blockage of the helix 5/5 gap by CE molecules are compatible with experimentally demonstrated product inhibition in CE-rich particles, inhibition that is directly correlated with particle size, i.e., the number of CE molecules per particle (29). Shih, et al (44) have also performed CGMD simulations of CE-rich HDL particles but did not include UC in their structure. They also do not report on the relationship of CE to the apoA-I sequence other than reporting that the CE core was central and dynamic.
The enzyme LCAT transacylates the acyl chains of PC to cholesterol, forming CE and lysolecithin (45). Although a model of the secondary structure of LCAT has been proposed (46), the lack of a detailed three-dimensional model of LCAT complicates any attempt to understand the molecular nature of the interactions of the hypothesized presentation tunnel with the enzyme.
As noted, one function of apoA-I is presentation of acyl chains of long chain PC to the hydrolytic active site of LCAT; this is the rate limiting step in the process. We hypothesize that the binding of the putative lipase ‘lid’ of LCAT to PL-rich HDL particles—probably in the vicinity of the helix 6/4 domain (47–50)—opens a passage from the bilayer interior to the active sites of LCAT (for a review of LCAT structure see (51)). Since we propose that apoA-I acts as an amphipathic conduit for movement of hydrophobic acyl chains and amphipathic UC to the active site of LCAT, apoA-I in effect becomes part of the LCAT molecule when substrate is in the tunnel and the tunnel connects to the active site of LCAT. Studies by John Parks (52) suggest that the poor LCAT reactivity for long-chain sn-2 acyl groups is multifactorial, involving turnover of substrate molecules at the active site of the enzyme and activation of the enzyme by its cofactor apoA-I. Our model is compatible with these observations since it suggests that turnover at the active site should be indistinguishable from activation, i.e., both would depend upon apoA-I structure.
The majority of studies—although not all (53)—designed to identify the unique structural elements of apoA-I involved in activation of LCAT have implicated the central antiparallel double belt domain of discoidal HDL (47–49, 54–63), specifically residues 99–164 (helixes 4–6). Although several publications conclude that helix 6 is the lone LCAT activating domain in apoA-I (29, 56, 60, 61, 63, 64), the rationale for this conclusion is open to question: i) The conclusion is based upon the use of an incorrect structural model for apoA-I on discoidal HDL, the so-called picket fence model (64). When the double belt model is considered, major structural mutations of helix 6, such as deletion, replacement with the helix 10 sequence or reversal of its sequence (56, 60, 61), are difficult to interpret, since these mutations would have negative effects on antiparallel double belt salt bridge interactions of the mutated helix 6 with its pairwise mate, helix 4. Since this pairwise interaction occurs on both sides of the helix 5/5 domain, these mutations also should have a major effect on the intervening antiparallel helix 5/5 structure. As an example of pairwise cross-talk, a rather subtle mutation affecting the solvent inaccessible salt bridge in helix 4, E110A/E111A, has been shown to significantly affect LCAT activation (49). Since this mutation would be expected to influence the structure of both helix 4/6 domains and the intervening helix 5/5 domain, these results cannot be used to prove that helix 4 is the lone LCAT activating domain in apoA-I. ii) Mutations of three Arg residues (R149, R153 and R160) in helix 6 (residues 143–164) appear to produce loss of LCAT activity, apparently not through loss of activation but through loss of binding between apoA-I and LCAT due to disruption of salt bridge interactions (47, 48). iii) Wu, et al (50) using hydrogen-deuterium exchange mass spectroscopy have suggested that LCAT associates with apoA-I near the helix 6/helix 7 junction. iv) Mutations in other portions of apoA-I, e.g. the N-terminal domain (53), appear to have significant effects on LCAT activation.
ApoA-I is not unique as an activator of LCAT. Other apolipoproteins also containing class A amphipathic helical domains, such as apoA-IV, apoC-I, and apoE (65), can activate LCAT but are less than 20% as effective as apoA-I (29, 64). Model synthetic class A amphipathic helical peptides have been shown also to activate LCAT, again less effectively than apoA-I (66, 67). It appears, therefore, that class A amphipathic helixes are necessary, but not sufficient, for the degree of activation of LCAT provided by apoA-I, suggesting that one or more unique structural elements of apoA-I are involved.
As we noted in our previous T-jump paper (17), the terminal domains (N- and C-termini) represent the most mobile and generally least stable regions of the apoA-I double belt. We showed in that paper, since confirmed for PL-rich particles containing full length apoA-I (Gu, Jones, Chen, Patterson, Catte, Jerome, Li, and Segrest, unpublished data), that the mobility of the terminal domains lead to exposure of large amounts of phospholipid acyl surfaces as the particles increase in size and through this mechanism the terminal domains act as polar lipid remodeling-switches, regulating polar lipid exchange and particle fusion.
Since apoA-I is not unique as an activator of LCAT, we suggest that the N- and C-terminal domains, and other parts of the double belt, by exposing acyl chains to solvent, have the ability to activate LCAT nonspecifically. We propose that the central helix 4-helix 5-helix 6 domain in apoA-I is the unique structural element that provides full specificity to apoA-I activation of LCAT through substantially more LCAT binding and substrate binding and presentation to LCAT thorough the presentation tunnel.
Our use of molecular dynamics to examine the role of apoA-I in activation of LCAT involved simulations of discoidal particles containing 160 POPC, 24 UC and two full length apoA-I. Due to current limitations on simulation times achievable with current computer architectures for such large structures—200,000–300,000 atoms—standard all-atom molecular dynamics was deemed inadequate. Therefore, we used a combination of all-atom MDSA developed by trial and error to bypass kinetic energy barriers to the global structures of these particles and CGMD to achieve 18–20 µsec simulations. Both methods, MDSA and CGMD, have strengths and weakness but not the same ones: i) The strength of MDSA is that it is all-atom; the principle weaknesses are use of high temperature jumps and a failure to robustly sample long simulation times. ii) The strength of CGMD is that it samples longer simulation times; the principle weakness is use of pseudoatoms to represent clusters of carbon atoms, more of a problem for proteins than for lipids.
Based upon RMSDs plotted in Fig. 1B–D, the protein portions of the four structures simulated by the MDSA protocol appear to have converged after a total simulation time of 30 ns. All four MDSA simulations resulted in a gap between the helix 5/ 5 domains (15) that exposed a small region of the underlying PC acyl chains to the solvent; we refer to this gap as the amphipathic presentation tunnel. In the last 20% of the three simulations in which the juxtaposed pairwise central residue, K133, forms intrahelical salt bridges with E136 (MDSA-1, MDSA-2 and MDSA-4), at least one terminal methyl is inserted into the amphipathic presentation tunnel to within 6 Å of G129 in 99% of the frames (Fig. 6A).
In our initial use of the particle reduction methodology (16) we showed that the voids created by POPC deletion—voids filled with water molecules by using the VMD solvation module—collapsed with water expulsion into bulk solvent within 250 ps. The lipids in these particles therefore equilibrate much more rapidly than the protein. The lipid bilayer does expand during the initial 500 K simulation but relaxes back to a minimal surface bilayer structure after cooling to 310 K (Gu, Jones, Chen, Patterson, Catte, Jerome, Li, and Segrest, manuscript submitted for publication), an important point suggesting equilibration. During surface expansion caused by the T-jump step, the curved palmitoyloleoylphosphatidylcholine (POPC) bilayer surfaces approach planarity. Relaxation back into saddle-shaped structures after cool-down and equilibration supports our published saddle-shaped particle model (16), a model recently confirmed for reconstituted particles by Miyazaki, et al (68). The average change in secondary structure during the 10 ns of 500 K simulation has been quantified (Gu, Jones, Chen, Patterson, Catte, Jerome, Li, and Segrest, manuscript submitted for publication). The average helicity for the four 160:24:2 particles decreased from about 92% to 73% during the heating phase. During the next 20 ns of cooldown and equilibration at the end of the MDSA protocol, the average helicity increased again to 77%. The apoA-I in experimentally reconstituted particles with a comparable composition and size has been show by circular dichroic spectroscopy to have a helicity that ranges from 72–80% (69, 70), a range quite compatible with our MDSA results.
While a CG model can be mapped to certain helical variations—the classic α helix, the 3/10 helix and the π helix—there is no present algorithm for the coarse graining of α11/3 helixes. However, the α11/3 helix represents one extreme of the spectrum of helix pitches available to α helixes (the pitch can range from 3.5 to 3.7 residues per turn, the α11/3 helix being 3.66667). As we noted in our initial paper describing the α11/3 helix (3) and elsewhere (71), the classic α18/5 helix produces a twist in the hydrophobic face of 20° for every 22-residue repeat; over three repeats the twist would be 60°. This twist might be a problem over the full length of apoA-I but, due to the ability of proline residues to partially readjust the twist, 20° per repeat is not a significant problem over three or so repeats. Since β turn and hairpin loop conformations are scattered throughout each of the final four MDSA structures (unpublished observations), the effects on lipid interactions of coarse graining the α11/3 helixes of the all-atom models to α18/5 helixes are minimized since the coarse grained α11/3 helical regions can adjust the angle of their hydrophobic faces to interact with the lipid. In fact, examination of the final structures of the three CGMD simulations of the PL-rich discs shows that the hydrophobic residues have indeed aligned themselves to face the lipid along the full length of each apoA-I forming the double belt.
There is a disparity in percent insertion of terminal methyls in MDSA versus CGMD that requires discussion. Due to differences in scale of all-atoms versus pseudoatoms, distance measurements between the two methods cannot be directly compared. Examination of Fig. 6 shows that terminal methyl counts go to zero at 2 Å for MDSA (Fig. 6A) but at 4 Å for CGMD. Assuming a disparity of 2 Å between the two methods, Fig. 6B shows that, in the CGMD simulations, at least one terminal methyl is inserted into the amphipathic presentation tunnel to within 8 Å of G129 (CG equivalent of 6 Å in MDSA) in only 71% of the frames compared to 99% for MDSA. Why is that?
Part of the reason, of course, is that UC is competing with the terminal methyls for insertion. Another likely reason is the CG terminal methyl pseudoatom represents the last four acyl carbons rather that the terminal methyl carbon per se. Differences in force fields and size between the two terminal methyl representations likely inhibit insertion of the pseudoatom representation into the amphipathic presentation tunnel compared to the all-atom representation.
A second issue requiring discussion is disparity in UC insertion between the two methods. Part of the reason may be that the relative inhibition of insertion of the pseudoatom representation of the terminal methyl into the amphipathic presentation tunnel compared to the all-atom representation allows an increased insertion of UC (e.g. the CG pseudoatom, representing as it does four carbons is larger and may be considerably more hydrophobic that the all-atom version that represents a single carbon). Another reason for the disparity may be differences in simulation times between the two methods. With longer simulation times, all-atom MD simulations might result in UC insertion. The fact that UC is inserted in one of four final MDSA structures for the CE-rich particles (Fig. 10B) supports this possibility; these particles have been subjected to two MDSA simulations, first as a PL-rich particle and then as a CE-rich particle, resulting in an increase in sampling of conformational space. In any case, the exact molar ratios between terminal methyls and UC at equilibrium in both PL-rich and CE-rich particles, and between those two molecules and CE in the CE-rich particles cannot be determined from the MD simulations reported here.
Parametization in CHARMM of the polar UC-OH group oxygen gives it a small negative charge of −0.23 (one electron = −1.0). Once an UC molecule approaches the helix 5/5 or helix 7/3 domains, it would be retained by forming a charge pair (72) with either one or both of the K133 or the single R177. This charge pair formation at least partly explains the long residence time of up to 6 µsec observed for a single UC in the amphipathic presentation tunnel during CGMD simulations (data not shown).
The concentration of UC in the vicinity of the helix 5/5 domain shown in Fig. 9D and E appears to be related in part to the formation of charge pairing and hydrogen bonding of UC-OHs with basic and polar residues, respectively, on the R-face of helix 5. Fig. S4A of the Supporting Information is a histogram plot of the average number of UC-OH oxygens in MDSA-1, MDSA-2 and MDSA-4 that are within 5Å of any given protein charge pair acceptor (Lys, Arg, Gln, Asn, Ser, or Thr); residues R131, Q128, S142, and Q127 in helix 5 are most closely associated with UC-OH. The locations of these four residues are shown in Fig. S5A of the Supporting Information. A similar histogram plot for MDSA-3 (Fig. S4B of the Supporting Information) shows significant bonding of UC-OH to only one residue in helix 5, K133, the residue that is twisted away from the wheel position 2 to the R-face (Fig. 2C). Plots of the average number of UC molecules within 3 Å of any given residue (data not shown) indicate that the UC molecules that form the charge pairs with the four acceptors have closest van der Waals contacts with three Leu, L126, L137 and L134, and R131( Fig. S4B of the Supporting Information) that cluster around K133. The individual UC that make contact with the helix 5/5 domain are denoted in Fig 2C for all four MDSA simulations.
In our first MD simulation of PL-rich HDL particles (30), using continuum electrostatic calculations we observed after a short 1 ns simulation at 300K that a positive electrostatic potential extended across the hydrocarbon region of the particle in the middle of the bilayer. A positive electrostatic potential in PL-rich HDL particles might induce a small net diffusion of UC toward the helix 5/5 and helix 7/3 sites. Further, both K133 in the helix 5/5 domain and the single R177 in each helix 7 are in wheel position 2, a solvent inaccessible position (3, 42). Amphipathic helix wheel position 2 places a positive charge near the center of the bilayer and, as we suggested in our previous publication (72), may contribute significantly to the positive electrostatic potential extending across the hydrocarbon region of the particles.
As noted, one major difference in the conformation of the helix 5/5 domain of the MDSA-3 simulation and the other three simulations is twisting of one K133 away from the center of the amphipathic presentation tunnel in MDSA-3 (Fig. 2C). To explore the possible role of a transparticle bilayer electrostatic potential in concentrating UC in the vicinity of the helix 5/5 and helix 7/3 sites shown in Figs. 9D and E, we calculated electrostatic potentials for the four MDSA simulations looking for differences between MDSA-3 and the other three all-atom simulations. Examination of electrostatic potentials of the helix 5/5 domains of apoA-I in the different particles suggests that the lipid-facing surface of the amphipathic presentation tunnel in MDSA-3 has a slightly lower positive potential than the comparable domains in the other three particles (Fig. S5A and B in the Supporting Information). However, although there is a pronounced positive transparticle bilayer potential in each, there are no clear electrostatic potential differences between MDSA-3 and the other three particles (Fig. S5C and D in the Supporting Information). We conclude that transparticle hydrocarbon positive electrostatic potentials are not globally affected by the rearrangement of the helix 5/5 domain that occurs in the MDSA-3 structure. A fuller understanding of the origin and/or function of the positive transbilayer hydrocarbon electrostatic potential of PL-rich HDL particles awaits future computational and in synthetico experiments on PL-rich particles containing apoA-I mutants.
It is conceivable that the conformational change that occurs in MDSA-3, swinging of one K133 away from the center of the amphipathic presentation tunnel to the R-interfacial edge (Fig. 2C), could represent a functional intermediate involved in a later step in LCAT activation. For example, it is tempting to speculate that a 180° swing in K133 might facilitate the movement of one complete acyl chain of POPC through the presentation tunnel into the interior of LCAT for hydrolysis to the fatty acid intermediate.
Since insertion of methyl groups to within 4–5 Å of G129 would minimize the surface area of acyl groups exposed to the solvent (methylene groups cannot be so tightly inserted), the major driving force for insertion and retention of acyl chain methyl groups in the amphipathic presentation tunnel of HDL particles (Movie S3 in the Supporting Information) should be the hydrophobic effect. A similar mechanism is likely for blockage of the amphipathic presentation tunnel by CE in CE-rich HDL (Fig. 10–Fig. 12). The major driving force for insertion and retention of UC molecules in the amphipathic presentation tunnel of HDL particles (Fig. 7, Fig. 10 and Fig. 11) is likely a combination of charge pair interactions of K133 with the UC-OH group and the hydrophobic effect acting on the steroid ring.
The driving force for flipping UC-OH from the bilayer surface to the interior of the amphipathic presentation tunnel is less clear. Perhaps some larger scale motion of the helix 5/5 domain is involved, such as the flip-flop of K133 seen in MDSA-3. Based upon the CGMD simulations of the PL-rich HDL particles, opening and closing of the amphipathic presentation tunnel is dynamic; inspection of the trajectory of one CGMD simulation suggests that the tunnel is closed approximately 1/3 of the time.
It has been shown that UC-OH shows preference for the interior of polyunsaturated lipid membranes as the result of an increased rate of flip-flop (73, 74). It is therefore interesting to speculate that the known atheroprotective effect of a diet rich in polyunsaturated fats (2) might perhaps be related to an enhanced affinity of UC-OH for the presentation tunnel of apoA-I in HDL enriched in polyunsaturated PL.
As noted, the enzyme LCAT transacylates sn-2 acyl chains of POPC to form CE with approximately 90% fidelity (25, 26). Over the last 2 ns of the four MDSA simulations on the PL-rich particles the average number of sn-2 terminal methyls within 6 Å of G129 is 7.3. The percentage of sn-2 versus sn-1 terminal methyls within 4–10 Å of G129 reaches a maximum of 88% between 6 and 7 Å and decreases to 85% and 83% of total methyls at 5 and 4 Å, respectively. These results suggest the possibility that apoA-I, rather than LCAT alone, may be at least partially responsible for the sn-2 specificity of LCAT. This is an interesting observation, but as it is based on only four simulations, must await a larger ensemble of simulated particles for statistical validation.
In conclusion, we propose the working hypothesis that the helix 5/5 domain on both PL-rich and CE-rich HDL particles diverge to create an amphipathic presentation tunnel that exposes acyl chain terminal methyls and UC-OH to solvent. After attachment of LCAT to HDL, the helix 5/5 domains in apoA-I form amphipathic tunnels for migration of hydrophobic acyl chains and amphipathic UC-OH from the bilayer to the phospholipase A2-like and esterification active sites of LCAT, respectively. We further suggest that the product of the reaction, CE, can inhibit LCAT by blocking the amphipathic presentation tunnel.
Since the gap in the helix 5/5 pair occurs in essentially every simulation to date, it seems reasonable to conclude that the gap is real and not an artifact of simulation. Because this gap conceivably could have a function unrelated to LCAT activation, the LCAT presentation tunnel hypothesis is currently being tested by site-directed mutagenesis and additional MD simulations.
We would like to thank UAB Information Technology and Department of Mechanical Engineering for use of the IBM Blue Gene/L rack that they jointly maintain.
†This work was supported in part by the National Institutes of Health Grant P01 HL-34343 (to JPS).
SUPPORTING INFORMATION AVAILABLE
Figs. S1–S5 detail additional dynamics of the amphipathic presentation tunnel including a movie (S3). Fig. S6 illustrates examples of electrostatic potential calculations. This material is available free of charge via the Internet at http://pubs.acs.org.