Modeling of the Triglyceride-Rich Core in Lipoprotein Particles - The

13772

J. Phys. Chem. B 2008, 112, 13772–13782

Modeling of the Triglyceride-Rich Core in Lipoprotein Particles Anette Hall Department of Physics, Tampere UniVersity of Technology, P.O. Box 692, FI-33101 Tampere, Finland

Jarmila Repakova Laboratory of Physics, Helsinki UniVersity of Technology, P.O. Box 1100, FI-02015 HUT, Finland

Ilpo Vattulainen* Department of Physics, Tampere UniVersity of Technology, P.O. Box 692, FI-33101 Tampere, Finland, Department of Applied Physics, Helsinki UniVersity of Technology, P.O. Box 1100, FI-02015 HUT, Finland, and MEMPHYSsCenter for Biomembrane Physics, UniVersity of Southern Denmark ReceiVed: May 5, 2008; ReVised Manuscript ReceiVed: August 18, 2008

Triglycerides are a major component of many important biological entities such as lipoproteins and lipid droplets. This work focuses on two common triglycerides, tripalmitin and triolein, which have been simulated through atomistic molecular dynamics at temperatures of 310 and 350 K for 300-700 ns. In these systems, both structural and dynamical properties have been characterized, paying particular attention to understanding the packing of triglyceride molecules and their molecular conformations. Additionally, we study the liquidto-crystalline phase transition of tripalmitin through a temperature quench from the high-temperature isotropic liquid phase to 310 K, corresponding to a polymorphic, crystalline-like phase. The transition is characterized in detail through density, average molecular shape, and, in particular, the relevant order parameter describing the transition. 1. Introduction Triglycerides (see Figure 1), triesters of glycerol with fatty acids, are neutral lipids insoluble in water, commonly known as fats and oils. They are an important energy source in a normal diet and utilized for storage in body fat cells, adipocytes, as lipid droplets.1 Inside of the body, triglycerides are transported in lipoproteins, carrier particles of cholesterol and fatty acids. The biological relevance of triglycerides is highlighted by the fact that they are significant components of lipoprotein particles. Lipoproteins are one of the most studied carrier particles in cell biology. This is due to the fact that dysfunctions in their metabolism lead to altered concentrations of these particles in the bloodstream, which in turn is associated with cardiovascular disease, a major health risk in modern society.2-4 However, little is still known and understood of the structures and functions of lipoproteins at a molecular level or how these might affect the formation of atherosclerotic lesions, known as arterial plaques, that may lead to strokes and infarcts as they are ruptured from their origin.2,5 Two lipoproteins in particular are the key players in the development of the disease. Elevated blood concentration levels of low-density lipoprotein (LDL), colloquially known as “bad cholesterol”, or depressed levels of high-density lipoprotein (HDL), the “good cholesterol”, contribute toward a high atherosclerotic risk. This is not the whole truth, however, as it has been found that small dense triglyceridepoor low density lipoprotein particles might be more lethal than normal sized ones.2,3 No single explanation for this has been presented. Thus, the importance of understanding the structural details of lipoproteins is evident. * To whom correspondence should be addressed. E-mail: Ilpo.Vattulainen@ tut.fi.

Despite an extensive amount of experimental work focused on lipoproteins, the understanding of their structures has remained limited. This is largely due to the fact that lipoproteins are soft matter readily influenced by thermal fluctuations, for which reason lipoproteins do not have stable fixed structures. Additionally, lipoproteins are nanosized vehicles as their sizes range from about 8 nm of HDL through 25 nm of LDL to 90 nm for large VLDL (very low density lipoprotein) particles,2 which poses major conditions on the resolution of experimental techniques used to gauge their structural features. Consequently, atomistic as well as coarse grained simulations have recently been carried out to complement experiments with atomistic and molecular resolution to elucidate the details of lipoprotein properties.6-11 These studies have provided a great deal of insight into the protein-lipid interactions and the structures of the protein component surrounding the lipid core of lipoproteins. However, in almost all cases simulated by far, the description of the lipid core has been simplified and based on a single species of phospholipids only, while in practice, the molecular composition of the lipid core is considerably more complex,5 including significant amounts of phosphatidylcholines (PCs), sphingolipids, cholesterol, cholesteryl esters, lyso-PCs, triglycerides, and smaller amounts of other lipids. To understand lipoprotein properties, the full molecular composition of the lipid droplet should also be taken into account; otherwise, the protein-lipid interactions and their influence on protein dynamics remain only as partially understood. A recent simulation study by Catte et al.11 of HDL highlighted, for example, the importance of interactions between the cholesterol moiety of cholesteryl oleates and the ApoA-I proteins that surround the lipid droplet.

10.1021/jp803950w CCC: $40.75  2008 American Chemical Society Published on Web 10/10/2008

Modeling of the Triglyceride-Rich Core in Lipoproteins

J. Phys. Chem. B, Vol. 112, No. 44, 2008 13773

Figure 1. Chemical structures of the simulated molecules, (A) triolein and (B) tripalmitin. The numbering used to distinguish specific carbon atoms is also shown.

less common cases, triglycerides also exist in the shape of a trident, where all three fatty acid chains are side by side, pointing in the same direction. There is evidence that the molecule takes this form in a monolayer at hydrophobic-hydrophilic interfaces, such as an oil-water interface or mica.18 As for the biologically more relevant liquid phases, transitions between them have been found in experiments but have not been studied to a great extent. Typically, liquid triglycerides are isotropic in structure, but liquid-crystal forms have also been identified.19 The structure of liquid-crystal triglycerides is believed to be smectic, although a nematic structure has also been suggested.17,18 As far as atomistic simulation studies of triglyceride systems are concerned,20-23 there is reason to highlight the recent work by Sum et al.,23 who established a model for the behavior of triglycerides and broadly studied their thermodynamic and transport properties, focusing especially on viscosity. The present work furthers previous considerations by focusing on a number of issues related to fluid triglyceride systems. First, we characterize their molecular conformations and structural properties and consider how triglycerides pack themselves. Second, we deal with some of the dynamic properties related to mixing and diffusion. Finally, we characterize the dynamics associated with a structural transition where a triglyceride system in the liquid phase goes through a phase transition in the physiological temperature of 310 K, corresponding to the polymorphic phase R. The observed results provide insight into the structure of a triglyceride-rich core inside of lipoprotein particles and allow the development of coarse-grained model descriptions for triglycerides. The latter is currently being used in large-scale studies of lipoproteins (work in progress).

Given the complexity of lipoproteins, atomistic simulations have approached this topic by increasing the complexity of the system step by step. In a previous study, Heikela et al. used atomistic simulations to elucidate the structure and dynamics of a fluid mixture of cholesteryl esters,12 the aim being to model the core region of lipoproteins rich in cholesteryl esters. Here, we use the same approach for triglycerides. We employ atomistic molecular dynamics simulations to characterize the structural and dynamic properties of triglycerides in the absence of solvent, the principal aim being to better understand the structural features of triglycerides residing in the hydrophobic core of lipoproteins. Previously, triglyceride systems have been characterized to some extent, though the focus has largely been on the crystalline phases instead of the fluid ones. In the case of crystalline triglycerides, the molecule usually adopts a shape where two of the chains point in the same direction and the third one in the opposite one.13,14 When triglyceride molecules crystallize, their sn-1 and sn-3 side chains pack side by side, while the sn-2 chain is alone on the other side. Depending on the fatty acid distributions, the sn-2 chain can also pack alongside one of the side chains, leaving the third chain to point in the opposite direction. The former of these cases is called the “tuning fork” and the latter the “chair” configuration.13 When these structures stack together, they form layers of double or triple chain length, which further can adopt a certain orientation with respect to one another. Moving on, for crystalline triglycerides, there are three major polymorphic forms, the R, β, and β′. The R form is the least stable and has the lowest melting point. The β forms both stack at an angle, but the β is more stable than the β′. The different forms can be distinguished by X-ray diffraction studies by identifying the spacing in lattice structure.13,15-17 In some

2. Methods 2.1. Model and Simulation Details. We studied two triglycerides, tripalmitin (TP) and triolein (TO). They were chosen as oleic and palmitic acid are the two most abundant fatty acids in the human adipose tissue.24 Tripalmitin is a saturated triglyceride, having three palmitoyl chains connected to the glycerol backbone. Triolein composed of three oleoyl chains has one double bond in each of its fatty acid chains, making the chains monounsaturated. Schematic pictures of the studied molecules are shown in Figure 1. The figure also shows the numbering used to distinguish certain atoms in the molecules, namely, the central and the carboxyl carbons. The parameters of the force field for bonded and nonbonded interactions in TP were taken from similar combinations of atom types in the dipalmitoylphosphatidylcholine (DPPC) lipid molecule (see http://moose.bio.ucalgary.ca/files/). The partial charge distribution for atoms was set up according to that for the sn-1 acyl chain of DPPC. The topology for the TO molecule was constructed similarly, except for the double bond regions where the parameters were taken from the oleoyl chain of palmitoyloleoylphosphatidylcholine (POPC).25 For the purpose of completeness, after a short 5 ns simulation of the system of 108 triolein molecules using the above force field, the partial charges were also computed by QM calculations from eight different conformations using the approach described in ref 26. The results were found to be consistent. The initial configurations for large systems were built from the conformations of four tripalmitin molecules, which were put together in a cell that, in turn, was cloned three times in every direction. This resulted in a simulation box of 108 molecules. The density of the box was scaled successively to 1030 kg/m3, which is in the range of the experimental density

13774 J. Phys. Chem. B, Vol. 112, No. 44, 2008 of a LDL core.27 Energy was minimized after each increase in density. The triolein system was built similarly, except for the initial structure, which was done by editing the tripalmitin structure; the double bond was added into each chain, followed by an addition of two extra CH2 groups. Further systems were constructed as follows. The largest system of 500 tripalmitin molecules was constructed from the final configuration of the 108 triglyceride system. Additionally, another configuration for a system comprised of 128 tripalmitin molecules was made by removing excess molecules from the 500 tripalmitin system. All systems were simulated at two different temperatures, 310 and 350 K. TO systems are expected to be in an isotropic liquid phase at both temperatures as the melting point of TO from the liquid to the first polymorphic, crystalline phase (β) is around 278 K.28 For TP, the corresponding transition temperature is about 339 K.28 Here, let us briefly comment on the commonly used nomenclature and its relation to the simulations that we have conducted. The high-temperature phase is obviously an isotropic liquidlike phase characterized by substantial entropy and disordered molecular conformations. In terms of concepts in polymer physics, this phase could be considered as a liquid polymer melt. Then, the low-temperature phases are known as polymorphic forms (R, β′, β). In the literature, these polymorphic forms are often referred to as crystals;29 thus, in this work, we follow the same practice. Yet, there is reason to stress that the time scales considered in our simulations are on the order of 700 ns or less, which, despite their considerable length, are definitely too short for the formation of true crystalline structures. Instead, we expect that the ordered structures observed in our simulations for TP at 310 K (see below) are crystalline-like, describing intermediate structures prior to the final crystalline phase. Anisotropic pressure coupling was used for the 128 tripalmitin systems to examine its effect on ordering and possible layer formation. Initially, both systems had a full six directional anisotropy, which unfortunately caused instability for the 350 K system. This simulation was restarted with diagonal anisotropy. The other systems were simulated with isotropic pressure coupling. The effects of different pressure coupling types on the systems can be understood through the densities discussed below in detail. Anisotropic pressure coupling has no effect on liquid systems which are isotropic by nature. The TP system at 350 K has exactly the same average density for both the 128 molecule diagonally anisotropic system and the 500 molecule isotropic system. The other analyzed properties were also mostly identical and, as they provide no further information, are not presented in the results. Then, the 500 TP at 310 K system was not found to crystallize during the simulation, partly due to the large size of the system which did not allow studies of times beyond 200 ns. Nonetheless, we expect that the isotropic pressure coupling used in this case played a role too since it likely does not allow the proper formation of a layered structure due to the inability of the system box vectors to freely fluctuate. Anisotropic pressure coupling was also tested on the system, but the simulation of 100 ns was not long enough for any significant change to happen, and therefore, no data of this is shown. For these reasons, we do not discuss the 500 tripalmitin system at 310 K any further in this paper. The phase transition associated with the TP system is discussed through the small 128 TP system. The simulations were performed with the GROMACS software package using versions 3.1 and 3.3.30-32 A time step

Hall et al. of 2 fs was used for all MD simulations, with the LINCS algorithm applied to constrain all of the bonds in the system. The van der Waals interactions were calculated up to a cutoff radius of 1 nm and the Particle Mesh Ewald (PME) technique33,34 was utilized for the long-range Coulombic forces, with the grid size of 0.12 nm and cubic interpolation. The simulations were conducted in the NpT ensemble. The pressure was kept constant at 1 bar by the Berendsen barostat35 with a time constant of τ ) 1 ps. The Berendsen temperature coupling35 was used for equilibration, after which it was changed to the Nose´-Hoover thermostat36,37 for the rest of the simulation. The time constant for the temperature coupling was τ ) 0.1 ps. The 108 molecule triolein systems were simulated for 300 ns at both temperatures, 310 and 350 K. The 128 tripalmitin system was simulated for 700 and 300 ns at 310 and 350 K, in respective order. In a similar manner, the system of 500 tripalmitins was simulated for 200 ns at 310 and 350 K. The equilibration period was approximately 50-100 ns for all of the systems, except for the system of 128 tripalmitins at 310 K, in which case the equilibration period was 550 ns. On the basis of system density and radial distribution functions for a number of atom-atom pairs considered, these time scales were sufficient to guarantee proper equilibration. Thus, unless mentioned otherwise in the below discussion, only the trajectories after the equilibration periods were used for analysis. The results for varying system sizes were compared to one another to study the effects of pressure coupling, periodic boundaries, and finite size. 2.2. Data Analysis Methods. As the simulations are performed in constant pressure, the density of the system is left to alter and fluctuate. These variations can be monitored to examine the extent of equilibration of the system. The fluctuations are rather large, but the running average shows the tendencies of the density to drift or maintain stability as the equilibrium is reached. As such, the density can be used as a measure for the systems’ degree of equilibration. The time evolution of other properties, such as the radial distribution function (RDF) for a given pair of interacting bodies, can also give insight into equilibration processes. The RDF is a pair correlation function between particles A and B. It characterizes the local structure of a fluid by describing the probability to locate a particle at a certain distance relative to a random distribution at equal average density. It has the form NA 〈FB(r)〉 1 gAB(r) ) ) 〈FB〉local 〈FB〉localNA i∈A

NB

∑∑ j∈B

δ(rij - r) 4πr2

(1)

in which 〈FB(r)〉 is the particle B average density at distance r around particles A, 〈FB〉local is the average density of particles B inside of the maximum radius, NA and NB are the total number of particles A and B, and δ(rij - r) is a delta function. In the calculation, r stands for a differential spherical slice, rather than an absolute value. In the RDFs calculated for this work, groups A and B are equal. The center of mass RDF is especially useful in characterizing phase transitions. For a perfect crystalline solid, the RDF would be simply a collection of delta peaks at radii corresponding to neighbors in the lattice structure. On the other hand, for an ideal liquid, the radial distribution would quickly rise to and stabilize at unity. Because of this, the RDF can be used to establish the degree of solidification or fluidity of a substance, especially near its phase boundary. A solidifying liquid exhibits deeper depres-



sions and sharper peaks than those in a more liquid state; even additional longer-range peaks, caused by remnants of the lattice structure, may be visible. Furthermore, the coordination number (calculated by integrating over the first peak of the RDF) is a property that describes the number of neighboring particles as given by the radial distribution function. The orientation ordering in the systems is examined via an order parameter that measures how strongly correlated some vector describing a molecular orientation is aligned as a function of intermolecular distance. This is defined as the average of the second-order Legendre polynomial

1 S ) 〈3 cos2 φ - 1〉 2

(2)

in which φ is the angle between two vectors (defined in the text below) separated by a distance r, each vector characterizing the orientation of a triglyceride molecule (or its acyl chain). This order parameter is frequently used to study the phase structure of liquid-crystal systems. To describe the shape and molecular structure of triglycerides, we have considered distributions for a number of angles between the different acyl chains in the same molecule, which also yield insight into the planarity of the molecules. By doing so, one is also given some view of the molecules’ conformations that can be utilized in comparison with and development of coarsegrained models for these molecules. The radius of gyration is a measure of the size and compactness of a molecule often used for polymers and other molecular structures. It is defined as the mean square distance of each particle in the structure to its center of mass rcm

RG ≡

N

∑

1 (b r -b r cm)2 N i)1 i

(3)

where b ri is the position of particle i ) 1, 2,..., N in a molecule comprised of N particles. The moment of inertia is the rotational analogue to mass, a property that resists movement around a certain axis. It is generally defined relative to an axis as N

Ixx )

∑ mixi2

(4)

i)1

where mi is the mass of an atom i and xi is distance from the axis in question. As the primary rotational axes for the molecules are not known, a tensor is used to determine the moments of inertia relative to the center of mass of a molecule.38 The offdiagonal elements of the tensor are obtained from N

Ixy ) -

∑ mixiyi ) Iyx

(5)

i)1

This symmetrical tensor can then be diagonalized to find the moments of inertia around the principal axes of the molecule. The components of the tensor are arranged so that I1 is the largest value and I3 the smallest. From the moments of inertia in the directions of the principal axes can be calculated the radii of gyration corresponding to primary axis j through

RGj )

∑

Ij N m i)1 i

(6)

Here, i runs over all of the atoms in the molecule. From these three radii of gyration, one can also calculate the total radius of gyration RG.39 As for dynamical properties, we consider diffusion that describes how rapidly molecules move in space resulting from thermal forces and molecular collisions; at long time scales, the motion of the molecules can be viewed as a random-walk through the system. The diffusion discussed here refers to tracer diffusion, the movement of a single tagged particle. Here, we calculate diffusion for a group of particles through the meansquared displacement (MSD), in which case the diffusion coefficient D is given by

〈[b(t r + t0) - b(t r 0)]2〉 D ) lim tf∞ 6t

(7)

The coordinate b r stands for the center of mass of a molecule, for which the time evolution is examined, and the brackets denote an ensemble average over all molecules and origins in time in the simulated trajectories. 3. Results and Discussion 3.1. Density and Equilibration. The densities of the systems are presented in Figure 2 as 1 ns running averages. The degree of the systems’ equilibration can be studied from the stability of the densities. All of the systems in Figure 2A turn out to be stable by 100 ns of simulation. Of particular interest is the sharp rise in the 310 K tripalmitin system (see Figure 2B), which is due to a phase transition from liquid-like to a more ordered liquid-crystalline-like state. We discuss this process in more detail below. The densities are evaluated from the most stable part of each system by calculating the running averages of 10 ns and then using these data to compute the averages and standard deviations. The results are shown in Table 1. For comparison, experiments give a value of 875.2 kg/m3 for tripalmitin at 347 K and 915 kg/m3 for triolein at 288 K.40 Sum et al.23 have presented density values of 901 kg/m3 for triolein at 313 K and 867 kg/m3 at 373 K. The molecular dynamics simulations reported in the same article give results that differ from the experimental values approximately by the same amount as the ones calculated here, albeit in different directions.23 Considering the differences in force fields, these changes are expected. It can be concluded that the values for the density obtained here are reasonably well in accordance with experimental results and previous simulation studies; the differences are only on the order of 2-3%. 3.2. Liquid-to-Crystalline Transition of Tripalmitin. The above results for density (Figure 2B) have already clearly indicated that there is a structural transition from the isotropic fluid to a liquid-crystalline phase for the tripalmitin system of 128 molecules at 310 K. This is expected because the melting point from the liquid to the ordered polymorphic phases of tripalmitin is about 339 K.28 Here, we characterize this transition in more detail through the time dependence of the order parameter associated with the transition. We define an order parameter that reflects the orientation order of triglyceride molecules with respect to the average direction of the chains in the layered structure. The average

13776 J. Phys. Chem. B, Vol. 112, No. 44, 2008

Hall et al.

Figure 3. The time evolution of the order parameter of the tripalmitin at 310 K system during its phase transition between 300 and 700 ns. The order parameter is calculated for the vectors formed by the snchains compared to the average direction of the chains; see text for details. Error bars calculated by the standard deviation of the mean are presented.

Figure 2. (A) Time dependence of the 1 ns running average density of the liquid triglyceride systems from 0 to 300 ns. (B) Time evolution of the 1 ns running average density of the tripalmitin system at 310 K from 0 to 700 ns.

TABLE 1: Densities of the Various Systems and the Time Interval from Which They Have Been Calculated system

density (kg/m3)

time interval

108 TO at 350 K 108 TO at 310 K 128 TP at 350 K 128 TP at 310 K 500 TP at 350 K 500 TP at 310 K

846.7 ( 0.3 876.9 ( 0.4 853.8 ( 0.3 952.0 ( 0.6 853.8 ( 0.1 885.4 ( 0.3

150-300 ns 100-300 ns 100-300 ns 550-700 ns 100-200 ns 50-200 ns

direction in the fully formed ordered structure, denoted by the normalized vector lbref, is therefore used as a reference and is computed from the end of the simulation where the structure is obvious. Three vectors (denoted by the normalized vectors lbi(t), i ) 1, 2, 3) are then formed from the central carbon (C0 in TP; see Figure 1) to each chain end. From these vectors, the time evolution of the order parameter, see eq 2, where φi is the angle between lbref and lbi(t) at time t, can be determined to characterize the transition. The results for tripalmitin at 310 K between 300 and 700 ns are shown in Figure 3. During the first 300 ns, the order parameter does not essentially change but remains at a small value, around 0.3-0.4. Then, rather suddenly, after 300 ns has passed, the order parameter rises from the average value of 0.4 to a highly ordered value of 0.9, where the system clearly stabilizes. The main steps of the ordering process take place in about 250 ns, the evolution toward full equilibrium requiring considerably longer times though. Similarities with the progression of density during the same time period are evident; see Figure 2B.

Figure 4 depicts examples of system structures during the evolution of the ordered phase. The formation of the crystallinelike phase is evident. There is reason to stress that the transition was found specifically for the small system, in which characteristic times for substantial structural changes are naturally much shorter compared to those in a considerably larger system such as the one composed of 500 TPs. The transition could also take place in the larger system, but due to the limited time scale of the simulation, it was not observed in the present work. The ordering process is rather rapid and can be decomposed into two parts. Right after the temperature quench, the system is in a metastable state and remains in the liquid phase. The time needed to initiate the formation of the crystalline phase may be profoundly long and was here shortened by considering a small system where minor perturbations are expected to readily drive the system toward the thermodynamically stable equilibrium state. Once the ordering process is initiated, we find the system to go through a phase transition in about 200 ns. Results of a similar nature have been observed by Marrink et al.,41 who employed coarse-grained model simulations to study the fluid-gel transition in a lipid bilayer and found that the gel formation process was preceded by a corresponding lag time of hundreds of nanoseconds, and the actual ordering of the gel phase also took place in a few hundred nanoseconds. The total ordering process occurred in about 1 µs. Molecular shapes and conformations reflecting the different thermodynamic conditions are discussed below. 3.3. Intermolecular Structure. 3.3.1. Radial Distribution Functions. Radial distribution functions are here exploited to gain insight into the intermolecular structure in the system. The most striking fact visible from the fluid-phase configurations (see snapshots in Figure 4) is the tendency of the oxygen in the carbon chains to group together forming a complex network or, in the case of a liquidcrystalline phase, a layered structure. To illustrate this behavior for the centers of the molecules to gather near each other and away from the chain ends, two RDFs are depicted, one for the central carbon (Figure 5A) and another for the center of mass of the molecules (Figure 5B). The radial distribution functions for the central carbon (C0 in Figure 1) show a very strong nearest-neighbor peak at 0.5 nm distance and the next nearest peak at 1.0 nm for all liquid



Figure 4. Snapshots of the TP system at 310 K during the evolution of the crystalline phase at different times: at (A) 300, (B) 360, (C) 420, and (D) 480 ns. Also shown are the last frames of the crystalline-like TP system (E) at 310 K at 700 ns and (F) at 350 K at 300 ns. The 310 K systems are presented with two simulation boxes to better illustrate the layer formation.

systems. Beyond this range, the distributions smoothly oscillate below and above the average value of one, slowly dampening. Overall, the systems display behavior characteristic of liquids characterized by short-range order. The only exception is the 128 molecule tripalmitin system at 310 K, whose RDF differs from that of the others due to the different phase. The crystallinelike phase for TP at 310 K is evident from the center of mass RDF. Here, by considering Figure 5B, one gets an impression that the RDF decays to one at a distance of about 2.5 nm and that there would not be further order beyond this range. This

interpretation is not correct, however. There is reason to stress that as the simulation box includes only one layer of tripalmitins, the RDF is computed only up to (roughly) 2.5 nm, which corresponds to a half of the linear size of the box. However, in practice, the layer repeats itself in the direction of the layer normal; thus, there will be a strong correlation with a correlation length determined by the interlayer distance, which, in our case, for crystalline TP, is about 4.5 nm (the system size in the layer normal direction). Due to the finite system size, this feature is not featured by the RDF in Figure 5, but it is nevertheless an


Hall et al.

Figure 6. The order parameter of the plane formed by the middle carbon (C0) and the sn-1 and sn-3 carboxyl carbons (C1 and C3) as a function of distance from the central carbons; see eq 2. Estimated error bars are shown for TP at 310 K.

Figure 5. (A) Comparison between the radial distribution functions of the central carbon for the systems of 108 and 128 molecules. (B) The radial distribution function of the center of masses of the triglyceride molecules for the 108 and 128 molecule systems. Estimated error bars for TP at 310 K are shown. For the liquid systems, the errors are significantly smaller due to the greater extent of equilibration.

inherent property of the system. Then, given the fact that the lamellar layer phase repeats itself in the layer normal direction, what would be the actual correlation length beyond which the correlations would die out? Unfortunately, our studies do not shed light on this question. It is evident that the correlations will not persist over longer distances due to thermal fluctuations that drive layers away from their ideal phase; thus, the ordering between the layers is expected to be of quasi-long-range. Data on the coordination numbers calculated from the radial distribution functions can be found from Supporting Information. 3.3.2. Order Parameter. Another way to examine intermolecular structure is to use the order parameter describing the orientation of triglyceride molecules with respect to each other. Here, it is calculated from a part of the glycerol backbone, namely, the plane defined by the central carbon and the sn-1 and sn-3 carboxyl carbons (C0, C1, and C3; see Figure 1). The studied angles are calculated from the normals of the plane. The distance is evaluated from the central carbon, so that comparison can be made with its radial distribution. Results for the order parameter are given in Figure 6. The tripalmitin system at 310 K exhibits strong ordering compared to that of the other systems, and it does not greatly fade at longer distances. For the fluid systems, some structure can be seen below the distance of 1.5 nm, after which the order parameter settles to zero. The liquid systems clearly exhibit a very similar ordering behavior to each other, while the crystalline-like system has a completely different pattern.

3.4. Intramolecular Structure. 3.4.1. Angle Distributions. The main focus for considering the intramolecular structure of the triglycerides is to obtain information that enables one to understand how these molecules in triglyceride-rich regions, such as in the core of low-density lipoproteins, pack themselves. Additionally, the data facilitates the construction of coarse-grained models for studies of related systems over scales much larger than those gauged in atomistic simulations. Figure 7 shows the distribution of angles between the chains in the glycerol region for the TP systems. These are defined through the vectors formed from the central carbon and the carboxyl carbons of the ester bond in each chain (C0 with C1, C2, or C3; see Figure 1). The side chains, sn-1 and sn-3, form a strong single peak with the maximum value at approximately 85°. Due to symmetry, distributions of the angles between the middle chain and either side chain are similar. There are two maxima at 120 and 159°. The TO systems have distributions very similar to the TP at 350 K. These are presented only as Supporting Information. When the different TO systems are compared with one another, one finds that the distributions broaden as the temperature increases, but the differences are minor. As for TP, the cases at 310 and 350 K are distinctly different in terms of quantitative behavior but closely related considering their qualitative features. Also in this case, the main difference lies in the width of the distributions; the lower the temperature, the more narrow the peaks. Figure 8 describes the planarity of the glycerol region. It compares the plane comprised of the central carbon and the two side chain carboxyl carbons to the vector from the central carbon to the middle sn-2 chain carboxyl carbon. It shows this area to be close to planar, with the most favored angle being 12°. The shape of the curve is dependent mostly on the temperature for the liquid systems. Here, the most distinct feature is the orientation of TP at 310 K, where the distribution freezes to a more well-defined structure characterized by a distribution that is more narrow than that for other systems. In addition to studying the angles of the glycerol region, the same analysis is made for the whole molecule. Figure 9 shows the angle distribution of the vectors formed from the central carbon with each chain end. The tripalmitin system at 310 K is strikingly different from the others. The liquid systems have broad and relatively even angle distributions, while the TP at


Figure 7. The distribution of angles formed by vectors from the central carbon to each of the carboxyl carbons; (A) TP at 310 K and (B) TP at 350 K. The angles are named according to the chains they represent.

Figure 8. The distribution of angles formed by the vector from the central carbon to the sn-2 carboxyl carbon and the plane formed by the central carbon and the sn-1 and sn-3 carboxyl carbons.

310 K has sharp peaks at 14°, and smaller ones at 170°, with no angles between 60 and 120°. This tells us that while liquid molecules can twist freely and form random orientations, the crystalline molecules are rigid and defined. Tripalmitin at 310 K also favors small angles between the side chains and large ones for side chain-middle chain distributions more strongly. The distributions for triolein can be found as Supporting Information. Analysis shows that in the crystalline-like phase (TP at 310 K), 55% of the molecules are in the tuning fork, 33% in the chair, and 9% in the trident formation. In the fluid


Figure 9. The distribution of angles of the vectors from the central carbon to the end of each chain; (A) TP at 310 K and (B) TP at 350 K. The angles are named according to the chains they represent.

systems, these conformations are not as common; 80% are in other forms. Figure 10 depicts examples of representative snapshots of the TP molecules in the fluid and crystalline phases. Similarly, Figure 11 shows the equivalent of the results in Figure 8 for the ends of the chain. Tripalmitin at 310 K exhibits clearly different behavior than the others. It has a higher affinity for smaller angles, resulting in more planar molecules. For the liquid systems, the distribution is constant for angles smaller than 13°, after which it begins to decline linearly. 3.4.2. Size and Shape. The radius of gyration gives an estimate of the size of a molecule. Figure 12 shows the distributions for the triglyceride systems. We find that that lowering the temperature enlarges the radius of gyration due to the chains becoming more rigid. Tripalmitin at 310 K has the sharpest and narrowest peaks, which is connected to the strongly ordered nature of the system. It is also the only one to have more than one peak in its distribution due to different specific molecular conformations. Generally, one finds that in the liquid phase, TP is smaller than TO. The moment of inertia is calculated with respect to the three principal axes, I1 being the largest and I3 the smallest. In Figure 13, the distribution of these moments is plotted separately for each system. Notably, the axis 3 has a moment of inertia much smaller than the two others for all systems. This would imply that the atoms of the molecule lie close to this axis, leading to a narrow elongated shape, such as a disk-like ellipsoid, like that in the tuning fork and chair conformations. The extra peaks in the tripalmitin system at 310 K on the other hand are connected to the trident formation.


Hall et al.

Figure 10. Examples of molecular conformations for tripalmitin: (A) tuning fork, (B) chair, (C) trident, and (D) liquid (random) conformation.

Figure 11. The distribution of angles of the vector from the central carbon to the end of the sn-2 chain and the plane formed from the central carbon and the sn-1 and sn-3 end of chain carbons.

Figure 12. Comparison between the distributions of the radii of gyration.

3.5. Diffusion. The mean-squared displacements (MSDs) were analyzed to quantify the rate of single-particle diffusion. The diffusion coefficient was extracted from the MSD; the attained results are presented in Table 2. Temperature enhances movement; thus, at 350 K, the systems have larger displacements. The tripalmitin molecules have shorter fatty acid chains, making it easier for them to move about as entanglements effects are weaker. This is the reason why triolein at 350 K has slower diffusion than tripalmitin, although unsaturation usually raises it.42 It can be also found that in the liquid phase, the diffusion coefficient becomes about 10 times larger as the temperature is raised from 310 to 350 K. The increase in the diffusion rate as the system crosses from a crystalline to a liquid phase is more substantial. For comparison, the experimental data for triolein gives values of 6.0 × 10-7 cm2/s at 350 K and 1.7 × 10-7 cm2/s at 310 K.42 These are somewhat larger than the simulated results but yield the correct qualitative trends. Direct quantitative comparison is difficult to realize, however, since the simulated time scale is short compared to the time scales followed in experiments.

To get some flavor to the time scales associated with the diffusion of triglycerides, let us consider TP and TO in the fluid phase and assume them to comprise the hydrophobic core of LDL, whose radius is roughly 15 nm. The mixing of different lipid components in the core is limited by the diffusion rate characterized by the diffusion time τD ) l 2/6D, where l is the distance diffused by a molecule and D is the corresponding diffusion coefficient. Given that D ≈ 3 × 10-7 cm2/s and l ≈ 15 nm, one finds τD ≈ 1 µs. Thus, the mixing of molecules inside of lipoproteins is a rather slow process. In practice, this estimate implies that simulation time scales should preferably be considerably larger than 1 µs to guarantee full mixing of different molecular components. Considering the substantial size of LDL, this is not feasible through atomistic simulations at the moment, but it is likely within reach through coarse-grained model simulations. 4. Concluding Remarks In this work, the structure, conformations, and dynamics of two common triglycerides have been characterized in detail



Figure 13. The distributions of moments of inertia about the principal axes: (A) TP at 310 K, (B) TP at 350 K, (C) TO at 310 K, and (D) TO at 350 K.

TABLE 2: Diffusion Coefficients for the Triglyceride Systemsa system

diffusion coefficient (cm2/s)

TO 350 K TO 310 K TP 350 K TP 310 K

2.2 × 10-7 2.4 × 10-8 3.7 × 10-7 1.0 × 10-9

a The uncertainty is estimated as about (25% using data from the three independent coordinates for diffusion.

using atomistic molecular dynamics simulations. The obtained results indicate that TP at 310 K, which is below its melting point, undergoes a phase transition to a crystalline-like state. This has allowed us to not only study the properties of liquid triglycerides but also to compare their differences to those of a crystalline-like state. The behavior of the liquid systems is clearly different from that found in the crystalline-like system. Additionally, information on the phase transition itself was obtained through the time development of the system density and the relevant order parameter. Let us here raise a few issues that are relevant in the context of lipoproteins. First, considering the core of lipoprotein particles and, in particular, the hydrophobic bulk region of low-density lipoproteins, it is largely comprised of triglycerides and cholesteryl esters. Recent atom-scale simulation studies have shown12 that the diffusion rate of cholesteryl esters is considerably smaller than the diffusion coefficient of triglycerides (as found in this work) when the systems are considered in the fluid phase. This is understandable since the steroid moiety of cholesteryl esters is expected to promote the packing of the molecules in the core and hence slow down the rate of diffusion.

Similarly, it has been found that cholesterol slows down diffusion in fluid lipid membranes;43 the larger the molar fraction of cholesterol, the smaller the diffusion coefficient. Second, assuming that we can again use the analogy to lipid membranes, it is worth pointing out that cholesterol has been found to eliminate the main phase transition between fluid and solidlike membranes when the amount of cholesterol in twocomponent mixtures of phospholipids and cholesterol is larger than about 30 mol %.44,45 That is, cholesterol promotes the order within a membrane but eliminates the formation of a crystalline state. In the present context, one is therefore tempted to propose that in the core of lipoproteins that consists of a mixture of triglycerides and cholesteryl esters, there is no true crystalline state under physiological conditions. Considering the above, detailed studies of mixtures of triglycerides and cholesteryl esters would be highly useful. Due to the long time scales needed for gauging ordering phenomena in a polymer-melt-like system, which is typically characterized by slow dynamics, it is unlikely that atomistic simulations would be feasible for this purpose. It is likely more appropriate to consider coarse-grained model simulations instead. To this end, the simulations presented in this work and those in ref 12 provide a realistic starting point for coarse graining in the bottom-up manner. Work in this direction is underway. Acknowledgment. We thank the financial support by the Glycoscience Graduate School (A.H.) and the Academy of Finland (A.H., I.V.). We also wish to thank the HorseShoe (DCSC) supercluster at the University of Southern Denmark and the Finnish IT Centre for Scientific Computing for computer resources.

13782 J. Phys. Chem. B, Vol. 112, No. 44, 2008 Supporting Information Available: Coordination numbers calculated from the radial distribution functions and the angle distributions for triolein. We also give the topology file for tripalmitin. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Martin, S.; Parton, R. G. Nat. ReV. 2006, 7, 373–377. (2) German, J. B.; Smilowitz, J. T.; Zivkovic, A. M. Curr. Opin. Colloid Interface Sci. 2006, 11, 171–183. (3) Tulenko, T. N.; Sumner, A. E. J. Nucl. Cardiol. 2002, 9, 638–649. (4) Jackson, R. J.; Morrisett, J. D.; Gotto, A. M. Physiol. ReV. 1976, 56, 259–316. (5) Hevonoja, T.; Pentikainen, M. O.; Hyvonen, M. T.; Kovanen, P. T.; Ala-Korpela, M. Biochim. Biophys. Acta 2000, 1488, 189–210. (6) Sheldahl, C.; Harvey, S. C. Biophys. J. 1999, 76, 1190–1198. (7) Klon, A. E.; Segrest, J. P.; Harvey, S. C. J. Mol. Biol. 2002, 324, 703–721. (8) Shih, A. Y.; Denisov, I. G.; Phillips, J. C.; Sligar, S. G.; Schulten, K. Biophys. J. 2005, 88, 548–556. (9) Shih, A. Y.; Freddolino, P. L.; Sligar, S. G.; Schulten, K. Nano Lett. 2007, 7, 1692–1696. (10) Catte, A.; Patterson, J. J.; Jones, M. J.; Jerome, W. G.; Bashtovyy, D.; Su, Z.; Gu, F.; Chen, J.; Aliste, M. P.; Harvey, S. C.; Li, L.; Weinstein, G.; Segrest, J. P. Biophys. J. 2006, 90, 4345–4360. (11) Catte, A.; Patterson, J. C.; Bashtovyy, D.; Jones, M. K.; Gu, F.; Li, L.; Rampioni, A.; Sengupta, D.; Vuorela, T.; Niemela, P.; Karttunen, M.; Vattulainen, I.; Segrest, J. P. Biophys. J. 2008, 94, 2306–2319. (12) Heikela, M.; Vattulainen, I.; Hyvonen, M. T. Biophys. J. 2006, 90, 2247–2257. (13) Himawan, C.; Starov, V. M.; Stapley, A. G. F. AdV. Colloid Interface Sci. 2006, 122, 3–33. (14) Walstra, P.; Klock, W.; van Vliet, T. Crystallization and Polymorphism of Fats and Fatty Acids; Sato, K., Garti, N., Eds.; Marcel Dekker: New York, 1988. (15) Ferguson, R. H.; Lutton, E. S. Chem. ReV. 1941, 29, 355–384. (16) Ferguson, R. H.; Lutton, E. S. J. Am. Chem. Soc. 1947, 68, 1445– 1448. (17) Ueno, S.; Minato, A.; Seto, H.; Amemiya, Y.; Sato, K. J. Phys. Chem. B 1997, 101, 6847–6854. (18) Zdravkova, A. N.; van der Eerden, J. P. J. M. J. Cryst. Growth 2006, 293, 528–540. (19) Callaghan, P. T.; Jolley, K. W. J. Chem. Phys. 1977, 67, 4773– 4774. (20) Yan, Z.; Huhn, S. D.; Klemann, L. P.; Otterburn, M. S. J. Agric. Food Chem. 1994, 42, 447–452.

Hall et al. (21) Engelsen, S. B.; Brady, J. W.; Sherbon, J. W. J. Agric. Food Chem. 1994, 42, 2099–2107. (22) Pascher, I. Curr. Opin. Struct. Biol. 1996, 6, 439–448. (23) Sum, A. K.; Biddy, M.; de Pablo, J. J.; Tupy, M. J. J. Phys. Chem. B 2003, 107, 14443–14451. (24) Christie, W. W. http://www.lipidlibrary.co.uk (2008). (25) Ollila, S.; Hyvonen, M. T.; Vattulainen, I. J. Phys. Chem. B 2007, 111, 3139–3150. (26) Repakova, J.; Capkova, P.; Holopainen, J.; Vattulainen, I. J. Phys. Chem. B 2004, 108, 13438–13448. (27) Hevonoja, T.; Pentikainen, M. O.; Hyvonen, M. T.; Kovanen, P. T.; Ala-Korpela, M. Biochem. Biophys. Acta 2000, 1488, 189–210. (28) Hagemann, J. W.; Tallent, W. H. J. Am. Oil Chem. Soc. 1972, 49, 118–123. (29) Kodali, D. R.; Atkinson, D.; Redgrave, T. G.; Small, D. M. J. Lipid Res. 1987, 28, 403–413. (30) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43–56. (31) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306–317. (32) van der Spoel, D.; Lindahl, E.; B.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1. (33) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089– 10092. (34) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577–8593. (35) Berendsen, H. J. C.; Postma, J. P. M.; Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (36) Nose, S. Mol. Phys. 1984, 52, 255–268. (37) Hoover, W. G. Phys. ReV. A 1985, 31, 1695–1697. (38) Arfken, G. B.; Weber, H. J. Mathematical Methods for Physicists, 5th ed.; Harcourt Academic Press: Orlando, FL, 2001. (39) Rubinstein, M.; Colby, R. H. Polymer Physics; Oxford University Press: New York, 2003. (40) Lide, D., Ed. CRC Handbook of Chemistry and Physics, 87th ed.; Taylor and Francis: Oxford, U.K., 2007. (41) Marrink, S. J.; Risselada, J.; Mark, A. E. Chem. Phys. Lipids 2005, 135, 223–244. (42) Callaghan, P. T.; Jolley, K. W. Chem. Phys. Lipids 1980, 27, 49– 56. (43) Falck, E.; Patra, M.; Karttunen, M.; Hyvonen, M. T.; Vattulainen, I. Biophys. J. 2004, 87, 1076–1091. (44) Sankaram, M. B.; Thompson, T. E. Biochemistry 1990, 29, 10670– 10675. (45) Vist, M. R.; Davis, J. H. Biochemistry 1990, 29, 451–464.

JP803950W

Modeling of the Triglyceride-Rich Core in Lipoprotein Particles - The

Recommend Documents