Aggregation of Lysozyme in the Presence of a Mixed Bilayer of POPC and POPG.

Understanding the molecular mechanisms by which amyloidogenic proteins interact with membranes is a challenging task. Amyloid accumulates from many human diseases have been observed to contain membrane lipids. In this work, coarse-grained molecular dynamics simulations have been used to inspect hen egg white lysozyme (HEWL) aggregation and membrane association in the presence of a pure POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) bilayer and a POPC and POPG (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylglycerol) mixed bilayer. It was observed that, in both cases, two HEWLs formed aggregates. In the presence of a mixed bilayer, after aggregation, the aggregated system started to interact with the membrane. It has been found that one of the lysozymes which came closer to the mixed bilayer unfolded more. The process of the initial insertion of an aggregated system in the mixed bilayer has been analyzed. The structural rearrangements of the protein and lipids were analyzed as well along the course of the simulation. Although with a pure POPC bilayer, aggregation was observed, the aggregated system moved away from the membrane. We believe that our study will provide considerable insights into lysozyme aggregation in the presence of a membrane environment.

Introduction

Biological membranes have a miraculous ability to regulate a wide range of physiological and pathological processes,[1,2] including signal transduction and cellular trafficking.[3] These processes involve diverse lipid–protein interactions.[4,5] Membranes interact not only with native proteins but also with aggregated and denatured proteins. Protein aggregation[6] has attracted significant scientific attention over the last few decades, as a number of fatal diseases, including neurological diseases,[7] type-II diabetes,[8] spongiform encephalopathy,[9] and so forth, have been shown to be involved in the formation of amyloidogenic aggregates. It has been perceived that cell membranes might be targeted by protein aggregates to cause physical changes in the membrane, leading to dysfunction and cell death.[10] It has been conjectured that given the marginal stability of the native proteins, the surroundings may influence the progression from a monomeric to an aggregated state.[11,12]

The formation of amyloidogenic aggregates involves not only the association of monomers but also denaturation.[13] Membranes have been suggested to have a compelling role in protein aggregation.[14] The composition of membranes also plays a crucial part in determining the nature of the protein–membrane interaction.[15] In the present study, we have used hen egg white lysozyme (HEWL) as a model protein to study the effect of membranes on its dimerization and the subsequent interaction with the bilayer. Lysozyme is an essential enzyme of the innate immune system with antimicrobial, antitumor, and immuno-modulatory activities that play a major role in lipid-binding properties.[16,17] It has also been reported that the binding of the lysozyme to the phospholipid bilayer can modulate the aggregation behavior of the lysozyme.[18] Due to its tendency to form aggregates in the presence of lipid membranes, HEWL has been used as a model to study the effect of membranes on the unfolding and aggregation of pathologically important peptides/proteins.[19] Apart from that, HEWL has a well-defined three-dimensional structure[20] and a relatively small size, which makes it a good candidate for computer simulations.[21] It is also homologous with the human lysozyme[22] and has been studied extensively for mechanisms of its folding.[23,24]

CG MD simulations[25,26] have been used to explore the aggregation pathway of HEWL in the presence of a pure POPC layer and a mixed POPC/POPG (7:3)[27] bilayer. The mixture of POPC to POPG lipids in a molar ratio of 7:3 mimics the bacterial inner membrane.[28] Our main objective was to explore the effect of a specific bilayer mixture of lipids in protein aggregation followed by the initial insertion of the aggregated system in the bilayer along with monitoring the conformation of proteins and the changes in the structure and dynamics of the lipid bilayer due to protein aggregation. The present simulation study focuses on the nature of the interactions of HEWL and membranes along with the mutual effects they have on each other’s structure and dynamics. We have also looked into the effects of the mixed bilayer on the protein aggregation.

Methods

The Martini model[29] has been used in this work because this can provide the needed reduced set of interactions while still providing the necessary chemistry of the entire system. The Martini model has been widely used for the past few years to provide a representation of force-field parameters in the case of proteins/lipid systems.[25,30]

System Setup

The system was prepared with the help of insane.py script.[31] A mixed bilayer was constructed using POPC and POPG molecules. The ratio of POPC and POPG was 7:3 in both upper and lower leaflets.[27] The input geometry for the heterogeneous membrane system was specified with lipid moieties as per the membrane system ratio in both upper and lower leaflets. The number of POPC was 252 per leaflet. The number of POPG was 108 per leaflet. The area per molecule of the membrane is as follows: for the first 1 μs, it was 0.620 ± 0.004, and for the last 1 μs, it was 0.621 ± 0.004 (Figure S1, Supporting Information). Two lysozymes were placed 60 Å away from the center of the mixed bilayer (shown by the solid red line in Figure ). The distance between the centers of mass of two lysozymes was 50 Å (shown by the solid black line in Figure ). A coarse-grained (CG) water bead model was used.[26] System neutralization was achieved with the addition of 0.15 mM NaCl salt.[32] The simulation box size was 114 Å × 114 Å × 114 Å.

Figure 1

Molecular details of the membrane.

Figure represents the molecular details of the membrane. Protein1, protein2, POPC, POPG, water, and ions were represented in orange, green, purple, cyan, blue, and red beads, respectively.

Simulation Parameters

The simulations described in this paper were performed with the GROMACS simulation package, version 5.0.7.[33] Simulations were performed using the Martini2.2 force field under isothermal isobaric conditions. The system temperature was maintained at 300 K.[34] The temperature was kept constant using the Berendsen temperature coupling algorithm with a time constant of 1 ps.[35] Semi-isotropic pressure coupling was applied using the Berendsen algorithm, with a pressure of 1 bar, independently in the plane of the membrane and perpendicular to the membrane.[17] A time constant of 1.0 ps and a compressibility of 3.0 × 10–4 bar–1 were used. The relative dielectric constant εr was set to 15 in water[36] due to the absence of partial charges in the standard Martini water model.[37] The system was equilibrated before the production run (equilibration curve of the system is shown in the Supporting Information, Figure S2). The system was simulated with an integration time step of 20 fs.[38] The trajectory analysis was performed with GROMACS,[33] and the snapshots were generated by VMD.[39]

Results and Discussion

Aggregated System in a Membrane Environment

Three independent simulations were run for 10 μs each at 300 K. A 7:3 POPC/POPG mixed bilayer was chosen as a mimic of the bacterial inner membrane.[28] It has been reported earlier that in the presence of liposomes of different molar ratios of the zwitterionic POPC and the negatively charged POPG, the native structure of HEWL got modified leading to elongated aggregates above a charge-density threshold.[15] It has been reported that by taking an appropriate ratio of phosphatidylserine and phosphatidylcholine, the significance of the negatively charged lipid toward the aggregation of the lysozyme was explained.[19] In our study, we have observed that the two lysozymes formed an aggregate in the presence of both the bilayers. While in the presence of the mixed bilayer, the aggregate of the proteins had shown membrane insertion, in the control run with only the pure POPC bilayer, we did not observe any initiation of insertion of the aggregated protein in the bilayer, which is in agreement with the past experimental results.[15] It was observed that the two proteins formed an adduct but drifted away from the POPC membrane (Figure S3 of the Supporting Information). For the mixed bilayer, two independent simulations were executed with different starting velocities. Figure shows the snapshots of the system at different time intervals. Similar observations were obtained for both simulations, and analysis was done for both independent simulations. These observations are in agreement with the experimental observations.[15] Here, the results of one of the trajectories (simulation 1) were reported. Some of the analysis of simulation 2 is given in the Supporting Information. During simulation, the following steps were observed.

Figure 2

Snapshots showing the aggregated system with respect to time. Protein1, protein2, POPC, and POPG are presented as orange, green, purple, and cyan van der Waals spheres, respectively. Image rendering was done with VMD.

Interaction of two lysozymes with the membrane independently:

At the beginning of the simulation, the two proteins were 60 Å away from the center of the mixed bilayer along the z axis. It was observed that after equilibration, the two proteins tried to come close to the membrane surface. Within 50 ns of simulation, both proteins reached the surface of the bilayer and started to interact with the membrane (second snapshot in Figure ).

two lysozymes:

After nearly 80 ns, an interesting observation was noticed, that is, both the proteins leaving the membrane surface. Again, the two proteins tried to come close to each other to form an aggregate.

The whole aggregated system again started to interact with the membrane:

When the two proteins assembled, the whole aggregated system again tried to come close to the surface of the membrane to interact with the mixed bilayer.

Lipid–Protein Interaction

Analysis of the distance between the center of mass of two proteins (Figure S4, Supporting Information) and the number of contacts between them with simulation time (Figure a) was done. It was found that, after 2 μs, the number of contacts and the distance between the center of mass of the two proteins became almost constant. As lipid–protein interactions are crucial for many cellular processes,[4,5] the interaction energies of individual proteins with lipids were determined and the partition of the total interaction energy in terms of electrostatic and van der Waals was done.[40] From the graph in Figure b, it can be noticed that within 1 μs, there was a significant interaction (both van der Waals and electrostatic interactions combined) between protein2 and lipid. It is also observed that the interaction between protein1 and protein2 stabilized significantly after 1 μs. Similarly, to calculate the contribution of the interaction energies of individual lipids POPC and POPG, the individual interaction energies of POPC and POPG with proteins were calculated (Figure S6, Supporting Information). It was observed that POPG interacted more strongly with proteins than that of POPC, which is also similar to the previous experimental results.[15] This graph even showed a significant change in the interaction energy after 1 μs for both lipids. Simulation 2 also provided the same results (Figure S6, Supporting Information).

Figure 3

Determination of the lipid–protein interaction: (a) number of contacts between two proteins and (b) interaction energy vs time.

The role of the lipid type (bilayer mixture of POPC and POPG) was examined in HEWL aggregation followed by the initial insertion of the aggregated system in the mixed bilayer. In an earlier report, the role of the bilayer mixture of POPC and POPG in the transition of HEWL into the aggregation-component conformation was explained by analyzing the change in the secondary structure of HEWL along with the interaction with the membrane.[17]Figure describes the stages of protein aggregation in the presence of a membrane. It has been reported earlier that membranes may help in the process of aggregation.[41,42] The increase in the protein–protein contact with the simulation time indicated the formation and stability of the aggregated system (Figure a). Figure b, shows that initially the lysozyme–lysozyme interaction energy was much more stabilizing than the lysozyme–membrane interaction energy. It was noticed that the lysozyme2–membrane interaction energy became more favorable after 1 μs. It is also seen from Figure b that negatively charged POPG has a stronger interaction than that of POPC. From the lipid–protein interaction energies, it can be concluded that both electrostatic and van der Waals interactions are crucial for the lipid–protein interaction, in agreement with previous reports.[11,18] It may be said that the lipid composition of the membrane, particularly its anionic phospholipid content, seems to play an important role in this process: the establishment of electrostatic interactions with acidic lipids favors membrane-induced protein misfolding and the subsequent aggregate nucleation.[17]

Protein Structure

Conformational analysis of the proteins is necessary to estimate the influence of the mixed bilayer on the structure and dynamics of proteins. This type of analysis has been done experimentally as well as theoretically.[15,17] In the present simulation, the changes in conformations of lysozymes were estimated by rmsd (root-mean-square deviation),[17]Rg (radius of gyration),[17,43] and SASA (solvent accessible surface area)[17,44] analyses. The details of the conformational drift were obtained by computing the rmsd for backbone beads of the proteins relative to the starting conformation.[45] For each protein, the rmsd increased initially and closely approached convergence within 1 μs (Figure a). It was noticed that, after 1 μs, the rmsd value of protein2, which was closer to the membrane (Figure ), had a higher rmsd than that of protein1. Protein2 showed an increase in Rg and SASA values after 1 μs (Figure b,c), compared to those of protein1. Rg (Table S1) and SASA (Table S2) values of the two proteins were calculated for the first and last 0.1 μs. It was found that averaged over the last 0.1 μs of CG simulations, the rmsd values of protein1 and protein2 were 0.71 ± 0.02 and 0.84 ± 0.03 nm, respectively, where it was observed that for the last 0.1 μs of the CG simulations, the average Rg values of protein1 and protein2 became 1.41 ± 0.01 and 1.43 ± 0.01 nm, respectively. For proteins in the pure POPC system, it was noticed that, for both proteins, the changes in rmsd, Rg, and SASA were almost the same (Figure S7). We have recently shown that, at 300 K, lysozymes can form aggregates, even in the absence of membranes.[46] When compared with the free simulations, it was observed that for protein2, which came closer to the mixed bilayer, the fluctuations in Rg and SASA were more than in the other systems. For example, the fluctuation in Rg for protein2 is in the range of ±0.35 nm compared to ∼±0.02 nm for the free system[46] and in the presence of the POPC bilayer. For SASA, the observed fluctuation for protein2 is ±3.19 nm2 in the mixed bilayer system and ±2.5 to 2.7 nm2 in the free system[46] and in the presence the POPC bilayer system.

Figure 4

Determination of the protein structure: (a) rmsd, (b) Rg, (c) SASA vs time, and (d) RMSF of the backbone bead of the protein.

The analysis of RMSF (root mean square fluctuation) of backbone beads of both proteins[47] (Figure d) was done. For both proteins, maximum fluctuations were observed for the beads from ∼90 to ∼110 and from ∼150 to ∼170, which were a part of the flexible loop. Even for simulation 2, we got similar observations (Figure S8, Supporting Information).

The higher values of rmsd, Rg, and SASA of lysozyme2 than those of lysozyme1 indicate the influence of the membrane on conformational changes in the lysozyme.[17] Up to 1 μs, the values of Rg, SASA, and rmsd of both lysozymes were the same, but after 1 μs, these value were higher for lysozyme2 because it was closer to the membrane (Figure ). The higher rmsd, Rg, and SASA values in the presence of membrane (Figure ) in comparison to the values of rmsd, Rg, and SASA in the absence of membrane (our previously published work) indicate the role of membranes in inducing conformational changes in proteins.[46] Similar to the experimental results, as in the POPC membrane, the aggregated system drifted away from the membrane and there was no initiation of insertion of the aggregated protein in the membrane;[15] for this reason, in this case, the changes in the values of rmsd, Rg, and SASA for both proteins were almost similar. This implies the importance of a specific type of bilayer in the case of protein aggregation in a membrane environment. The RMSF plot clearly gives information about the influence of the membrane on unfolding of both lysozymes.[48] It has been generally accepted that the main driving force for protein aggregation is the partial unfolding of the native state of the polypeptide chain into an aggregation-prone conformation with exposed hydrophobic regions.[49,50] Ample evidence from both theoretical and experimental studies suggests that the protein aggregation potential is substantially enhanced in a membrane environment.[11] The present simulation results have demonstrated that the aggregate–membrane association was accompanied by the destabilization of the lysozyme structure and exposure of the flexible part of the protein chain.

Effect of the Aggregated Protein on the Membrane Structure and Dynamics

A detailed understanding of changes in each lipid of the mixed bilayer is necessary to gain an insight into the influence of the proteins on the membrane due to the aggregation and insertion of proteins in the mixed bilayer, as no such simulation has been done previously by using the CG molecular dynamics (MD) simulation method. An earlier MD simulation method reported about the role of the POPC/POPG mixture in the transition of HEWL into the aggregation-component conformation for a short time scale only.[17] The bilayer structural properties (conformational ordering, head-group orientation, and lipid in-plane distribution), as well as the dynamic changes in the local lipids (if the phosphorus bead of the lipid was within 10 Å of the protein, it was referred as the local lipid and the rest were referred to as bulk lipids) have been highlighted.

To better understand the effect of protein aggregation on the lipid structure, changes in the number of local lipids (Figure ) and bulk lipids with simulation time (Figure S9, Supporting Information) was first examined. A striking feature was found that the number of local lipids for both POPG and POPC increased up to 1 μs, and after 1 μs, the values became constant for both POPC and POPG. Similarly, it was noticed that the numbers of bulk lipids decreased up to 1 μs, and after 1 μs, these values became almost constant (Figure S9). It is also to be noted that initially the ratio of local POPC to local POPG was 7:3, but at 10 μs this ratio became 1:1 (the ratio of the PO4 beads of POPC and POPG), indicating the accumulation of more POPG near the protein, whereas, for simulation 2, we obtained a ratio of 1:1.17.

Figure 5

Number of local lipids vs time.

To analyze the effect of the aggregated protein on bilayer properties, two-dimensional plots for bilayer thickness and lipid-order parameters (Figure S10, Supporting Information) over the entire simulation time period were generated. The bilayer thickness was estimated as the average distance between the PO4 beads of the two leaflets using previously developed analytical tools.[51] The obtained value was ∼4.23 ± 0.14 nm. The one-dimensional lipid-order parameter (P) for POPC and POPG in the bilayer systems for the first and last 100 ns was determined from the equation: P = 0.5 (3⟨cos2(θ)⟩ – 1). It was calculated for each consecutive bond in the CG lipids,[52] where θ denotes the angle between the bond vector connecting the consecutive beads and the bilayer normal. P values of 1, −0.5, and 0 indicate the perfect alignment, perfect anti-alignment, and random orientation, respectively. However, from this analysis, no significant changes were found, rather a stable conformation was observed.

The lateral mobility of the lipid molecules can be obtained with the help of their translational diffusion coefficients.[53] The diffusion coefficients of the upper leaflet lipids were estimated for the first and last 100 ns of simulation for both local and bulk POPG, and POPC lipids are given in Table . It was observed that the change in the lateral diffusion coefficient was significantly more for the local lipids than that of the bulk lipids. For simulation 2, we get a similar range of values (Table S3, Supporting Information). This is in agreement with the analysis of the density map[54] of local and bulk POPG and POPC for the first and last 100 ns of simulation, as given in Figure . The density map for the bulk lipid is given in the Supporting Information, Figure S11. To calculate the density map, first, local lipids were considered which were within 10 Å of the protein and the rest were considered to be bulk lipids. To see the changes in density, the densities of the local and bulk lipids for the first 100 ns and the last 100 ns of the trajectory for both POPG and POPC were calculated. Moreover, when the diffusion coefficients and densities of local POPG and POPC were compared, it was noticed that, in the case of POPG, the change in the values for the first and last 100 ns was more.

Table 1

Diffusion Coefficient (10–5 cm2/s)

lipid (upper leaflet)	first 100 ns	last 100 ns
POPG (local)	0.027 ± 0.011	0.003 ± 0.007
POPC (local)	0.046 ± 0.008	0.014 ± 0.004
POPG (bulk)	0.037 ± 0.002	0.037 ± 0.003
POPC (bulk)	0.038 ± 0.005	0.041 ± 0.006

Figure 6

Density map: (a) POPG (local) for first 100 ns, (b) POPG (local) for last 100 ns, (c) POPC (local) for first 100 ns, and (d) POPC (local) for last 100 ns.

One can calculate the radial distribution function (RDF), g, describing the probability of finding a particle of one type i, at a distance from another type j. RDFs of PO4 beads[55] of local POPG and POPC (Figure a,b) were calculated. Simulation 2 also provided a similar observation (Figure S12, Supporting Information). It can be seen from the figures that the local POPG lipids have aggregated more than that of local POPC. This indicates that the binding of the aggregated protein caused the disruption of the host POPG/POPC bilayer, as also has been reported earlier about the disarrangement of the phospholipid bilayer by the aggregated protein.[18]

Figure 7

RDF: (a) POPG local and (b) POPC local.

Through the analysis of bilayer thickness, the order parameter suggested that the structural properties of the mixed bilayer did not change to a large extent, and the number of local lipids with the simulation time plot (Figure ) suggested that the lysozyme pair had induced the local accumulation of lipids. Further analysis of the diffusion coefficient of both lipids (Table ) and density maps (Figure ) revealed that the interaction of the aggregated protein with local POPG was more than the others because for local POPG, we get more changes in the diffusion coefficient as well as in the density map.[4,56] Another interesting observation was the increase of the RDF of PO4[57] of local POPC and POPG of the mixed bilayer (Figure ), where it was observed that the close proximity was more for the last 1 μs than for the first 1 μs, which was also supported by the initiation of protein insertion. It provided a reliable insight about the aggregated protein and more probability distribution of PO4 beads of POPG than that of POPC around the protein and confirmed the higher affinity of POPG to the aggregated protein.

Initiation of the Insertion of Aggregated Proteins

The initiation of insertion is certainly an interesting topic to study with the help of MD as experimental tools to study the detailed motion and insertion of large biological molecules along with the aggregation phenomenon within the mixed bilayer interior are not readily available. Therefore, to get a clear description of the initiation of insertion of lysozymes along with the aggregation of lysozymes, one should consider about the interaction of the protein and the membrane across the whole bilayer. It was observed that the two lysozymes individually did not try to get inserted into the membrane, rather the whole system tried to interact with the membrane. During the initiation of insertion, the phospholipid group and the basic part of proteins came close to each other, which are also clearly observed from the VMD snapshots (Figure ).

Figure 8

Snapshot showing the electrostatic interaction and the initiation of insertion of the aggregated system with respect to time. Protein1, protein2, POPC, POPG, phosphate bead of lipids, and basic and acidic beads of the protein are presented as orange, green, purple, cyan, ocher, red, and blue van der Waals spheres, respectively. Image rendering was done with VMD.

From the information about the initiation of insertion of aggregated proteins, it may be concluded that there was a strong electrostatic interaction between the aggregated protein and the membrane. Figure also shows that at the interface of the mixed bilayer and the aggregated system, the basic residues of the protein and phosphate groups of the lipid were more close to each other than the acidic residues of the protein. This interaction helped the whole aggregated system for permeation in the mixed bilayer. This part may also be supported by another observation that, at the beginning of simulation, the lipid ratio of POPC and POPG was 7:3 (through the mixed bilayer) but, at 10 μs, the ratio became 1:1 (ratio of PO4 beads of POPC and POPG, which were within 10 Å of protein).

Conclusions

In this paper, we have used CG MD simulations to investigate the pathway of lysozyme aggregation in the presence of membrane and to monitor the changes in the properties of the mixed bilayer and lysozyme due to the interaction between them. Two simulations of 10 μs each have shown that the lysozymes first formed a pair and then approached the membrane as an aggregated entity. Though the length of the simulation is limited, it can provide a possible mechanism of membrane-induced denaturation-aggregation of lysozymes and the initiation of the process of membrane type binding by the aggregated system. The higher values of rmsd, Rg, and SASA of lysozyme2 than those of lysozyme1 indicate the significance of membrane on unfolding of proteins.

The structural/dynamic changes of the bilayer were recognized by analyzing different properties of the membrane, counting the number of local and bulk lipids, determining the diffusion coefficients of the local and upper leaflets of bulk lipids, plotting a density map of local and bulk lipids, and from the radial distribution of the phosphate group of the upper leaflet of local and bulk lipids. Moreover, initially, the ratio of local POPC and local POPG was 7:3 but, at 10 μs, this ratio became 1:1 (Figure ); all these observations indicated that POPG had a strong influence on lysozyme aggregation. At the same time, the initiation of lysozyme insertion into the mixed bilayer affected lipid ordering by generating a structural rearrangement of the bilayer. The combined results provide new biophysical insights into lysozyme aggregation in the presence of a mixed bilayer along with the effect of binding of the aggregated protein on the membrane structure and dynamics.