Dissipative particle dynamics simulation of critical pore size in a lipid bilayer membrane

We investigate with computer simulations the critical radius of pores in a lipid bilayer membrane. Ilton et al. (Ilton et al. 2016 Phys. Rev. Lett. 117, 257801 (doi:10.1103/PhysRevLett.117.257801)) recently showed that nucleated pores in a homopolymer film can increase or decrease in size, depending on whether they are larger or smaller than a critical size which scales linearly with film thickness. Using dissipative particle dynamics, a particle-based simulation method, we investigate the same scenario for a lipid bilayer membrane whose structure is determined by lipid–water interactions. We simulate a perforated membrane in which holes larger than a critical radius grow, while holes smaller than the critical radius close, as in the experiment of Ilton et al. (Ilton et al. 2016 Phys. Rev. Lett. 117, 257801 (doi:10.1103/PhysRevLett.117.257801)). By altering key system parameters such as the number of particles per lipid and the periodicity, we also describe scenarios in which pores of any initial size can seal or even remain stable, showing a fundamental difference in the behaviour of lipid membranes from polymer films.


Introduction
A recent article by Ilton et al. [1] examined the evolution of pores in a model membrane constructed from polystyrene in a water bath. A tightly focused laser was used to create a temperature gradient in the film, decreasing the local surface tension and driving the formation of a hole whose size was determined by the power and exposure time of the laser. It was found that holes below a critical size r c sealed, while larger holes began to increase in size, a behaviour which had previously been observed in solid-state membranes [2].
The liquid polymer films of [1] are easily observable via traditional microscopy techniques due to their large size (the thinnest homopolymer film studied was approximately 100 nm thick, with most films of the order of 1 mm). By contrast, a typical lipid membrane, as might be found in a biological cell, is much smaller (approx. 10 nm thick) [3,4]. Pore formation and evolution in a lipid membrane thus have a number of experimental difficulties owing to the small spatial and temporal scales [5]. Previous work has modelled membrane pore dynamics with Monte Carlo mesh simulation [6], continuum elasticity and energy arguments [7][8][9] and small-scale particle simulations [7,10,11]; in this paper, we simulate a porated lipid membrane with a coarse-grained method which allows feasible computation at large scales. The simulated system will be shown to reproduce many of the results observed in [1].
Simulating dynamics on the scale of individual atoms and molecules (such as lipids) is often done with molecular dynamics, a class of commonly used simulation methods which operate by numerically solving Newton's equations of motion for each particle [12]. Because particles are simulated individually, the computational cost of such simulations becomes infeasible at even moderate scales. This is particularly evident when simulating particles or structures immersed in fluids (e.g. membranes), whose characteristic time of motion often differs significantly from that of the solvent.
To explicitly simulate a lipid bilayer membrane, we employ dissipative particle dynamics (DPD), a modification of molecular dynamics which reduces computational complexity by aggregating small groups of like atoms or molecules into single 'dissipative particles'. Extensive comparisons with molecular dynamics and Navier-Stokes simulations have shown that, with the proper choice of intermolecular forces, DPD simulations maintain the correct hydrodynamic behaviour across a wide range of spatial and temporal scales [13][14][15]. The result is a method which explicitly models the solvent (via coarse-grained dissipative particles) while being computationally feasible at the scale of large biological structures [12,16].
We begin by describing the explicit formulation of DPD in §2. The construction of the membrane and measurement techniques are detailed in §3. Sections 4 and 5 present the simulation results. Finally, the results are contextualized and discussed in §6.

The DPD simulation method
Let the mass, velocity and position of dissipative particle i be given by m i , v i and r i , respectively. The DPD equation of motion for particle i comprises three pairwise contributions: F C ij is a conservative force deriving from a potential exerted on particle i by particle j, similar to the usual pairwise forces implemented in molecular dynamics schemes. Here, we adopt the common explicit form ( in terms of a conservative coefficient a ij , the inter-particle displacement r ij ¼ r i 2 r j (with magnitude r ij and unit vectorr ij ) and a cutoff radius r c . The DPD conservative force is a soft potential (i.e. does not diverge as r ij ! 0), and so particles can overlap or even occupy the same point in space, corresponding to the notion of DPD particles as coarse-grained clusters of smaller atoms or molecules [14,17]. The dissipative force F ij D and random force F ij R function as a thermostat and are given by where v D ( Á ) and v R ( Á ) are position-dependent weight functions, v ij ¼ v i 2 v j , and g ij and s ij are the dissipative and random strengths, respectively. Noise is introduced by the Gaussian white-noise term u ij , which satisfies the stochastic conditions In order to ensure conservation of momentum, it is assumed that the noise terms are symmetric in i and j, i.e. u ij ¼ u ji . Españ ol & Warren [13] provided an additional condition in order to preserve the invariant distribution of the system with conservative forces alone, namely, the fluctuation-dissipation relation 2 and where k B is the Boltzmann constant and T the equilibrium temperature. In total, the DPD equations of motion can then be written as the following set of coupled stochastic differential equations: To proceed with the DPD simulation, the system is then stepped forward with a numerical solver. Here, we use the DPD velocity-Verlet scheme: given positions r n and velocities v n at step n and a timestep Dt, compute the half-step velocities then calculate the next step as This scheme is very similar to traditional velocity-Verlet, with the exception that the dissipative force term F D must be calculated twice per step since it is both position-and velocity-dependent [12,18].

Lipid bilayer simulation
Before carrying out the simulations, it remains to specify the coefficients of equation (2.1). Because there is a choice of scale for the coarse-graining, DPD coefficients are usually specified for the nondimensionalized system. In their foundational paper on DPD, Groot & Warren [19] found that the dimensionless compressibility of water could be matched by using the conservative coefficient a ij ¼ 25.0, dissipative coefficient g ij ¼ 4.5, cutoff radius r c ¼ 1.0 and numerical density r ¼ 3 for unitmass DPD particles. To model a lipid bilayer membrane, we additionally introduce particles with modified coefficients to model the lipids, here referred to as 'head' and 'tail' particles. The pairwise a ij and g ij , chosen to be similar to the existing literature on DPD membranes [20 -22] and to reproduce mesoscopic properties such as lateral fluidity and a stable bilayer structure [23], are shown in table 1.
Head and tail particles are connected as in figure 1c by harmonic bonds with dimensionless energy E ij ¼ 64(r ij 2 0.5) 2 , so that the resting length of a bond is 0.5. Along the tails, three-body potentials E ijk ¼ 20(1 þ cosf ijk ) are introduced to provide stiffness, where f ijk is the angle between bond ij and bond jk; the resulting hydrophobic chains are thus governed by a bending energy with zero curvature. Copies of the lipids are then placed in two layers at the vertices of a square lattice, with the tails facing inward.  The outside of the membrane is initialized with water particles placed on a cubic lattice with numerical density r ¼ 3. Periodic boundary conditions are imposed around the simulation box. For the first simulation, the periodic simulation box of 144 Â 144 Â 40 was filled entirely in the first two dimensions with a bilayer membrane comprising 28 700 lipids with three head particles and two tails of six particles each, so that the side length of the square lattice on which lipids were initialized was approximately 1.202. The number of lipids was chosen experimentally to achieve stability in the membrane. To simulate abrupt pore formation by laser ablation, all lipids intersecting an orthogonal cylinder of the fixed radius were removed, and the cylinder was included in the initialization region for fluid particles. The resulting initial state can be seen in figure 1b.
In each simulation, the time evolution of the pore was observed using DPD velocity-Verlet with a temperature T ¼ 1.0 and timestep Dt ¼ 0.005. To track the size of the pore, the locations of all lipids were recorded every 1000 steps. Images of the system were created by rendering spheres of radius 0.5 at the location of every head and tail particle, then projecting the three-dimensional system into a two-dimensional plane tangent to the membrane surface. The resulting images were then smoothed and thresholded, resulting in binary matrices M i with value 1 within the pore and 0 without, as in figure 1d. Finally, the two-dimensional centre of mass position x i was computed, and the sum P x M i (x)jx À x i j was calculated. For a circle of radius R, the result should be approximated by the integral and so the approximate radius of the pore is given by : To compute the membrane width, the three-dimensional representation was instead projected into a plane bisecting the membrane (i.e. an orthographic side view) as in the side view of figure 1b. For each row in the projected image, the number of non-background pixels was summed, approximating the width of the membrane at that location; to reduce the effect of the natural thermal fluctuations in the membrane surface, the minimum of all row widths was used as the membrane width for that image. royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181657 The resulting value was time-averaged over 20 images (one every 1000 steps) to obtain a reference value h for the membrane width.

Existence of a critical radius
Results for pores of various initial radii can be seen in figure 2. The behaviour of the pore is a function of its radius: sufficiently small pores begin to shrink and eventually seal entirely, while sufficiently large pores begin expanding, eventually severing the membrane. The time-evolving radius of each pore, as measured by equation (3.1), is shown in figure 3. The growth of supercritical pores proceeds in a roughly exponential fashion, in agreement with the findings for homopolymer films in [1]. For very large times (t . 1000), the pore approaches the order of the simulation box, and so begins interacting with itself across the periodic boundary, resulting in slower growth. The subcritical pore, which closed around t ¼ 400, yielded a stable membrane for the remainder of the simulation. The observed critical ratio 1.34 , r c /h , 1.47 is larger than the ratio r c /h ¼ p/4 derived by Ilton et al. [1] for a homopolymer film, suggesting an additional free energy cost for pore formation due to the lipid structure. Note that we define the ratio in terms of the initial size r 0 , meaning it is possible a pore which evolves to be smaller than the critical radius may still expand as t ! 1.
A second set of simulations examined the effect of increasing internal pressure in the membrane. The number of lipids in the periodic simulation box was increased by 4:4% to 29 970, resulting in a stable membrane with a lipid excess. For short times, a pore opened in the modified membrane begins to close regardless of pore size as the membrane relieves internal pressure. As seen in figure 4, the existence of a critical ratio remains in this scenario, as sufficiently large pores reverse the initial collapse and expand exponentially as in the first simulation. This behaviour is governed by equilibrium size, rather than initial size; the pore with r 0 /h ¼ 1.88, well above the critical ratio, closes regardless in this new scenario despite having expanded in the initial case of figure 3.

Membranes of varied thickness
The curvature of the membrane around the lip of the pore decreases with the membrane width h; increasing the thickness should decrease the energy cost of pore formation, thereby increasing the critical radius r c . It was shown in [1] that for a homogeneous film, where the free energy cost of pore formation is derived entirely from edge tension, the scaling should be linear as r c ¼ hp/4. To examine this scaling for the simulated lipid membrane, we changed the membrane thickness by altering the number of particles per lipid tail.
Initially, the tail length was increased from six to seven particles. The resulting membrane was found to be stable and at equilibrium with 29 200 lipids in the same simulation box, i.e. a lipid density increase of 1.74%. The resulting membrane had a width of approximately h ¼ 8.63 (8.15% thicker than the original membrane).
Results for this set are shown in figure 5. The single additional tail particle induced an upward shift in the critical ratio, to 1.48 , r c /h , 1.59, an increase of between 0.68 and 18.66% from the six-particle case. In addition, there was a notable change in the time scale on which pores evolved away from the critical region. In particular, the pore initialized at r 0 /h ¼ 1.59 was nearly stable for the entire duration of the previous simulations (t % 2000) before eventually opening up into supercritical growth. Figure 5 also demonstrates the stochastic nature of the dynamics at this spatial scale: in these realizations, the pore with initial size r 0 /h ¼ 1.48 closed faster than an initially smaller pore with radius r 0 /h ¼ 1.36.
Next, the tail length was increased further to 12 particles. The resulting time series can be seen in figure 6. For very small pores, the behaviour was unchanged, with the hole quickly being sealed. For pores of moderate size, the doubling of lipid tail length afforded significantly increased stability-above some size threshold, all pores in the membrane are stable indefinitely, showing a marked contrast with the pore radius/membrane width Tail lengths of eight and nine lipids were also considered but were found to produce nearly identical results to the case of 12 lipids per tail. The simulated bilayer membrane thus exhibits a sharp change in stability as the number of lipids per tail increases from seven to eight.

Pores in finite membranes
To further understand this change in stability, we finally considered the case of a finite membrane, i.e. a free-floating square patch of the membrane of finite size. To simulate such a membrane, lipids placed on a lattice in the initialization phase were truncated a fixed distance of 5 units away from the edge of the simulation box. The empty space left by truncating the membrane was included in the region of initialization for fluid particles.
All simulated finite membrane pores invariably sealed, regardless of initial radius. The membrane with six particles per lipid tail, which formerly exhibited a critical pore radius, transitioned through a metastable torus configuration to a layered cluster (see figure 7). The membrane with 12 particles per lipid tail, whose pores were stable above a small threshold radius, sealed its pore but remained stable in a finite bilayer disc. We hypothesize that the existence of a critical radius of pores in the periodic case corresponds directly to the stability of the finite membrane; a stable finite membrane prevents pores from expanding regardless of their initial size.

Discussion
Many of the experimental findings about the polymer films of [1] were also observed in our dissipative particle dynamics simulation of a bilayer membrane. In particular, the existence of a critical pore radius was observed in the periodic simulations when membrane tails comprised six to seven monomers. Ilton et al. [1] explain the mechanism for such a phenomenon by writing the energy cost of pore formation DG(r) as a function of the pore radius: (A 2 2pr 2 )g, in terms of the surface area A of the pore edge and the per-area surface tension of the film g. Modelling the pore edge as the inner half-surface of a regular torus with diameter h (the membrane width), the edge surface area A is given by p 2 hr 2 ph 2 . The resulting cost DG(r) is a concave function maximized when @ @r DG(r) ¼ (p 2 h À 4pr)g ¼ 0, yielding a critical radius r c ¼ hp/4 which scales linearly with the thickness h of the film. The cost DG(r) is a barrier for pore formation: once the critical size is reached, further expansion of the pore begins to reduce the free energy.
This expression ignores any energy cost associated with the molecular structure of the membrane; the authors also describe a modified argument for a diblock film, which has been used as a simple model of a lipid bilayer membrane in theoretical work [24]. Owing to the additional cost of rearranging molecules on the curved surface around the pore, they predict a critical radius larger than for a homopolymer film by a factor proportional to the non-dimensional curvature L/h, where L is the equilibrium thickness of the lamellar layers. Our simulation results agree in this respect: critical radii of the lipid membranes with six to seven monomer tails were found to be 1.34h , r c , 1.47h and 1.48h , r c , 1.59h, respectively, compared to the homopolymer r c ¼ hp/4 % 0.79h.
There also exist significant differences between our simulations and the work of Ilton et al., most notably in the case of the thicker membrane. As the diblock copolymer correction term scales with the non-dimensional curvature, its influence should decay as the membrane width increases, reaching the same limit of r c /h ¼ p/4 as h ! 1. By contrast, our simulations of a thicker lipid bilayer membrane showed pores above a certain size to be stable indefinitely, i.e. a critical radius above which pores expanded no longer existed. This suggests that the structure of a lipid bilayer membrane is fundamentally different to the structure of a diblock copolymer. Unlike the polystyrene film and PS-b-PMMA diblock copolymer of [1], whose critical radius was a continuous function of membrane width, the simulated bilayer membrane exhibits a phase transition from a regime where critical radius relates to thickness (less than or equal to seven monomer lipid tails) to a regime where arbitrarily large pores are stable (greater than or equal to eight monomer lipid tails).
Our simulations additionally provide insight into the time scale on which growth occurs. Although pores above the critical radius (in simulations where a critical radius existed) were found to increase in size roughly exponentially, the time scale of the exponential growth was markedly different between the membranes with tail lengths of six and seven particles. To compare these time scales, the time series for time evolution of perforated finite membranes . The original membrane seals its pore but simultaneously transitions into a layered cluster. Conversely, the thicker membrane seals its pore but remains in a stable finite disc indefinitely. The interior of the membrane in (b) comprises only tail particles and is not shown.
royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181657 the radius of the smallest supercritical pore was fit to an exponential r 0 e t/t in terms of the characteristic time t. The thinner membrane ( §4.1) was found to have a characteristic time of t 1 % 900, while the slightly thicker membrane ( §4.2) was found to have a characteristic time of t 2 % 2000. The time scale of pore evolution for lipid bilayer membranes is thus significantly affected by the structure/width, potentially in addition to chemical properties (in this context, the force coefficients for particle interaction, which were not changed between simulations). Since the DPD method explicitly models the solvent and does not make equilibrium assumptions, it is also suitable for examining the behaviour of lipid membrane pores in the presence of fluid flows or pressure gradients, such as those observed in biological cells during movement. Future work can examine the stability and dynamics of such pores in a variety of contexts of interest in cell biology.
Data accessibility. Code used to generate the visualizations and data in this study can be found in a public repository at https://github.com/clark-bowman/LAMMPS-PoratedMembrane.