Final Project: Molecular Simulations of Lipid Membranes
James Pelletier
MAS.864 Nature of Mathematical Modeling
Question: How to insert big molecules into biological cells?
This figure shows an image from my research project  these are bacteria about 1 micron in diameter. The fluorescence signal shows whole chromosomes, initially outside the cells. Our goal is to push the chromosomes into the bacteria. At the start, the chromosomes appear like rings because they coat the surfaces of the spherical cells, and through the microscope we view a projection of the system.
I have several interests in lipid membranes, the dynamic structures that contain and subdivide biological cells. Lipid membranes are semipermeable: many small molecules can diffuse through membranes; however, highly charged molecules such as DNA cannot pass through the hydrophobic barriers. In order to install DNA in cells, we must transiently open holes in the membranes. Transient electric fields and osmotic shock can permeabilize cells. Chemical treatment with divalent ions followed by osmotic shock can even push enormous pieces of DNA (megabases) into bacteria (about 1 micron in diameter), or fuse bacteria together. I would like to gain some intuition for the mechanisms that govern hole opening/closing in membranes and transport of materials through the holes.
In particular, how do osmotic pressures rupture membranes? How large are the holes? For how long do the holes live?
Prior art
Due to the molecular complexity of the biomolecules, atomistic simulations are limited to 10s of nanometers and a few nanoseconds integration time. Furthermore, membrane structure is a function of solvent quality  and explicit inclusion of all the solvent molecules costs a lot of computational effort. Therefore, I have been reading about coarsegrained, solventfree molecular dynamics simulations. "Coarsegrained" means we simulate minimalist molecular models to reduce the number of degrees of freedom and focus in on the core physical phenomena (Figure 1).
"Solventfree" means we do not explicitly include the solvent molecules; rather, we implicitly include the solvent molecules by adjusting the form of the intermolecular interactions between the coarsegrained lipids. I would like to learn about Molecular Dynamics/Conformations simulations, I think this would be a great system to explore.
These image come from Lin, C.M., Wu, D. T., Tsao, H.K. & Sheng, Y.J. Membrane properties of swollen vesicles: growth, rupture, and fusion. Soft Matter 8, 6139 (2012).
This image comes from Deserno, M. Mesoscopic Membrane Physics: Concepts, Simulations, and Selected Applications. Macromol. Rapid Commun. 30, 752–771 (2009).
My MD simulations
I combined the MD simulation method from Lin et al. and the coarsegrained lipid structure from Deserno et al.
Here is my code.
Dissipative Particle Dynamics (DPD)
 mesoscopic simulation method
 link between simulation parameters and FloryHuggins $\chi$parameters
 DPD fluid has quadratic equation of state
 simulate fragments to derive $\chi$parameters
 possible problem with timescale separation between particle diffusion and momentum diffusion (viscosity)
 soft condensed matter systems have relevant length scale between atomistic and macroscopic scales
 polymer gels have length scale set by distance between crosslinks
 simple bead and spring model predicts universal linear viscoelastic behavior
 indicates lifetimes/structures of crosslinks matters more than chemical details
 surfactant mesophases have length scale set by wavelength of fluctuations
 diverse phase behavior: lamellar structures, gels, micelles, networks
 lattice techniques not suited to branched or multiblock copolymers
 macroscopic simulations of associative liquids do not describe molecular details or hydrodynamic interactions
 [11, 12] Hoogerbrugge and Koelman introduced dissipative particle dynamics (DPD)
 soft spheres with collision rules
 bead and spring particles to simulate polymers
 [15] Español and Warren studied fluctuationdissipation theorem in DPD
 three forces: conservative, dissipative, random
 for statistical mechanics to correspond to canonical ensemble, dissipative and random forces must satisfy a relation
 check statistical mechanical validity as a function of the time step size
 simple scaling relation, to relate model parameters to FloryHuggins polymer theory
 can calculate solubility parameters from atomistic simulations
DPD Simulation Method
 time evolution from Newton's equations of motion
 set masses to unity
 sum forces with cutoff radius $r_c$
 set cutoff radius to unity
 conservative force, soft repulsion along line of centers
$$F_{ij}^C = \left\{ \begin{array}{c} a_{ij} \left(1  r_{ij}\right) \hat r_{ij}, r_{ij} < 1 \\ 0, r_{ij} \ge 1 \end{array} \right.$$
 dissipative force: acts opposite relative velocity between particles
$$F_{ij}^D =  \gamma w_D \left( r_{ij} \right) \left(\hat r_{ij} \cdot \vec v_{ij} \right) \hat r_{ij}$$
 random force: to conserve momentum, should use same value of random variable for both $F_{ij}^R$ and $F_{ji}^R$
$$F_{ij}^R = \sigma w_R \left( r_{ij} \right) \theta_{ij} \hat r_{ij}$$
 $w_D$ and $w_R$ are weight functions that vanish beyond the cutoff radius $r_c$
 $\theta_{ij}$ is a random variable with Gaussian statistics
$$\left< \theta_{ij} (t) \right> = 0, \left< \theta_{ij} (t) \theta_{kl} (t^\prime) \right> = \left( \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} \right) \delta (t  t^\prime)$$
 both dissipative and random forces act along line of centers, conserve linear and angular momentum
 independent random function for each pair (emphasized, see note above) of particles
 review fluctuationdissipation theorem
$$w_D(r) = \left( w_R(r)\right)^2 = (1  r)^2 r < 1, \sigma^2 = 2 \gamma k_B T$$
 chose $w_R(r)$ same functional form as conservative force
 FokkerPlanck equation governs time evolution of probability distribution
$$\frac{\partial \rho}{\partial t} = \mathcal{L}^C \rho + \mathcal{L}^D \rho$$
 $\mathcal{L}^C$ the Liouville operator of Hamiltonian system with conserved forces $F^C$
 $\mathcal{L}^D$ contains dissipative and noise terms
 without dissipation and noise, GibbsBoltzmann distribution a solution of $\partial \rho / \partial t = \mathcal{L}^C \rho = 0$
 with dissipation and noise, must still satisfy GibbsBoltzmann distribution $\mathcal{L}^D \rho = 0$
 review why dissipation and noise should not violate GibbsBoltzmann
 if new distribution matches GibbsBoltzmann, can use standard thermodynamic relations
 NVT molecular dynamics can simulate Hamiltonian system in canonical ensemble, why DPD?
 can view DPD as novel thermostat
 and DPD preserves hydrodynamics
 annealing defects in ordered mesophases
 dynamic density functional theory or Monte Carlo methods do not preserve hydrodynamics
 soft interaction potential another advantage over simulations with LennardJones atoms
 allows much larger time step
 to obtain soft potentials, average molecular field over rapid atomic motions
 does not diverge at $r = 0$ and lacks characteristic minimum of LJ potential
 needs complicated integrations, similar to virial coefficients in low density expansion
 so no analytical expression for effective interaction potential
 neglects correlation in local density fluctuations
 Groot and Warren suggest a pragmatic alternative
 top down rather than bottom up, match thermodynamics
 set $k_B T = 1$ specifies time, as rms velocity $\sqrt{3}$ from MB distribution
 how to choose time step
 compromise between fast simulation and equilibrium condition
 monitor system temperature
 no statistical difference between uniform and Gaussian random numbers
 uniform random numbers take less CPU time than Gaussian numbers
 artificial temperature increase influences time step choice
 compared to Euler, Verlet with acceleration term gives 50 fold improvement (both needed)
 equilibration time decreases as noise/dissipation amplitude increases
 temperature increases with noise amplitude
 how to choose repulsion parameter
 must choose $a$ to describe fluctuations in the liquid, related to compressibility
 analogous to WeeksChandlerAnderson perturbation theory of liquids
$$\kappa^{1} = \frac{1}{n k_B T \kappa_T} = \beta \left. \frac{\partial \rho}{\partial n} \right_T$$
 $\kappa_T$ is the isothermal compressibility
 obtain pressure via virial theorem as
$$p = \rho k_B T + \frac{1}{3V} \left< \sum_{j > i} \left( \vec r_i  \vec r_j \right) \cdot \vec f_i \right>$$
 radial distribution function: how density varies with respect to distance from particle
 surface tension: integral over stress across an interface
 excess pressure scales linearly with repulsion parameter, proportional to $\rho^2$
 CPU time increases with square of density
 choose lowest density that satisfies scaling relation
 $a = 75 k_B T / \rho$
 soft repulsion, apparent absence of $\rho^3$ term at high densities
 cannot produce liquidvapor coexistence
 review details
 steep repulsion at small $r$ spoils temperature control
$$\frac{F}{k_B T} = \frac{\phi_A}{N_A} \ln \phi_A + \frac{\phi_B}{N_B} \ln \phi_B + \chi \phi_A \phi_B$$
 do not favor contact $\chi > 0$

 to advance state
 simple Euler algorithm
 must include a factor of $\Delta t^{1/2}$ in random force
 integrate stochastic differential equations
 interpret random force as Wiener process
 let random force have zero mean and variance $\sigma^2$
 if variance has no dependence on time step, unphysical behavior
 mean square value of time integral proportional to the square distance i.e. to the time
 mean diffusion must have finite value independent of integration step size
 spread of random force increases as we subdivide the time interval
 modified version of velocity Verlet algorithm
 prediction for new velocity corrected in last step
 force updated once per iteration
 temperature updated after last step
 stochasticity makes order of algorithm unclear
 lambda = 1/2
Models
Per molecule, we have two nonpolar monomers and one polar monomer  this approximates the aspect ratios of the lipid constituents, important to the morphologies achieved by the lipids. Furthermore, the size of the head group compared to the size of the two tail groups can dramatically influence the equilibrium phase behavior.
 Monomer size/penetration: WeeksChandlerAndersen potential (WCA) $$V_{\rm rep}(r, b) = \left\{ \begin{eqnarray} 4 \epsilon \left[ \left( \frac{b}{r} \right)^{12}  \left( \frac{b}{r} \right)^6 + \frac{1}{4} \right] &,& \,\,\,\, r \le r_c \\ 0 &,& \,\,\,\, r > r_c \end{eqnarray} \right.$$
 Bonds between connected monomers in lipids: Finite extensible nonlinear elastic (FENE) $$V_{\rm bond} =  \frac{1}{2} k_{\rm bond} r_\infty^2 \ln \left[ 1  \left( \frac{r}{r_\infty} \right)^2 \right]$$
 Bending rigidity of lipid: harmonic spring $$V_{\rm bend} = \frac{1}{2} k_{\rm bend} (r  4 \sigma)^2$$
 Interaction between nonpolar tail groups: Need to be careful with this one in order to achieve fluid properties (if interaction too short range, then only achieve gas and gel phases)$$V_{\rm attr}(r) = \left\{ \begin{eqnarray}  \epsilon &,& \,\,\,\, r < r_c \\  \epsilon \cos^2 \frac{\pi \left( r  r_c \right)}{2 w_c} &,& \,\,\,\, r_c \le r \le r_c + w_c \\ 0 &,& \,\,\,\, r > r_c + w_c \end{eqnarray} \right.$$
Questions
I seek to understand membrane rupture. I will take this strategy
 First model vesicles.
 Next, add inert spheres inside/outside the vesicles.
 At what ramp rates/concentrations do the vesicles implode?
 How big are holes that form?
Approach
 1. To find equilibrium starting configurations, identify local minima in energy landscapes, using material from the search and optimization sections. To decrease time required to evaluate energy of system, truncate range of potential.
 2. To obtain some dynamical information, solve Newton's equations directly. Use an algorithm such as the Verlet algorithm. This will use information from the random systems and partial differential equations section.
 3. I will compete this against a package for molecular dynamics simulations, an MD environment called ESPResSo (Extensible Simulation Package for RESearch on SOft Matter), developed for simulation of soft matter systems with electrostatic interactions. I downloaded and complied Tcl/Tk and ran a few Espresso tutorials. For visualization, I have downloaded and used plugins to export Espresso output into VMD (Visual Molecular Dynamics).
Data structure
 from Numerical Simulation in Molecular Dynamics: Numerics, Algorithms, Parallelization, Applications
 by Michael Griebel, Stephan Knapek, and Gerhard Zumbusch
 our DPD model has shortrange interactions, need only consider neighbors
 in 2D, a 4Ndimensional phase space
 assume a square or rectangular simulation domain
 periodic boundary conditions (PBC) to compensate for limited size of simulation domain
 for LennardJones potential, cutoff radius typically about $2.5 \sigma$
 in our case, we define cutoff radius in DPD interactions
Linked Cell Method
 split physical simulation domain into uniform subdomains
 choose cell size equal to or greater than cutoff radius
 then particles in a cell interact only with particles in the cell and in neighboring cells
 compute force on particle as sum over particles within sum over neighboring cells
 complexity reduced from $\mathcal{O}\left(N^2 \right)$ to $\mathcal{O}\left(N \right)$
 need appropriate data structures
 in 2D, assign indices $\left(c_x, c_y \right)$ to cells; each cell has 8 neighbors
 particles located in one cell stored in a linked list associated with the cell
 need a dynamic data structure for the particles of each cell
 [559, 560] linked list references
 each element consists of the particle and a pointer next
to the next element of the list
 remove from one list, and insert to the beginning of another list
 i>next = *root_list;
 i is a pointer to the inserted particle
 sets the next pointer to point to the current first particle
 *root_list = i;
 i is a pointer to the inserted particle
 sets the pointer to the new first particle i
 with removal operation, must change next of previous element
 in a simple list, no direct way to find previous element
 easy to implement if pointers in both directions
 if inflow or outflow conditions, should remove element and perhaps release memory
 with PBC, need not worry can just move elements from one list to another
 cell represented as root pointer of the list
 size and position of cell determined from global parameters
 define vector grid of cells of type cell (root_list pointers)
 take 2d and vectorize
 i>p points to the first particle in the linked list
 i = i>next to advance from one element to the next
 outer loop over all cells, inner loop over all particles
 can put the simulation into the idle loop
 decided to make one structure, with an integer that flags the molecule as either a lipid or a water, so the linked list just points to a single structure
 had trouble making an array of pointers
 fixed by using the ampersand operator
 figured out how to change the view angle in OpenGL
 can use the glOrtho to do 3d slices through the 3d volume as soon as I set that up
 finished initialization of the particles
 GOT TO FIGURE OUT pointers and linked lists!
 remember to put ; after structure definitions!
 got the insert function working! woo hoo!
 decided to make a doubly linked list in order to facilitate the remove function
 memory allocation problems..
 got the linked list working! :) woo hoo!
Implementation
 linked cell method https://www10.cs.fau.de/Teaching/Courses/SS2009/AlgCSE/presentation_aprelkin.pdf
 linked list associated with each cell
 cell represented by root pointer to the cell
 grid cells stored in a vector
 can decrease number of forces by summing over pairs/forces rather than particles
 C tutorials by Alex Allain
 linked list: http://www.cprogramming.com/tutorial/c/lesson15.html
Rendering
 thought to use VMD
 realized just plotting spheres, did not need all the complexity
 decided to practice with OpenGL
 found a tutorial here http://www.swiftless.com/tutorials/opengl/collision.html
 had to remember how to compile the openGL code on Mac
 openGL has a stack of matrices
 push copies the top matrix and adds it to the top of the stack
 pop deletes the top matrix from the top of the stack
 initially, each of the stacks contains one matrix, the identity
 on a Mac, need to include #include header file
 to compile need to include framework flags
 gcc gl_spheres.c o gl_spheres framework OpenGL framework GLUT
 can make a bunch of spheres by doing
glPushMatrix();
position = (float)x / 20;
glTranslatef(position, position, position);
glutSolidSphere(0.05, 40, 40);
glPopMatrix();
 tried to make spheres at random locations, will need randomness in force expression
 decided to write LFSR in c, had written in python for a previous assignment
 learned about number representations in c
 the "0b" prefix did not work right away
 read about hexadecimal representation, starts with "0x"
 from http://www.cprogramming.com/tutorial/bitwise_operators.html
 the byte is lowest level at which we can access data
 most significant bit to the left
 need to be careful about storage, can end up with modular arithmetic
 bitwise AND: ampersand &, can use to find the bit in the nth position
 can use bitwise leftshift to generate the mask
 bitwise OR: line 
 bitwise complement: ~
 bitwise XOR: ^
 tried the bitwise complement and encountered twos complement
 struggled with adding numbers together  e.g., tried to add 1 and 8 to make 9, and instead got 32776, I think it has something to do with adding the numbers as binary variables rather than as integers, and with the signed integers, twos complement thing.
 NO! just had to do with initial value of variable when declared was 2^15..!?
 ran into trouble passing variables into main  seemed fixed after I put the function arguments into main rather than passing them in (function main takes zero or two arguments)
 Charles taught me about textmate, a convenient editor that lets me run the code without switching back and forth between the editor and the terminal to compile then run the code also saved at the same time
 finished implementation of the LFSR  can treat numbers as binaries and act on them directly.
 wait a sec  worked with M = 4, but did not immediately work with M = 5
 was indexing in the wrong direction, and had gotten lucky because before had worked because M = 4 takes the first and last one, whereas M = 5 the first and second to last, so the second to last didn't catch the first and all the bits got pushed off the register. note that all the indices in the register have the first bit in the tap, otherwise 1 would not cycle back
 M = 23 somehow looked much worse than M = 24
 no made a mistake (assigned 2 and 3 to 0 and 1), but still looked random
 great openGL tutorials at http://www.swiftless.com/tutorials/opengl/color.html
 just needed glEnable(GL_COLOR_MATERIAL) and then could set color with glColor
 decided to start in 2d and then graduate to 3d
 as not to waste too much time with the rendering etc.