Modification of lipid membrane compressibility induced by an electric field

K. R. Prathyusha ,1,*,† Ignacio Pagonabarraga ,2,3,4 and P. B. Sunil Kumar1,5,‡ 1Department of Physics, Indian Institute of Technology Madras, Chennai, India 2CECAM, Centre Européen de Calcul Atomique et Moléculaire, École Polytechnique Fédérale de Laussane (EPFL), Batochime, Avenue Forel 2, 1015 Lausanne, Switzerland 3Departament de Física de la Matèria Condensada, Universitat de Barcelona, C. Martí Franquès 1, 08028 8 Barcelona, Spain 4University of Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, 08028 10 Barcelona, Spain 5Indian Institute of Technology Palakkad, Ahalia Integrated Campus, Kozhippara, Palakkad 678557, Kerala, India


I. INTRODUCTION
Lipid molecules spontaneously self-assemble in an aqueous environment to form a bilayer membrane [1,2]. As a result of the charges on their polar head groups and the low permeability of the hydrophobic tail to the solvent, lipid membranes are strongly influenced by the action of the external electric field. Two widely used techniques involving electric fields are electroporation and electroformation [3,4]. In electroporation, the large applied electric field induces pores in the bilayer which can be used to inject macromolecules in vesicles or cells and has enormous applications in biology, biotechnology, and medicine [4][5][6][7]. In electroformation, the application of a low-frequency AC electric field induces peeling of membranes from multilamellar stacks on a conducting substrate, leading to the formation of giant unilamellar vesicles [8].
Recently, electroformation technique is employed to construct lipid tubes [9] when the applied electric field is parallel to the bilayer interface and to produce double vesicles [10] when a combination of sinusoidal and amplitude-modulated electric field is applied. As both techniques are linked to various applications in cell biology, biotechnology, and pharmacology, it is crucial to understand the response of bilayer membranes to an external electric field.
In the past few decades, many aspects of the response of the bilayer to electric fields have been investigated us-ing experimental, theoretical, and computational tools. On the experimental side, extensive efforts have been made to understand the influence of the electric field on the bilayer. It has been observed that the electric field affects the rate of swelling of the bilayer significantly. [3,11,12]. Studies on giant vesicles have demonstrated that membrane deformation and poration depend crucially on the electric pulse duration [13,14]. X-ray scattering measurements on a supporting bilayer have revealed that an alternating electric field could increase the bending stiffness due to interaction between charges inside the electric double layer and can decrease the surface tension, caused by the amplification of charge fluctuations [15].
The influence of an electric field on the biological membrane has been studied theoretically [16][17][18][19][20][21][22][23][24][25]. Instability of the membrane can occur due to thickness fluctuations and bending modes [6,17]. Sens and Isambert [17] demonstrated that negative surface tension and ionic current near the membrane interface induce the undulation of the membrane. The applied electric field increases the membrane bending rigidity and decreases the membrane tension, at sufficiently high potential, stretching instability induced by negative tension results in the formation of pores [19]. A model that included the explicit coupling between the orientation of the dipolar lipid head groups and the membrane shape explained that undulation instability arises through the buckling modes involving the thickness fluctuation, leading to the localized membrane breakdown and pore formation [18].
Molecular dynamics simulations have shown that the applied electric field enhances the membrane permeability due to pore formation which assists the transport of ions or small molecules [26][27][28][29][30][31][32]. Coarse-grained simulations have verified that an applied electric field would give rise to additional stress in the membrane [33] or a spontaneous curvature of biologically relevant size for the case of a floating bilayer [34].
In this paper, we investigate the influence of the electric field on the mechanical properties of a membrane, in the absence of explicit ions, using molecular dynamics simulations. The length scale of interest here, the size of the vesicles formed by electroformation, is of the order of micrometers and the time scale is in the range of a few microseconds to milliseconds [13]. Since an all-atom simulation is inadequate to cover this range of length and time scales, the coarse-grained MARTINI force field [35] is used to model the membrane lipids and the solvent. For simplicity, we assume the solvent to be one with constant dielectric permittivity and consider lipids with fixed dipoles. In our model, we ignore the effect of mobile charges and any hydrodynamic flows. These are nontrivial extensions and will be studied in the future.
We have simulated membrane patches up to a size of 0.4 μm for approximately 4 μs to investigate the influence of the electric field on the membrane mean area, the compressibility, and the fluctuation spectrum. In Sec. II, we describe the computational method in detail which is employed to simulate an equilibrium tensionless membrane. Section III contains our main results and identifies the critical electric field beyond which membrane undulations of the scale of the patch set in. We finish with a discussion of the main results obtained.

A. Model
We use the MARTINI force field (MARTINI 2.0) [36], a coarse-grained force field for biomolecular simulations, to study the membrane dynamics. In this model, four atoms of a given species are equivalent to one representative element. For example, four water molecules are coarsegrained to form a single water bead (P 4 ). We consider single component membranes composed of zwitterionic dipalmitoylphosphatidylcholine (DPPC), which is modeled using 12 coarse-grained (CG) beads [35]. In the coarse-grained representation of DPPC, the zwitterionic phosphatidylcholine (PC) head group is modeled by positively charged particles (choline group) labeled as Q 0 and a negatively charged particle (phosphate group) labeled as Q a . The choline group is deprived of hydrogen bonding capability, and the phosphate group can act as a hydrogen bonding acceptor. The glycerol-ester group in DPPC is modeled as two apolar (N a ) particles. Each tail has four hydrophobic particles modeled as nonpolar (C 1 ) beads, representing 16 methyl and methylene groups.
The MARTINI force field considers that all nonbonded particles interact through a shifted Lennard-Jones (LJ) potential energy function given by where σ i j characterizes the size of the particles, r is the separation between them, and i j is the strength of their interaction. The σ i j is set as 0.67 nm for Q-and C 1 -type particles, and for all other pairs, σ i j is set to 0.47 nm. The cutoff distance r cut is set to 1.2 nm. A standard Gromacs shift function [37] is used to shift the LJ interaction smoothly between r shift = 0.9 TABLE I. The MARTINI interaction matrix for CG water and DPPC system used in the paper. Different levels of interaction: O, = 5.6 kJ/mol; I, = 5.0 kJ/mol; II, = 4.5 kJ/mol; III, = 4.0 kJ/mol; IV, = 3.5 kJ/mol; VI, = 2.7 kJ/mol; VIII, = 2.0 kJ/mol; and IX, = 2.0 kJ/mol, σ = 0.67 nm. II  I  III  IX  N a  III  III  III  III  VI  C 1  VIII  IX  IX  VI  IV nm and r cut . The interaction parameters used in this model are given in Table I. The complete interaction matrix for this model follows the one provided in Ref. [36].
The electrostatic interaction between all charged groups is through a truncated and shifted Coulomb potential given by where r = 15 is used for explicit screening. This potential is also shifted smoothly from r shift = 0.0 nm to r cut = 1.2 nm using the same Gromacs shift function [37]. The bonded particles interact only via a harmonic potential, and a three-body interaction for bending stiffness, The equilibrium distance r bond = 0.47 nm and the force constant K bond = 1250 kJ mol −1 nm −2 . The bending stiffness parameter K angle is set to 25 kJ mol −1 and the equilibrium bond angle θ 0 is set to 180 • .

B. Simulation details
Preassembled bilayers are used as initial configurations. Bilayers are constructed using 6400 and 1600 coarse-grained DPPC molecules. The bilayers are sandwiched by a solution containing 269 660 CG water beads with 6400 lipids and 67 415 CG water beads with 1600 lipids. Unless otherwise specified, all results presented in this paper correspond to simulations with 6400 lipid molecules. The equilibrium simulations run for about 160 ns.
Simulations are performed at a constant temperature of T = 323 K and a constant pressure of P = 1 bar, using the Nose-Hoover method [38] with lateral and normal pressure independently coupled to the barostat to obtain an equilibrium tensionless membrane. Periodic boundary conditions are employed in all three directions. An integration scheme, developed by Tuckerman et al. [39], is adapted for the time integration of non-Hamiltonian equations of motion with an integration time step of t = 40 fs. The molecular dynamics package LAMMPS [40] is used for all simulations.

C. Membrane in the presence of the electric field
To simulate the effect of an external electric field of magnitude E , applied along e z which is chosen as the direction perpendicular to the initial bilayer surface, a force F i = q i E is applied on particle i bearing charge q i [41,42]. We explore a wide range of electric field magnitudes, from 0.01 to 0.4 V/nm, to study the various deformation regimes of the membrane. The magnitude of the electric field, required to deform the membrane appreciably, is much larger than the magnitudes used in electroformation experiments. The primary reason for this difference could be the inability of the model used in the simulation to reproduce the dielectric properties of water quantitatively. However, the strength of the applied electric fields is comparable to the magnitude required to form pores in previous simulation studies [43,44].

D. Calculation of pressure profile
For planar membranes, the common practice to obtain pressure profile across the membrane is to divide the simulation box in small slabs along the z direction and measure the pressure in each slab, assuming that the normal of the membrane corresponds to the initial orientation of the bilayer (in our studies along the z axis) [45]. The method breaks down at high electric fields due to large membrane deformations.
To determine the local normals to the membrane surface, we marked the membrane-water interfaces using the head group positions of the lipids and then employed the Delaunay triangulation technique [46] to recreate the interface. The surface normal is then determined at each triangle. For a good balance between resolution and accuracy we averaged the position of three head groups to mark a vertex, and each bilayer is treated separately. The pressure (P) profile as a function of the local normal is now obtained by averaging over a small volume around the center of each triangle. We compute the pressure profiles using the normals of the top layer and the bottom layer independently. Once we have computed the pressure tensor at each position along the membrane, we must transform it to local coordinates, P = H T PH, by using the House-holder matrix H = I − 2ww T , where w =n −ẑ |n−ẑ| andn is the unit normal of the triangle [47]. P zz and a 2 × 2 minor of P give the stress along the normal and the tangential plane, respectively.

III. RESULTS AND DISCUSSION
The zero-field membrane thickness is measured as the distance between the average head-group position in the upper and the lower leaflet of the membrane. This value is found to be 4.01 ± 0.04 nm [35], close to the experimentally determined thickness of 3.85 nm for the lamellar liquid crystalline phase of DPPC [48].
We find that when subject to an electric field, the response of the membrane can be linear, nonlinear, or unstable, depending on the strength of the electric field. As seen in Figs. 1(a) and 1(b) the membrane deformations are negligible for low fields, E < 0.04 V/nm. At intermediate electric field strengths, 0.04 V/nm E < 0.3 V/nm, large amplitude undulations of the membrane develop, as shown in Figs. 1(c) and 1(d). Finally, for stronger fields, E 0.3 V/nm the membrane becomes unstable and the water-membrane interface re-orients parallel to the electric field. We see that the membrane undulations with wavelengths comparable to the system size dominate. We cannot rule out the existence of a characteristic wavelength larger than the system sizes we have been able to reach. Sens and Isambert [17] showed that, for a freely floating membrane having a bending modulus of κ = 5 × 10 −20 J and a thickness of 5 nm, in a solvent of viscosity η = 10 −3 Pa s, the wavelength of the fast-growing undulating mode is ∼0.5 μm, larger than the system sizes of ∼43 nm in our present study. However, whether this is the relevant length scale for uncharged membranes remains an open question.
After this qualitative description of the impact that an applied electric field has on the membrane morphology, we analyze now in more detail the change in mechanical properties of the membrane.

A. Structure factor
Membrane fluctuations can be used to extract information about the membrane mechanical properties, in particular its rigidity and tension. From the local height of the membrane in the plane corresponding to the initial orientation of the membrane, in our case h(x, y), we can compute the circularly averaged structure factor, S(q) = |h(q)| 2 , where q = √ q 2 x +q 2 y , and q x and q y are Fourier transform variables of x and y, respectively. This structure factor quantifies the membrane undulations, and according to the Helfrich Hamiltonian, 2 ], and applying the equipartition theorem, the structure factor reads where L x and L y are the box size in the x and y directions, respectively. This expression relates the membrane height fluctuations to the mechanical properties of the membrane through γ , the membrane surface tension, and κ, its bending modulus. Figure 2 displays k B T /[Aq 2 S(q)] as a function of q 2 for different values of the applied electric field, where A = L x L y . The intercept and the slope of this curve correspond to γ Sq and κ, respectively. The main panel of Fig. 2 shows that the intercept is close to zero and is insensitive to the electric field, indicating that the surface tension is always negligible even when undulations are large. It is clear from the inset of Fig. 2 that, though there are large undulations in the membrane, the bending modulus is unaffected by the electric field and stays close to the experimental value for DPPC [49].

B. Lateral pressure profile across the membrane
Since the surface tension measurement is noisy, to verify whether the membrane remains tensionless at large electric fields of E > 0.04, we use the local anisotropy in the pressure profile along and perpendicular to the bilayer normal.   Using the Irving-Kirkwood formula [50], for a planar lipid bilayer in the x-y plane with its surface normal oriented along the z axis, the components of the pressure profile, parallel and perpendicular to the membrane, can be expressed as P L (z) = 1 2 [P xx (z) + P yy (z)] and P N (z) = P zz (z) (For membranes undergoing large deformations, these quantities are calculated using the local normal to the membrane). The surface tension is obtained from the average pressure difference between both components, γ p f = dzp z , where p z (z) = P L (z) − P N (z). Figure 3 shows the anisotropic pressure profile across the bilayer. The middle peak in the anisotropic excess pressure profile corresponds to the repulsion between the hydrocarbon chains, while the side peaks account for the head-head repulsion and the troughs are due to the hydrophobic effect [2]. In the zero-field pressure profile, the heights of the interfacial and head group peaks are comparable with all-atom simulations [36,51]. From Fig. 3, it is clear that the overall pressure profile and the peak positions are insensitive to the applied electric field. As the electric field increases, there is a systematic decrease in the peak and troughs, due to the reduction in intermolecular interactions.
This decline suggests the average separation between lipids is increasing, leading to the increment in the total area of the membrane (see Sec. III D). Above the field required to set the undulation amplitude to be equal to the bilayer thickness, E c = 0.04 V/nm, the depth of the troughs decreases (see inset of Fig. 3). This increase in membrane area with the increase in the electric field strength causes the "membrane stretching." Figure 4 displays the surface tension extracted from the lateral pressure profile and the structure factor. It is evident from the figure that the two methods give consistent results and there is no systematic change in the value of surface tension with the electric field. Instead, we see a significant membrane lateral expansion in response to the electric fields, as discussed in detail in Sec. III D. This is in contrast to the results obtained from simulations with a fixed projected area where electric-field-induced surface tension is observed [33].

C. Membrane polarizability
In equilibrium, the balance between the lipid dipolar interaction in the two membrane leaflets, thermal fluctuations, and molecular shape determine the dipolar orientation of the lipids. As a result, lipids exhibit a preferential orientation, characterized by the equilibrium angle θ eq that they make with respect to the outward membrane normal of two leaflets. In our system, we observe, θ eq 60 • for both upper and lower leaflets, comparable to the value for a membrane made of palmitoyl-oleoyl-phosphatidylcholine (POPC) lipids [28].
The electric field aligns the molecular dipoles along with it. We observe that, as we increase the strength of the electric field, θ eq decreases for the lipids in the upper leaflet and it increases for the dipoles in the bottom leaflet. Such a shift in θ eq (θ eq 38 • for the upper leaflet and θ eq 83 • for the lower leaflet) with the electric field has been previously reported for all-atom simulations [28,44]. This change in the dipolar angular distribution implies that the net leaflet dipole moment will vary asymmetrically in the two leaflets that in turn affect the membrane polarizability. We can quantify the membrane polarizability by measuring the dipolar membrane moment, where p i = qd i and d is the separation vector between positive and negative charges in each lipid molecule and N lipids stands for the total number of lipids. We now analyze P mem,z , the component of the polarization normal to the membrane. Figure 5(a) displays the total membrane dipole moment along with the total dipole moment of the two leaflets, quantifying the asymmetric response to the applied electric field. The comparison between the three dipole moments also indicates that the bottom leaflet is more sensitive to the applied field and is the main contributor to the overall membrane dipole moment. Additional insight into how each leaflet contributes to the overall membrane dipole moment is provided by subtracting its zero-field value. Figure 5(b) demonstrates that there is a significant contribution from the dipoles in the lower leaflet. From the figure, it is clear that, at the low electric field, the response in both leaflets is essentially the same. As the electric field is increased, the bottom layer deviates more strongly from its zero-field value. This asymmetry in the dipole moment between the two leaflets indicates that the membrane polarizability will have a strong nonlinear dependence and could be one of the reasons for the instability observed in our simulations. Our results suggest that, when the dipoles tend to align with the field, the lipid molecules change their orientation locally, to minimize the energy, promoting membrane undulations. A similar trend in the dipole moment with electric fields is presented in a previous work; however, one cannot exclude the additional contributions coming from the presence of mobile ions in their simulation. An earlier the-  6. Comparison of the averaged projected area A p (scaled by A p0 ) and the averaged total surface area A s (scaled by A s0 ) as a function of E for two difference system sizes. The projected area decreases while the membrane surface area increases with the electric field. The membrane area is measured in units of nm 2 .
oretical study [18] shows that coupling between dipole and membrane orientation destabilizes a planar membrane; a more thorough comparison to identify the microscopic origin of these phenomenological couplings will be useful to clarify the role that the asymmetric dipolar response to the applied electric field plays in membrane rupture.

D. Area of the membrane and compressibility modulus
The membrane "stretching" as a function of the electric field strength can be substantiated by measuring the in-plane membrane surface area A s . We have computed A s from the area of the triangles obtained using the Delaunay triangulation technique described in Sec. II D. We have computed the membrane surface area from the upper and lower leaflets and have found the difference to be negligible in all regimes. Figure 6 depicts the average membrane surface area (A s ) and membrane projected area (A p ), scaled by their zero-field values (A s0 and A p0 ), as a function of the electric field for both 1600 and 6400 lipids. Both A p and A s show a plateau below the critical electric field E c . The observed values of the critical electric field E c , which characterize the low-electricfield regime plateau, are 0.04 and 0.1 V/nm for 6400 and 1600 lipids, respectively. Above E c static undulations become prominent and A p decreases while, correspondingly, A s increases.
The changes in the total and projected membrane surface areas impact the membrane compressibility. The area compressibility moduli for both quantities are defined as where i should be replaced by p or s, when the compressibility modulus is measured using projected area A p or the interfacial surface area A s , respectively. The compressibility moduli for both the projected area A p and the total surface area A s are shown in Fig. 7(a). In the low-electric-field regime, for E = 0.01 V/nm the value of both compressibility moduli are found to be ≈240 mN/m, consistent with the compressibility modulus, K A = 260 ± 40 mN/m, computed for a membrane patch containing 6400 lipids [35] (both K Ap and K As ), in the absence of an external field. The values and the system size dependence of area compressibility moduli [see Fig. 7(a)] at the low electric field are consistent with previous studies [52][53][54].
When the electric field is increased up to E c , both compressibility moduli decrease due to the induced membrane undulatory modes. The deviations between K Ap and K As are small when E < E c . For E > E c , when the static undulations set in, K Ap increases while K As decreases. Further increasing the electric field does not change the projected area, and the corresponding compressibility modulus shows a monotonic increase. However, the interfacial area continues to increase in this regime, leading to a saturation of the corresponding compressibility. Figure 7(a) shows that while the system size is not relevant for strong electric fields, it affects the onset of membrane undulation. Larger system size simulations will be required to clarify the implications of the reduction of the critical electric field with system size; it remains unclear if there will exist a minimum finite E c or whether the nature of the unstable modes become singular for sufficiently large system sizes. In order to see if there is a generic trend, we have scaled the membrane compressibility by their zero-field value and the applied electric field by the critical value E c when nonlinear undulations develop. The resulting graph, shown in Fig. 7(b), identifies the weak and strong electric field regimes, although we do not see any scaling regime for the different membrane compressibilities.

E. Instability of the membrane at higher electric fields
For strong enough electric fields, E 0.3 V/nm, the continuous increase in the undulation amplitude and the decrease in membrane density lead to pore formation; pores appear at the highly curved regions of the membrane. Figure 8 shows the time evolution of a pore on the membrane. One can see that once a pore formation is initiated, the diameter of the pore increases until it becomes comparable to the system size. This instability reorients the water-bilayer interface parallel to the z axis. In this new configuration, dipoles orient perpendicularly to the electric field. In an earlier simulation with only 126 lipids, complete rotation of the membrane was seen at high hydration [55]. Our simulations show that for the larger system size membrane undulations precede poration and eventually reorientation of the interface. The driving mechanism for this instability is related to the well-known fact that, when an electric field is applied to a system containing two dielectric media, the low-energy configuration corresponds to the one where the interface aligns parallel to the applied field [56]. Our simulations represent an infinite system, in the plane perpendicular to the electric field, through periodic boundary conditions. Since such an infinite system cannot reorient, an alignment of the interface parallel to the field is possible only through a breakup of the membrane. In Figs. 8(g) and 8(h), it is clear that the interface is parallel to the external electric field. However, we emphasize that, except for limiting the largest length scale possible for the undulatory mode, usage of the periodic boundary condition does not affect our results. Therefore, the nature of the instability we have observed is qualitatively different from previous instabilities described for charged membranes [17,57] and is controlled exclusively by the dielectric properties of the membrane and the medium it is embedded in.

F. Summary
We have carried out molecular dynamics simulations of a lipid bilayer membrane in the presence of an external electric field in the isothermal-isobaric (NPT) ensemble. We have observed the onset of static undulations above a critical field, with the wavelength of undulations set by the system size. Above the critical field, we observe an increase in the inplane area of the membrane, while the projected area remains constant. The undulation amplitude increases with the electric field, and beyond a system-dependent threshold, E = 0.32 V/nm for 6400 lipids, pores form at the crests or troughs of the undulating membrane, leading to membrane breakup characterized by the tendency of the water-bilayer interface to align parallel to the electric field due to the dielectric mismatch between the membrane and the solvent medium. We have computed the pressure profiles, the height fluctuation spectrum, the membrane dipole moment, and the area compression moduli to understand the nature of the morphological changes of the membrane. It is observed that the membrane mechanical parameters such as surface tension and bending rigidity are unaffected by the electric field. There is a systematic decrease in the lateral pressure profile at the critical electric field for the onset of static deformations. These results, along with the fact that the membrane thickness is unchanged when forced by the electric field, suggest that the membrane responds by increasing the interlipid distance leading to membrane stretching without modifying its surface tension. The reason for this increase has not been completely identified, although the membrane asymmetric dipolar reorientation under the action of the electric field favors membrane distortion and promotes membrane undulations. The change in compressibility observed in the simulations, and its understanding in terms of membrane deformation and stretching, has not been reported earlier.
We point out that, in our current study, we excluded the contributions arising from the motion of explicit ions or fluid flows. To fully model the experimental situation, one should take into account the accumulation of charges near the surface of the membrane. However, our study helps to segregate and understand the effects of the electric field on the membrane due to the coupling between the electric field and the lipid dipole. We claim that interaction with dipoles in the head group and the external electric field is enough to produce curvature and undulation in the membrane. Undulations at the high electric field can lead to pore formation and instability. We believe that our findings can benefit new electroformation protocols when sucrose and deionized water is used with stainless steel electrodes and designing biophysical applications of the membrane when one has to screen bioactive agents or when having salt is a disadvantage [58]. The effect of ionic impurities in the system could be crucial, as they contribute to electrokinetic effects, and future studies will incorporate them. In order to have a meaningful comparison with the experimental results, additional simulations are required by increasing the system size and modifying the parameters to understand whether any characteristic length is associated with the membrane instability and pore formation. The model introduced can be extended further to explore the effect that the architecture of lipids and cholesterol have on the instability, in the presence of the electric field.
The MARTINI model used here has several advantages and shortcomings when employed for extensive area simulations [59]. As the model uses short-range forces, and it has good scalability with computer clusters, one can employ parallel computational approaches and simulate large systems. A ma-jor computational bottleneck is the large number of water particles needed to stabilize a large membrane patch despite the four to one mapping. Electrostatic interaction is expressed using a shifting potential, and the electrostatics is screened with a relative dielectric constant, which can lead to certain artifacts. To incorporate the correct dielectric properties of water, one can choose a polarizable MARTINI water model [60]. In addition to the implicit screening in the model, the neglect of long-range electrostatic interaction is another drawback. More detailed limitations of the MARTINI model can be found in Ref. [59]. In a real system, the electric field can enhance the generation of reactive oxygen species and alter the composition of lipids leading to lipid peroxidation. Various studies show that unsaturated lipids are easily susceptible to peroxidation, while saturated lipids do not undergo any oxidation [61,62]. In our simulations, we consider that the bilayer is only made of saturated lipid DPPC; therefore, we ignore the effect of peroxidation. To capture the lipid peroxidation process, one certainly need to do an atomistic scale simulation. While the reactive contribution must be taken in a biological membrane simulation, wherein both saturated and unsaturated lipids are present, our simplified model illustrates that, even in the absence of these processes, pores can be produced. This insight is relevant to understand the different physical mechanisms that lead to pore formation.