Simulations reveal that antimicrobial BP100 induces local membrane thinning, slows lipid dynamics and favors water penetration.

BP100, a short antimicrobial peptide, produces membrane perturbations that depend on lipid structure and charge, salts presence, and peptide/lipid molar ratios. As membrane perturbation mechanisms are not fully understood, the atomic scale nature of peptide/membrane interactions requires a close-up view analysis. Molecular Dynamics (MD) simulations are valuable tools for describing molecular interactions at the atomic level. Here, we use MD simulations to investigate alterations in membrane properties consequent to BP100 binding to zwitterionic and anionic model membranes. We focused on membrane property changes upon peptide binding, namely membrane thickness, order parameters, surface curvature, lipid lateral diffusion and membrane hydration. In agreement with experimental results, our simulations showed that, when buried into the membrane, BP100 causes a decrease in lipid lateral diffusion and lipid acyl-chain order parameters and sharp local membrane thinning. These effects were most pronounced on the closest lipids in direct contact with the membrane-bound peptide. In DPPG and anionic-aggregate-containing DPPC/DPPG membranes, peptide flip (rotation of its non-polar facet towards the membrane interior) induced marked negative membrane curvature and enhanced the water residence half-life time in the lipid hydrophobic core and transmembrane water transport in the direction of the peptide. These results further elucidate the consequences of the initial interaction of cationic alpha-helical antimicrobial peptides with membranes. This journal is © The Royal Society of Chemistry.

Introduction

According to the World Health Organization, “antibiotic resistance is one of the biggest threats to global health, food security, and development today”.[1] Particularly, Gram-negative bacteria are threatening, as they are naturally more resistant to antibiotics, mainly due to the extra protection provided by an outer-membrane formed mainly by negatively charged lipids and lipopolysaccharides.[2,3] Antimicrobial peptides (AMPs) are active against antibiotic-resistant bacteria, even in biofilms.[4,5] AMPs are a broad class of molecules that destroy or inhibit microbes' growth by different mechanisms, particularly by disrupting the bacterial outer membrane.[4,6,7] This action mechanism is exciting and is effective against a large variety of Gram-negative bacteria.[8-10] Membrane disruption can occur by different pathways (barrel-stave, carpet, toroidal pore),[6] which may involve multiple peptides in cooperative processes.[11] Although AMPs have been widely studied, their therapeutic utilization is hindered by the lack of detailed understanding of their mechanism of action on membranes, a fundamental step in drug development.[4] FDA approved therapeutical use, mainly as topical medications, only a few of the more than three thousand known AMPs.[12]

BP100 (H-KKLFKKILKYL-NH2) is a short hybrid cationic alpha-helical AMP (CHAMP) designed by combining cecropin A and melittin, two antimicrobial peptides.[13,14] BP100 is highly selective towards Gram-negative bacteria, displays low minimal inhibitory concentrations and low cytotoxicity, making this peptide a potential candidate for drug development.[13-15] BP100 acts on the membrane, and its binding depends on the ratio of anion/zwitterionic lipid ratio of the membranes, being anionic lipid-rich membranes more prone to disruption than charge-neutral membranes.[16] Zeta potential measurements show that the increase of peptide/lipid ratio leads to the neutralization of the interfacial charge of large unilamellar vesicles and that the mechanism of membrane action also depends on the peptide/lipid ratio.[16] While BP100 is a random coil in bulk solution,[17,18] it folds into an amphipathic alpha-helix when bound to membranes, and the extent of alpha-helix formation when bound depends on the anionic lipid fraction in the vesicle.[16-18] Finally, BP100 can promote inner content leakage in unilamellar vesicles.[16,18] However, the peptide's action mechanism on membranes seems to differ in low and high peptide-content or anionic lipid-content regimes: a gradual leakage is observed below a certain peptide/lipid ratio threshold, or below a certain anionic lipid content, and an instantaneous leakage above these ratios.[16] Atomistic molecular dynamics simulations (MD) in water and membranes show that this peptide is a random coil in an aqueous solution but maintains an alpha-helix conformation in membranes containing negatively charged lipids.[19] The hydrophobic facet of the alpha helix is buried inside the membrane hydrophobic core, while the hydrophilic facet of the helix is exposed to the bulk solution.[19] These findings correspond with the reported experimental results, even though such systems are heavily dependent on force field choices.[20]

Simpler computational models can be used to evaluate the mechanism of membrane disruption by AMPs,[21-24] due to its reduced computational cost.[25] For instance, Brownian dynamics simulations of coarse-grained peptides in implicit membrane models[21] suggest that BP100 disrupt membranes via the carpet model. Although these models can describe the membrane's effect on the structure and dynamics of the AMPs adsorbed at the interface, a few exceptions[23] only captures the effect of AMPs on the membranes. However, due to the cooperative nature of the peptide action on membranes,[11,26,27] it is essential to correctly depict a monomer's effect on the bilayer, as the binding of one peptide will disturb the membrane, altering the adsorption of subsequent peptides.

Thus, the understanding of isolated AMPs effects on membranes is necessary before describing its mechanisms of action at higher peptide concentrations correctly.

Here we present the findings from atomistic MD simulations of a single BP100 in membranes composed of anionic, zwitterionic, and mixtures of both lipids. We employed a distance-based scheme for a layered analysis of the outcomes of BP100 binding on bilayers. Membrane thickness, lateral lipid diffusion, local membrane curvature, membrane hydration, and order parameters were investigated, and our results show more dramatic peptide effects on closer lipids from both monolayers. Also, these effects were dependent on the local composition of the membrane, i.e., more severe effects of BP100 on bilayers were observed when its local composition was richer in anionic lipids. This detailed description of BP100 effects in membranes of varied composition may provide the necessary information for future models to adequately describe the mechanism of action of BP100 on membranes at higher concentrations.

Methods

Simulation set-ups and details

We studied initially folded α-BP100 in lipid systems: in pure membranes (DPPC and DPPG), and mixed membranes DPPC : DPPG (50 : 50). All membranes contained 64 lipids on each leaflet. Two symmetrical mixed bilayers were studied: one with lipids randomly distributed in each monolayer (PCPG) and one with an initial aggregate of 16 DPPG lipids in the center of the monolayer with the remaining 16 DPPG and 32 DPPC lipids randomly distributed (PCPG*) (Table S1†).

All-atom force fields Amberff99sb-ildn-NMR[28] and SLipids[29,30] were chosen for peptide and lipids were respectively after testing and validation.[19] Pure membranes topologies (DPPC and DPPG) were from the SLipids developer group website,[31] and the mixed bilayers were assembled with PACKMOL[32] software. TIP3P water molecules solvated our systems reaching an approximately 53 water/lipid ratio. A single BP100 was initially positioned parallel to the membrane surface approximately 2 nm away. Aqvist[33] Na+ ions and Cl− ions (Dang[34]) were introduced to counter-balance charges from BP100 and DPPG.

To integrate the equations of motion, we used a Leap-Frog integrator[35] with a 2 fs time-step. All simulations were performed in the NPT ensemble in which the temperature was kept constant at 323 K with the V-rescale thermostat[36] coupling separately peptide, lipids, water, and counter-ions with a coupling constant of 0.1 ps. The Berendsen barostat[37] was used to keep the pressure constant with semi-isotropic pressure coupling at 1 bar with a coupling constant of 1 ps and compressibility of 4.5 × 10−5 bar. LINCS algorithm[38] was used to constrain all bonds. For long-range electrostatic correction Particle–mesh–Ewald method,[39] with a real-space cut-off of 1.5 nm, was applied. van der Waals interactions were truncated at a distance of 1.5 nm with a switch function from 1.4 nm. Each set-up was simulated for a total time of τtot = 2 μs, after an initial equilibration step.

To analyze the influence of BP100 on local membrane properties, we performed our lipid analysis in layers, taking into account the gradual distance to the peptide. Lipids were ranked according to the distance between their phosphorus atoms and the peptide atoms, using GROMACS trjorder command. For membrane thickness (DHH) and surface curvature angle distribution (θdistr), the first 10, 30 lipids closest to BP100 and the entire monolayer were investigated (L10, L30, L64, Fig. 1C). For deuterium order parameters (SCD) and lipid lateral diffusion coefficients (DL) L10, L20, which excludes L10, and L64 were considered (Fig. 1D). Each monolayer was analyzed separately (Fig. 1), and for comparison, pure bilayer simulations were used as control. Due to peptide and lateral lipid movement during the simulation, lipid group indexing was updated every 10 ns (Fig. 2). This interval was chosen after simulation trajectories visual inspection.

Fig. 1

Lipid indexing for local and overall membrane analysis. (a) BP100 molecular structure. (b) Lipid indexing in relation to its phosphorus distance to BP100. Both monolayers were analyzed separately (c). For all analysis, L10 represents the first 10 closest lipids to BP100, and L64, the entire monolayer. For membrane-shape-dependent properties, such as membrane thickness and surface curvature angle distribution, SuAVE analysis package was utilized and L30 represents the first 30 closest lipids to BP100 (d). For deuterium order-parameters and lipid lateral diffusion coefficients, L20 groups the first 20 closest lipids to BP100 excluding the first 10 lipids (e).

Fig. 2

Lipid indexing update and calculation of membrane properties. Lipid indexing for lipid groups (L10, L20, and L30) were updated on every 10 ns to ensure accurate lipid selection. All membrane properties (SCD, DHH, θdistr, and DL) were calculated in 10 ns windows. Then, 20 average and standard error values were obtained over blocs of 100 ns each.

Two hundred windows of 10 ns each was analyzed, considering the total simulation time τtot of 2000 ns. Then, they were averaged in blocks of 10 windows for all calculated lipid properties, resulting in 20 average values of 100 ns each with their respective standard error. This scheme is illustrated in Fig. 2.

Membrane thickness (DHH) and surface curvature angle distribution (θdistr)

Accurate calculation of membrane thickness (DHH) and surface curvature (θdistr) requires proper consideration of the bilayer shape. This calculation was achieved by using SuAVE software package,[40] which fits a rectangular grid mesh (Fig. 3) based on the location of a selection of atoms (i.e., phosphorus from the lipids phosphate group) and then calculates the average distance between the upper and lower grids for DHH (Fig. 3A) and the distribution of the angles (θdistr, Fig. 3B) between the z-axis and the normal vector of each rectangular grid for surface curvature analysis. Higher angles indicate the membrane is curved.

Fig. 3

Membrane thickness and surface curvature angle distribution analysis using SuAVE. SuAVE software generates a grid (a) onto the position of selected atoms of the membrane (in our case, phosphorus atoms, in cyan) considering a specific number of rectangular partition bins along the x and y axes, chosen by the user. Surface curvature is analyzed by measuring the angle between the z-axis and the normal vector of each surface rectangular grid (b) and an angle distribution is generated. Higher angles indicate membrane curvature. Membrane thickness (b) was obtained by calculating the average distance between both upper and lower grids.

Using SuAVE,[40] membrane surface grids for L10, L30, and L64 were obtained using phosphorus atoms positions at every 10 ns, and their respective DHH and θdistr values were computed.

Deuterium order parameters (SCD) and lipid lateral diffusion coefficient (DL)

Deuterium order-parameters (SCD) and lipid lateral diffusion coefficients (DL) were obtained using GROMACS 5.0.2.[41-43] The local order of a lipid bilayer can be described by 2H-NMR deuterium order parameters (SCD) as SCD = (1/2)〈3 cos2 θ − 1〉, where θ is the angle between vector carbon–deuterium (actually C–H bond vectors in our simulations) and the membrane normal, whereas the brackets indicate the average during the simulation.

|SCD| values for the sn1 acyl chain were calculated for all simulations and lipid groups (L10, L20, L64) in both leaflets. The sn1 acyl chain was selected as a representative sample of lipid order parameters.

The lipid lateral diffusion coefficient (DL) can be obtained from MSD = 4DLt, where MSD is the mean square displacement of an atom as a function of time, t. MSD of the phosphorus atoms from each lipid group were calculated at every 10 ns. Then, DL was obtained from the slope of the least square fitting between 2 to 5 ns of each 10 ns-time-window.

Membrane hydration

Simulated membranes were sectioned along the z-axis in 3 regions (R1, R2 and R3) having the center of the bilayer as the origin. R1 is the central region in the hydrophobic core of the bilayer delimited by planes z = +0.3 and −0.3 nm while R2 is delimited by z = +0.5 and −0.5 nm, which defines the membrane hydrophobic core, or tail-to-tail region. R3 encloses the region between the upper and lower surface grids generated by the SuAVE package, adjusted over the phosphorus atom positions on both leaflets (Fig. 3B). Therefore, this R3 region defines the membrane/water interface. Additionally, to better identify the peptide effects in the membranes, these 3 regions were divided in sub-regions: (i) upper that has positive z-coordinates and is closer to the peptide and (ii) lower that has negative z-coordinates. The overall membrane hydration was investigated by computing the water density profile across the z-axis over the last 100 ns of the simulation trajectories and by computing the average number of water molecules in R3. Water penetration into the membrane hydrophobic core was investigated by computing the number of water molecules in R2 and R1 regions, their z-coordinate probability distribution and the residence time also in R2 and R1 regions. In addition, the average number of water molecules, 〈Nw〉, and the probability of finding water molecules, PW, inside the membrane hydrophobic core were computed using a time resolution of Δtr = 10 ps. Therefore, considering the whole simulation time τtot, the total number of trajectory frames analyzed was NTF = τtot/Δtr, in which only some frames have N water molecules inside the hydrophobic core, NWF(N), where N ranges from 1 to 4. Then, the total number of frames with at least one water molecule was obtained as NWF = ΣNWF(N). Considering only these NWF frames, the average number of water molecules was calculated as 〈Nw〉 = Σ(N × NWF(N))/NWF. However, considering all NTF frames, the frequency or the probability of finding N water molecules in the hydrophobic core was calculated as PW(N) = NWF(N)/NTF and the probability of finding any amount of water molecules in the hydrophobic core was calculated as PW = NWF/NTF.

Results

We analyzed four simulation sets that run for τtot = 2 μs each. BP100 was initially folded as an alpha-helix and positioned approximately 2 nm away from the membranes (DPPC, DPPG, PCPG and PCPG*). Control simulations of pure membranes gave the expected area per lipid and membrane thickness.[44-48] Peptide secondary structure analysis and the average number of lipids in contact with BP100 were investigated previously.[19] Moreover, in the peptide/membranes simulations, BP100 approached the membrane with its positively charged facet and rotated, facing its non-polar residues to the membrane core, while maintaining its alpha-helical structure (Fig. S2B†). This motion, which we previously defined as peptide flip (Fig. S3†), is accompanied by peptide dehydration and an increase in non-polar contacts between peptide and membrane.[19] However, in DPPC, although the peptide flipped, its alpha-helical structure was lost (Table S4†), and therefore we will refer to it as a semi peptide-flip. The following sections present our focus on the peptide effects on the membrane.

Membrane thickness

All trajectories showed BP100 adsorbed onto the membranes and continued embedded throughout the simulations. We observed peptide flip in the BP100 on DPPC (∼700 ns), DPPG (∼1500 ns), and PCPG* (∼390 ns) simulations.[19]Table 1 presents the values of membrane thickness averaged over the entire simulation, and control shows DHH values for membrane-only simulations. In all simulations, compared to the overall membrane thickness (L64), L10 shows lower values, indicating bilayer thinning where BP100 was bound. The overall membrane thickness values (L64) take into account all monolayer lipids; thus, the actual thickness difference between L10 and the rest of the membrane was higher in all cases. L10 thickness reduction was detected for all simulations, and it was significantly higher in membranes with an anionic lipid aggregate (Table 1 and Fig. S5†). L30, L64, and control membranes shared similar thicknesses, indicating BP100 thinning activity was limited to the first near lipids, close to the peptide.

BP100 induces membrane thinning in anionic bilayers. Average membrane thickness (DHH) and standard deviation obtained for all simulations

	Membrane thickness (nm)
	L10	L30	L64	Pure membrane
BP100 in DPPC	3.71 ± 0.08	3.75 ± 0.04	3.77 ± 0.03	3.79 ± 0.03
BP100 in DPPG	3.50 ± 0.12	3.67 ± 0.05	3.72 ± 0.03	3.70 ± 0.05
BP100 in PCPG*	3.57 ± 0.10	3.77 ± 0.05	3.83 ± 0.04	3.82 ± 0.04
BP100 in PCPG	3.72 ± 0.10	3.82 ± 0.05	3.84 ± 0.02	3.83 ± 0.04

Simulations with peptide flip and anionic membranes, namely BP100 in DPPG and PCPG*, showed the lowest membrane thicknesses in L10, with averaged values of 3.50 ± 0.12 and 3.57 ± 0.10, respectively (Table 1). In the BP100 in DPPC simulation, despite peptide flip, the peptide lost some of the alpha-helical conformation (45% of helicity, Table S4†).

An anionic lipid aggregate, either as a membrane patch (PCPG*) or in the whole membrane (DPPG), seemed essential for the peptide-induced membrane thinning (Table 1). This can also be visualized by comparing the obtained thicknesses in L10 for BP100 in PCPG and PCPG*. While both simulations share the same membrane composition and only differ in their lipid distribution, PCPG* shows an average thickness of 3.57 ± 0.10 nm in L10 and PCPG, 3.72 ± 0.10 nm (Table 1).

According to the peptide distance, the gradual membrane thinning can be visualized through number density graphs (Fig. S6 and S7†). Number density graphs of before and after peptide flip revealed the distance shortening between upper-L10 and lower-L10, particularly in BP100 in DPPG (Fig. S6B†) and BP100 in PCPG* (Fig. S7A†) simulations. Although the upper-L10 position was affected by BP100 the most, shown by the curve broadening of upper-L10 signal, lower-L10 is also pulled towards the peptide (Fig. S6B† and S7A), demonstrating that BP100 binding effect has an impact on the opposite bilayer leaflet. Fig. 4 shows membrane thickness and peptide density in 2D mapping obtained from the BP100 in DPPG simulation, where membrane thinning was the most significant (Table 1). A simulation snapshot at 1600 ns is shown in Fig. 4C. A negative curvature was visible by overlaying the BP100 position with SuAVE grids generated based on the DPPG phosphorus atoms positions (Fig. 4C). The overlay of 2D membrane thickness and peptide density mappings gives information on the peptide position and its outcomes on membrane thickness. In DPPC, semi peptide flip showed mild effect (Fig. S6A†) on membrane thickness, while in BP100 in DPPG and PCPG* simulations, a valley on the membrane matching the peptide position is evident (Fig. 4 and S8C†). We confirmed that peptide flip increased membrane thinning by analyzing the data before and after the flip in DPPG (Fig. S8B†). The membrane was affected by peptide adsorption (Fig. S8B† – before), but peptide flip made the membrane thinner, decreasing the membrane thickness of the first surrounding lipids from 3.2 nm to 3.0 nm, on average (Fig. S8B† – after). In PCPG*, we observed similar results (Fig. S8C†) as in DPPG, with membrane thickness decrease after the flip. In PCPG (Fig. S8D†), where no flip occurred, no change was detected, comparing the simulation beginning and end.

Fig. 4

BP100 induces membrane thinning in anionic bilayers. Membrane thickness map for BP100 in DPPG (a) and peptide density map (b), averaged from 1600 to 1700 ns of simulation, after peptide flip. BP100 induced membrane thinning is clearly visible by the superposition of both maps. In (c), a snapshot at 1600 ns; BP100 positively charged residues are colored in blue and the non-polar residues in green. The upper and lower grids were generated by SuAVE taking into account the positions of the DPPG phosphorus atoms. BP100 thinning activity is stronger on anionic-lipid-aggregate-containing membranes (PCPG* and DPPG).

Surface curvature

2D membrane thickness mappings indicated that the membrane surface was altered by BP100 binding and flip. To investigate these effects quantitatively, we calculated the distribution of the surface curvature angle θ, defined as the angle between the normal vector of the surface rectangular grid partitions and the z-axis. Larger θ values indicate greater membrane surface curvature. Fig. 5 shows the distribution of θ in different upper monolayer regions for BP100 in PCPG* simulation.

Fig. 5

Peptide flip induces negative membrane curvature. Membrane surface curvature angle (θdistr) distribution before (a) and after (b) the flip in the BP100 in PCPG* simulation. Higher angles indicate increased membrane curvature. Control curves represent θdistr from pure membrane simulations.

Comparison of the θ distribution for all upper lipid groups (L10, L30, L64) against control showed an increase in the percentage of lipids populating higher θ (20°–50°) (Fig. 5). These data show that the binding of a single peptide changed the overall membrane topology. The depression shown in 2D membrane thickness mappings (Fig. S8†) indicated BP100 caused negative curvature on anionic bilayers.

For BP100 in DPPG and PCPG* simulations, the surface angle distribution for upper-L10 and upper-L30 showed a most significant portion of θ populating lower values before the flip (Fig. 5A, S9B and S10A†). After peptide flip, the surface of upper-L10 showed a broader and flattened distribution curve, with an increase in the population of higher angle values (30°–50°), compared to the overall upper monolayer (upper-L64), and lower monolayer (lower-L64), demonstrating peptide flip caused local membrane curvature (Fig. 5B, S9A and S10B†).

In the BP100 in DPPC (Fig. S9A†) and BP100 in PCPG (Fig. S10B†), no remarkable shifts in angle distribution were detected. Semi-peptide flip was insufficient to alter upper-L10 surface angle distribution in the BP100 in DPPC simulation (Fig. S9A†). These data substantiate our findings in that not only anionic lipid content is essential for peptide activity on membranes, but its distribution is crucial as well.[19]

Lipid order parameters

Lipid chain order parameters were evaluated to investigate whether BP100 binding influenced lipid hydrophobic chain ordering and dynamics. Sn1 carbon chain of lipids (Fig. S11†) was taken as a representative of the overall lipid acyl chains. Fig. 6 shows lipid sn1 acyl chain order parameters of distinct membrane regions from BP100 in DPPG simulations before and after the simulation.

Fig. 6

BP100 decreases local lipid chain order parameter. Order parameters of sn1 DPPG acyl chain in the BP100 in DPPG simulation before (a) and after (b) the flip. The sn1 lipid chains were used for calculating lipid order parameters of the 10 and 20 closest lipids to BP100 the entire monolayer on both monolayers. Control sn1 SCD data was taken averaging over an entire peptide-free membrane simulation.

Even in earlier stages of peptide binding to the membrane (Fig. 6A), upper-L10 and lower-L10 lipids showed more disordered chains than control and other membrane regions. Nevertheless, with peptide flip, a sharp decrease in SCD was observed for upper-L10, showing the BP100 effect on its immediate neighboring lipids after flipping (Fig. 6B).

In all simulations, the SCD profiles of the upper-L10 and lower-L10 acyl chains were significantly lower compared to other regions of the upper monolayer (Fig. 6 and S12†). This behavior was observed in simulations with negatively charged bilayers (Fig. 6 and S12†), with both upper-L10 and lower-L10 having higher conformational freedom due to BP100 binding.

In simulations where peptide flip was observed, the decrease in upper-L10 order before and after the flip was apparent for all acyl chain carbons (Fig. 6 and S12†), even for BP100 in DPPC with a semi-flip (Fig. S12A†). Table 2 shows the averaged order parameters of all sn1 chain carbons in upper-L10. In simulations where peptide flip occurred, namely DPPC, DPPG, and PCPG*, a decrease of 9%, 23%, and 17% in 〈|SCD|〉 is observed after the flip, compared to control. In contrast, in PCPG albeit the presence of anionic lipids, we found a reduction of 5% in 〈|SCD|〉 compared to control (Table 2). The decrease in overall lipid acyl chain order parameters in upper-L10 can be explained by the insertion of BP100 hydrophobic residues into the bilayer core, increasing the acyl chain(s) freedom.

Averaged lipid sn1-chain order parameter, 〈|SCD|〉, for lipids in upper-L10, obtained from the initial and last 300 ns of simulation. Peptide flip was observed in DPPC, DPPG and PCPG* simulations, thus their 〈|SCD|〉 values from the last 300 ns of simulation are post-peptide flip. Control values are presented for comparison

	〈|SCD|〉 for upper-L10
	Initial 300 ns	Last 300 ns	Pure membrane
BP100 in DPPC	0.172 ± 0.019	0.159 ± 0.018	0.174 ± 0.038
BP100 in DPPG	0.151 ± 0.021	0.133 ± 0.021	0.172 ± 0.040
BP100 in PCPG*	0.162 ± 0.024	0.151 ± 0.018	0.182 ± 0.041
BP100 in PCPG	0.166 ± 0.024	0.170 ± 0.021	0.185 ± 0.040

In the BP100/DPPC simulations, semi peptide flip lowered the order parameters of the closest lipids in contact (upper-L10, Fig. S12A† and Table 2). However, the loss of the alpha-helical conformation led to less peptide penetration (Fig. S6A†) and limited insertion of hydrophobic residues into the bilayer. Lower peptide penetration explains the discrepancy between SCD profiles of other lipid groups and upper-L10 after the flip (Fig. S12A†). The BP100/PCPG order parameter profile (Fig. S12B†) showed upper-L10 and lower-L10 decrease in chain order late in the simulation. Presumably, the binding of BP100 produced a weak effect on the acyl-chain order parameters, as evidenced by the higher |SCD| values (Table 2) for all carbons of the acyl chain compared to simulations where peptide flip took place (DPPC, DPPG, and PCPG*). These data highlight the importance of anionic lipid aggregates for the peptide to maintain a helix conformation and, later, a full flip.

Lipid lateral diffusion

We investigated alterations in lipid dynamics upon peptide binding in the near and far vicinity of BP100. Lipid lateral diffusion coefficients for L10, L20, and L64 in both leaflets were calculated from the slopes of the mean-square displacements (MSD) in the xy-plane (Table 3).

Lipid lateral diffusion values of different lipid regions for all simulations (Fig. 1) and available experimental and theoretical data are also presented for comparison. Our simulation data reveal lipid lateral mobility is approximately 50% lower in the first surrounding of BP100 (upper-L10) compared to a peptide absence scenario (lower-L64) and membranes with higher anionic lipid content are more affected due to stronger electrostatic interactions and flip

	Lipid lateral diffusion (×10−7 cm2 s−1)
	Upper-L10	Upper-L20	Upper-L64	Pure membrane	Literature
BP100 in DPPC	0.9 ± 0.2	1.3 ± 0.2	1.6 ± 0.2	1.6 ± 0.3	1.8,[44]1.78,[47] ∼1.5 [45]
BP100 in DPPG	0.6 ± 0.2	1.1 ± 0.2	1.2 ± 0.2	1.3 ± 0.2	0.9 [48]
BP100 in PCPG*	0.6 ± 0.2	1.0 ± 0.2	1.2 ± 0.2	1.2 ± 0.2	—
BP100 in PCPG	0.7 ± 0.3	1.1 ± 0.3	1.2 ± 0.2	1.3 ± 0.2	—

The calculated lateral diffusion coefficients in upper-64 for DPPC and DPPG were 1.6 ± 0.2 and 1.2 ± 0.2, respectively. These values were well within experimental results and similar to our controls (Table 3). Systematic analysis of lipid dynamics considering the proximity to BP100 reveals an approximately 50% decrease of lipid lateral diffusion in upper-L10 compared to pure membrane systems, for all simulations (Table 3).

Similar to changes in membrane thickness, order parameters, and membrane curvature, BP100 effect on lipid lateral diffusion is more substantial on upper-L10, while upper-L20 and lower-L64 show similar values (Table 3). We also calculated the diffusion for the 5 closest lipids to BP100 and the results found are similar to upper-L10, showing that 1 : 10 peptide/lipid ratio is where BP100 shows the most capacity for membrane disturbing effects (Table S13†).

Alterations in membrane thickness, lipid tail order parameters, and lipid lateral diffusion caused by peptide binding could lead to changes in membrane hydration. And analyzing the behavior of water molecules in the membrane could point to the mechanism of action of CHAMPs in membranes. Membrane hydration was investigated by computing the water density profile along the z-axis (Fig. 7A and S14†) and water penetration into the inner membrane/water R3 region (Fig. 7B and S15†). Additionally, in the membrane hydrophobic core, it was investigated the amount of water molecules (Fig. 7C, D, S16 and S17†), their z-coordinate probability distribution (Fig. 8A, B, D–F and S18†) and their residence time (Fig. 8C, F and S18†).

Fig. 7

Water distribution in BP100/DPPG. (a) Number density profiles of BP100 in DPPG after peptide flip (solid lines) and pure DPPG (dashed lines), over 100 ns of simulation. The R1, R2 and R3 regions are also illustrated showing the upper and lower sub-regions. (b) Percentage of water content in upper and lower R3 sub-regions and percentage of peptide volume inserted in R3/upper sub-region. Number of water molecules, N, detected in (c) R2/upper and (d) R2/lower sub-regions before and after the peptide flip. Analyses performed using a time resolution of 10 ps.

Fig. 8

Membrane core hydration (for BP100/DPPG, BP100/PCPG* and controls) and average residence time of water molecules within R2 region. (a, b, d and e) z-coordinate histograms of the center of mass of the water molecules that have entered into the R2 region coming from the upper monolayer (in blue) or from the lower monolayer (in red), in the simulations of (a) BP100 in DPPG, (b) pure DPPG, (d) BP100 in PCPG* and (e) pure PCPG*. The insert in (a) and (b) illustrate the trajectory of one water molecule that crosses the membrane. In (c and f) the distribution of the residence time (symbols) of water molecules within R2 region (upper, lower and control) are shown for the same systems and the values of the residence half-life time (τ1/2) obtained from exponential fits (solid lines). Analyzes performed over the last 500 ns of simulations using a time resolution of 10 ps.

Fig. 7 shows data obtained from the BP100 in DPPG simulation. Similar water density profiles were obtained for both peptide membrane and pure membrane systems (Fig. 7A and S14†), which indicates that a single-peptide binding did not affect the membrane overall hydration. The average number of water molecules accessing the R3 region is similar for BP100/membranes (393 ± 31 for DPPC, 332 ± 18 for PCPG, 334 ± 19 for PCPG* and 309 ± 18 for DPPG) and its associated controls (390 ± 28, 341 ± 20, 336 ± 17 and 307 ± 18, respectively) (Fig. S15†). The water molecules in R3 region are symmetrically distributed in R3/upper (∼50%) and R3/lower sub-regions (Fig. 7B). However, significant differences between peptide/membrane and control simulations are seen in the water positioning across the leaflets and its dynamics in the R2 and R1 regions.

When computing the amount of water molecules within R2 (and R1) region for both peptide/membranes and pure membranes simulations, the large majority of trajectory frames have no water inside (N = 0), but for those frames that have water the most common is only one water (N = 1) and occasionally N = 2, 3 or 4 (Fig. 7C, D, S16, S17 and Tables S19–S21†). For example, for BP100/DPPG, the whole simulation of 2000 ns generates a total of NTF = 2 × 105 frames in intervals of 10 ps, but only a few frames have water molecules inside the hydrophobic core (in R2 region): NWF = 1215, that was obtained as a sum of all frames with N water molecules, NWF(1) = 1102, NWF(2) = 96, NWF(3) = 17 and NWF(>3) = 0 (Table S20†). Then, considering only the NWF frames, that have at least one water molecule, the average number of water was calculated as 〈Nw〉 = (1 × 1102 + 2 × 96 + 3 × 17)/1215 = 1.11. However, considering all NTF frames, the probability of finding water molecules in the hydrophobic core was calculated as PW = 1215/2 × 105 = 60.8 × 10−4 (Table 4). Therefore, finding water molecules in the hydrophobic core of the membranes is a rare event that happens a few times at the nanosecond time scale.

Probability of finding water molecules inside the hydrophobic core in the R2 upper and lower sub-regions, PW, in BP100/membranes and pure membranes. The modules of difference between upper and lower values is represented by ΔPW. The analyzes were performed: (i) over the whole simulations time for pure membrane (900 ns) and peptide/membrane systems (2000 ns), (ii) for the first 300 ns in parenthesis, and (iii) for the last 300 ns in brackets, for the peptide/membrane simulations, using a time resolution of 10 ps. For DPPG and PCPG* the first and the last 300 ns describe the situation before and after the peptide flip, respectively

		P w (×10−4)
		DPPC	DPPG	PCPG*	PCPG
Pure membrane	R2	33.7 ± 1.8	43.4 ± 1.7	25.9 ± 1.6	31.4 ± 1.6
	R2/upper	18.2 ± 1.8	21.1 ± 1.3	12.8 ± 1.1	15.9 ± 1.6
	R2/lower	15.4 ± 1.6	22.3 ± 1.0	13.1 ± 1.0	15.6 ± 1.1
	ΔPW	2.8 ± 3.4	1.2 ± 2.3	0.3 ± 2.1	0.3 ± 2.7
BP100 + membrane	R2	38.3 ± 1.4	60.8 ± 1.8	51.5 ± 2.3	32.7 ± 1.7
	R2/upper	20.1 ± 1.0	42.7 ± 1.3	36.2 ± 2.0	15.5 ± 1.0
	(First 300 ns)	(22.3 ± 1.1)	(21.7 ± 1.3)	(14.7 ± 1.8)	(19.3 ± 1.1)
	[Last 300 ns]	[24.3 ± 1.2]	[87.3 ± 2.2]	[50.7 ± 1.9]	[11.0 ± 0.7]
	R2/lower	18.2 ± 0.9	18.1 ± 0.8	15.3 ± 0.9	17.2 ± 1.0
	(First 300 ns)	(22.7 ± 1.0)	(23.7 ± 0.9)	(9.7 ± 0.7)	(24.7 ± 0.9)
	[Last 300 ns]	[20.7 ± 0.8]	[12.7 ± 0.8]	[18.0 ± 1.3]	[12.0 ± 1.1]
	ΔPW	1.9 ± 1.9	24.6 ± 2.1	20.9 ± 2.9	1.7 ± 2.0
	(First 300 ns)	(0.4 ± 2.1)	(2.0 ± 2.2)	(5.0 ± 2.5)	(5.4 ± 2.0)
	[Last 300 ns]	[3.6 ± 2.0]	[74.6 ± 3.0]	[32.7 ± 3.2]	[1.0 ± 1.8]

For all systems, in the average approximately one single water molecule penetrates in the hydrophobic core (R2 region) per water containing frame, 〈Nw〉 = 1.06 ± 0.03, 1.11 ± 0.04, 1.12 ± 0.05, 1.04 ± 0.04 water for BP100 in DPPC, DPPG, PCPG* and PCPG membranes, respectively, and 1.05 ± 0.03, 1.10 ± 0.04, 1.06 ± 0.03 and 1.03 ± 0.04 water for the respective pure membranes. Therefore, the total amount of water molecules within R2 (and R1) region is barely affected by peptide binding, but the frequency or probability in which water is detected changes with the peptide presence. Table 4 shows the values obtained for the probability of finding water molecules, PW, in upper and lower R2 sub-regions for both peptide/membrane (whole simulation of 2000 ns, before and after the peptide flip, i.e. the first and last 300 ns, respectively) and control systems. For the pure membranes, the DPPG membrane has a slightly larger probability of finding water in its hydrophobic core than the others, (43.4 ± 1.7) × 10−4. The difference between the lower and upper sub-regions is small (ΔPW < 2.8 × 10−4) and is in the same order of the uncertainty (±2.1–3.4 × 10−4), showing a symmetry or isotropy between both sides of the membrane, as expected. However for the peptide/membrane systems, a significant difference between the PW upper and lower was observed for DPPG, ΔPW = (74.6 ± 3.0) × 10−4, and for PCPG*, ΔPW = (32.7 ± 3.2) × 10−4, after the peptide flip, but not before the peptide-flip, (2.0 ± 2.2) and (5.0 ± 2.5) × 10−4, respectively. These results show an increase of approximately 4 times in the probability of finding water in the hydrophobic core of DPPG and PCPG* close to the peptide side (upper sub-region) after the peptide flip compared to pure membranes (from 21.1 to 87.3 × 10−4 for DPPG and 12.8 to 50.7 × 10−4 to for PCPG*). This increase is not a consequence of a greater number of water molecules entering into the membrane hydrophobic core, because most of the time only 1 water molecule was detected in the MD trajectory frames (Fig. 7C, D, S16, S17, Tables S20 and S21†). This higher probability is related to the higher residence time of water molecules in R2/upper region (Fig. 8).

The trajectory of each water that penetrates the hydrophobic core of the membranes was analyzed separately for molecules that come from the upper monolayer (where the peptide is bound) and from the lower monolayer, and the z-coordinate distributions of these water molecules were obtained and compared between the two sides where water enters. For simplification and easier comparison, the water inlet side was always plotted at the right side of the graphics in Fig. 8A, B, D and E for DPPG and PCPG* membranes and Fig. S18A, B, D and E† for DPPC and PCPG membranes. Then, the left side of these graphics quantifies the transmembrane water transport. As expected, for control systems (Fig. 8B and E), the water distributions are equivalent regardless of the entrance side, but they present a higher number of occurrences near this side due to some molecules that get in and get out of the hydrophobic region through the same side. The peptide interaction with the upper region of the membrane induces an increase in the distribution of water molecules in the hydrophobic core by increasing the probability of water detection in the upper region (Fig. 8A and D in blue). This effect is much larger for DPPG and PCPG* membranes. Additionally, the transmembrane water transport was analyzed by summing the distributions from the center of the membrane (z = 0) to the water outlet side (left side of the graphics, i.e. z = −0.5 nm for the blue distribution and z = 0.5 nm for the red distribution) (Fig. 8A and D and S18A–D†). For DPPG and PCPG*, the transport from lower to upper R2 region (Fig. 8A and D in red) is around 35% higher than those in the opposite direction (Fig. 8A and D in blue), and for DPPC and PCPG, no significant change was identified (Fig. S18†). This result reveals that the peptide binding and flip creates a small water flux through the membrane hydrophobic core in the direction of the peptide, going from lower to upper leaflet.

Fig. 8C, F and S18† show the water residence time distributions inside the membrane hydrophobic core, R2 region. The values of the half-life residence time (τ1/2) were obtained using best fit of these distributions with double exponential decay functions. For pure membranes, τ1/2 presents a small variation, 21.6 ps for DPPG, 23.7 ps for PCPG*, 28.0 ps for PCPG and 32.1 ps for DPPC. But for peptide/membranes, τ1/2 are much larger for the R2/upper than R2/lower after the peptide flip, 49.9 ps and 8.9 ps for DPPG and 36.3 ps and 10.3 ps for PCPG*, respectively. However, for PCPG and DPPC, the non-flipping or semi-flipping peptide does not induce significant change, τ1/2 are 23.5 and 27.6 ps, respectively. This result combined with the water distribution (Fig. 8A and D) indicates that the water molecules entering in R2/upper region reside longer due to peptide influence, increasing the probability of finding those single water molecules in the region. The BP100 in DPPG simulation shows the most expressive τ1/2 compared to the control (Fig. 8C, F and S18†), and the maximum water residence time in R2 region changes from approximately 350 ps in pure DPPG to more than 1000 ps in BP100/DPPG. These results refer to a single-peptide binding and higher peptide concentrations should show more pronounced effects in the membrane hydration.

Discussion

In our simulations, BP100 promoted membrane thinning, depending on the membrane's composition (Table 1, Fig. 4 and S8†). While marginal thinning was observed in membranes without PG agglomerates, more severe effects were observed in simulations where the peptide was in close contact with PG-enriched regions. Furthermore, our layered analysis showed that the peptide effect was more substantial to lipids closer to the peptide and that the effect vanished after a couple of lipid layers (after L20). These results show, simultaneously, the effect of BP100 on membranes and its dependence on membrane composition. Peptide-induced membrane thinning has been observed both experimentally in microscopy experiments,[49] X-ray difraction,[50] circular dichroism,[51,52] and NMR.[53-55] Molecular simulations of membrane/CHAMPs systems[56-58] also support the experimental observations.

Additionally, BP100 is known for its activity-dependence on membrane charge, and its disrupting effect increases with negatively charged membranes,[16] due to the electrostatic interactions with negatively charged groups of the membrane.[19] Finally, the considerable observed difference between the PC/PG mixture membranes with and without PG agglomerate points to the fundamental importance of lipid organization close to the BP100, i.e., the effect of the peptide on the membrane does not solely depend on the amount of PG content in the membrane, but also the local distribution of PG in the BP100 vicinities. We emphasize here that lipid clustering is reported in some CHAMPs/membrane systems.[8,59,60]

Using SuAVE analysis software,[40] which fits a surface at the membrane interface, a local and negative Gaussian membrane curvature was observed, induced by BP100 (Fig. 5). Once more, the lipids closer to the peptide were more disturbed by the peptide, as the angle distribution function for the L10 group is skewed towards higher angles, compared with other lipids in the membrane. Also, the peptide effect increased when peptide flip occurred. The L20 was also disturbed but not to the same extent observed for the L10 group. Similar to the thinning effect, no significant peptide effect was observed in DPPC membrane (Fig. S9A†) or PC/PG membrane without PG agglomerate (Fig. S10B†), highlighting the importance of electrostatic interactions between the positively charged groups of BP100 and the negatively charged groups of PG. CHAMPs-induced membrane curvature is already reported, both by experimental[61-63] and computational methods,[64-66] and it seems to be a necessary condition for membrane pore formation,[61,67,68] usually a collective process that involves several peptides. Our simulations show that a single BP100 can induce high local negative curvature, changing the local lipid packing and, consequently, the binding of additional peptides in this peptide-perturbed region of the membrane.

Along with membrane thinning and local negative curvature, BP100 promoted a significant decrease in the calculated order parameter, SCD (Fig. 6), throughout the lipids acyl chain (sn1), which corresponds to an increase in conformational freedom of the hydrophobic chains upon BP100 binding. In line with previous properties, this SCD decrease depended on the proximity between the lipid and the peptide, having closer lipids to BP100 (L10 group) higher conformational freedom than the other lipids in the membrane. Local lipid distribution modulated the effect of BP100 on the SCD of the lipids, as the observed reduction was greater in the simulation with (than without) the PG agglomerate. These findings corroborate the trends observed in experimental measurements of SCD of BP100-containing systems.[53,55] Our results showed that the orientation of the peptide with respect to the membrane changes the computed SCD values, as in DPPG and PCPG*, where a full flip took place, SCD profiles of upper-L10 and lower-L10 (Fig. S12† and Table 2) were significantly lower than other lipid groups, their respective controls, and other simulations (Table 2).

Lateral diffusion of lipids was analyzed from our simulations, as it constitutes an essential component of the membrane dynamics and cell functions. Lipid lateral diffusion alteration upon peptide binding has been studied by several peptides experimentally[69-73] and theoretically.[74,75] For example, the neurotoxic Alzheimer's disease peptide amyloid-beta increases lipid diffusion as shown by neutron scattering experiments.[71,73] Alpha helical and/or cationic peptides such as gramicidin, mellitin, and cWFW slow down lipids.[69,70,76] In our simulations, BP100 decreased the upper-L10 lipid lateral diffusion coefficients (DL) by up to 50% (Table 3) compared with the lipid diffusions in membranes without the peptide, regardless of the composition of the membrane. We note that we also computed the DL of upper L5 (Table S13†), and the results were similar to those of upper L10, suggesting that the ca. 10 lipids closest to BP100 [19] are more affected than non-neighboring lipids. No significant effects of peptide flip on DL were found in our simulations. Similarly, to previous computed lipid/membrane properties, this decrease was most significant in upper-L10 than in other lipid groups, highlighting BP100 capacity of disturbing nearby lipids in the membrane. Although lateral diffusion of lipids in membranes is a more complex process than a simple 2D motion of lipids,[77] our results showed that BP100 reduces local lipid mobility. This reduction is compatible with recent reports on both in vivo and in vitro effect of a CHAMP on membrane fluidity.[76]

Lastly, the effects of the peptide binding in the membrane hydration were investigated. Previous experimental[78-83] and theoretical[84-89] studies have reported the hydration of pure lipid bilayers and water penetration in the hydrophobic core. In our simulations, we have shown that a single unit of BP100 does not affect the average water density in membranes, even with peptide flip (Fig. 7 and S14†). Due to the larger volume of lipids in the membrane hydrophobic core compared to the small amount of water molecules that enter into this region in the nanosecond scale, the water number density profile analysis is not suit for accurate water detection and its dynamics (Fig. 7A).

However, the computed probability of finding a single water molecule in the membrane hydrophobic core was higher in those where full peptide flip took place (Table 4, S19 and S20†). Such deep water penetration, although thermodynamically unfavorable in bilayers, has been reported experimentally.[78-81] Recently, it was shown that water concentration in the hydrophobic core of POPC bilayers should be less than one water molecule per lipid,[82] and water permeability was observed in the hydrophobic core of oleate vesicles.[83] The process of water permeability through a membrane is thought as a diffusive process with a fairly constant free energy barrier, which depends on the breaking of hydrogen bonding between water molecules that have already penetrated the membrane and on the energy gain from random collisions.[84] As our results show, peptide binding followed by peptide flip can lead to an increase in the probability of finding a single water molecule in the hydrophobic core of the membrane due to the increasing water residence half-life time in the peptide-bound leaflet. Additionally, the increase of 35% in the transmembrane water transport from lower to upper R2 region compared with the opposite direction, creates a water flux through the membrane hydrophobic core in the peptide direction. Peptide flip decreases the lipid tail order parameters (Fig. 6) which would give rise to free volume pockets large enough to accommodate water molecules, and the low viscosity in the core region would favor water higher diffusivity compared to bulk water.[84,90,91]

Our results could possibly explain the earlier stages of water permeability findings of antimicrobial peptides that favor water penetration in bacterial membranes even at very low peptide concentrations, without destroying the membrane integrity.[92-95] However, simulations with higher peptide concentrations and larger membranes should be performed to confirm the positive correlation between peptide concentration and water flow and concentration within membranes.

Conclusion

Our results demonstrated that a single BP100 peptide can largely affect anionic or zwitterionic/anionic lipid membrane local properties. Membrane thinning was observed, while membranes acquired a curved conformation, concomitant with increased lipids chain conformational freedom, reduced lipid lateral diffusion and favored the water residence half-life time in the membrane hydrophobic core, and transmembrane water transport into the peptide direction. These effects depended upon the lipids neighboring BP100: a minimum amount of anionic lipid is necessary for BP100 to disturb the membrane. Thus, our results confirmed the need for correct representation of single peptide–membrane interaction, as several membrane properties can be (locally) affected.

Experimental results (circular dichroism, NMR, dynamic light scattering, zeta potential, electrophoretic mobility and leakage of dyes with fluorescence) of BP100 action upon membranes of various anionic lipid content suggest that the clustering of negatively charged lipids was required for membrane disruption.[16] Taken together with previous reports that alpha-helical BP100 inserts its hydrophobic residues inside membrane hydrophobic region while maintaining the electrostatic interactions of charged residues with lipid headgroups,[19] the present results provide a more clear picture of the early stages of BP100-anionic membrane interaction while providing insights on BP100 mechanism of action for membrane disruption.

The electrostatic long-range peptide–lipids interaction must drive BP100 towards the membrane when adsorption occurs. The highly positively charged peptide may promote anionic lipid clustering, locally disturbing membrane while acquiring an alpha-helix conformation. Such events would trigger a cooperative process, with anionic lipid clustering increasing the probability of subsequent BP100 binding, in line with the mechanism of action of CHAMPs upon membranes proposed by Epand and Epand.[96,97]

Authors contributions

P. P., F. S. L., H. C., L. R. F. and K. C. designed the simulations, and P. P. and L. R. F. performed simulations, analyzed and wrote the main manuscript text with tables and figures. L. R. F. wrote the main analysis scripts. P. P., L. R. F., F. S. L, H. C., I. M. C. and K. C. analyzed the data and reviewed the manuscript.

Conflicts of interest

The authors declare no competing interests.