Identification of electroporation sites in the complex lipid organization of the plasma membrane

  1. Lea Rems  Is a corresponding author
  2. Xinru Tang
  3. Fangwei Zhao
  4. Sergio Pérez-Conesa
  5. Ilaria Testa
  6. Lucie Delemotte  Is a corresponding author
  1. KTH Royal Institute of Technology, Dept. Applied Physics, Science for Life Laboratory, Sweden
  2. University of Ljubljana, Faculty of Electrical Engineering, Slovenia
  3. University of Chinese Academy of Sciences, China

Abstract

The plasma membrane of a biological cell is a complex assembly of lipids and membrane proteins, which tightly regulate transmembrane transport. When a cell is exposed to strong electric field, the membrane integrity becomes transiently disrupted by formation of transmembrane pores. This phenomenon termed electroporation is already utilized in many rapidly developing applications in medicine including gene therapy, cancer treatment, and treatment of cardiac arrhythmias. However, the molecular mechanisms of electroporation are not yet sufficiently well understood; in particular, it is unclear where exactly pores form in the complex organization of the plasma membrane. In this study, we combine coarse-grained molecular dynamics simulations, machine learning methods, and Bayesian survival analysis to identify how formation of pores depends on the local lipid organization. We show that pores do not form homogeneously across the membrane, but colocalize with domains that have specific features, the most important being high density of polyunsaturated lipids. We further show that knowing the lipid organization is sufficient to reliably predict poration sites with machine learning. Additionally, by analysing poration kinetics with Bayesian survival analysis we show that poration does not depend solely on local lipid arrangement, but also on membrane mechanical properties and the polarity of the electric field. Finally, we discuss how the combination of atomistic and coarse-grained molecular dynamics simulations, machine learning methods, and Bayesian survival analysis can guide the design of future experiments and help us to develop an accurate description of plasma membrane electroporation on the whole-cell level. Achieving this will allow us to shift the optimization of electroporation applications from blind trial-and-error approaches to mechanistic-driven design.

Introduction

The plasma membrane of a cell is a complex assembly of hundreds of different types of lipids and membrane proteins, which tightly regulate transmembrane trafficking and participate in cell signalling (Krapf, 2018; van Meer et al., 2008). The molecular organization of the plasma membrane and its integrity are essential for the life of the cell. However, when the cell is exposed to external forces, the membrane integrity can become transiently disrupted by formation of transmembrane pores. Such disruption can be useful in many clinical applications, for example when nucleic acids need to be delivered across the plasma membrane into the cell interior, where they can carry out their tasks (Gary and Weiner, 2020; Glass et al., 2018). Various physical methods can induce transmembrane pores, including ultrasound, light, electric field, and mechanical deformation (stretching/squeezing) (Gurtovenko et al., 2010; Ding et al., 2017; Yang et al., 2020; Schneckenburger, 2019). In terms of clinical applications, poration by the application of electric fields or electroporation is the most widely used method. It is approved for treatment of solid tumours, and it is being tested in clinical trials for gene therapy, vaccination against cancer and infectious diseases, and for cardiac ablation (Campana et al., 2019; Geboers et al., 2020; Algazi et al., 2020; McBride et al., 2021).

One of the bottlenecks of electroporation is that its protocol, especially the parameters of the applied electric pulses (amplitude, duration, number, repetition rate), needs to be optimized for each specific application and also for each specific cell/tissue type (Cemazar et al., 2009; Rols and Teissié, 1998; Hunter et al., 2021). When using electroporation for intracellular delivery of nucleic acids, the cells need to survive the treatment to be able to express the transgene. On the contrary, when using electroporation as an ablation modality, the cells need to die. Any electroporation-based treatment thus needs to be designed to either avoid or reach the point of no return leading to cell death. At present, the understanding of this point of no return is very limited (Batista Napotnik et al., 2021). One of the main reasons for this is our insufficient understanding of the molecular mechanisms that govern the increased cell membrane permeability induced by the applied electric field (Kotnik et al., 2019). If we can identify the molecular alterations of the cell membrane, we can begin to connect them to the biological response of the cells to pulsed electric fields.

The most accepted models, that describe electroporation on the whole-cell level, consider that pores can form only in the lipid domains of the plasma membrane and that all pores exhibit a similar kinetic behaviour (Krassowska and Filev, 2007; Li and Lin, 2011; Gowrishankar et al., 2013). However, accumulating evidence from experiments and simulations on model systems speaks against these assumptions. Poration kinetics in pure lipid bilayers has been shown to depend on the type of lipids and their phase state (Perrier et al., 2017; Sengel and Wallace, 2016). Since the lipids in the plasma membrane organize in domains (Lu and Fairn, 2018; Levental et al., 2020), there must exist locations which are more and less prone to poration. Moreover, our research suggests that pores can nucleate within some membrane proteins, causing protein denaturation and lipid rearrangement (Rems et al., 2020). Such lipid/protein pores can be more stable than pure lipid pores and are more likely to explain the persistent increase in plasma membrane permeability following exposure to electric pulses. Studies have further shown that pore formation and/or expansion is affected by the actin cytoskeleton, either via actin’s influence on lipid organization or the mechanical properties of the membrane (Muralidharan et al., 2021; Perrier et al., 2019). The current challenge is to gather this ensemble of findings into a coherent and predictive mathematical model describing electroporation of the living cell’s plasma membrane. In a living cell’s plasma membrane, pores cannot form anywhere: as soon as a sufficient number of pores are formed, the transmembrane voltage drops, preventing formation of new pores (DeBruin and Krassowska, 1999; Smith et al., 2014). In other words, pores will form preferentially in specific sites with the highest poration propensity. However, it remains to be elucidated which are the properties of these sites.

The challenge of studying pores in the plasma membrane experimentally is that pores are nanometre-sized and open only transiently, whereby most of them appear to rapidly close (ns–μs range) after turning off the electric field (Melikov et al., 2001; Bennett et al., 2014; Sözer et al., 2020). Pores have been imaged in giant unilamellar vesicles (Riske and Dimova, 2005; Lira et al., 2021); however, these pores have reached sizes on the order of 1 μm, which have not been observed in cell plasma membranes, likely because the actin cytoskeleton limits pore expansion (Perrier et al., 2019). Pores have also been imaged in real time in droplet interface bilayers with TIRF (total internal reflection fluorescence) microscopy (Sengel and Wallace, 2016); however, the membranes were exposed to seconds-long electric pulses, which are much longer than pulses used in electroporation applications (ns−ms range), and which would likely not be tolerated by living cells. A few attempts have been made to visualize pores in cells using electron microscopy (Chang and Reese, 1990; Lee et al., 2012); however, the observed pores were suggested to be artefacts of sample preparation (Teissie et al., 2005). Overall, the current state of experimental methods does not appear to be at a stage where it would provide the spatiotemporal resolution required to understand the molecular mechanisms of plasma membrane electroporation in its entirety.

In this study, we thus resort to molecular modelling methods to investigate plasma membrane electroporation. In particular, we use coarse-grained molecular dynamics simulations, building on their success in studying membrane lateral organization and dynamic behaviour (Duncan et al., 2017; Marrink et al., 2019; Khalid and Rouse, 2020). By running electroporation simulations on lipid membranes mimicking the realistic composition of plasma membranes, we confirm that pores do not form homogeneously across the membrane, but colocalize with domains that have specific features, particularly high content of polyunsaturated (PU) lipids. By training machine learning algorithms, we further demonstrate that knowing the local lipid distribution is sufficient to predict with ~80–90% accuracy the locations, which are most likely to be porated. Additionally, by analysing poration kinetics with Bayesian survival analysis we show that poration does not depend solely on local lipid arrangement, but also on membrane mechanical properties and the polarity of the electric field. Finally, we discuss how atomistic and coarse-grained molecular dynamics simulations, machine learning methods, and Bayesian survival analysis combined can help us develop more accurate cell-level models, which are required to foster new and better electroporation-based applications.

Results

To study plasma membrane electroporation we have used coarse-grained membranes consisting of >60 different lipid types parametrized with the Martini force field (Marrink et al., 2007; de Jong et al., 2013). The membranes mimic the composition of either an idealized average mammalian plasma membrane (APM) or a human brain plasma membrane (BPM) and have been developed and equilibrated in earlier work (Ingólfsson et al., 2014; Ingólfsson et al., 2017). The lipid composition of both APM and BPM is asymmetric, the lipids in the outer leaflet being different from those of the inner leaflet. Both compositions contain similar lipid types but differ in their fractions (Figure 1, Figure 1—figure supplement 1, Tables 1 and 2).

Figure 1 with 1 supplement see all
The average mammalian and human brain plasma membranes.

Two equilibrated membranes taken from the study of Ingólfsson et al., 2017 were cut into four 30 nm x 30 nm large pieces each. The pie charts show the fraction of lipid subgroups in the inner and outer leaflets of the membranes. CHL, cholesterol; PC, phosphatidylcholine; PE, phosphatidylethanolamine; SM, sphingomyelin; GM, gangliosides; PS, phosphatidylserine; FS, fully saturated lipids; MU, monounsaturated lipids; PU, polyunsaturated lipids. Figure inspired by Ingólfsson et al., 2017.

Table 1
Separation of lipids into groups based on their headgroup type or tail saturation.

The lipids are grouped into subtypes, depending on their headgroup type and tail saturation. Lipids are considered to be fully saturated, if they contain no double bonds in any of the tails. Lipids are considered to be monounsaturated, if they contain exactly one double bond in one or both of the tails. Lipids are considered to be polyunsaturated, if they contain at least two double bonds in one or both of the tails*.

Group nameAbbrev.Martini lipids
PhosphatidylcholinesPCDAPC, DOPC, DPPC, OIPC, OUPC, PAPC, PEPC, PFPC, PIPC, POPC, PUPC
PhosphatidylethanolaminesPEDAPE, DOPE, DUPE, OAPE, OIPE, OUPE, PAPE, PIPE, POPE, PQPE, PUPE
SphingomyelinsSMBNSM, DBSM, DPSM, DXSM, PBSM, PGSM, PNSM, POSM, XNSM
GangliosidesGMDBG1, DPG1, DXG1, PNG1, POG1, XNG1, DBG3, DPG3, DXG3, PNG3, POG3, XNG3, DBGS, DPGS, PNGS, POGS
CeramidesCEDBCE, DPCE, DXCE, PNCE, POCE, XNCE
LysolipidsLPCAPC, IPC, OPC, PPC, UPC, IPE, PPE
DiglyceridesDAGPODG, PIDG, PADG, PUDG
PhosphatidylserinesPSDAPS, DOPS, DPPS, DUPS, OUPS, PAPS, PIPS, POPS, PQPS, PUPS
PhosphatidylinositolsPIPOPI, PIPI, PAPI, PUPI
Phosphatic acidsPAPOPA, PIPA, PAPA, PUPA
Phosphatidylinositol phosphatesPIPPAP1, PAP2, PAP3, POP1, POP2, POP3
CholesterolCHOLCHOL
Fully saturated tailsFSDPPC, DBSM, DPSM, DXSM, PBSM, DPPS, DBCE, DPCE, DXCE, PPC, PPE, DBG1, DPG1, DXG1, DBG3, DPG3, DXG3, DBGS, DPGS
Monounsaturated tailsMUDOPC, POPC, DOPE, POPE, BNSM, PGSM, PNSM, POSM, XNSM, DOPS, POPS, POPI, POP1, POP2, POP3, POPA, PODG, PNCE, POCE, XNCE, OPC, PNG1, POG1, XNG1, PNG3, POG3, XNG3, PNGS, POGS
Polyunsaturated tailsPUOIPC, OUPC, PAPC, PEPC, PFPC, PIPC, PUPC, OAPE, OIPE, OUPE, PAPE, PIPE, PQPE, PUPE, OUPS, PAPS, PIPS, PQPS, PUPS, PAPI, PIPI, PUPI, PAP1, PAP2, PAP3, PAPA, PIPA, PUPA, PADG, PIDG, PUDG, APC, IPC, UPC, IPE, DAPC, DUPE, DAPE, DAPS, DUPS
  1. *

    The way in which lipids are grouped by tail saturation is different than in Ingólfsson et al., 2017, where the grouping is based on the total number of double bonds in both lipid tails. Here, the grouping is motivated by the role of lipid oxidation in electroporation, whereby lipid tails containing two or more double bonds are considerably more prone to oxidative damage than tails containing a single double bond. This is because bis-allylic hydrogens are much more easily abstracted by free radicals compared to allylic hydrogens (Reis and Spickett, 2012). Furthermore, membranes made of polyunsaturated lipids (by our definition) were found to be considerably more prone to poration/rupture by mechanical stretching compared to membranes made of lipids containing a single bond in one or both lipid tails (Olbrich et al., 2000). Thus, we consider that a lipid is polyunsaturated only if it contains at least one polyunsaturated tail.

Table 2
The total number of lipids and the percentage of individual lipid groups in each of the eight membranes.

The numbers correspond to inner/outer leaflet.

Average plasma membraneBrain plasma membrane
mem #1mem #2mem #3mem #4mem #1mem #2mem #3mem #4
Total number of lipids in each leaflet
in / out1503/16381564/16631571/16501553/16631728/18911740/18471770/18861764/1906
Percentage of lipid groups
CHOL27.7/31.730.4/32.130.8/32.129.6/33.044.5/45.745.1/46.145.4/46.344.4/46.4
PC17.9/36.418.9/34.315.9/34.717.8/34.312.4/23.012.4/23.614.8/25.814.2/21.9
PE24.5/6.621.9/4.624.8/5.325.2/6.222.8/9.721.8/11.719.6/9.220.9/10.5
SM10.6/18.78.9/19.59.7/19.68.6/19.12.1/8.92.1/9.02.3/9.22.8/8.5
GM0.0/3.70.0/7.00.0/5.20.0/4.70.0/11.30.0/8.60.0/8.00.0/11.5
PS11.3/0.010.7/0.09.6/0.010.6/0.09.8/0.011.0/0.09.1/0.08.7/0.0
PI3.7/0.04.5/0.05.0/0.05.0/0.05.4/0.04.9/0.05.0/0.05.1/0.0
PIP2.0/0.01.7/0.01.1/0.01.5/0.01.4/0.00.9/0.01.3/0.01.7/0.0
PA1.6/0.02.1/0.02.0/0.00.9/0.00.3/0.00.6/0.00.5/0.00.2/0.0
LPC0.0/1.50.0/0.80.0/1.40.0/1.00.5/0.30.2/0.50.5/0.60.3/0.3
DAG0.5/1.00.3/0.50.4/0.80.3/1.00.4/0.20.3/0.10.8/0.31.0/0.1
CE0.3/0.40.6/1.10.6/0.90.6/0.70.2/0.70.6/0.30.7/0.50.7/0.7
FS6.6/12.24.4/14.85.3/14.74.8/14.05.1/22.45.2/19.36.1/19.85.8/20.6
MU24.6/26.725.4/29.524.4/27.125.4/26.613.2/15.512.9/15.914.1/16.615.2/17.2
PU41.1/29.439.8/23.639.5/26.240.2/26.537.2/16.336.8/18.734.4/17.234.5/15.8

The following subsections present the results and analysis as follows. First, we present the results from electroporation simulations and demonstrate that all membranes exhibit preferential poration sites. Then we determine local membrane properties and examine which of them increases/decreases the poration propensity. We further investigate the importance of local membrane properties by training machine learning algorithms and predicting the sites, which are most likely to be porated. Finally, we apply Bayesian survival analysis to investigate how membrane properties influence the poration kinetics and to develop an underlying kinetic model.

Membranes exhibit preferential poration sites

The original APM and BPM membranes from Ingolfsson et al. (Ingólfsson et al., 2014; Ingólfsson et al., 2017) were cut into four pieces each (Figure 1) to increase sampling, since the analysis was focused on the first poration event. After a short re-equilibration, each of the eight membranes was subject to electroporation simulations, 60 simulations under hyperpolarizing field and 60 simulations under depolarizing electric field. Both polarities of the electric field were considered, since during exposure of a cell to an electric field, the plasma membrane becomes hyperpolarized on the side facing the positive electrode and depolarized on the side facing the negative electrode. The electric field was set to ±127.7 mV/nm, which was high enough to induce a pore within ~15 ns. Being able to observe poration over short time scales was important to minimize lateral lipid diffusion and provide a reliable mapping of local membrane features before electric field application to the likelihood of a poration event. Note that the value of the electric field strength imposed in simulations is not directly comparable to the electric field strengths reported in experimental studies (see section Molecular dynamics simulations for further explanation). After identifying when and where the pores formed in each simulation, we observed that pore locations do not distribute homogeneously along the membrane surface but often cluster together (Figure 2). An individual membrane exhibits one or more of such clusters. The location of pores is similar albeit not completely identical under depolarizing and hyperpolarizing electric field.

Location and kinetics of pore formation in each membrane upon depolarization and hyperpolarization.

The pore locations are expressed relative to the dimensions of the simulation box at the poration time, to correct for the natural fluctuations in the membrane area during the simulation. All data points correspond to the first poration event. The colour of the circle codes for poration time according to the colour bar. (A) Average plasma membranes. (B) Brain plasma membranes.

Poration sites colocalize with domains with specific features, in particular with a high density of PU lipids

We hypothesized that pores preferably form in nanodomains with specific features. We used the recently developed tool MemSurfer (Bhatia et al., 2019) to extract from each membrane the local area per lipid (APL), membrane thickness, mean curvature, cosine of the dipole angle (cos θdip), charge, and lipid tail order parameter (Figure 3—figure supplement 1). We also determined the local density of individual groups of lipids, grouping the lipids either according to their head architecture or their tail saturation. These membrane features were extracted from 101 frames of a 10-ns long trajectory before electric field application, from each leaflet separately and in points corresponding to locations where a pore has or has not formed in any of the electroporation simulations. The first and second group of points are labelled as ‘porated’ and ‘non-porated’ locations, respectively (Figure 3A, left). We estimated the distribution of the values of the various features in non-porated and porated locations by constructing histograms (Supplementary file 1) and quantified the difference between probability density estimates by means of the Kullback–Leibler divergence (Figure 3B). By far the most significant feature to distinguish the locations where pore formed from locations where they did not was the presence of PU lipids, in both the APM and BPM membranes. This finding is corroborated by visualizing the colocalization of pore clusters with nanodomains enriched with PU lipids (Figure 3A). However, the analysis further showed that both the lipid head and tail architecture influence poration, with pores being favoured in regions with higher density of phosphatidylcholine (PC) lipids and lower densities of cholesterol (CHOL), gangliosides (GM), and fully saturated (FS) lipids. In BPM, the pores are also favoured in regions with higher content of phosphatidylethanolamine (PE) lipids, where most PE lipids in BPM are polyunsaturated (Figure 1—figure supplement 1). In addition, the analysis showed that pores in all membranes are favoured in regions with greater area per lipid and lower lipid order. None of the analysed features appeared markedly dependent on the electric field polarity, even when we contrasted them with 10-ns-long trajectories obtained under non-porating electric field (Figure 3—figure supplement 2).

Figure 3 with 2 supplements see all
Membrane features that favour and disfavour poration.

(A) (Left) Definition of porated and non-porated locations in one of the membranes; (middle) most pore locations colocalize with increased polyunsaturated lipids density; (right) the corresponding histogram of polyunsaturated lipids density values in porated and non-porated locations. (B) Distances (KL divergence) between probability density estimates of individual features in porated and non-porated locations. The higher the bar, the more the feature influences poration. Positive and negative bars show whether a feature favours or disfavours poration, respectively. APL, area per lipid; PC, phosphatidylcholine; PE, phosphatidylethanolamine; SM, sphingomyelin; GM, gangliosides; PS, phosphatidylserine; PI, phosphatidylinositol; PA, phosphatic acid; PIP, phosphatidylinositol phosphate; LPC, lysolipids; CE, ceramides; DAG, diglycerol; CHL, cholesterol; FS, fully saturated lipids; MU, monounsaturated lipids; PU, polyunsaturated lipids.

Knowing the lipid distribution is sufficient for machine learning models to reliably predict poration sites

The local poration propensity in both APM and BPM and under both electric field polarities is governed by similar features, for example, high density of PU lipids, as suggested by Kullback–Leibler divergence (Figure 3B). To corroborate this finding further, we trained three machine learning models, namely random forest, support vector machine (SVM), and multilayer perceptron neural network (Fleetwood et al., 2020), on selected subsets of data, that is, using features from APM or BPM and/or features obtained under hyperpolarization or depolarization. The accuracy of predicting poration in another/different subset of data typically surpassed 80% in all three models tested, with random forest exhibiting slightly superior performance (Table 3). Similar accuracy was obtained regardless of whether we trained and tested the models on (1) datasets from different membrane pairs of the same composition and the same electric field polarity (rows 1–6 in Table 3); (2) membranes of the same composition but different electric field polarity (rows 7–8 in Table 3); or (3) membranes of different composition (rows 9–11 in Table 3). Visual inspection of the prediction showed that locations corresponding to pore clusters are reliably predicted, whereas those that are scattered away from the clusters tend not to be predicted as accurately (Figure 4A).

Table 3
Prediction accuracy by machine learning models, reported for the training and test datasets; dep, depolarization; hyp, hyperpolarization.
#Training datasetTest datasetRandom forest all featuresSVM all featuresNeural network all featuresRandom forestfour features
TrainTestTrainTestTrainTestTrainTest
Same membrane composition, same E-field polarity; train on dataset from two out of four membranes
1APM-dep, mem 1 and 2APM-dep, mem 3 and 4100%83%97%83%97%81%100%82%
2APM-dep, mem 1 and 3APM-dep, mem 2 and 4100%80%96%80%94%77%100%76%
3APM-dep, mem 1 and 4APM-dep, mem 2 and 3100%85%96%85%91%81%100%82%
4BPM-dep, mem 1 and 2BPM-dep, mem 3 and 4100%87%97%86%94%86%100%84%
5BPM-dep, mem 1 and 3BPM-dep, mem 2 and 4100%82%97%84%97%79%100%85%
6BPM-dep, mem 1 and 4BPM-dep, mem 2 and 3100%75%97%83%96%77%100%82%
Different membrane composition, same E-field polarity
7APM-depBPM-dep100%83%94%82%93%83%100%83%
8APM-hypBPM-hyp100%83%94%82%95%83%100%83%
Same membrane composition, different E-field polarity
9APM-depAPM-hyp100%86%95%83%94%81%100%84%
10BPM-depBPM-hyp100%90%95%86%95%84%100%88%
11APM and BPM-depAPM and BPM-hyp100%88%94%84%93%82%100%86%
Same membrane composition, same E-field polarity; train on 60% dataset from all four membranes
12APM-dep (60%)APM-dep (40%)100%99%94%91%91%91%100%92%
13APM-hyp (60%)APM-hyp (40%)100%99%94%90%92%92%100%91%
14BPM-dep (60%)BPM-dep (40%)100%99%95%93%94%94%100%93%
15BPM-hyp (60%)BPM-hyp (40%)100%99%94%91%95%94%100%93%
16APM and BPM-dep (60%)APM and BPM-dep (40%)100%99%94%92%92%92%100%92%
17APM and BPM-hyp (60%)APM and BPM-hyp (40%)100%99%92%92%91%91%100%92%
Random forest output.

(A) Comparison between the predicted poration sites (red dots) and real poration sites (yellow dots). Only locations that were classified as porated in >90% of the trajectory frames are shown. (B) Feature importance score for average plasma membrane (APM) and brain plasma membrane (BPM) obtained under hyperpolarization and depolarization. The feature importance score quantifies how much each feature is used in each tree of the random forest.

The prediction accuracy became considerably improved when a single dataset was randomly split into 60% training/40% test dataset (rows 12–17 in Table 3). In such case, random forest reached the highest accuracy of 99%. The improvement in accuracy suggests that there is something unique about poration of each membrane, even when membranes have practically identical overall composition, which the machine learning models can capture. However, this high accuracy could be partially biased by oversampling, which was performed to balance the starting number of porated locations (~60) and the number of non-porated locations (300), see Machine learning methods. When randomly choosing 60% of a given dataset for training, the values of features from all (not just 60%) actual porated locations are effectively taken into account. Therefore, we made additional tests, where we trained random forest on the first 32 porated and 150 non-porated locations for each of the APM membrane under depolarization. The accuracy of the prediction for the rest of the porated and non-porated locations was 92%, 97%, 94%, 95% for membranes 1–4, respectively. The accuracy was lower than 99%, but still considerably higher than when training/testing on pairs of different APM membranes (rows 1–3 in Table 3). Consequently, this exercise suggests that prediction accuracy might be improved by finetuning machine learning models on data from more different membranes or by adding additional features to the analysis in future studies.

From the datasets for both APM and BPM, we determined the feature importance score, which confirmed that the density of PU lipids is the most important for poration, followed by density of GM, CHOL, and FS lipids (Figure 4B). Furthermore, the feature importance score suggested that knowing the distribution of these lipid groups is sufficient to predict poration sites. The densities of other lipids and other membrane properties including the area per lipid, thickness, …, and lipid tail order can be practically neglected. We confirmed this by training random forest using four most important features, that is the density of PU, GM, CHOL, and FS lipids. The accuracy of the prediction was indeed very similar (≥76%) as when using all features (Table 3, last column).

Electroporation kinetics depends on membrane composition and electric field polarity but can be described with a universal model

The analysis presented so far demonstrates that the pore location mainly depends on the local lipid arrangement, particularly the density of PU lipids, regardless of the type of membrane and electric field polarity. However, this does not necessarily mean that poration is equally likely or equally fast in all membranes. To investigate whether the membranes differ in poration kinetics we compared the distributions of first poration times for each membrane. These distributions were indeed found to depend both on membrane composition and electric field polarity (Figure 5A). APM membranes tend to porate faster than BPM, and hyperpolarizing fields tend to porate faster than depolarizing fields.

Figure 5 with 6 supplements see all
Kinetics of membrane poration.

(A) Probability distribution of the first poration times, sampled from 60 simulations, and approximated with kernel density estimation. (B) The time course of λ0(t). (C) Credible intervals of exp(βi) for average plasma membrane (APM) and brain plasma membrane (BPM) under depolarization and hyperpolarization. (D, E) Correlation between βi and the steady-state relative change in membrane thickness Δd/d0 and projected area ΔA/A0 normalized by E2, extracted from least-square fits presented in Figure 5—figure supplements 5 and 6. Square symbols show the average over all four membranes, whereas the error bars show the range of values in four membranes. The ratio ΔA/A0 / E2 is inversely proportional to the area compressibility (or stretching) modulus (Needham and Hochmuth, 1989).

A Bayesian survival model was then trained (1) to quantify and study the statistical significance of the differences in poration rates, (2) to investigate if an underlying universal model can be used to describe the poration kinetics, and (3) to parametrize the kinetic parameters. Cox’s proportional hazards survival model (Cox, 1972) is a suitable model for this purpose, since it allows for an arbitrary time dependence of the poration rate, and it allows to then compare different systems or situations. This model assumes that for a given system, i, the instantaneous probability of an event happening, λi(t) (the rate pore formation), is a function of time:

(1) λi(t)=λ0(t)exp (βi)

where λ0(t) is a time-dependent and system-independent baseline rate of pore formation and βi are regression coefficients accounting for the difference in poration kinetics between systems (APM-depolarization, APM-hyperpolarization, BPM-depolarization, BPM-hyperpolarization). The baseline rate captures the common time dependency whereas the regression coefficient captures the system specificity. The inferred model was able to successfully reproduce the measured data (Figure 5—figure supplement 1) highlighting its robustness and suggesting the existence of a common kinetic model for electroporation. To further test the model validity, we trained a more complex model which allowed the βi coefficients to vary with time, allowing for a system-specific time variability of λi(t). The βi(t) obtained were found to slightly drop in the interval 0 < t < 3 ns until becoming practically time independent (Figure 5—figure supplement 2). Given that the βi(t) for all systems follow the same trend and that the time-dependent model only adds information in the initial stage of electroporation and has a high uncertainty, the time-independent model was selected for its easier interpretability and lower complexity.

The time-independent Bayesian survival model shows that the baseline poration rate λ0(t) is initially zero and then increases, reaching a constant steady-state value after approximately 3 ns (Figure 5B). To determine the mechanistic basis for the initial transient kinetic regime, we performed additional simulations at non-porating electric fields. These simulations revealed that the electric field E causes the membrane to thin and expand its area, whereby the relative change in membrane thickness and area are proportional to E2 (Figure 5—figure supplement 5 and Figure 5—figure supplement 6), as expected for Maxwell stress (Lewis, 2003). The membrane thickness and area reach a steady state ~3 ns after electric field onset (Figure 5—figure supplement 4), suggesting that the time scale of the initial transient kinetic regime is related to the progressive thinning and area expansion caused by Maxwell stress.

The calculated βi credible intervals (Figure 5C, Figure 5—figure supplement 3) confirm that the differences in poration kinetics between different systems are statistically significant. APM membranes porate ~5 and ~2 times faster than BPM for hyperpolarizing and depolarizing fields, respectively. Hyperpolarizing fields porate membranes ~3 and ~2 times faster than depolarizing fields for APM and BPM membranes, respectively. We found that the values of βi correlate with the ability of a membrane to thin and expand its area under electric field, that is, they corelate inversely with the membrane area compressibility (or stretching) modulus (Figure 5D, E).

Discussion

In this study, we focused on lipid pores formed by electric field (electroporation) and we asked three main questions: Where do pores form in membranes with realistic plasma membrane lipid composition? Which membrane features/properties govern the most likely poration sites? Which membrane features/properties govern the poration kinetics?

Features important for membrane poration

Research conducted so far has mainly focused on studying formation of pores in model membrane systems with up to three different lipid types. In membranes containing a single type of fluid-phase lipids, poration was found to depend on both the head and tail architecture. Longer tail length, larger headgroups (e.g. heads containing sugar moieties), and stronger intramolecular interactions between both lipid tails and heads generally reduce the propensity for poration (Ziegler and Vernier, 2008; Piggot et al., 2011; Polak et al., 2013; Polak et al., 2014; Gurtovenko and Lyulina, 2014). Poration was further found to be more difficult in membranes with gel-phase lipids compared fluid-phase lipids (Knorr et al., 2010; Majhi et al., 2015; Liu et al., 2014). In binary mixtures of fluid- and gel-phase lipids, experiments suggested that pores form in the fluid domains (Perrier et al., 2018). Mixing a lipid with cholesterol was found to either decrease or increase the poration propensity, whereby this appears to depend on the cholesterol concentration and the architecture of the lipid (Portet and Dimova, 2010; Fernández et al., 2010; Mauroy et al., 2015; Casciola et al., 2014; Perrier et al., 2017). Studies were also done on ternary mixtures containing cholesterol, where the lipids organize in the liquid-ordered and -disordered domains. Both molecular dynamics simulations and experiments showed that pores preferentially form in liquid-disordered domains (Reigada, 2014; Sengel and Wallace, 2016). Overall, studies on simple model systems showed that there are many different factors influencing poration and it is quite impossible to predict which ones will be most important in a complex mixture such as the plasma membrane.

Our findings from membranes composed of >60 different lipid types mimicking the realistic composition of plasma membranes are fully in agreement with findings from simpler model systems. We observed that pores form preferentially in locations with lower lipid tail order, and that pores avoid regions enriched with gangliosides containing large sugar headgroups, regions enriched with cholesterol and fully saturated lipids. A new finding from our study is that, in plasma membranes, the most important factor governing poration propensity is the local density of polyunsaturated lipids, regardless of the membrane type (average or brain plasma membrane) and polarity of the electric field. Note that the majority of polyunsaturated lipids in both average and brain plasma membranes have either PC or PE headgroup (>75% of all polyunsaturated lipids; Figure 1—figure supplement 1), whereby PC and PE headgroups also contribute favourably to poration (Figure 3B).

Polyunsaturated lipids and lipid oxidation

The finding that poration is strongly facilitated in the presence of polyunsaturated lipids is very important from the standpoint of lipid oxidation. Polyunsaturated lipids are highly prone for oxidation by free radicals, whereby upon oxidation a part of the lipid tail becomes hydrophilic making the membrane considerably more permeable to ions and hydrophilic molecules (Yin et al., 2011; Boonnoy et al., 2015; Rems et al., 2019). Several experimental studies have shown that electroporation is associated with oxidative lipid damage, and it is believed that lipid oxidation plays an important role in increased membrane permeability after electric field exposure (Maccarrone et al., 1995; Breton and Mir, 2018). Our study suggests that most pores form in domains enriched with polyunsaturated lipids. If pores somehow act as precursors for lipid oxidation, such as by improving the access for free radicals to lipid tails, this could provide a clue on how electroporation is related to lipid oxidation.

Electroporation kinetics

Bayesian survival analysis demonstrated that all membranes are characterized by two distinct kinetic regimes. Initially, no poration occurs, but progressively as the membrane area expands due to Maxwell stress, poration becomes easier and the poration rate increases, plateauing at a maximum rate when the membrane has reached its steady-state area. The steady-state poration rate was found to be inversely correlated with membrane area compressibility modulus – a measure of how resistant the membrane is to compression or expansion (Figure 5E). Inverse correlation between electroporation propensity and membrane area compressibility modulus has also been found experimentally (Needham and Hochmuth, 1989). The increase in membrane area is directly related to a decrease in membrane thickness (Figure 5—figure supplement 4). Membrane thinning is expected to facilitate pore formation, as water molecules need to travel shorter distance when bridging the membrane. Note that when comparing properties in porated and non-porated locations (Figure 3B) the variation in local membrane thickness was not large enough to discriminate porated locations. However, pore formation was indeed found to be favoured in regions with larger area per lipid.

In the steady-state kinetic regime, the poration rate becomes constant and the time distribution of the poration events becomes exponential. A constant poration rate and an exponential distribution of first poration times are characteristic of a Poisson behaviour of the number of pores formed per unit time. Indeed, the most accepted kinetic models of electroporation assume a constant poration rate at a fixed value of the transmembrane voltage (or the electric field within the membrane) (Neu and Krassowska, 1999; Vasilkoski et al., 2006), which is consistent with our model in the steady-state regime. Observing the non-steady-state transient kinetic regimes is beyond the time resolution of most experiments. As such, our model has the advantage of characterizing the transient initial kinetic regime, which should be important when exposing cells to increasingly used pulses with duration in the (sub)nanosecond range (Pakhomov and Pakhomova, 2020; Neuber et al., 2019).

Electroporation rate on the whole-cell level is typically modelled as (Neu and Krassowska, 1999; Vasilkoski et al., 2006; Glaser et al., 1988):

(2) λ=A exp (δBΔΨ2kT)

where ΔΨ is the transmembrane voltage, δ is the energy barrier for pore formation at ΔΨ = 0 V, B is a proportionality constant, k is the Boltzmann constant, and T is the temperature. A is a prefactor that is proportional to the number of possible pore nucleation sites and the frequency of lateral lipid fluctuations. Comparing the ratios of the pore formation rates of two systems i and j in Equations (1-2), assuming AiAj and a steady-state kinetic regime we obtain:

(3) βi-βj-δi-BiΨ2kT+δj-BjΨ2kT

Therefore, the difference between βi of two systems is approximately equal to the negative difference of their steady-state poration-free energy barrier in kT units. This relates the βi regression coefficients to their physical interpretation. Assuming approximately equal prefactors AiAj is reasonable for our membranes because all membranes have the same total area. In addition, the sites with the highest poration propensity are comprised of similar lipids, so we can safely assume that those regions have similar lipid fluctuation frequency. Note that the pore formation barrier δ can be independently determined by free energy methods (Hu et al., 2015), whereas the parameter B can be inferred from simulating electroporation at different values of the transmembrane voltage; however, both these approaches are computationally more demanding. Bayesian survival analysis thus offers an alternative and simpler way to obtain parameters of Equation 2 for different systems.

Kinetic differences between the average and brain plasma membranes

The average plasma membranes (APMs) exhibit shorter poration times compared to the brain plasma membranes (BPMs). Bayesian survival analysis confirmed that the poration kinetics in APMs and BPMs is statistically different. Compared with BPMs, APMs contain a smaller fraction of cholesterol and fully saturated lipids and a greater fraction of polyunsaturated lipids, all favouring poration and likely increasing the poration rate. Note that the polyunsaturated lipids in BPMs have on average greater number of double bonds than in APMs (Ingólfsson et al., 2017) however, this is not enough to balance the poration rate with APMs. In experimental studies, cells of different types, or even cells in different phases of the cell cycle, have been found to exhibit different electroporation thresholds (i.e., different electroporation propensities) (Cemazar et al., 2009; Towhidi et al., 2008; Golzio et al., 2002). According to our results, this difference in thresholds is, at least in part, related to the fact that cells of different types can have considerable differences in their lipid composition, which further changes along the cell cycle (Zhang et al., 2017; Atilla-Gokcumen et al., 2014).

Kinetic differences between hyperpolarizing and depolarizing electric field

When a cell is exposed to an electric field, its membrane becomes hyperpolarized on the side facing the positive electrode (anode) and depolarized on the side facing the negative electrode (cathode). Experiments have shown that plasma membranes of different types of cells can become more permeabilized either on the anodic or cathodic side, whereby this asymmetry in permeabilization is still not completely understood (Mehrle et al., 1985; Kinosita et al., 1992; Tekle et al., 1994; Sözer et al., 2017). Our results suggest that the asymmetric lipid composition present in all mammalian plasma membranes favours pore formation on the anodic side, which is hyperpolarized. We found that hyperpolarization induces more profound membrane thinning compared to depolarization, consequently increasing the poration rate. Greater membrane thinning could be associated with electrophoretic drag of negatively charged lipids, which are mainly present in the inner leaflet and are pulled towards the membrane interior by hyperpolarizing electric field. Moreover, in our previous study, in which we characterized pores that formed in voltage sensors of sodium voltage gated channels (Rems et al., 2020), we also observed that such pores are more easily formed under hyperpolarization, albeit this was associated with asymmetric distribution of charged protein residues. Nevertheless, there are numerous types of membrane proteins in the plasma membrane, some of which might become denaturated more easily under depolarization. Whether the cell is permeabilized more on the anodic or cathodic side might therefore depend on the plasma membrane’s lipid–protein content and the preferential sites of pore formation. This exemplifies the need to identify preferential poration sites in membranes with complex organization, which we discuss further in Towards building accurate cell-level models of electroporation.

Towards building accurate cell-level models of electroporation

The cell membrane is a complex organization of lipids and proteins, whereby insights from atomistic molecular dynamics simulations suggest that the electric field can form pores both in the lipid domains and within some membrane proteins. The limitation of molecular dynamics simulations is that they are only able to model a small part of the membrane, and that they are not able to take into account the dynamic changes in the transmembrane voltage, which are present when electroporating whole cells (Hibino et al., 1993; Frey et al., 2006). Namely, on the whole-cell level, the induced transmembrane voltage varies with position on the membrane and depends on the spatiotemporal profile of the membrane conductivity. As sufficient number of pores form in the membrane, they increase the membrane conductivity to the extent that starts decreasing the transmembrane voltage and prevents formation of additional pores. In other words, in molecular dynamics simulations it is possible to observe poration of practically any membrane; however, in a real cell membrane only the sites with the highest poration propensity can be porated.

To understand electroporation of living cells, we need to be able to develop equations describing the poration rate of different types of pores (pores in the lipid domain, different membrane proteins) and embed them into a system of ordinary and partial differential equations that describe electroporation on the whole-cell level (DeBruin and Krassowska, 1999; Smith et al., 2014). Such models can then be used to study the increase in membrane conductivity, transmembrane molecular transport of different types of molecules, and changes in the membrane resting potential and/or action potentials induced by different parameters of electric pulses. Electroporation can be caried out with different pulse waveforms, where the duration of individual pulses can range from a few 100 picoseconds to tens of milliseconds. Exploring the pulse-parameter space in silico instead of through trial-and-error experimental approaches will facilitate optimization of electroporation-based applications in vitro and in vivo (Rems and Miklavčič, 2016).

We envision that by combining atomistic and coarse-grained molecular dynamics simulations, machine learning methods and Bayesian survival analysis we will be able to improve existing kinetic models describing electroporation on the whole-cell level. Coarse-grained simulations, together with more mesoscopic models and experiments, are in the future anticipated to enable modelling of whole plasma membranes providing their detailed molecular organization (Pezeshkian and Marrink, 2021). An exciting finding from our study is that knowing the lipid distribution is sufficient for identifying the most likely poration sites in lipid domains by using machine learning methods. As such, we anticipate that we will be able to use machine learning to estimate the membrane area, which is most amendable to poration, and hence estimate the prefactors in Equation 2. By performing electroporation simulations on selected membrane regions and applying Bayesian survival analysis, we can characterize the poration kinetics and simplify the determination of the corresponding kinetic parameters. For example, by following these approaches using coarse-grained membranes associated with actin filaments (Schroer et al., 2020) we can investigate and quantify how the presence of actin cytoskeleton influences poration kinetics, either by affecting the local distribution of lipids or by influencing the mechanical properties of the membrane or both. To enable finer decomposition of the relative importance of local vs. global membrane properties, future studies could combine local lipid neighborhood analysis with pore initiation rate analysis.

Coarse-grained simulations have their disadvantages. At present, coarse-grained simulations cannot be used to study poration of membrane proteins, as the protein secondary structure typically needs to be constrained for the protein to remain stable (Periole et al., 2009). Nevertheless, poration of membrane proteins and its corresponding kinetics can be inferred from atomistic molecular dynamics simulations (Rems et al., 2020). Furthermore, coarse-grained lipid bilayers are known to be more difficult to porate than corresponding atomistic bilayers, likely due to coarse-graining of multiple water molecules into a single particle (Hu et al., 2015). Despite being more difficult to porate, coarse-grained systems are able to represent the differences in the energy barriers for pore formation in bilayers composed of different lipid types as well as the influence of the membrane mechanical properties on poration (Hu et al., 2015). This confirms that we can use coarse-grained simulations to identify the lipid domains, which are the most likely to be porated. By backmapping (Wassenaar et al., 2014), these regions to atomistic representation, we should be able to obtain a more accurate estimation of the kinetic parameters of Equation 2. Nevertheless, we hope our study will motivate further exploration on the validity of coarse-grained membrane models for studying membrane electroporation, both in comparison with corresponding atomistic computational membrane models and experimental model membrane systems.

Conclusions

Pores in the plasma membrane can be formed under diverse physicochemical conditions. They can be formed in various physiological processes by pore-forming proteins, and when the membrane is subject to external mechanical or electromagnetic forces (Gilbert et al., 2014; Sun et al., 2020). In this study, we investigated pores formed by electric field. Electroporation simulations of coarse-grained membranes mimicking realistic lipid composition of plasma membranes showed that pores preferentially form in domains enriched with polyunsaturated lipids and that pores avoid domains enriched with gangliosides, cholesterol, and fully saturated lipids. The density of polyunsaturated lipids is the most important feature governing the preferential pore location, regardless of the overall membrane composition and electric field polarity, as corroborated by machine learning methods. However, the poration kinetics does depend significantly on membrane composition and electric field polarity, as demonstrated by Bayesian survival analysis. The poration rate is higher under hyperpolarizing compared to depolarizing electric field and correlates with the ability of a membrane to expand its area under electric field. We envision that by combining atomistic and coarse-grained molecular dynamics simulations, Bayesian survival analysis and machine learning models, we will be able to improve existing kinetic models describing electroporation on the whole-cell level. Although we have focused on pores induced by electric fields, the findings are likely to be applicable to other ways of membrane poration, as we found that the lipid organization is much more important for poration than the electrical features of the membrane (charge, dipole angle).

Materials and methods

Molecular dynamics simulations

System preparation

Request a detailed protocol

The starting point for our systems was the topology files for the coarse-grained membranes used in the study of Ingólfsson et al., 2017, available at https://bbs.llnl.gov/neuronal-membrane-data.html. The membranes are parametrized with the Martini 2.2 force field (Marrink et al., 2007; de Jong et al., 2013). We took the frames extracted after 80 μs of the simulation (the files confout-80us.gro). The original membranes were about 70 nm × 70 nm large. We cut four 30 nm × 30 nm pieces from the original membranes. After cutting, we removed an appropriate number of Na or Cl ions, such that the final system had zero net charge. The NaCl concentration was ~150 mM. We replaced the non-polarizable water model with its polarizable version to have a more accurate system representation for electroporation studies (Yesylevskyy et al., 2010). We also added more water to each system, such that the simulation box in z direction was ~19 nm after equilibration. This procedure was done using functions from VMD (Visual Molecular Dynamics) (Humphrey et al., 1996), Gromacs (Abraham et al., 2015), and custom scripts. In the end we had eight membrane systems, four with lipid composition corresponding to the APM and four corresponding to the BPM composition.

Equilibration and production runs

Request a detailed protocol

After preparing the new systems and topology files, we made two steps of system minimization and four steps of system equilibration following the equilibration protocol for Martini membranes from charmm-gui.org (Qi et al., 2015). We then ran a short 500 ns simulation using the reaction-field electrostatics followed by 50 ns equilibration using Particle Mesh Ewald (PME) electrostatics (Darden et al., 1993). Some of the PC lipid’s heads were constrained in their z position to reduce bilayer undulations, as done in the original publication (Ingólfsson et al., 2017). Other MD parameters were equal to the default parameters for simulations with Martini membranes from charmm-gui.org in April 2020: leap-frog integration of the equations of motion using 20 fs time step; plain cutoff of van der Waals and Coulomb interactions at distance 1.1 nm; relative dielectric constant ɛr = 2.5; temperature coupling using velocity rescaling with a stochastic term (Bussi et al., 2007) with time constant 1 ps and temperature of 310 K; semi-isotropic pressure coupling using Parrinello–Rahman barostat (Parrinello and Rahman, 1981) with time constant 12 ps and reference pressure of 1 bar, with the compressibility of the system set to 3e−4 bar–1. The pressure coupling allows the size of the simulation box to adjusts to the changes in the membrane area. The trajectory was saved every 0.1 ns. All simulations were carried out with Gromacs 2019.4 (Abraham et al., 2015).

Electroporation simulations

Request a detailed protocol

After equilibration, and for each of the eight membranes, we ran multiple replicas of electroporation simulations where the membranes were exposed to an electric field of +127.7 mV/nm (60 replicas) or −127.7 mV/nm (60 replicas). No positional restraints were imposed in these simulations. The electric field was chosen such that we observed the formation of at least one pore within ~15 ns. Having a short poration time was important, because we aimed to map the local membrane features before poration to the likelihood of a poration event, which means that we needed to avoid considerable lateral lipid diffusion. The value of the electric field E imposed in simulations, and reported above, cannot be directly compared to that reported in experiments. The reason for this is the manner in which E is implemented in molecular dynamics. Namely, each particle carrying a charge qi is assigned an additional force F = qiE. As such, E corresponds to the electric field strength that would exist in vacuum. However, in our molecular dynamics systems there are many electric dipoles, including water molecules and headgroups of zwitterionic lipids, which respond to E by changing their average orientation, that is, they polarize. The ability of a material to polarize under electric field is characterized by relative dielectric permittivity, which corresponds to the factor by which the electric dipoles within the material reduce the electric field that would exist in vacuum. Since the relative dielectric permittivity of water is around 80 (both experimentally and for the polarizable MARTINI water model Yesylevskyy et al., 2010), the macroscopic electric field that establishes in the aqueous compartment of our molecular dynamics systems is about two orders of magnitude lower compared to the imposed electric field E (Vernier et al., 2013), that is, on the order of 1 mV/nm = 10 kV/cm. This macroscopic electric field strength is typically reported in experiments after being determined, for example, as the voltage-to-distance ratio, if placing the sample between a pair of parallel-plate electrodes. The macroscopic electric field strengths in our simulations are well within the range of those used in experiments, when exposing cells to submicrosecond pulses (order of 1–100 kV/cm). Furthermore, the probability of pore formation depends exponentially on the magnitude of the applied electric field (Equation 2; Böckmann et al., 2008). The electric field strength required for detectable electroporation thus reduces with increasing the duration of the applied electric pulse, and consequently macroscopic electric field strengths of the order of 0.1–1 kV/cm are typically sufficient for electroporation when exposing cells to conventional microsecond and millisecond electric pulses. While the conditions in our study strictly correspond to nanoseconds-long exposures to electric field, we expect similar behaviour also for longer electric pulses. Nevertheless, on a longer time scale, electrodeformation of the cell membrane and/or increases in local membrane curvature caused by electric field might play additional roles in the pore formation process (Perrier et al., 2017; Riske and Dimova, 2005).

Additional simulations under non-porating electric field

Request a detailed protocol

For each membrane, we also ran five 10-ns-long simulations under lower electric field strengths of 0, ±42.5, ±63.8, ±85.1, and ±106.3 mV/nm. No positional restraints were imposed in these simulations. The initial coordinates (initial gro file) were the same as for electroporation simulations. With the exception of a few simulations at ±106.3 mV/nm, these electric field strengths were too low to induce poration within 10 ns.

Pore localization

To determine the poration time and pore location in the electroporation trajectories we used a custom semi-automatic procedure with Python and MDAnalysis (https://www.mdanalysis.org/; Michaud-Agrawal et al., 2011). The procedure consisted of three main steps, as described below. For analysis, we generally considered only the first poration event. Additional pores could form after the first one and all pores eventually expanded until destroying the membrane, as is usual for simulations under constant electric field (Fernández et al., 2012; Delemotte and Tarek, 2012). We focused on the first poration event, as the local electric field changes after poration, changing the transmembrane voltage. In up to ~20% simulations, two or more pores formed practically simultaneously. In such a case, we considered all of these pores for analysis.

Step 1 – coarse search

Request a detailed protocol

Start from frame 12 (since no pore formed before time 1.2 ns), search every six frames (i.e. 12, 18, 24, …). In each frame, divide the membrane into 100 small pieces. In each small piece, look for water molecules around the centre of the lipid bilayer. If there are more than 12 water molecules, stop and write down the frame (t) and the region (i,j), otherwise, continue the loop. This step is used to speed up and improve the accuracy of the algorithm. The condition of 12 water molecules ensured that only fully formed pores were considered for the second step.

Step 2 – precise Search

Request a detailed protocol

Start from frame t − 5, search every frame until t (i.e. t − 5, t − 4, t − 3, …). In each frame look for water molecules within the centre of the lipid bilayer in the regions (i − 2, j − 2), (i − 2, j − 1), …, (i + 2, j + 1), (i + 2, j + 2). If there are more than four (APM) or six (BPM) water molecules, check the following two frames. If there are water molecules also in the following two frames, return the poration time and the centre of mass of these water molecules as the pore location. Otherwise, continue the loop.

Step 3 – manual check and correction

Request a detailed protocol

The procedure failed for BPM, when the pore formed later than ~7 ns after the onset of the electric field, mainly because of the increase in BPM curvature with time. In this case, the poration time and location were determined by visualization of the trajectory in VMD.

All poration times and pore locations were manually verified by extracting the porated frames from all trajectories and translating the systems in (x,y) such that the pore location moved to a predefined position. These extracted frames were then visualized in VMD.

Analysis of membrane properties

Request a detailed protocol

To determine local membrane properties we used the MemSurfer tool (Bhatia et al., 2019), available at https://github.com/LLNL/MemSurfer, (Lawrence Livermore National Laboratory, 2022; copy archived at swh:1:rev:51c7c0534f7b73e74e7233e44c08a7fb269a0188) MemSurfer is a 3D membrane surface fitting tool, which fits a surface to the inner and outer membrane leaflet and uses Delaunay triangulations and surface parameterizations to compute membrane properties of interest (Figure 3—figure supplement 1). Each vertex is mapped to the position of a lipid headgroup. The triangulated surface is then further smoothed (Figure 3—figure supplement 1D).

MemSurfer can be used to measure local membrane thickness, area per lipid, and mean curvature, among others. We added functions that analyse the dipole angle of zwitterionic lipids, the lipid tail order, and the headgroup charge. The dipole angle for each zwitterionic lipid was computed as the angle between the local membrane normal (determined by MemSurfer) and the vector connecting two beads corresponding to the phosphate group and the choline or amine group (PO4 and NC3 or NH3). To determine lipid tail order, we first computed the angles between the local membrane normal and the vectors connecting all adjacent beads in the lipid tails. Then we computed the average of the cosine of all angles and determined the average order parameter as

(4) P2=123cosθi2-1

where the average goes over all bond angles in both lipid tails. This definition of the lipid order turned out to be more sensitive for separating values in porated and non-porated regions, compared to the more common definition P2=123cos2θi-1 (Ingólfsson et al., 2017). We also extracted the type of lipid at a given vertex and the charge of this lipid. The scripts for running MemSurfer in this study, as well as the following analysis carried out in Matlab R2021a and described below, are available at https://github.com/learems/Electroporation-CGmem-MemSurfer, (copy archived at swh:1:rev:3e7d2d393b8ec98e04aa1e8915b55ee024840b97; Rems, 2021a).

MemSurfer returns the values of the above-listed properties/features at vertices, which correspond to individual lipids’ headgroups. For a given trajectory frame, the values of most of the properties of interest like area per lipid, lipid order parameter, charge, … can exhibit large variations among the adjacent vertices. Therefore, to extract the values at a selected (xp,yp) position, we used a Gaussian smoothing kernel, which determined a weighted average of the values Vi at all vertices within a radius 3rsmooth from the (xp,yp) position

(5) Vsmoothxp,yp=iViwiiwi

where the weights for a given adjacent vertex i are

(6) wi=exp(12(xixp)2r2smooth12(yiyp)2r2smooth)

When determining the presence/density of lipid groups (PC, PE, SM, etc.), we assigned a value of 0 or 1 if a given lipid group was located at a given vertex. Gaussian smoothing was carried out in the same way as for other features. The same smoothing procedure was also done for plotting the surface plots in Figure 5. We chose rsmooth = 1.5 nm, because this resulted in the best separation of the feature values in non-porated and porated regions.

The porated locations were defined as the locations at which a pore formed in any of the 60 electroporation simulations. For non-porated locations, we first divided the membrane into a grid with 31 by 31 points. We excluded all points which were within ~2 nm (6.7% of the box size) of any porated location. We also randomly excluded excess points such that the final number of non-porated locations was 300 (Figure 3A).

The distances between histograms of feature values in porated and non-porated locations (Figure 3B) were quantified by symmetrized version of the Kullback–Leibler divergence (Fleetwood et al., 2020) after performing a kernel smoothing probability density estimate:

(7) distance=±12x1x2pxlogp(x)q(x) +qxlogq(x)p(x)dx

where p(x) and q(x) are the distributions of the data in porated and non-porated locations, respectively. The limits x1 and x2 were defined as the lowest and highest value of x, where either of the distributions fell under 1% of their peak value. In Figure 3B, the distances correspond to probability density estimates, which were obtained after averaging the value of a given feature and at a given point over both membrane leaflets. Distances computed for each leaflet separately are shown in Figure 3—figure supplement 2.

Machine learning methods

Request a detailed protocol

The values of all the features extracted in non-porated and porated locations from 101 frames of a 10-ns-long trajectory before electroporation were used to train three machine learning models with Python and scikit-learn library (Pedregosa et al., 2011): random forest, support vector machine and multilayer perceptron neural network. The input data contained the (x,y) locations on the membrane and the extracted membrane properties at those locations (i.e., features, denoted as X). The input data were separated into two classes with label Y = 1 for porated locations and Y = 0 for non-porated locations. As the number of non-porated locations (300) and porated locations (~60) was imbalanced, we used SMOTE in Python to randomly oversample the data in porated locations and balance them with the data in non-porated locations. After oversampling the data, we trained the above-listed machine learning models and used them to predict the pore locations. The accuracy of the prediction was evaluated as:

(8) accuracy=no. correctly predicted points (porated and nonporated) in all framesno. all points × no. all frames

The feature importance score was determined using the attribute feature_importances_ of the scikit-learn random forest model. The codes used for this analysis are available at https://github.com/learems/Electroporation-CGmem-MachineLearning, (copy archived at swh:1:rev:1e32cf42c04af4543fb101ece2c92f35d793d62e; Rems, 2021b).

Bayesian survival analysis

Request a detailed protocol

To model the poration rates, we used Bayesian survival analysis. Given the non-exponential shape of the first poration times and the heterogeneity of poration rates between systems, it was crucial to use a model that allowed for an arbitrary time dependence of the rate and that was able to account for the system dependence. Cox’s proportional hazards model has these features and defines the event rate as

(9) λ(t)=λ0texpiβixi

λ0(t) is constructed as a piecewise constant function defined in time intervals with endpoints: 0 ≤ s0≤…≤ sN such that λ0(t) = λj for sjt < sj+1. xi are binary categorical variables such that if the poration time describes a system k, then xii=k = 1 and xi≠kk0, that is, one-hot-encoding representations of the system category. A good signal-to-noise ratio was found with 1.5 ns time intervals. λj were given independent priors in the form of gamma distributions with a shape parameter of 50 and a scale parameter of 10 ns–1. βi were given normally distributed priors centred at the origin and with a standard deviation of 100. Both priors are fairly uninformative.

The Bayesian inference of the model parameters {λi, βi } was done using as input features the first poration times and the system category (APM-dep, APM-hyp, BPM-dep, BPM-hyp) encoded in xi. In this way, the model was inferred on the electroporation events of all systems and a posteriori the model for a particular system i, λ0(t) exp(βj) was computed. Posterior predictive checks validated the quality of our models: the model can generate data that reproduce accurately the observed data (Figure 5—figure supplement 1).

To test the universality of λ0(t) and check the time independence of βi, we inferred another model allowing βi to be time dependent:

(10) λ(t)=λ0texpiβitxi

In this alternative model, βi(t) is also a piecewise constant function analogously to λ0(t). The heights of the steps of βi,j for sjt < sj+1 were given gaussian random walk priors: βi,j is given a gaussian prior with standard deviation of 1 but centred on their previous βi,j−1.

Since we are using Bayesian statistics, the parameters of the model {λj, βi} are treated as random variables whose distributions are inferred as a result. We can calculate their 94% credible intervals, that is, the interval containing 94% of the probability density of the variable. Therefore, uncertainties are presented as the median of the variable ± the distance to the credible interval limits.

The models were inferred using PYMC3 (Salvatier et al., 2016) and follow a similar methodology to one of the library’s case studies (available here). The full details of the implementation are available at https://github.com/sperezconesa/electroporation_modeling, (copy archived at swh:1:rev:1cc26d28bd76b5c254b4cbbf85703b2394aff775; Pérez-Conesa, 2022).

Data availability

The current manuscript is a computational study. Data required to reproduce the results is available at https://osf.io/fv98a/. Analysis codes are available at https://github.com/learems/Electroporation-CGmem-PoreLocalization (copy archived at swh:1:rev:af74428dc722ea4ecdb3e1f1c62178ea3f91b65b) , https://github.com/learems/Electroporation-CGmem-MemSurfer (copy archived at swh:1:rev:3e7d2d393b8ec98e04aa1e8915b55ee024840b97), https://github.com/learems/Electroporation-CGmem-MachineLearning, (copy archived at swh:1:rev:1e32cf42c04af4543fb101ece2c92f35d793d62e), and https://github.com/sperezconesa/electroporation_modeling, (copy archived at swh:1:rev:1cc26d28bd76b5c254b4cbbf85703b2394aff775).

The following data sets were generated
    1. Rems L
    2. Perez-Conesa S
    3. Delemotte L
    (2021) Open Science Framework
    ID fv98a. Electroporation through CG-MD, machine learning and bayesian survival analysis.

References

  1. Book
    1. Kinosita K
    2. Hibino M
    3. Itoh H
    4. Shigemori M
    5. Hirano K
    6. Kirino Y
    7. Hayakawa T
    (1992)
    Events of membrane electroporation visualized on a time scale from microsecond to seconds
    In: Chang DC, Chassy BM, Saunders JA, Sowers AE, editors. Guide to Electroporation and Electrofusion. San Diego: Academic Press. pp. 29–46.
    1. Pedregosa F
    2. Varoquaux G
    3. Gramfort A
    4. Michel V
    5. Thirion B
    6. Grisel O
    7. Blondel M
    8. Prettenhofer P
    9. Weiss R
    10. Dubourg V
    11. Vanderplas J
    12. Passos A
    13. Cournapeau D
    14. Brucher M
    15. Perrot M
    16. Duchesnay É
    (2011)
    Scikit-learn: Machine Learning in Python
    Journal of Machine Learning Research 12:2825–2830.

Decision letter

  1. Qiang Cui
    Reviewing Editor; Boston University, United States
  2. José D Faraldo-Gómez
    Senior Editor; National Heart, Lung and Blood Institute, National Institutes of Health, United States
  3. Helgi I Ingolfsson
    Reviewer
  4. Rumiana Dimova
    Reviewer; Max Planck Institute of Colloids and Interfaces, Germany

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Identification of electroporation sites in the complex lipid organization of the plasma membrane" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by Qiang Cui as Reviewing Editor and José Faraldo-Gómez as Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Helgi I. Ingolfsson (Reviewer #2); Rumiana Dimova (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) Additional discussion on the magnitude of voltage/electric field relative to typical experimental values and implications on the simulation results.

2) Analysis of simpler membranes to explicitly establish key factors that dictate poration kinetics and mechanism.

The reviewers also raise a number of issues related to the validation of the coarse-grained model, presentation of technical details and wording in certain discussions.

Reviewer #2 (Recommendations for the authors):

I realize validation of the coarse-grained results are likely out of scope for this already quite extensive manuscript (i.e. currently a number of the lipids in the complex membranes used don't have AA force field parameters and a in-between complex mixture would have to be used, etc, etc), therefore, I tried to suggest the validation as future directions. I do have a number of more specific questions and comments listed below:

One perception and/or wording comment. The brain is commonly thought to be enriched in polyunsaturated lipids. But you write that APM has more PU than the BPM:

L100 "Compared with BPM, APM has lower fraction of cholesterol, higher fraction of phosphatidylcholine (PC) lipids and higher fraction of monounsaturated (MU) and polyunsaturated (PU) lipids in the top leaflet." and

L391 "Compared with BPMs, APMs contain a smaller fraction of cholesterol and fully saturated lipids and a greater fraction of polyunsaturated lipids,"

This looks to be correct from how you defined the PU group (including also double unsaturated tails), compared to the 0, 1, 2, 3+ grouping in Figure 1A in Ingolfsson et al., 2017 (ref 41). Also, same paper Table 1, shows that on average (for the lipid fraction) there are more unsaturations per tail in the BPM than in the APM. Maybe some explanation or redefinition could prevent confusion for future readers.

Somewhat ill-defined wording comment. When first reading the manuscript, I got the sense that the Bayesina analysis was giving different result than the lipid feature/ML analysis. After reading it more in-depth I understand that they show different results on different "things" propensity/time vs spatial location and therefore are complementary. Maybe some rewording in the manuscript can make this clearer. I think my biggest stumbling block was "however" instead of something like additionally.

P1,L22, "learning. However, by" However ….

P2 "However, by analysing poration kinetics with Bayesian survival analysis we then show that"

I really like the analysis in section 2.4, but confused by one thing. It looks like form Figure 5 – sup. 1 and Figure 5B that the model is unable to capture the initial delay in pore formation time, the initial ~2.5 ns discussed later due to the bilayer expansion. Is this my misunderstanding or due to granularity of the binned data and/or lack of x-d offset term in β?

P16, L354 "Since membranes are practically volume-incompressible, the increase in area is directly related to a decrease in thickness.". Did you check this in your case? i.e. this is true for actual bilaye area, and projected area for very small flat bilayers but not larger undulating bilayers. I know both the APM and BPM were equilibrated with Z restrains to make them mostly flat but Ingolfsson et al., 2017 Figure 4A shows they are not totally flat and the BPM has more local undulations.

P16, L355 "Membrane thinning is expected to facilitate pore formation, as water molecules need to travel shorter distance when bridging the membrane.". Yes, but is it possible to utilize the power of the simulations and check for water bridging more directly? Here I realize this might be outside the scope of this manuscript, but have you considered looking at local bilayer properties such as "defects" (water penetration to the tails) or local min water distance across the bilayer (water dependent and not lipids dependent like the very related bilayer thickness)?

P9, L215+. I was quite confused reading that paragraph (L215 to L226), granted I am not an ML expert, but maybe some further explanation is in order. Also, could the questions of overfitting vs need for more data for higher accuracy be solved by increasing the number of depolarizable/hyperpolarizable simulation repeats? and/or leaving out some data for validation checks?

Reviewer #3 (Recommendations for the authors):

1. One aspect that remains unclear in this work relates to how universal the reported findings are and whether the large fraction of species present in the simulations is essential to the reported outcome, especially regarding the dependence on membrane elasticity. Indeed, to establish this correlation, it would have been interesting, using the selected approaches to first explore simpler membranes, and in particular, directly compare the poration probability and kinetics as a function of composition and membrane leaflet asymmetry. Membrane compositions with only a few representatives of the major classes of components would have been a good start as this would allow direct comparison to experimental results and would be helpful to resolve the role of the main actors involved.

2. The introduction (paragraph starting at line 70) explains that experimentally, pores in lipid membranes are difficult to assess. Work on giant unilamellar vesicles showing direct imaging of pores in lipid membranes suggests the opposite and should be referenced here (see e.g. DOI: 10.1529/biophysj.104.050310 and DOI: 10.1002/advs.202004068).

3. It will be good to explicitly discuss the following: The poration of the membranes in the simulations occurs within the first 15 nanoseconds, implying that the findings apply to pulses of high amplitude and nanosecond duration, but not necessarily to pulses of lower field strength and in the micro- or millisecond range duration (i.e. including conditions used in medical applications). A comment about this is due.

4. The authors indicate that such short poration times are needed to minimize lateral diffusion and allow mapping of local membrane features. Presumably, lipid diffusion and mixing which is substantial on longer times scales, might jeopardize the validity of the findings for longer pulse duration. Could the authors introduce a discussion also about this aspect?

5. The authors should explicitly specify the solution in which the membranes are simulated and the ionic strength.

6. The finding that gangliosides are quite important for poration (Figure 4B) is very interesting. The authors should discuss possible reasons and potential implications. Along these lines (even though not directly comparable), there has been a recent report on the poration of GM1-doped vesicles which exhibit much longer pore lifetimes compared to PC membranes (DOI: 10.1073/pnas.1722320115).

7. The authors should explicitly specify whether the size of the simulation box adjusts to accommodate the area of the pores in the membrane and whether more than one pore are simultaneously detected in one membrane patch. In line 379, the authors state that one can assume that Ai ≈Aj because all membranes have the same total area. Do they refer to the area before or after the application of the field? Could the authors further clarify the connection between the number of possible pore nucleation sites (as suggested on line 373) and the total area of the membranes?

8. Similarly, how does the pore size compare to the size of the data points displayed in Figure 2. In the caption of Figure 2, to avoid confusion, the authors should specify that all points correspond to the first poration event.

9. Figure 3 —figure supplement 3: The authors state that these quantities are computed from equation (1), however the caption above the graphs indicates that these two plots correspond to data where either no field was applied or a non-porating field was applied (i.e., situations where no pore is supposed to form in the membrane). How were these two diagrams plotted as they compare non porated with porated locations? Please clarify this in the caption.

10. The data for change in the membrane area as a function of the square of the nonporating field strength (Figure 5D and supplement 4) should be discussed in terms of the stretching elasticity modulus of the membranes. Typical stretching elasticity moduli of single-lipid membranes lie around 250 mN/m (see e.g. DOI: 10.1016/S0006-3495(00)76295-3). To claim the validity of the reported correlation between poration and mechanics, it has to be clarified, whether the applied simulation force fields can be used to correctly reflect the membrane elasticity. To address this (if not already reported for the selected force fields), the authors should measure the stretching elasticity for single-lipid membranes and compare with experimental values.

11. Line 270: "is fairly time independent" should be reformulated and justified by further discussion. As the authors refer to the values of β_i of Figure 5 – Supplement figure 2, the trend in this figure shows two main problems: (i) the values are positive for the whole time interval which is inconsistent with Figure 5 showing only negative values, and (ii) the values look actually time dependent as for instance APM-hyp shows a drop of 25% of its initial value on the time interval where its related probability density displayed Figure 5 A is non zero. Please discuss this issue.

https://doi.org/10.7554/eLife.74773.sa1

Author response

Essential revisions:

1) Additional discussion on the magnitude of voltage/electric field relative to typical experimental values and implications on the simulation results.

This has been addressed, as described in response to comment (R1.1).

2) Analysis of simpler membranes to explicitly establish key factors that dictate poration kinetics and mechanism.

We would like to argue that key factors that dictate poration kinetics and mechanism have already been established on simpler membranes. However, from studies on simple systems it is difficult to predict what will be the key factors dictating poration in a complex lipid mixture such as the plasma membrane. Our study is specifically designed with the aim to explore poration in a complex plasma membrane mimic. In addition, we develop a methodology that will allow us to establish quantitative models for the poration kinetics in plasma membranes at the whole-cell level.

The results presented in our study (especially Figure 3B, Figure 4, and Figure 5D,E) show that the following factors influence poration:

1. Lipid headgroup architecture (PC lipids are most prone for poration and gangliosides least prone for poration)

2. Lipid tail architecture (polyunsaturated lipids are more prone for poration than mono- and fully saturated lipids)

3. Cholesterol content and lipid tail order (less cholesterol and less order more likely poration)

4. Membrane thickness (thinner membrane more likely poration)

5. Membrane area compressibility or stretching modulus KA (lower KA more likely poration)

Below we provide a brief literature review summarizing previous findings on simpler lipid systems, using molecular dynamics (MD) simulations, experiments on giant unilamellar vesicles (GUVs), or experiments on planar lipid bilayers. This literature review shows that our findings are fully in agreement with previous findings on simpler systems. As Reviewer 1 correctly pointed out, our findings of what influences poration are all somewhat expected qualitatively. However, the literature review shows that there are many different factors influencing poration and it is quite impossible to predict which ones will be most important in a complex mixture. Consequently, our study fills this knowledge gap.

Influence of lipid headgroup and tail architecture

Previous studies showed that poration propensity in membranes containing fluid phase lipids depends on the architecture of both the lipid head groups and tails. Polak et al., studied electroporation of monocomponent bilayers consisting of lipids with different tails and headgroups using atomistic MD simulations (Polak et al., 2014, 2013). They demonstrated a considerate effect of methyl branches in the lipid tails, as well as the type of linkage between the head group and the carbonyl region. They observed that poration propensity decreases, respectively, in linear-chained DPPC lipids, methyl-branched DPhPC with ester linkages, and DPhPC with ether linkages, all in the fluid phase. Based on their analysis, they proposed that the presence of methyl branches could reduce the mobility of water molecules in the hydrophobic core and hence reduce the poration propensity. Note that DPhPC lipids are not present in our average and brain plasma membranes; however, knowing the characteristics of this lipid is important for understanding certain results referenced later. Polak et al., further studied electroporation of archaeal lipids, which have the same tail structure as DPhPC-ether lipids, whereas the archaeal head groups are formed by large sugar moieties (Polak et al., 2014). Compared with DPhPC-ether, archaeal lipids are more difficult to porate, associated with stronger interactions between the archaeal head groups. Similarly, Piggot et al., found that outer membranes of Gram-negative bacteria, which contain lipopolysaccharides, are more difficult to porate compared to membranes of Gram-positive bacteria, which are devoid of lipopolysaccharides (Piggot et al., 2011). These findings are consistent with our finding that gangliosides, which contain large sugar moieties, are more difficult to porate compared to other lipids in our membranes.

Atomistic MD studies of Gurtovenko and Lyulina showed that POPE bilayer is more difficult to porate than POPC (Gurtovenko and Lyulina, 2014). Lower poration propensity has been attributed to the primary amines in the POPE head groups capable of intra and intermolecular hydrogen bonding, in contrast to the choline moieties in the POPC head groups. POPE lipids are thus packed more densely than the POPC lipids, which hinders the penetration of water molecules in the bilayer and slows down the reorientation of the lipid head groups into the pore. Mixing these two lipids in an asymmetric bilayer (POPE in one and POPC in the other leaflet) results in an intermediate poration propensity with respect to pure POPC and POPE. Our average and brain plasma membranes contain an asymmetric mixture of PC and PE (as well as other) lipids; therefore, our analysis shows comparable poration propensity of PC and PE lipids.

Mixing different lipids does not always influence the poration propensity. Through atomistic MD simulations, Levine and Vernier showed that 30% of POPS lipids in one leaflet of an otherwise POPC bilayer does not have a significant effect on poration propensity compared to pure POPC (Levine and Vernier, 2012). Consistent with their findings, we did not see considerable influence of PS lipids on poration propensity.

Hu et al., determined the free energy for pore formation in 18 monocomponent coarse-grained Martini bilayers (Hu et al., 2015). They studied fully saturated and monounsaturated lipids with different chain lengths containing PC, PE, PG or PE headgroup. Their analysis confirmed that poration propensity depends on the architecture of both lipid headgroups and tails in coarse-grained membranes, similarly as in atomistic membranes. They found slightly lower poration propensity in monounsaturated lipids compared to fully saturated lipids. However, Breton and Mir showed that GUVs from polyunsaturated DLiPC lipids are considerably more prone to poration compared to monounsaturated DOPC lipids (Breton and Mir, 2018). The finding of Hu et al., and Breton and Mir are consistent with our finding that pores strongly favor polyunsaturated lipids compared to both fully saturated and monounsaturated lipids.

Interesting are also the measurements of van Uitert et al., on planar lipid bilayers, which shows that mixing polyunsaturated PE or π lipids with DPhPC lipids increases the poration propensity compared to pure DPhPC membranes (van Uitert et al., 2010). This shows how the architecture of lipid tails can become more important for poration than the architecture of lipid headgroups (as mentioned above, PE and π lipids are less prone for poration than PC lipids containing the same tail). These results support our finding that pores preferentially form in regions enriched with polyunsaturated lipids.

Influence of cholesterol and lipid order

As cholesterol is added to the system, the lipids organize in the liquid-ordered phase. The cholesterol organizes itself in the hydrophobic core of the bilayer, where it can condense the lipids and it can alter the mechanical properties of the membrane, such as the thickness, the bending stiffness and the fluidity. However, the addition of cholesterol does not always lead to the same results. Depending on the concentration of the cholesterol and the architecture of the lipid, cholesterol can either decrease or raise the electroporation threshold (Portet and Dimova, 2010). Also, mechanical studies on bilayers have shown the non-universal and lipid-specific effect of cholesterol (Gracià et al., 2010; Pan et al., 2008). Studies of Mauroy et al., have shown that an increasing concentration of cholesterol on POPC decreases the poration propensity, whereas this increased cholesterol shows no significant influence poration propensity of Egg PC (Mauroy et al., 2015). Similar results for Egg PC have been shown before by Portet and Dimova (Portet and Dimova, 2010). In addition, they reported that increasing cholesterol could increase poration propensity for DOPC vesicles. Surprisingly, the experimental results on the effect of cholesterol on poration propensity of different lipid bilayers have not been fully supported by MD simulations. Simulations on the effect of cholesterol on poration propensity of POPC show similar results as found experimentally on GUVs (Casciola et al., 2014). Nevertheless, MD simulations of Fernandez et al., on DOPC showed an decrease of poration propensity when adding cholesterol (Fernández et al., 2010), which is in disagreement with the experimental results on GUVs (Portet and Dimova, 2010). Overall, the influence of cholesterol on poration propensity of a lipid bilayer is non-universal and strongly dependent on the architecture of the lipids.

By mixing two different lipids together with cholesterol, coexisting liquid-ordered and liquid-disordered phases can occur in the membrane. Liquid-ordered domains contain saturated lipids and cholesterol and liquid-disordered domains contain unsaturated lipids and possibly a low level of cholesterol. Van Uitert et al., studied this effect of cholesterol on poration propensity in planar bilayers made from binary lipid mixtures (van Uitert et al., 2010). They observed that the effect of cholesterol on poration propensity depends on the cholesterol percentage. At low percentages, poration propensity increased slightly with respect to pure binary mixture without cholesterol. However, above a certain threshold percentage, the poration propensity decreased together with the increase in cholesterol. From the experimental results it is difficult to interpret the molecular mechanisms of this biphasic influence of cholesterol percentage on poration propensity. With MD simulations on heterogeneous membranes made from ternary mixtures, Reigada showed that the probability of pore formation is highest in the middle of the liquid disordered phase (Reigada, 2014). Reigada’s observations have been corroborated by Sengel and Wallace by visualizing pores in droplet-interface bilayers containing liquid-order and liquid-disordered domains (Sengel and Wallace, 2016). In agreement with the findings described above, we found that in the average and brain plasma membranes, pores tend to avoid regions enriched with cholesterol and prefer locations with lower lipid order.

Influence of membrane thickness and membrane area compressibility modulus

Atomistic MD simulations by Ziegler and Vernier indicated that for fluid-phase fully-saturated or monounsaturated lipids with a PC head group, poration propensity tends to decrease with membrane thickness and thus with increasing the chain length of the lipid tails (Ziegler and Vernier, 2008). Similarly, Hu et al., determined the free energy for pore formation in 18 monocomponent coarse-grained Martini bilayers (Hu et al., 2015) and showed correlation between the probability for poration and membrane thickness, with thinner membranes tending to be more easily porated (Author response image 1) . Note, however, that around a given value of thickness the spread of energy barriers for poration is considerable (Author response image 1) , showing that thickness is not the only parameter influencing poration, as the lipid headgroup and tail architecture have an influence as well. Consistent with the thickness not being the only factor governing poration, we see greater poration of average plasma membranes, which are slightly thicker than brain plasma membranes (~3.8 nm and ~3.5 nm, respectively). Similarly, we found that the local membrane thickness is not a good predictor of pore locations (note also that the variation of the local thickness in our membranes is small, only about ±0.25 nm). However, we found that for a given type of membrane, the poration propensity is greater for membranes that can be thinned more by the electric field. The change in membrane thickness caused by electric field is associated with an increase in the membrane area. The ability of a membrane to change its area under electric field is inversely proportional to the membrane area compressibility (or stretching) modulus. Needham and Hochmuth studied GUVs made from SOPC lipids, SOPC:cholesterol mixture, and lipids extracted from red blood cells. They showed that the electroporation propensity increases with decreasing membrane area compressibility modulus, consistent with our results (Needham and Hochmuth, 1989).

Author response image 1
Correlation between membrane thickness and the free energy barrier for creating a 3-nm-radius pore in a lipid bilayer.

Each dot corresponds to one of 18 monocomponent coarse-grained Martini membranes, made from fully saturated and monounsatured lipids with different chain lengths and different headgroups (PC, PE, PS, and PG). Reproduced from data tabulated in (Hu et al., 2015).

References

Breton M, Mir LM. 2018. Investigation of the chemical mechanisms involved in the electropulsation of membranes at the molecular level. Bioelectrochemistry 119:76–83. doi:10.1016/j.bioelechem.2017.09.005

Casciola M, Bonhenry D, Liberti M, Apollonio F, Tarek M. 2014. A molecular dynamic study of cholesterol rich lipid membranes: comparison of electroporation protocols. Bioelectrochemistry, Bio-Electroporation organised by COST TD1104 100:11–17. doi:10.1016/j.bioelechem.2014.03.009

Fernández ML, Marshall G, Sagués F, Reigada R. 2010. Structural and Kinetic Molecular Dynamics Study of Electroporation in Cholesterol-Containing Bilayers. J Phys Chem B 114:6855–6865. doi:10.1021/jp911605b

Gracià RS, Bezlyepkina N, Knorr RL, Lipowsky R, Dimova R. 2010. Effect of cholesterol on the rigidity of saturated and unsaturated membranes: fluctuation and electrodeformation analysis of giant vesicles. Soft Matter 6:1472–1482. doi:10.1039/B920629A

Gurtovenko AA, Lyulina AS. 2014. Electroporation of asymmetric phospholipid membranes. J Phys Chem B 118:9909–9918. doi:10.1021/jp5028355

Hu Y, Sinha SK, Patel S. 2015. Investigating hydrophilic pores in model lipid bilayers using molecular simulations: Correlating bilayer properties with pore-formation thermodynamics. Langmuir 31:6615–6631. doi:10.1021/la504049q

Levine ZA, Vernier PT. 2012. Calcium and Phosphatidylserine Inhibit Lipid Electropore Formation and Reduce Pore Lifetime. J Membrane Biol 245:599–610. doi:10.1007/s00232-012-9471-1

Mauroy C, Rico-Lattes I, Teissié J, Rols M-P. 2015. Electric Destabilization of Supramolecular Lipid Vesicles Subjected to Fast Electric Pulses. Langmuir 31:12215–12222. doi:10.1021/acs.langmuir.5b03090

Needham D, Hochmuth RM. 1989. Electro-mechanical permeabilization of lipid vesicles. Role of membrane tension and compressibility. Biophys J 55:1001–1009.

Pan J, Mills TT, Tristram-Nagle S, Nagle JF. 2008. Cholesterol Perturbs Lipid Bilayers Nonuniversally. Phys Rev Lett 100:198103. doi:10.1103/PhysRevLett.100.198103

Piggot TJ, Holdbrook DA, Khalid S. 2011. Electroporation of the E. coli and S. aureus membranes: Molecular dynamics simulations of complex bacterial membranes. J Phys Chem B 115:13381–13388. doi:10.1021/jp207013v

Polak A, Bonhenry D, Dehez F, Kramar P, Miklavčič D, Tarek M. 2013. On the Electroporation Thresholds of Lipid Bilayers: Molecular Dynamics Simulation Investigations. J Membrane Biol 246:843–850. doi:10.1007/s00232-013-9570-7

Polak A, Tarek M, Tomšič M, Valant J, Ulrih NP, Jamnik A, Kramar P, Miklavčič D. 2014. Electroporation of archaeal lipid membranes using MD simulations. Bioelectrochemistry, Bio-Electroporation organised by COST TD1104 100:18–26. doi:10.1016/j.bioelechem.2013.12.006

Portet T, Dimova R. 2010. A New Method for Measuring Edge Tensions and Stability of Lipid Bilayers: Effect of Membrane Composition. Biophysical Journal 99:3264–3273. doi:10.1016/j.bpj.2010.09.032

Reigada R. 2014. Electroporation of heterogeneous lipid membranes. Biochimica et Biophysica Acta (BBA) – Biomembranes 1838:814–821. doi:10.1016/j.bbamem.2013.10.008

Sengel JT, Wallace MI. 2016. Imaging the dynamics of individual electropores. PNAS 113:5281–5286. doi:10.1073/pnas.1517437113

van Uitert I, Le Gac S, van den Berg A. 2010. The influence of different membrane components on the electrical stability of bilayer lipid membranes. Biochimica et Biophysica Acta (BBA) – Biomembranes 1798:21–31. doi:10.1016/j.bbamem.2009.10.003

Vernier PT, Levine ZA, Gundersen MA. 2013. Water bridges in electropermeabilized phospholipid bilayers. Proceedings of the IEEE 101:494–504. doi:10.1109/JPROC.2012.2222011

Ziegler MJ, Vernier PT. 2008. Interface Water Dynamics and Porating Electric Fields for Phospholipid Bilayers. J Phys Chem B 112:13588–13596. doi:10.1021/jp8027726

The reviewers also raise a number of issues related to the validation of the coarse-grained model, presentation of technical details and wording in certain discussions.

We have addressed all Reviewers’ comments, as described on the following pages.

Reviewer #2 (Recommendations for the authors):

I realize validation of the coarse-grained results are likely out of scope for this already quite extensive manuscript (i.e. currently a number of the lipids in the complex membranes used don't have AA force field parameters and a in-between complex mixture would have to be used, etc, etc), therefore, I tried to suggest the validation as future directions. I do have a number of more specific questions and comments listed below:

One perception and/or wording comment. The brain is commonly thought to be enriched in polyunsaturated lipids. But you write that APM has more PU than the BPM:

L100 "Compared with BPM, APM has lower fraction of cholesterol, higher fraction of phosphatidylcholine (PC) lipids and higher fraction of monounsaturated (MU) and polyunsaturated (PU) lipids in the top leaflet." and

L391 "Compared with BPMs, APMs contain a smaller fraction of cholesterol and fully saturated lipids and a greater fraction of polyunsaturated lipids,"

This looks to be correct from how you defined the PU group (including also double unsaturated tails), compared to the 0, 1, 2, 3+ grouping in Figure 1A in Ingolfsson et al., 2017 (ref 41). Also, same paper Table 1, shows that on average (for the lipid fraction) there are more unsaturations per tail in the BPM than in the APM. Maybe some explanation or redefinition could prevent confusion for future readers.

Indeed, we grouped lipids based on unsaturation in a different was as Ingolfsson et al., 2017. The reason for this is related to lipid oxidation and its role in electroporation, as we discuss in Section 3.2. We added the following explanation to Table 1:

* The way in which lipids are grouped by tail saturation is different than in (Ingólfsson et al., 2017), where the grouping is based on the total number of double bonds in both lipid tails. Here, the grouping is motivated by the role of lipid oxidation in electroporation, whereby lipids tails containing two or more double bonds are considerably more prone to oxidative damage than tails containing a single double bond. This is because bis-allylic hydrogens are much more easily abstracted by free radicals compared to allylic hydrogens (Reis and Spickett, 2012). Furthermore, membranes made of polyunsaturated lipids (by our definition) were found to be considerably more prone to poration/rupture by mechanical stretching compared to membranes made of lipids containing a single bond in one or both lipid tails (Olbrich et al., 2000). Thus, we consider that a lipid is polyunsaturated only if it contains at least one polyunsaturated tail.

From the first paragraph in Results we deleted the following sentence, which was in any case incomplete, because there are further differences in the compositions, visible in Figure 1 but not mentioned in the sentence:

“Compared with BPM, APM has lower fraction of cholesterol, higher fraction of phosphatidylcholine (PC) lipids and higher fraction of monounsaturated (MU) and polyunsaturated (PU) lipids in the outer leaflet.”

We also added the following sentence to Section 3.4:

“Note that the polyunsaturated lipids in BPMs have on average greater number of double bonds than in APMs (Ingólfsson et al., 2017); however, this is not enough to balance the poration rate with APMs.”

Somewhat ill-defined wording comment. When first reading the manuscript, I got the sense that the Bayesina analysis was giving different result than the lipid feature/ML analysis. After reading it more in-depth I understand that they show different results on different “things” propensity/time vs spatial location and therefore are complementary. Maybe some rewording in the manuscript can make this clearer. I think my biggest stumbling block was “however” instead of something like additionally.

P1,L22, “learning. However, by” However ….

P2 “However, by analysing poration kinetics with Bayesian survival analysis we then show that”

As suggested, we changed the word “however” to “additionally”.

I really like the analysis in section 2.4, but confused by one thing. It looks like form Figure 5 – sup. 1 and Figure 5B that the model is unable to capture the initial delay in pore formation time, the initial ~2.5 ns discussed later due to the bilayer expansion. Is this my misunderstanding or due to granularity of the binned data and/or lack of x-d offset term in β?

We thank the reviewer for the positive feedback. The delay time the reviewer is referring to is possibly the short initial time interval where there are no poration events (roughly the first nanosecond). This time is not captured by our model since it is smaller than our discretization grid. In contrast, ~3 ns discussed in the manuscript is the time the system takes to reach the steady state kinetics, the time at which l0(t) plateaus. During the pre-steady state regime poration does occur but at a lower rate.

Please note that we found an error in the script for plotting l0(t), which has now been corrected.

P16, L354 “Since membranes are practically volume-incompressible, the increase in area is directly related to a decrease in thickness.”. Did you check this in your case? i.e. this is true for actual bilaye area, and projected area for very small flat bilayers but not larger undulating bilayers. I know both the APM and BPM were equilibrated with Z restrains to make them mostly flat but Ingolfsson et al., 2017 Figure 4A shows they are not totally flat and the BPM has more local undulations.

We verified that the increase in the projected membrane area is correlated with a decrease in membrane thickness, whereby we determined the average membrane thickness by averaging the thickness values computed at all lipid positions by MemSurfer. We also estimated the change in membrane volume, by summing the products of area per lipid a and half of the membrane thickness d, computed at each lipid position by MemSurfer

ViNlipidsaidi2

where Nlipids is the number of all lipids in the membrane. The relative changes in projected area, average thickness, and volume for different magnitudes of nonporating electric fields are shown in Figure 5 —figure supplement 4. The relative change in membrane volume is indeed nonzero; however, it is much lower than the relative change in projected area or membrane thickness.

Figure 5 —figure supplement 5 shows the steady state values of the relative change in membrane thickness and how they depend on the electric field E, similarly as shown for projected membrane area in Figure 5 —figure supplement 6. The relative change in membrane thickness decreases linearly with E2.

Finally, we verified that the factors bi are inversely correlated with the relative change in membrane thickness, as shown in Figure 5D. Therefore, we can confirm that a greater poration rate is associated with membrane thinning.

Based on the comment, we made the following modifications to the manuscript:

– We added Figure 5 —figure supplement 4.

– We added Figure 5 —figure supplement 5.

– We added Figure 5D.

– We made slight modifications to the text to accomodate the additional results.

– We changed the problematic sentence in Section 3.3 from "Since membranes are practically volume-incompressible, the increase in area is directly related to a decrease in thickness.” to ”The increase in membrane area is directly related to a decrease in membrane thickness (Figure 5 —figure supplement 4).”

P16, L355 "Membrane thinning is expected to facilitate pore formation, as water molecules need to travel shorter distance when bridging the membrane.". Yes, but is it possible to utilize the power of the simulations and check for water bridging more directly? Here I realize this might be outside the scope of this manuscript, but have you considered looking at local bilayer properties such as "defects" (water penetration to the tails) or local min water distance across the bilayer (water dependent and not lipids dependent like the very related bilayer thickness)?

This is an excellent suggestion. We have thought about looking at water defects, but we first wanted to focus on the relationship between bilayer properties in the absence of an electric field and the probability of pore formation when exposing the membrane to an electric field. This is what we would in general like to do – predict the most likely poration sites in a membrane with a new composition without additional extensive simulations and analyses. As we obtained reliable prediction of porated and nonporated locations with machine learning models by just looking at local lipid densities, we decided to leave it at that for the present study.

We expect that the rate of formation of water defects depends on the electric field strength, similarly as the rate of formation of transmembrane water bridges and pores. Therefore, monitoring water defects in the absence of an electric field might turn out uninformative. Another variable that is presumably related to occurrence of water defects and worth further investigation is the local lateral pressure in the headgroup and tail region, which has already been correlated with poration rate in simpler atomistic bilayer systems (Polak et al., 2013, http://dx.doi.org/10.1016/j.bioelechem.2013.12.006; Gurtovenko and Lyulina, 2014, http://dx.doi.org/10.1021/jp5028355). We think this will be very interesting to study in the future, especially in combination with simulations on corresponding atomistic systems obtained after backmapping regions of interests from larger coarse-grained systems.

P9, L215+. I was quite confused reading that paragraph (L215 to L226), granted I am not an ML expert, but maybe some further explanation is in order. Also, could the questions of overfitting vs need for more data for higher accuracy be solved by increasing the number of depolarizable/hyperpolarizable simulation repeats? and/or leaving out some data for validation checks?

We have addressed comments R2.6 and R2.12 together. Please refer to R2.12 for our response.

Reviewer #3 (Recommendations for the authors):

1. One aspect that remains unclear in this work relates to how universal the reported findings are and whether the large fraction of species present in the simulations is essential to the reported outcome, especially regarding the dependence on membrane elasticity. Indeed, to establish this correlation, it would have been interesting, using the selected approaches to first explore simpler membranes, and in particular, directly compare the poration probability and kinetics as a function of composition and membrane leaflet asymmetry. Membrane compositions with only a few representatives of the major classes of components would have been a good start as this would allow direct comparison to experimental results and would be helpful to resolve the role of the main actors involved.

We thank the reviewer for the comment, and we fully agree that it will be very interesting to compare coarse-grained simulations on simpler systems with experiments on membranes with corresponding lipid composition, especially with the aim of carefully validating coarse-grained simulations for assessing electroporation.

However, probing what governs electroporation in simpler systems has already been performed, leading us to focus our efforts on complex systems. The literature review presented in R0.1 shows that there are many different factors influencing poration in simpler membranes and it its quite impossible to predict which ones will be most important in a complex mixture. Consequently, our study provides a way to fill this gap. The literature review, and our comments along, further show that our findings are fully in agreement with previous findings on simpler system.

To clarify this, we added the following text to Section 3.1:

“Overall, studies on simple model systems showed that there are many different factors influencing poration and it is quite impossible to predict which ones will be most important in a complex mixture such as the plasma membrane.

Our findings from membranes composed of >60 different lipid types mimicking the realistic composition of plasma membranes are fully in agreement with findings from simpler model systems.”

Regarding the dependence of poration probability on membrane elasticity: Our results show that the electroporation propensity increases with decreasing membrane area compressibility (or stretching) modulus. This is in agreement with experimental results of Needham and Hochmuth 1989 (https://dx.doi.org/10.1016%2FS0006-3495(89)82898-X), who studied electroporation thresholds in GUVs made from SOPC lipids, SOPC:cholesterol mixture, and lipids extracted from red blood cells.

In the revised manuscript, we have added the reference to the paper of Needham and Hochmuth in Section 3.3:

“Inverse correlation between electroporation propensity and membrane area compressibility modulus has also been found experimentally (Needham and Hochmuth, 1989).”

To motivate further studies combining simulations and experiments, in accordance with suggestions made by Reviewer #2, we added the following sentence at the end of Section 3.6.

“Nevertheless, we hope our study will motivate further exploration on the validity of coarse-grained membrane models for studying membrane electroporation, both in comparison with corresponding atomistic computational membrane models and experimental model membrane systems.”

2. The introduction (paragraph starting at line 70) explains that experimentally, pores in lipid membranes are difficult to assess. Work on giant unilamellar vesicles showing direct imaging of pores in lipid membranes suggests the opposite and should be referenced here (see e.g. DOI: 10.1529/biophysj.104.050310 and DOI: 10.1002/advs.202004068).

We agree. The following sentence, together with suggested references, has been added:

Pores have been imaged in giant unilamellar vesicles (Lira et al., 2021; Riske and Dimova, 2005); however, these pores have reached sizes on the order of 1 mm, which have not been observed in cell plasma membranes, likely because the actin cytoskeleton limits pore expansion (Perrier et al., 2019).

3. It will be good to explicitly discuss the following: The poration of the membranes in the simulations occurs within the first 15 nanoseconds, implying that the findings apply to pulses of high amplitude and nanosecond duration, but not necessarily to pulses of lower field strength and in the micro- or millisecond range duration (i.e. including conditions used in medical applications). A comment about this is due.

As a response to comment (R1.1) made by Reviewer 1 we now explain the relationship between the electric field used in simulations and experiments in Section 5.1, among which we also address comment R3.3:

“While the conditions in our study strictly correspond to nanoseconds-long exposures to electric field, we expect similar behaviour also for longer electric pulses. Nevertheless, on a longer time scale, electrodeformation of the cell membrane and/or increases in local membrane curvature caused by electric field might play additional roles in the pore formation process (Dimova et al., 2009; Perrier et al., 2017).”

4. The authors indicate that such short poration times are needed to minimize lateral diffusion and allow mapping of local membrane features. Presumably, lipid diffusion and mixing which is substantial on longer times scales, might jeopardize the validity of the findings for longer pulse duration. Could the authors introduce a discussion also about this aspect?

The membranes that we used in our study have already been equilibrated by Ingolfsson et al., in an 80 μs-long simulation, therefore the lipid mixing is at equilibrium when we start the simulation. Of course, if we kept the simulation of a chosen membrane, e.g. APM #1, running for tens of microseconds before applying an electric field, the local lipid organization would be somewhat different since domains laterally move along the membrane. However, our results from 4 APM membranes, all containing slightly different lipid organizations, are consistent, indicating that our results are still valid even if we consider some additional waiting time before pore creation starts.

The situation with the waiting time mentioned above resembles the situation when a cell is exposed to an electric field. From an electrical point of view, a cell behaves like a capacitor. Therefore, the transmembrane voltage gradually increases with time until it becomes sufficiently high to trigger pore creation. Author response image 2 shows the time course of electroporation at the pole of a cell (θ = 0), predicted by a model developed by DeBruin and Krassowska in 1999, doi: 10.1016/S0006-3495(99)76973-0. As soon as the transmembrane voltage exceeds certain value (~1.3 V in the model), pores are rapidly created, i.e. the pore density exhibits a jump. This is because pore creation rate depends exponentially on the transmembrane voltage (or the local electric field strength within the membrane). With increasing number of pores, the membrane conductivity also increases, which in turns starts decreasing the transmembrane voltage and stops pore creation. This model therefore shows that most pores are created on a nanosecond time scale, no matter how long the pulse is.

Author response image 2
Time course of the induced transmembrane voltage, pore density, and membrane conductivity at the cell pole (q = 0) upon exposure to 1-ms-long electric pulse.

Only the first 10 ms after the onset of the electric pulse are shown for clarity. The results are reproduced from the paper by DeBruin and Krassowska 1999, doi: 10.1016/S0006-3495(99)76973-0.

Consequently, we have no reason to think that lipid diffusion and mixing on longer times scales would jeopardize the validity of the findings for longer pulse duration.

5. The authors should explicitly specify the solution in which the membranes are simulated and the ionic strength.

Thank you for noticing that we forgot to report the NaCl concentration. We added a sentence to Section 5.1:

“The NaCl concentration was ~150 mM.”

6. The finding that gangliosides are quite important for poration (Figure 4B) is very interesting. The authors should discuss possible reasons and potential implications. Along these lines (even though not directly comparable), there has been a recent report on the poration of GM1-doped vesicles which exhibit much longer pore lifetimes compared to PC membranes (DOI: 10.1073/pnas.1722320115).

Our study shows that pores avoid domains enriched with gangliosides. This is expected according to previous atomistic simulations on phospholipid and archaeal lipid membranes by Polak et al., 2014 (Ref. 54), which showed that large sugar moieties in the lipid headgroups make pore formation more difficult, which we already discussed in Section 3.1. We unfortunately cannot comment on pore lifetimes, as we have not investigated the kinetics of pore closure, nor have we investigated bilayers without gangliosides.

7. The authors should explicitly specify whether the size of the simulation box adjusts to accommodate the area of the pores in the membrane and whether more than one pore are simultaneously detected in one membrane patch. In line 379, the authors state that one can assume that AiAj because all membranes have the same total area. Do they refer to the area before or after the application of the field? Could the authors further clarify the connection between the number of possible pore nucleation sites (as suggested on line 373) and the total area of the membranes?

The following sentence about adjustments of the box size was added to Section 5.1:

“The pressure coupling allows the size of the simulation box to adjusts to the changes in the membrane area.”

The information about multiple simultaneous poration events was already in Section 5.2:

“In up to ~20% simulations, two or more pores formed practically simultaneously. In such a case, we considered all of these pores for analysis.”

When assuming AiAj we consider the area before poration. Nevertheless, due to the exponential dependence of the pore formation rate on the energy barrier for pore formation, small variations in the membrane area are not very influential.

The theoretical description of the pore creation rate is based on ideas from the classical nucleation theory, where a pore is considered as the nucleus of a new phase. In oversimplified terms, each lipid is considered as a possible nucleation site. We refer to the following references, which provide further useful information:

Glaser, R W, S L Leikin, L V Chernomordik, V F Pastushenko, and A I Sokirko. ‘Reversible Electrical Breakdown of Lipid Bilayers: Formation and Evolution of Pores’. Biochimica et Biophysica Acta 940, no. 2 (24 May 1988): 275–87. https://doi.org/10.1016/0005-2736(88)90202-7.

Evans, Evan, and Volkmar Heinrich. ‘Dynamic Strength of Fluid Membranes’. Comptes Rendus Physique 4, no. 2 (March 2003): 265–74. https://doi.org/10.1016/S1631-0705(03)00044-6. (See their Section 3)

Vasilkoski, Zlatko, Axel T. Esser, T. R. Gowrishankar, and James C. Weaver. ‘Membrane Electroporation: The Absolute Rate Equation and Nanosecond Time Scale Pore Creation’. Physical Review E 74, no. 2 (3 August 2006): 021904. https://doi.org/10.1103/PhysRevE.74.021904.

8. Similarly, how does the pore size compare to the size of the data points displayed in Figure 2. In the caption of Figure 2, to avoid confusion, the authors should specify that all points correspond to the first poration event.

We did not evaluate the pore size, as all pores expand upon creation until the system explodes, which is typical for simulations at constant electric field strength and is related to the relatively small system size. The size of points in Figure 2 was also coding for poration time, similarly as the color of the points. To avoid confusion, we corrected the figure and made all points equal size.

We also added the following sentence to Figure 2 caption:

“All data points correspond to the first poration event.”

9. Figure 3 —figure supplement 3: The authors state that these quantities are computed from equation (1), however the caption above the graphs indicates that these two plots correspond to data where either no field was applied or a non-porating field was applied (i.e., situations where no pore is supposed to form in the membrane). How were these two diagrams plotted as they compare non porated with porated locations? Please clarify this in the caption.

Thank you for the comment. We had mistakenly referred to equation (2) instead of equation (7), which probably led to the confusion. This has now been corrected.

The figure caption already gives information on how the distances were obtained:

Distances between probability density estimates of individual features in porated and nonporated locations, calculated according to equation (7). … The graphs on the left side correspond to values extracted from a 10-ns-long trajectory before electric field application. The graphs on the right side correspond to values extracted from a 10-ns-long trajectory, where a non-porating electric field of +106.3 mV/nm (analysis for depolarization) or -106.3 mV/nm (analysis for hyperpolarization).

10. The data for change in the membrane area as a function of the square of the nonporating field strength (Figure 5D and supplement 4) should be discussed in terms of the stretching elasticity modulus of the membranes. Typical stretching elasticity moduli of single-lipid membranes lie around 250 mN/m (see e.g. DOI: 10.1016/S0006-3495(00)76295-3). To claim the validity of the reported correlation between poration and mechanics, it has to be clarified, whether the applied simulation force fields can be used to correctly reflect the membrane elasticity. To address this (if not already reported for the selected force fields), the authors should measure the stretching elasticity for single-lipid membranes and compare with experimental values.

In the manuscript we referred to the membrane area compressibility modulus in the last paragraph of Section 2.4. We now also added the term “stretching” modulus:

We found that the values of βi correlate with the ability of a membrane to thin and expand its area under electric field, i.e., they corelate inversely with the membrane area compressibility (or stretching) modulus (Figure 5D,E).

The membrane area compressibility modulus KA has been evaluated for coarse-grained Martini bilayers previously. For DPPC bilayer in fluid phase with size comparable to our membranes, Marrink et al., 2004 (https://doi.org/10.1021/jp036508g) obtained KA = 260 ± 40 mN/m, which compares well with the experimental estimate for DPPC KA = 250 of Nagle et al., 2000 (https://doi.org/10.1016/S0304-4157(00)00016-2). For POPC, Chacón et al., 2015 (http://dx.doi.org/10.1063/1.4926938) obtained a value of KA = 310 ± 20 mN/m, which is close to 272 mN/m reported by Janosi and Gorfe 2010 (https://doi.org/10.1021/ct100381g) for atomistic POPC bilayers described with Charmm force-field, and K = 277 ± 10 mN/m reported by Braun et al., 2013 (https://doi.org/10.1021/jp401718k) for a similar DOPC lipid described by united atom force field. It is also close to the experimental value of 265 ± 18 mN/m measured in DOPC giant unilamellar vesicles by Rawicz et al., 2000 (https://doi.org/10.1016/S0006-3495(00)76295-3).

11. Line 270: "is fairly time independent" should be reformulated and justified by further discussion. As the authors refer to the values of β_i of Figure 5 – Supplement figure 2, the trend in this figure shows two main problems: (i) the values are positive for the whole time interval which is inconsistent with Figure 5 showing only negative values, and (ii) the values look actually time dependent as for instance APM-hyp shows a drop of 25% of its initial value on the time interval where its related probability density displayed Figure 5 A is non zero. Please discuss this issue.

We thank the Reviewer for the valuable comments.

(i) During the revision we found a bug in the plotting code used to make Figure 5 —figure supplement 2. The corrected plot is shown in Author response image 1. After correcting the code, the βi(t) values of the time-dependent model are closer to those of the time-independent model. Nevertheless, the values of βi in the time-independent and time-dependent model cannot be directly compared.

We have added the following explanation to Figure 5 —figure supplement 2 caption:

“Note that the value of βi(t) are different from the ones obtained in the time-independent model in Figure 5. Since all coefficients are parametrized simultaneously with all the data, different combinations of magnitudes of λ0(t) and βi can result in an equivalent λi(t). Therefore, in order to compare how rates increase or decrease in two different systems i and j, only the relative change can be studied λi(t)/λj(t) = exp(βi)/exp(βj).”

For completeness, we also added the plot of λ0(t) for the time-dependent β model.

(ii) We agree with the reviewer that this should have been explained in more detail. We have included the following explanations to Section 2.4:

“The βi(t) obtained were found to slightly drop in the interval 0 < t < 3 ns until becoming practically time-independent (Figure 5 —figure supplement 2). Given that the βi(t) for all systems follow the same trend and that the time-dependent model only adds information in the initial stage of electroporation and has higher uncertainty, the time-independent model was selected for its easier interpretability and lower complexity.”

https://doi.org/10.7554/eLife.74773.sa2

Article and author information

Author details

  1. Lea Rems

    1. KTH Royal Institute of Technology, Dept. Applied Physics, Science for Life Laboratory, Solna, Sweden
    2. University of Ljubljana, Faculty of Electrical Engineering, Ljubljana, Slovenia
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft, Writing - review and editing
    For correspondence
    lea.rems@fe.uni-lj.si
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7470-4367
  2. Xinru Tang

    1. KTH Royal Institute of Technology, Dept. Applied Physics, Science for Life Laboratory, Solna, Sweden
    2. University of Chinese Academy of Sciences, Beijing, China
    Contribution
    Data curation, Formal analysis, Investigation, Visualization, Writing – original draft
    Contributed equally with
    Fangwei Zhao
    Competing interests
    No competing interests declared
  3. Fangwei Zhao

    1. KTH Royal Institute of Technology, Dept. Applied Physics, Science for Life Laboratory, Solna, Sweden
    2. University of Chinese Academy of Sciences, Beijing, China
    Contribution
    Data curation, Formal analysis, Investigation, Visualization, Writing – original draft
    Contributed equally with
    Xinru Tang
    Competing interests
    No competing interests declared
  4. Sergio Pérez-Conesa

    KTH Royal Institute of Technology, Dept. Applied Physics, Science for Life Laboratory, Solna, Sweden
    Contribution
    Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing - review and editing
    Competing interests
    No competing interests declared
  5. Ilaria Testa

    KTH Royal Institute of Technology, Dept. Applied Physics, Science for Life Laboratory, Solna, Sweden
    Contribution
    Conceptualization, Funding acquisition, Investigation, Supervision, Validation, Writing - review and editing
    Competing interests
    Reviewing editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4005-4997
  6. Lucie Delemotte

    KTH Royal Institute of Technology, Dept. Applied Physics, Science for Life Laboratory, Solna, Sweden
    Contribution
    Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Validation, Writing - review and editing
    For correspondence
    lucie.delemotte@scilifelab.se
    Competing interests
    Reviewing editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0828-3899

Funding

Science for Life Laboratory

  • Ilaria Testa
  • Lucie Delemotte

Vetenskapsrådet (2018-04905)

  • Lucie Delemotte

Gustafssons Stiftelse

  • Lucie Delemotte

European Commission (H2020-MSCA-IF 893077)

  • Lea Rems

Slovenian Research Agency (J2-2503)

  • Lea Rems

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by grants from the Science for Life Laboratory, from the Gustaffsson foundation from Swedish Research Council (VR 2018-04905) to L.D., and by a synergy postdoc grant to L.D. and I.T. from KTH Royal Institute of Technology. The work was also supported by funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 893077 to L.R., and by funding from Slovenian Research Agency (ARRS) under project no. J2-2503. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing at parallelldatorcentrum (PDC) Centre for High Performance Computing and at High Performance Computing Center North (HPC2N).

Senior Editor

  1. José D Faraldo-Gómez, National Heart, Lung and Blood Institute, National Institutes of Health, United States

Reviewing Editor

  1. Qiang Cui, Boston University, United States

Reviewers

  1. Helgi I Ingolfsson
  2. Rumiana Dimova, Max Planck Institute of Colloids and Interfaces, Germany

Publication history

  1. Preprint posted: October 16, 2021 (view preprint)
  2. Received: October 16, 2021
  3. Accepted: February 22, 2022
  4. Accepted Manuscript published: February 23, 2022 (version 1)
  5. Version of Record published: March 10, 2022 (version 2)

Copyright

© 2022, Rems et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 672
    Page views
  • 128
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Lea Rems
  2. Xinru Tang
  3. Fangwei Zhao
  4. Sergio Pérez-Conesa
  5. Ilaria Testa
  6. Lucie Delemotte
(2022)
Identification of electroporation sites in the complex lipid organization of the plasma membrane
eLife 11:e74773.
https://doi.org/10.7554/eLife.74773

Further reading

    1. Structural Biology and Molecular Biophysics
    Maicon Landim-Vieira et al.
    Research Article Updated

    Phosphorylation and acetylation of sarcomeric proteins are important for fine-tuning myocardial contractility. Here, we used bottom-up proteomics and label-free quantification to identify novel post-translational modifications (PTMs) on β-myosin heavy chain (β-MHC) in normal and failing human heart tissues. We report six acetylated lysines and two phosphorylated residues: K34-Ac, K58-Ac, S210-P, K213-Ac, T215-P, K429-Ac, K951-Ac, and K1195-Ac. K951-Ac was significantly reduced in both ischemic and nonischemic failing hearts compared to nondiseased hearts. Molecular dynamics (MD) simulations show that K951-Ac may impact stability of thick filament tail interactions and ultimately myosin head positioning. K58-Ac altered the solvent-exposed SH3 domain surface – known for protein–protein interactions – but did not appreciably change motor domain conformation or dynamics under conditions studied. Together, K213-Ac/T215-P altered loop 1’s structure and dynamics – known to regulate ADP-release, ATPase activity, and sliding velocity. Our study suggests that β-MHC acetylation levels may be influenced more by the PTM location than the type of heart disease since less protected acetylation sites are reduced in both heart failure groups. Additionally, these PTMs have potential to modulate interactions between β-MHC and other regulatory sarcomeric proteins, ADP-release rate of myosin, flexibility of the S2 region, and cardiac myofilament contractility in normal and failing hearts.

    1. Biochemistry and Chemical Biology
    2. Structural Biology and Molecular Biophysics
    Atefeh Rafiei et al.
    Research Article Updated

    Doublecortin (DCX) is a microtubule (MT)-associated protein that regulates MT structure and function during neuronal development and mutations in DCX lead to a spectrum of neurological disorders. The structural properties of MT-bound DCX that explain these disorders are incompletely determined. Here, we describe the molecular architecture of the DCX–MT complex through an integrative modeling approach that combines data from X-ray crystallography, cryo-electron microscopy, and a high-fidelity chemical crosslinking method. We demonstrate that DCX interacts with MTs through its N-terminal domain and induces a lattice-dependent self-association involving the C-terminal structured domain and its disordered tail, in a conformation that favors an open, domain-swapped state. The networked state can accommodate multiple different attachment points on the MT lattice, all of which orient the C-terminal tails away from the lattice. As numerous disease mutations cluster in the C-terminus, and regulatory phosphorylations cluster in its tail, our study shows that lattice-driven self-assembly is an important property of DCX.