The effect of orientation dynamics in melittin as antimicrobial peptide in lipid bilayer calculated by free energy method

Melittin is a widely studied antimicrobial peptide which has gained interest because of its potential in therapeutic use. To perform its antimicrobial function, the melittin prefer to be in a certain state of membrane. In this study, we simulate the melittin in horizontal orientation towards the surface of membrane and in vertical orientation inside transmembrane state. The free energies are then calculated to asses the most favorable orientation, indicating its stable structure, of this peptide interacting with the lipid bilayer. Umbrella sampling with and without pulling simulations is carried out to obtain the sampling windows which used in observing the free energy changes along the reaction coordinate. From our simulations, we found that melittin in vertical orientation in transmembrane state has smaller free energy value than the horizontal one indicating that melittin is more favorable in this configuration. Our results show that melittin undergoes reorientation process from horizontal to vertical to perform its antimicrobial function.


Introduction
Antimicrobial peptides (AMPs) are important and effective components in innate host defenses against infectious pathogens. Simply described, innate immunity refers to all aspects of protection against infection that are genetically predetermined and do not require prior exposure to the foreign potential pathogens [1]. These peptides are classified into two types correspond to their cell selectivity: those that are toxic for bacterial cells but not mammalian cells and those that are toxic to both bacterial and mammalian cells.
Melittin is one of the most extensively studied of antimicrobial peptide (AMP). According to Huang [2], it is believed that the antimicrobial action of melittin involves forming a toroidal pore in the bacterial membrane that eventually leads to the cell death. This action can be occurred influenced by several factors, such as the concentration, the orientation of the peptides towards the surface of the lipid bilayer and the membrane state in which the peptide presence. In addition, the material properties of the bilayer is also expected to play a significant role in this formation process [3].
For melittin to insert into a pore and become a part of this pore, it has to be in transmembrane state [4]. Irudayam in [5] investigated the forming of this pore by monitoring the structural changes in the POPC bilayer and analyzed its free energy. It is observed that reorientation of melittin structure has to occurred from the surface state to a transmembrane state in order to create a pore in the lipid bilayer. However, the mechanism of how melittin reaches this state remains ambiguous. Thus, many studies have been carried out to investigate the interaction of melittin in membrane system. The most frequently used method to study this interaction is through the use of molecular dynamics (MD) simulation.
Here, we attempt to analyze which conformation of melittin showing the most stable structure when interacts with lipid bilayer. Thus, in this study, we perform MD simulation of melittin peptide in the mixed 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC):1-palmitoyl-2oleoyl-sn-glycero-3-phosphoethanolamine (POPE) lipid bilayer, in regard of different orientation towards membrane state. We then calculate the free energy of these conformations and analyze their stable structure based on these values.

Methods
We carried out MD simulations of melittin in different orientation towards lipid bilayer with a certain membrane state. There are two main models of melittin have been constructed; I) the structure in horizontal orientation is placed close to the surface of bilayer, and II) in vertical orientation towards the bilayer embedded with its C-terminus residues anchors to the headgroup and N-terminus is inserted into the core of bilayer.
I a) Melittin in horizontal orientation on the surface of bilayer: umbrella sampling with pulling simulation. To estimate free energy changes of melittin from the water then adsorbed into the surface of lipid bilayer, umbrella sampling with pulling simulation was performed to get is conformation samplings. Firstly, crystal structure of melittin was placed 35Å away from the bilayer center, along the z direction. We pulled the structure as the function of the distance from the center of melittin to the center of bilayer. A 500 ps pulling simulation was run with the force constant for umbrella potential of 200 kJ/mol. The interval distance for each window sampling was about ∼0.2 nm. Next, a 10 ns umbrella simulation with protein restraint was carried out for each window samplings. I b) Regular MD simulation of melittin in horizontal orientation on the surface of bilayer. To validate the simulation performed using umbrella sampling with pulling simulation, we conducted a regular MD simulation by placing the crystal structure of melittin 35Å away from the bilayer center on the upper leaflet, along the z direction. The conformation was oriented horizontal towards the x -axis and rotated by 60 • along the x -axis. Such setting is considered in order to set the hydrophobic side faces the water and the hydrophilic side faces the lipid headgroup region, as suggested by Irudayam [4]. We then run the simulation of 10 ns for this system under NPzT configuration. The structure of the last 10 ns simulation was taken and further equilibrated for 1 ns exploring its conformations for umbrella sampling simulation.
II a) The melittin in vertical orientation inside the transmembrane state: umbrella sampling with pulling simulation. In this simulation, free energy changes were observed as the function of the distance between the center of mass of the first three residues of the N-terminus and the center of bilayer. During pulling simulation, these first three residues of the N-terminus was acting as pulling group, while the C-terminus being restrained. The pulling simulation was performed for 500 ps with the force constant of 200 kJ/mol. The interval distance of 0.2 nm was applied. However, the number of sampling window which fulfilled this interval was only three coordinates, not a sufficient sampling to calculate the free energy changes . Hence, we took the additional coordinates for every 50 ps. Afterwards, umbrella simulation was carried out for each window sampling for 10 ns.
II b) Regular MD simulation of melittin in vertical orientation embedded in transmembrane state. The peptide was oriented along the z -axis and rotated for 90 • on the y-axis. This setting was applied such that the C-terminus of melittin anchors to the headgroup atoms and the N- terminus points to the core of bilayer. This orientation would induce much less perturbation as the C-terminus has a charge of +4 and the N-terminus has a charge of +1, therefore, it requires less energy [4]. This simulation were then run for 10 ns. Also, the last coordinate of this simulation which showed the stability of the peptide was taken and further equilibrated to obtain the sampling window for umbrella sampling. Simulation details: The crystal structure of melittin was retrieved from NMR structure from protein data bank (PDB) with PDB ID of 2MLT [6]. The coordinate of chain A of the crystal structure was used in construction all models. The lipid bilayer consisted of 64 POPC and 64 POPE lipid molecules randomly placed on the upper and lower leaflets. This lipid bilayer structure was taken from the pre-equilibrated of 50 ns performed by Hub [7], while the bilayer topologies was downloaded from Tieleman's site [8]. The GROMOS96 43a1 [9] and Berger [10] lipid were used as force field parameters . This peptide has a net charge of +6 because of its Nterminus is protonated by amine group and hence 6 Cl -was added to neutralize the system. The SPC [11] water model was adopted to describe water molecules. The temperature and pressure were maintained at 310 K and 1 bar using Nose-Hoover thermostat [12,13] and semi-isotropic pressure coupling with Parinello-Rahman barostat [14]. The temperature and pressure coupling constants were set to 0.5 ps. The simulation time step was 2 fs. The electrostatic interaction was calculated using Particle Mesh Ewald [15] method. We used LINCS [16] algorithm to constrain covalent bond lengths. All of simulations were performed using GROMACS [17] package: version 4.6.5. To extract the potentials of mean force (PMF) values, the weighted histogram analysis method (WHAM) which included in GROMACS as g wham [18] utility was used. The root mean square deviation and fluctuation were calculated using g rms and g rmsf scripts of GROMACS.

Results and Discussion
We discuss the results of our study follow the same sequence of simulations as in Methods section.
I a) To estimate the free energy changes of melittin in horizontal orientation from the water adsorbed into the surface of lipid bilayer on upper leaflet, we performed umbrella sampling with pulling simulation. This simulation was independently carried out from the regular melittin simulation. Here, as in the regular simulation of melittin in horizontal orientation, the peptide was placed 35Å away from bilayer center. During the pulling simulation of 500 ps, we obtained the initial coordinates for umbrella simulation with the interval of ∼0.2 nm between each of coordinate sampling.   Figure 1 shows the position of melittin at difference distances from the bilayer. The peptide retains its helical structure above ∼2 nm while below this distance it tends to adsorb into the surface of bilayer and adopted a "U-shaped" conformations. We also discover that at the distance of 1.2 nm, the three residues around the center (LEU9, THR10 and THR11) having structural changes from the helix to a turn forms. The free energy change of melittin in horizontal orientation is presented in PMF curve in figure  2. The PMF curve is plotted after discarding the first 3 ns simulation of umbrella samplings as we started the simulation from a non-equilibrium structures. From figure 2, we can see that the melittin reaches its minimum energy around the distance of 3 nm with the value of -7.78 kcal/mol. When the peptide starts to attach to the bilayer surface (∼2.4 nm) even adsorbed into it (less than ∼2 nm), the free energies tend to increase and move toward 0 kcal/mol. While, if the peptide was located near to the surface, around 2.5-3.0 nm away from the bilayer center, the PMF curve shows the free energies tend to reach the local minimum. These results indicate that the favorable position of melittin in horizontal orientation is close to the surface of bilayer.
Irudayam and Berkowitz [4] conducted the umbrella sampling with pulling simulation to estimate the free energy changes of melittin binding from the water to the bilayer surface in parallel orientation. They obtained the PMF value of -10.66 kcal/mol, showing the strong affinity of peptide to the lipid bilayer.
. In this simulation, we placed the melittin at the distance of 35Å from the bilayer center on the upper leaflet, along the z direction. The hydrophobic residues of the peptide were facing the water molecules and the hydrophilic residues faced the lipid headgroup. The snapshot of the system taken at 10 ns simulation are shown in figure 3.   As shown in figure 3, the four residues of C-terminus and two residues of N-terminus start to unfold and extend on the surface of bilayer as the peptide getting closer to the lipid head. This structural changes are detected on the two residues of C-terminus (ARG24 and GLN25) within the first 150 ps simulation, followed by the N-terminus around 600 ps. Almost of all residues maintain their structure in alpha-helix during simulation. Also, about 5 residues having structural changes by forming coil, bend, and turn structures.
We further analyze the fluctuation of each residue by calculating the RMSF value after the system become stable about after 4ns (see figure 4(a)). Based on the RMSF value, GLY1, ILE2, TRP19, LYS23 and ARG24 have fluctuation more than 1.5Å. Indicating that these residues are more fluctuate than the others during simulation. We observe that as our melittin conformation was placed such that the hydrophobic side faces the water molecule while its hydrophilic faced the lipid headgroup, it prefers to keep this conformation setting as effect of hydrophobic and hydrophilic interaction. The distance of melittin towards the center of bilayer is 27.68Å at the end of simulation. This structural changes is also analyzed by using RMSD calculation. For this simulation, average of RMSD is 2.34Å, showing no large structural changes of melittin occurred in the system. The RMSD curve is shown in figure 4(b). Figure 4(c) shows the PMF curve of the free energy changes for melittin in horizontal orientation. As the used conformations in umbrella sampling are taken from an equilibrated structures, all of trajectories can be used in calculation of free energies. As seen from this figure, the melittin reaches the stability around 2.73 nm away from the bilayer center with the value of -0.56 kcal/mol. Although this result is quite larger than the one of umbrella sampling with pulling simulation, the range of distance in which the peptide obtains the energy minimum is in agreement, which is 2.6-2.85 nm.
II a) We also performed umbrella sampling with pulling simulation to approximate the free energy changes of melittin in transmembrane state. The peptide in vertical orientation with its C-terminus residues points to the lipid headgroup while the N-terminus faces to the core of bilayer. In this simulation, we observe the free energy changes as the function of the distance of the first three residues of N-terminus to the bilayer center. Figure 5. Snapshots of melittin at different distance of the N-terminus residues from the bilayer center in pulling simulation. The peptide is colored by the residue type, as in Figure 1.
The snapshots of melittin in vertical orientation is displayed in figure 5. In this 500 ps of pulling simulation, the three residues near the center of melittin (LEU9, THR10 and THR11) are having the structural changes almost over the simulation time. The THR11 residue starts to unfold and extends to form a turn at 50 ps. Around 100 ps, it folds again into a helix structure. Starting from this 100 ps, the peptide conformations begin to slightly moving along y-axis from these three residues to the N-terminus yielding a "U-shaped" conformations. As these residues tend to bend near the center of melittin, the N-terminus residues are therefore being pulled to the center of bilayer. We also discover that these three residues frequently transform into a coil, a turn and helix shapes in a short range of simulation time.  Snapshots of melittin at different distance of the N-terminus residues from the bilayer center in pulling simulation. The peptide is colored by the residue type, as in Figure 1. Figure 6 shows PMF curve of melittin in vertical orientation in transmembrane state. The trajectories of the first 3 ns simulation are discarded as equilibration process. According to this figure, melittin obtains its minimum energy at the distance of 1.29 nm with the value of -11.38 kcal/mol, showing a strong affinity of the peptide to the lipid bilayer. When the distance of the first three residues of the N-terminus get closer to bilayer center (less than ∼0.9 nm), the free energies tend to increase and move toward 0 kcal/mol. Whereas, the free energies tend to stable if the distance around 0.9-1.3 nm from the bilayer center. By observing these results, we can state that as the distance of the first three residues of N-terminus getting closer to the bilayer center affected by the bending of the three residues near the center of peptide ("Ushaped" conformations) is not the favorable position for the melittin in transmembrane state. In contrast, the peptide finds its preferable position when the structures are more like a "bent-rod" shape with the bending angle about 140 • -160 • .
II b). The second simulation is melittin in vertical orientation inserted into the core of bilayer. Here, the C-terminus was placed points to the headgroup of bilayer. During the 10 ns simulation, we observe that there is no large structural changes developed on the melittin. Instead of clearly unfolded like in previous simulation, the helix structure of melittin is retained for most of simulation. There are three residues having seen transform into coil and turn forms. We expect that these three residue are LEU9, THR10 and THR11 which located near to the center of bilayer. These structural changes thus validate the previous pulling simulation which also shows similar transformation. Moreover, as these residues are hydrophilic, its interaction with the hydrophobic lipid tails causes the shifting of its position to the inward and outward the center of bilayer, therefore affects the bending angle. The interactions of residues on the upper and lower leaflets also influence the bending near the bilayer center. The unfolded on N-terminus starts to developed around 5 ns. The surface of lipid bilayer on the upper leaflet is pulled by the C-terminus residues such that the N-terminus getting closer the lipid headgroup in lower leaflet (see figure 7). The fluctuation of each residue is shown in figure 8(a). We calculate the RMSF value after the peptide reach the stability after about 4 ns. From this figure, the ILE2 and GLN26 have more fluctuation than the other residues, which is more than 1.5Å. But overall residues are less fluctuate during simulation. We also calculate the RMSD value for this simulation which is shown in figure 8(b), with the average of RMSD of 1.44Å indicating that the overall dynamics structure are close to the native structure without large structural changes.

Conclusions
We have performed MD simulations of melittin in different orientation towards lipid bilayer. The melittin is placed in horizontal orientation close to the surface of bilayer and in vertical orientation in transmembrane state. Two types simulations, umbrella sampling with and without pulling simulations, are conducted to observe the favorable orientation, indicating the stable structure of melittin interacting with lipid bilayer. We found that overall residues of melittin are less fluctuate, thus leading to no large structural changes occurred during simulations. We also calculate the free energy changes of peptide by integrating the PMF. We discovered that melittin in vertical orientation in transmembrane state has smaller value of free energy than the one in horizontal orientation on the surface state, both in umbrella sampling with and without pulling simulations. These results mean that the melittin is more preferable in vertical orientation in transmembrane state. We therefore can conclude that the melittin undergoes the rorientation process from horizontal orientation on the surface to the vertical one within the membrane to perform its antimicrobial action.