Analysis of Effects of Shockwaves on Cholesterol-Containing Lipid Bilayer Models

In the following work, the effect of shock impulses (2.5–10.0 mPa * s) on lipid bilayer membranes consisting of POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phoshocholine) and varying amounts of cholesterol (10% and 50%) has been studied using molecular dynamics. Both models were equilibrated for 7 nanoseconds in an NPT ensemble at 310 K and 1 atm. After equilibration, the 50% cholesterol model showed a larger average membrane thickness and an increased average deuterium order parameter (SCD) of carbons in lipid tails when compared to the 10% model. In addition, values such as the total area per lipid and area per POPC were reduced in the 50% model when compared to the 10% cholesterol model. After application of similar shock impulses to the 50% and 10% cholesterol models, the structural effects on each respective bilayer model were compared. The change in bilayer thickness due to shock impulse was measured for both models at all impulses tested. These values were used to calculate the percent decrease in bilayer thickness of both models. On average, the 50% cholesterol model showed a 1-2 % lesser reduction in bilayer thickness then the 10% cholesterol model at the shock impulses tested. After application of shock impulses to both models, the temporal change of the instantaneous averaged order parameter SCDi of the upper and lower monolayers of both models was calculated at shock impulses of 7.5 mPa * s and 10. mPa. * s. In both models, regardless of shock impulse, the SCDi value of the upper monolayer decreased followed by a decrease in the SCDi. parameter of the lower monolayer. These effects can be linked to the collapse phase of the lipid bilayer due to shockwave. Once the shockwave has fully traveled through the lipid bilayer, the SCDi parameter of the lower monolayer begins to increase which indicates beginning of the rebound phase. Finally, the increased fluidity of the lipid bilayer models was examined through calculation of the lateral displacement of the lipids in both models during the shockwave impulse ranges tested. The 50% model showed a slightly higher later displacement at the shock impulse ranges tested but not significantly higher to make any conclusions. In addition, the efficiency of the shock impulse based on the distance of the shock water slab from the lipid bilayer was analyzed. The right boundary of the shock slab was placed 1 Ǻ and 10 Ǻ away from the lipid bilayer and similar shock impulses of 2.5, 5.0, 7.5 and 10 mPa * s were applied. It was shown that placement of the shock slab 10 Ǻ away from the lipid bilayer causes a greater decrease in average bilayer thickness and slight increases in kinetic energy, temperature and pressure of the system. This could be owed to better momentum transfer among molecules achieved when the water slab is placed 10 Ǻ away rather than 1 Ǻ away.


Abstract
In the following work, the effect of shock impulses (2.5-10.0 mPa * s) on lipid bilayer membranes consisting of POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phoshocholine) and varying amounts of cholesterol (10% and 50%) has been studied using molecular dynamics. Both models were equilibrated for 7 nanoseconds in an NPT ensemble at 310 K and 1 atm. After equilibration, the 50% cholesterol model showed a larger average membrane thickness and an increased average deuterium order parameter (S CD ) of carbons in lipid tails when compared to the 10% model. In addition, values such as the total area per lipid and area per POPC were reduced in the 50% model when compared to the 10% cholesterol model.
After application of similar shock impulses to the 50% and 10% cholesterol models, the structural effects on each respective bilayer model were compared. The change in bilayer thickness due to shock impulse was measured for both models at all impulses tested. These values were used to calculate the percent decrease in bilayer thickness of both models. On average, the 50% cholesterol model showed a 1-2 % lesser reduction in bilayer thickness then the 10% cholesterol model at the shock impulses tested.
After application of shock impulses to both models, the temporal change of the instantaneous averaged order parameter S CDi of the upper and lower monolayers of both models was calculated at shock impulses of 7.5 mPa * s and 10. mPa. * s. In both models, regardless of shock impulse, the S CDi value of the upper monolayer decreased followed by a decrease in the S CDi. parameter of the lower monolayer. These effects can be linked to the collapse phase of the lipid bilayer due to shockwave. Once the shockwave has fully traveled through the lipid bilayer, the S CDi parameter of the lower monolayer begins to increase which indicates beginning of the rebound phase.
Finally, the increased fluidity of the lipid bilayer models was examined through calculation of the lateral displacement of the lipids in both models during the shockwave impulse ranges tested. The 50% model showed a slightly higher later displacement at the shock impulse ranges tested but not significantly higher to make any conclusions.
In addition, the efficiency of the shock impulse based on the distance of the shock water slab from the lipid bilayer was analyzed. The right boundary of the shock slab was placed 1 Ǻ and 10 Ǻ away from the lipid bilayer and similar shock impulses of 2.5, 5.0, 7.5 and 10 mPa * s were applied. It was shown that placement of the shock slab 10 Ǻ away from the lipid bilayer causes a greater decrease in average bilayer thickness and slight increases in kinetic energy, temperature and pressure of the system. This could be owed to better momentum transfer among molecules achieved when the water slab is placed 10 Ǻ away rather than 1 Ǻ away. Depth of potential well in Lennard-Jones potential function σ Distance where inter-particle potential (U LJ ) becomes zero r ij Distance between particle i and particle j U Coulomb Coulomb potential energy r eq Equilibrium distance between two bonded atoms θijk Bend angle between successive bonded atoms i, j,k  Figure 6.2: Change in bilayer thickness of 9:1 POPC/CHL model at varying shock impulses with a 10 Ǻ distance of water slab from bilayer……………………..…65 Figure 6.3: The evolution of bilayer thickness of the 1:1 POPC/CHL model at a 5 mPa * s shock impulse compared for different water slab distances tested (1 Ǻ vs 10 Ǻ)……………………………………………………………………………….....66 Figure 6.4: The evolution of bilayer thickness of the 1:1 POPC/CHL bilayer at a 10 mPa * s shock impulse compared for different water slab distances tested (1 Ǻ vs 10 Ǻ)…………………………………………………………………………….........67 Figure 6.5: The evolution of bilayer thickness of the 9:1 POPC/CHL bilayer at a 5 mPa * s shock impulse compared for different water slab distances tested (1 Ǻ vs 10Ǻ)……………………………………………………………………………….68 Figure 6.6: The evolution of bilayer thickness of the 9:1 POPC/CHL bilayer at a 10 mPa * s shock impulse compared for different water slab distances tested (1 Ǻ vs 10 Ǻ)…………………………………………………………………………….……68 Figure 6.7: Rise in kinetic energy of the 1:1 POPC/CHL system due to a 7.5 mPa shock impulse applied at different water slab distances (1 Ǻ vs 10 Ǻ)..………….70 Figure 6.8: Rise in kinetic energy of the 1:1 POPC/CHL system due to a 10 mPa shock impulse applied at varying water slab distances (1 Ǻ vs 10 Ǻ)……………70 Figure 6.9: Rise in kinetic energy of the 9:1 POPC/CHL system due to a 7.5 mPa shock impulse applied at different water slab distances (1 Ǻ vs 10 Ǻ)………..…71 Figure 6.10: Rise in kinetic energy of the 9:1 POPC/CHL system due to a 10.0 mPa shock impulse applied at different water slab distances (1 Ǻ vs 10 Ǻ)…………..71 Figure 6.11: Rise in pressure of the 1:1 POPC/CHL system due to a 10 mPa * s shock impulse applied at different shock slab distances (1 Ǻ vs 10 Ǻ)……………..…..73 Figure 6.12: Rise in pressure of the 1:1 POPC/CHL system due to a 10 mPa * s shock impulse applied at different shock slab distances (1 Ǻ vs 10 Ǻ)……………...….73 Figure 6.13: Rise in pressure of the 1:1 POPC/CHL system due to a 10 mPa * s shock impulse applied at different shock slab distances (1 Ǻ vs 10 Ǻ)………………....74 Figure 6.14: Rise in pressure of the 1:1 POPC/CHL system due to a 10 mPa * s shock impulse applied at different shock slab distances (1 Ǻ vs 10 Ǻ)………………....75 xiv List of Tables   TABLE  PAGE   Table 3.1: Atom composition and number of atoms in the 9:1 POPC/CHL and 1:1 POPC/CHL lipid bilayer systems…………………………………………………….21 Table 3.2 Initial Properties of the 9:1 POPC/CHL and 1:1 POPC/CHL lipid bilayer systems. ………………………………………………………………………………22 Table 4.1: Shockwave Simulation Setups for lipid bilayer models at varying shock impulses and water slab distances…………………………………………………….31  (Khan et al. 2013). Also present in the lipid bilayer are membrane proteins which can either serve as channel proteins such as the potassium ion channel or integral proteins, such as integrin. In addition, the lipid bilayer contains peripheral proteins which can either be attached on the outside of the membrane or on the cytosolic portion (Singer & Nicolson 1972). Overall, the lipid bilayer structure (figure 1) is composed of many different structures which all serve specific purposes both in regulation of the lipid bilayer as well as maintaining structural integrity. The Amphiphilic Nature of Phospholipids The ability of the lipid bilayer to self-assemble when exposed to water into its final structure is due to the amphiphilic nature of the lipid molecules present. The lipid bilayer is composed of phospholipids which have a hydrophilic phosphate head group and two hydrophobic fatty acid lipid tails. When exposed to water, the lipid molecules orient themselves as to maximize the interaction between water and their phosphate head groups and simultaneously minimize the exposure of the hydrophobic lipid tails to water (Alberts et al. 2002). This mechanism of assembly is also known as the hydrophobic effect. In addition, the favoring interactions between the lipid tails of each phospholipid keep water out of the hydrophobic region. This creates a sheet of lipid bilayer composed of two monolayers which is in a state of lowest free energy. In between the monolayers is the hydrophobic region which contains no water molecules.

Composition of Lipid Membranes
The mammalian lipid membrane is made up of almost 50% lipids while the remainder is proteins such as integral membrane proteins or peripheral proteins. The specific composition of the lipid bilayer is of heavy interest because it dictates properties of the membrane such as phase transition and fluidity. The amphiphilic lipid molecules contained in the mammalian lipid bilayer vary in head and tail structure but are generally limited to four different lipid types; the first three, phosphatidylethanolamine, phosphatidylserine, phosphatidylcholine, belong to the glycerophospholipids group whose structure consists of a sn-glycero-3-phosphate core that is esterified at its C1 and C2 positions into fatty acids tails (Alberts et al. 2002).
The fourth, sphingomyelin, belongs to the sphingolipid group and the structure consists of a choline molecule attached to a hydroxyl group from ceramide (Khan et al. 2013). The structures of these four phospholipids are shown in figure 1.3. In addition to structural difference that exists between these phospholipid types, the specific location of each phospholipid in the lipid bilayer can vary as well. Some phospholipids exist mostly in the cytosolic monolayer of a bilayer membrane while others are predominantly in the outer monolayer of the lipid bilayer (Janmey & Kinnunen 2006). For instance, glycolipids tend to be located in the outer monolayer and are thought to play a role in formation of lipid rafts. In mammalian red blood cells, choline-containing phospholipids such as phosphatidylcholine and sphingomyelin are primarily located in the outer monolayer while phosphatidylserines and phosphatidylethanolamines exist in the inner monolayer. Due to such a variety in lipid types, the lipid bilayer is asymmetrical and its structure and properties can vary even between lipid bilayers of the same cell types.

-Phosphotidylcholines and POPC
The most abundant phospholipid in mammalian cell membranes is phosphatidylcholine (Alberts et al. 2002). It is a constituent of human lung surfactants and serum lipoproteins and is the most frequently studied phospholipid in experimental lipid bilayer work. Therefore, the phase behavior and transition in the presence of water is of utmost interest to researchers in the field.
The phosphatidylcholine structure consists of a choline head group linked to a glycerophosphoric acid group attached to two fatty acid tails which can vary in structure and length (Alberts et al. 2002). Generally, phosphatidylcholines bear one saturated and one unsaturated fatty acid tail. An example is POPC (1-palmitoyl-2oleoyl-sn-glycero-3-phosphocholine), which contains a saturated (palmitoyl) and an unsaturated (oleoyl) fatty acid tail in its structure depicted in figure 1.4. Another key characteristic of lipid membranes is the temperature at which the lipid bilayer transitions from the solid to liquid phase, also known as the melting temperature (T m) (Berkowitz 2009). The T m is heavily affected by the fluidity and structural components of the lipid bilayer. For instance, the structure of the fatty acid tails of lipids in the membrane can either increase or decrease the T m . Longer fatty acid tails result in more van der WAALS attractive forces and increase the melting temperature while shorter chains decrease it (Alberts et al. 2002;Siegel et al. 1999). In addition to fatty acid tail length, the presence of double bonds in the structure creates kinks in the hydrocarbon chains and increases the mobility of lipids in the membrane.
The result of this is that lipids cannot pack as tightly together then if there were no double bonds (saturated) and therefore the T m is decreased. In general, the higher the degree of unsaturation (more double bonds), the lower the transition temperature of the lipid bilayer. Other factors that may influence the T m include the pH of the environment and the size of the bilayer. Overall, the combination of different lipid structures leads to mixed phase behaviors of the lipid bilayers which depend highly on which lipid types are contained in it and their orientation in the membrane.

-Cholesterol
In addition to phospholipids, another important lipid found in membranes is cholesterol, belonging to the sterol group (modified steroid). Its structure, shown in

The Condensation Effect of Cholesterol
The presence of cholesterol in lipid membranes, specifically in mammalian cells heavily influences the properties of the lipid bilayer. Typically, mammalian cells can contain a range of 20-30 % mol cholesterol in their membranes and up to 50% mol in red blood cells (Róg et al. 2009) and therefore the presence of cholesterol highly affects both the structural and thermodynamic properties of the lipid bilayer.
Structurally, cholesterol has been proven to have a "condensing effect" on lipid bilayers (Hung et al. 2007) (Daly et al. 2011). This condensing effect on the lipid bilayer has been shown in various experimental models (Mitchell & Litman 1998) (Oldfield et al. 1978) and in-silico models (Robinson et al. 1995) (Hung et al. 2007) (Martinez-Seara et al. 2008) . This is due to the structure of cholesterol which contains a polar hydroxyl head group which binds closely to the head groups of phospholipids around it, increasing the interaction strength between adjacent phospholipid molecules and ultimately leading to a more tightly packed lipid bilayer (Hung et al. 2007). This is the reason cells regulate the levels of cholesterol in their membranes as a way of controlling the fluidity and adjusting to different environmental conditions such as a higher or lower temperature. Overall, this condensing effect impacts structural properties of the lipid bilayer such as the thickness of the membrane, area per lipid, and the deuterium order parameter of the lipid tails (S CD ). Lipid bilayers with higher levels of cholesterol have been shown to have thicker bilayers on average, a reduced total area per lipid, and an increased order parameter of carbons in acyl chains. (Mannock et al. 2010).
In addition to structural effects, thermodynamic properties of lipid membranes are also altered by the level of cholesterol present. Specifically, the melting temperature is increased in bilayers with higher levels of cholesterol (Alberts et al. 2002). With more cholesterol, the lipid bilayer becomes more tightly packed during the gel phase and therefore requires more thermal energy to phase transition into the liquid-crystalline phase. It becomes clear that the presence of cholesterol greatly affects the properties of the lipid bilayer, one of the driving forces behind this research.

-Theory Behind Molecular Dynamics
In 1957, a group of scientists named Berni Alder and Thomas Wainwright published breakthrough research describing a numerical method to study the physical movements of atoms and molecules in an N-body simulation (Allen 2004). This research was the foundations of molecular dynamics computer simulations. Molecular dynamics simulations are commonly used to study specific phenomena at the microscopic level. Experimental work is important but is limited to studying macroscopic properties of a system in most cases. Therefore, molecular dynamics simulation can provide clues to experimentalists as to exactly what is happening to a system at the atomistic level. Many times, collaborations between experimentalist and computer simulation experts are key for unveiling the fundamentals behind important phenomena.

Foundation of Molecular Dynamics
Molecular dynamics simulation consists of numerically solving Newton's 2nd law of the classical equations of motion for a set of atoms in the system (2.1).

= ̈
(2.1) To perform this task, the force acting on every atom F i must be obtained and is derived from a potential energy function U (r N ) where r N represents the complete set of 3N atomic coordinates of the system studied (equation 2.2) (Allen 2004) .

= − (2.2)
This is because the force can be calculated as the negative gradient of the potential energy. Before the force is obtained, the potential energy acting on every atom in the system must be calculated. To calculate the potential energy on a given atom, there are several atomic interactions which must be considered for any given system.

Non-Bonded Interactions
The first type of interactions considered is one that occurs between a pair of nonbonded atoms. The most commonly used potential to describe non-bonded interactions is the Lennard-Jones potential (U LJ ) shown below, where ɛ is the depth of the potential well, σ is the distance where the inter-particle potential reaches 0 and r ij is the distance between the particles (Tieleman et al. 1997).
The r -12 term is the repulsive term and is based on the Pauli repulsion occurring between nearby atoms due to overlapping electron orbitals. Meanwhile the r -6 attractive term describes the attraction between two atoms due to long range forces such as van der Waals (Allen 2004). If charges are present between the set of interacting atoms, the Coulomb potential is used (equation 2.4), where Q 1 and Q 2 are the charges of the atoms, ε 0 represents the permittivity of free space and r ij is distance between the particles. This potential is repulsive for atoms with similar atomic charges and attractive for atoms with opposite electric charges.
Traditionally, calculation of non-bonded forces between atoms in a system is the most computationally expensive part of a simulation. Generally, during simulations a cutoff distance is provided which eliminates the calculation of non-bonded forces that are separated by a certain distance threshold. This cutoff distance saves calculation and computational time during simulations.

Bonded-Interactions
The intramolecular potential between a set of bonded atoms must be calculated as well. This consists of the sum of three separate potentials which consider bond distance, bend angles and torsion angles between a set of atoms (2.5a-2.5c). (2.5b) The first bonding potential considers the separation distance between a pair of bonded atoms where r ij = |r i -r j |. A harmonic form is utilized with a specified equilibrium distance r eq . The second term considers the bend angles, Θ ijk , between successive bond vectors of three atoms. Finally, the torsion angle term, Φ ijkl , describes the torsion angles between three successive bonds of four atoms. Traditionally, force-field simulations packages such as CHARMM , will provide the accurate forms of equations as well as the strength parameters, k, and other constants such as r eq and θeq .
Generally, these force-fields are parameterized by using a combination of empirical techniques as well as quantum-mechanics calculations. The accuracy of these forcefields is tested by the ability to reproduce structural, dynamic, bulk and thermodynamic properties of small molecules whose properties have been well described experimentally.
Once the potential energy function has been calculated, the forces acting on atoms in the system can be calculated (formula 2.2). Once the forces are calculated, a numerical integration algorithm, such as the velocity-Verlet algorithm (Swope et al. 1982) can be used to calculate the positions and velocities of atoms at time (t +Δt) given the initial coordinates and velocities of the system studied. A trajectory of the system is then produced describing the dynamic evolution of a system over an elapsed time length.

Thermodynamic Ensembles
In molecular dynamics, the choice of which statistical thermodynamic ensemble employed to study a specific model of interest depends on which properties the researcher is interested in. Generally, there are three common thermodynamic ensembles employed during molecular dynamics simulation, the microcanonical ensemble, the canonical ensemble, and the isothermal-isobaric ensemble. In the microcanonical ensemble (NVE), the system is represented as having a fixed number of atoms (N), a constant volume (V), and a constant energy (E) (Frenkel & Smit 2002). Generally, this statistical ensemble is utilized when the energy of a system is known but no information about the internal states is available. This is an unrealistic ensemble because in experimental setups, energy is generally a varying property.
In the canonical ensemble (NVT), the system is represented as having a fixed number of atoms (N), constant volume (V), and constant temperature (T) (Frenkel & Smit 2002). This ensemble accurately represents a system that is in immediate contact with a heat bath. Finally, the isothermal-isobaric (NPT) ensemble represents a system with a fixed number of atoms (N), constant pressure (P) and constant temperature (T) (Frenkel & Smit 2002). This thermodynamic ensemble is considered the more experimentally realistic ensemble due to temperature and pressure being relatively constant during experiments. It represents an experiment being carried out in an open flask exposed to environmental temperature and pressure. Regardless, the choice of ensemble depends highly on the in-silico model to be simulated. Certain ensembles are better suited at reproducing properties for one molecule but produce errors in properties in other models placed in these ensembles.

Molecular Dynamics Time step
The vital component of any molecular dynamics simulation is the time step chosen to perform numerical integration on the system modeled. To accurately reproduce the dynamics of an all-atom system and avoid discretization errors, the time step chosen must be smaller than the fastest vibrational frequency in the system modeled. Since the vibrational frequency of hydrogen molecule is 10 -14 s, a time step of 10 -15 s or 1 femtosecond (fs) must be employed (Allen 2004). At such a small time step, simulating phenomena of a few microseconds becomes quite impossible without the use of supercomputers. One way to increase this time step is to constrain the bond lengths of hydrogen-containing bonds during simulation using algorithms such as SHAKE or RATTLE which allows the use of a 2 fs time step. Overall, care must be taken in choosing which phenomena to model using molecular dynamics due to the time step constraint and the time scales that can be reached.

Shockwaves and Simulation Methods
The choice of what molecular event to study using molecular dynamics is highly dependent on the time lengths which can be reached through computer simulation. In an all-atom molecular dynamics system, it is highly unlikely that a researcher will want to observe an event which occurs in the time scale of milliseconds (10 -3 s) because simulating this event would years; depending on the number of atoms in the system and the number of CPUs (Central Processing Units) utilized to carry out the simulation. Even by coarse-graining particles in the system, these time scales are still out of reach. Realistically, events studied should occur within the nanosecond (10 -9 ) and at most microseconds (10 -6 s) time scale.
The high-speed phenomena that occurs during a shockwave event, is what makes this process so applicable for molecular dynamics simulation. Shockwaves occur when a wave exceeds the local speed of sound in a fluid and is characterized by a discontinuous change in pressure, temperature, and density of the medium it is occurring. The speed of sound in water is roughly 1432 m/s (Sliozberg 2013) and therefore since shockwaves travel above the speed of sound, this high-speed phenomenon is easily reachable in the time scales of molecular dynamics

-Importance of Studying Shockwaves
Drug delivery and gene therapy methods heavily utilize the application of shockwaves to manipulate the permeability of the lipid bilayer contained in cells of interest. Application of shockwaves has been proven to cause temporary permeability of the lipid bilayer allowing for delivery of drugs or oligonucleotides into a patient's cells. In cancer therapy, methods such as high intensity focused ultrasound (HIFU) which utilize shockwaves are used to treat cancerous tissue and cause destruction of cancerous cells through conversion of mechanical energy into heat (Espinosa et al. 2014) . In addition, shockwave pulses are used to treat patients suffering from kidney stones in a treatment called shockwave lithotripsy (Williams et al. 1999). Although the treatment is effective, kidney tissue can be damaged due to effects of the shockwave application. Regardless of application, the study of shockwaves and the effects they can have on our cells is of great interest.
At the macroscopic level, the effect of shockwaves on the lipid bilayer has been studied intensively. Microscopically, much work remains to be done to understand the effects on the lipid bilayer structure at the atomic level. This knowledge can provide the medical field with ways to make improvements to currently implemented methods which utilize shockwaves.
Past research has probed the effects of shockwaves on lipid bilayer structural properties in-silico through molecular dynamics research. Specifically, the Koshiyama group (Koshiyama et al. 2006a) modeled the effects of shockwaves on DPPC (1,2dipalmitoyl-sn-glycero-3-phosphocholine)-containing lipid bilayer models. They noted changes in the bilayer structure due to shockwaves such as the rebound and collapse of the lipid bilayer and increased water permeation. Properties such as the accumulated lateral displacement and the temporal evolution of the instantaneous averaged deuterium order parameter of lipid tails were measured. Two years later (Koshiyama et al. 2008), this same group tested the effects of shock impulses impacting at various incident angles to DPPC lipid bilayer models. In addition, several groups tested effects of shockwaves on coarse-grained lipid bilayers using different coarse grained models as well as different force fields; the first using the DPDE force field (Ganzenmüller et al. 2011) and the second group using the MARTINI force field (Santo & Berkowitz 2014). Although this research was important, the all-atom in silico lipid bilayer models tested were very simplistic, only containing one phospholipid type and did not contain cholesterol, a vital component in mammalian cell membranes. While in contrast, the coarse-grained models lacked both atomic detail and cholesterol in their membranes modeled. Regardless, the addition of cholesterol to lipid membranes alters structural properties such as the average membrane thickness, the area per lipid and the deuterium order parameter of acyl chains. The models previously tested, have different structural properties since they lack cholesterol. In addition, cholesterol-free lipid mammalian membranes hardly exist in nature and therefore, this research will provide a more realistic analysis of the effects of the lipid bilayer due to shockwave than prior research has.
The first goal of this research is to apply shock impulses of varying magnitude (2.5 -10 mPa * s) on two different in-silico lipid bilayer models. The first model contains a 9 to 1 ratio of POPC to cholesterol (10% cholesterol). This model will be referred to as the 9:1 POPC/CHL model. The second model consists of a 1 to 1 ratio of POPC to cholesterol (50% cholesterol). This model will be referred to as the 1:1 POPC/CHL model throughout this thesis. Structural effects such as the change in bilayer thickness, the temporal change in instantaneous averagde deuterium order parameter S CDi of lipid tails, and the lateral displacement of lipids will be quantified at the range of shock impulses tested for both models.
The second goal of this thesis research is to analyze the effect of the distance of the water slab from the lipid bilayer on the model systems. This will be tested by placing the shock slab 1 and 10 Ǻ away from the lipid bilayer models studied and measuring the effect on the change in bilayer thickness, kinetic energy and pressure on the systems studied. This will elucidate which mechanism is more effective in this setup for producing the strongest shock impulse.

Setup of Initial Lipid Bilayer Model Structures Containing POPC and Cholesterol
The membrane models tested in these simulations were created using the online software, CHARMM-GUI (Lee et al. 2016). This program allows the user to

Solvation and Ionization
Once the initial bilayer models were created, they were solvated explicitly using the TIP3P (Jorgensen et al. 1983) model for water molecules. A total of 23 nanometers (nm) of solvent was added in the z dimension. More specifically, 20 nm of solvent was added to the right of the bilayer model (z+) and 3 nm to the left (z-). The excess solvation in the z+ side of the lipid bilayer is needed in order to allow shockwave propagation to be modeled accurately, as previously noted by the Koshiyama group (Koshiyama et al. 2006a). Solvent was not added in the x and y dimension due to the lack of water molecules in the hydrophobic region of the bilayer.
Once solvated, sodium and chlorine ions were added to both systems until a concentration of .15M was reached. The solvation and ionization was performed using the program VMD (Visual Molecular Dynamics) (Humphrey et al. 1996).   In addition, snapshots of the final systems for both models are shown (

Equilibration of Lipid Bilayer Models
An initial 10,000 steps of minimization was carried out for both models using the conjugative gradient method (Fletcher & Reeves 1964) in NAMD (Malterer et al. 2005). Following minimization, equilibration of models was carried out in an NPT ensemble at 310 K and 1 atm using a time step of 2 fs. To use this timestep, all hydrogen-containing bonds were constrained using the SHAKE algorithm described here (Ryckaert et al. 1977). Constant temperature was maintained using Langevin dynamics with a damping coefficient of 1/ps. Constant pressure was controlled using a Nosé-Hoover Langevin piston algorithm (Feller et al. 1995 Long range electrostatic were treated using the Particle Mesh Ewald method (Essmann et al. 1995) which has been shown to be a good method in accurately calculating longrange forces in membrane simulations (Vermeer et al. 2007). The CHARMM 36 force field for lipids was used and is described here (Venable et al. 2010). An example equilibration configuration file for the 1:1 POPC/CHL model can be found in the appendix (APPENDIX II) provided. This configuration file is similar to the 9:1 model, varying only in file names.
The equilibration runs were set to 20 ns but after 7 ns of equilibration, values such as average membrane thickness, average area per lipid, temperature and pressure had converged for each system and therefore the equilibration runs were terminated. The final restart velocities and coordinates of these equilibrations were used in the following shockwave simulations

-Shockwave Implementation
The previously established methods of Koshiyama et al. (Koshiyama et al. 2006a) were used to model the shockwaves in the following simulations. The impulse in this work is defined as the impulse per unit area (Î), generally called the specific impulse.
The specific impulse can be defined as the change of momentum(M) over the area (A) where the pressure is exerted as shown in the formula below, Based on previous experiments by Kodama et al. (Kodama et al. 2000) that applied varying shock impulses to leukemia cells in water, it was shown that at the initiation of a shockwave, the shockwave front had not reached water and cells ahead of it and therefore the momentum of water molecules at shock initiation was 0. After a small distance traveled at time (t + ), the momentum of water molecules was shown to be the impulse applied times the area (I x A). From this finding, the shockwaves in the simulations were modeled as a change of momentum of water molecules in a slab of volume L z * A directly adjacent to the lipid bilayer; where L z represents the thickness of the water slab and A represents the lateral area of the water slab in the x and y plane. The water slab will be termed "shock slab" from this point on in the thesis. To  cholesterol (green). The blue slab represents the shock slab whose right boundary is placed 1 Ǻ away from the lipid bilayer while the white slab is placed 10 Ǻ away.
Other atoms are removed for visualization purposes.
Finally, to add the average velocity calculated, a MATLAB code was written to calculate the average velocity added to each water molecule in the slab by

-Shockwave Simulation Setup
Shock impulses of 2.5, 5, 7.5, 10 mPa * s were applied to each model. The shockwave simulations were carried out in an NVE (constant number of atoms, N, volume, V and total energy, E) ensemble in order to analyze non-equilibrium dynamics. An example configuration file for the shock simulations is provided in the appendix (APPENDIX III).
To determine the time length of each shockwave simulation at varying shock impulse, the approach by the Koshiyama group (Koshiyama et al. 2006) was used for a system simulated using periodic boundary conditions. The distance from the right boundary of the shock slab to the opposite boundary of the simulation box was divided by the average velocity added to determine length of simulation.
For instance, if the length from shock initiation to the opposite boundary of the simulation box is 200 Å and the average velocity added is 200 Å/ps, then the shockwave simulation will be carried out for a total of 1 picosecond.

Average Membrane Thickness
The Membrane Thickness tool contained in this software calculates the average membrane thickness over a trajectory by measuring the distance between the two density peaks of user-selected atoms (phosphorous) belonging to head group of POPC as well as the location of the middle point between them defined as the first and second central moment corresponding to the mass density profile of the phosphorous atoms.

% Decrease in Average Membrane Thickness
In addition, the percent decrease in the average membrane thickness was calculated for both models using the formula below Where Th t is the average bilayer thickess at time t and Th t+1 is the average bilayer thickness at time t+1.

Deuterium Order parameter S CD of Lipid Tails
The Deuterium Order Parameter S CD tool measures the deuterium order parameter of acyl chains in phospholipids. This parameter is typically derived through NMR experiments (Vermeer et al. 2007) and represents the orientational mobility of the carbon-hydrogen methylene bonds along the lipid tails of the phospholipid model (Heller et al. 1993). The formula below is used to measure the S CD parameter of lipid tails in this software, S CD = (-3/2(cos 2 Θ CD )-1)

Area per Lipid
The Area per Lipid tool calculates the total average area per lipids for all lipid species as well as the average area per individual lipid species such as POPC or cholesterol using a user-defined selection of atoms. The x and y coordinates of the selected species are projected onto a plane which is delimited by the simulation box. This plane is then subsequently divided into polygons and the area of each polygon is calculated.

Lateral Displacement
The two-dimensional (x and y plane) lateral square displacement LSD is calculated using formula 4.5, where N L is the number of lipids, and r i (t) the position of lipid I at time t.

-Equilibration Analysis
For any good in-silico bilayer model, equilibrium properties must be measured during equilibration to monitor how accurately the properties of these models, such as the area per lipid, match experimental values of similar models. Therefore, the model systems were equilibrated for a total of 7 ns and these equilibration trajectories used to calculate properties such as the average membrane thickness, average area per lipid, and average deuterium order parameter S CD for individual carbons in all lipid tails (palmitoyl + oleoyl) combined or for each lipid tail separately. In addition, the deuterium order parameter was calculated separately for carbons of lipid tails only in either the upper or lower monolayer.

Equilibration Analysis: Average Membrane Thickness
The average bilayer thickness calculated over the 7 ns equilibration trajectory is  Overall, the 50% cholesterol model showed an increase in average bilayer thickness.
Through experiments, it has been shown that the addition of higher levels of cholesterol will result in a thicker lipid bilayer than in models with a lower or no amount of cholesterol. Therefore, it is expected that the model 1:1 model with 50% cholesterol, should have a thicker membrane than the 9:1 model containing only 10% cholesterol.

Equilibration Analysis: Area per Lipid
During equilibration, the average area per lipid for both models was calculated for all lipid species (POPC +CHOL), as well per each individual lipid (POPC or CHL) specie. In figure 5.2, the average area per lipid for all lipid species in the system is computed and compared among the two models. Owing to cholesterol's condensation effect, the system with 1:1 POPC/CHOL ratio shows the lowest total area per lipid for all lipid species in the system.
In addition to the total area per lipid, the area per POPC molecule was calculated throughout the course of equilibration as well. In figure 5.3, the average area per POPC for each model is compared.    (Venable et al. 2010) . This value has also been validated in computational models of POPC-only membranes (Kučerka et al. 2006 (Hung et al. 2007) (Daly et al. 2011) and has been replicated in computer models (Smondyrev & Berkowitz 1999) (Zhang 2009) (Forrest & Sansom 2000) as well.

Equilibration Analysis: S CD (carbon-deuterium order parameter) of Lipid Tails
The deuterium order parameter (S CD ) of lipid tails in the lipid bilayer is an important parameter when it comes to analysis of lipid bilayer models. The deuterium order parameter is a measurement of how "ordered" lipid tails of phospholipids in the lipid bilayer are (Tieleman et al. 1997). It is expected that the addition of cholesterol will lead to a higher average deuterium order parameter of the carbons in the lipid tails of the POPC lipid species in the system due to cholesterol's condensation effect. This condensation effect leads to a tighter packed lipid bilayer and lipid tails which are more ordered "straight" as demonstrated both experimentally and in-silico (Róg et al. 2009). In figure 5.5, average order parameter S CD through equilibration time is calculated for all individual carbons in all lipid tails (oleoyl + palmitoyl).

-Shockwave Simulation Analysis
The implementation of a shockwave impulse using the current formulated method of Koshiyama et al. (Koshiyama et al. 2006a) leads to drastic changes on the lipid bilayer systems tested. Addition of a calculated average velocity to water molecules in the water slab leads to rises in the kinetic energy, temperature and pressure of each system. Therefore, to validate that the shockwave impulse is being added correctly, these variables must be measured and should show an initial increase due to the added velocity added followed by a slow decrease until an equilibrium value is reached for the system tested. This is only the case in an NVE ensemble because temperature and pressure are not controlled in contrast with an NPT ensemble where pressure and temperature are kept constant by thermostats and barostats and therefore the added kinetic energy is dissipated through the thermostat while the added pressure is dissipated through the barostat. The following set of analysis is carried out for both the 1:1 POPC/CHL and 9:1 POPC/CHL models at a range of shock impulses (2.5, 5.0, 7.5, 10 mPa * s). In addition, the utmost right boundary of the water slab was placed 1 Ǻ away from the lipid bilayer in each respective system.

Rise in Kinetic Energy
The addition of an average velocity to water molecules in the chosen water slab adjacent to the bilayer induces a rise in the kinetic energy of the system. With higher applied shockwave impulses, a higher rise in kinetic energy should be observed. The rise in kinetic energy due to shockwave impulses applied to each system is shown separately. Figure 5.8 plots the kinetic energy for the range of shockwave impulses tested on the 1:1 POPC/CHL model. Figure 5.9 represents the rise in kinetic energy for the range of impulses tested on the 9:1 POPC/CHL model.
These plots are simply a validation that the added velocity due to shockwave impulses is being applied correctly.

Rise in Pressure Due to Shockwave
Similar to kinetic energy, the application of an added velocity to water atoms in the system induces a rise in pressure as well. Figures 5.10 and 5.11 plot the respective rise in pressures induced by shockwave impulses for both models respectively. Once again, this is simply another validation that the added velocity is being applied correctly to both systems. With increase in shockwave impulse, an increase in the initial rise in pressure is seen as well when compared to lower impulses. Figure 5.10: Rise in Pressure of 1:1 POPC/CHL system at varying shock impulses. Figure 5.11: Rise in pressure of 9:1 POPC/CHL system at varying shock impulses.

Rise in Temperature Due to Shockwave
Finally, in addition to an increase in kinetic energy and pressure, the temperature in the system is increased rapidly as well due to shockwave impulses. A rise in kinetic energy leads to a rise in the temperature of the systems studied. Figures   5.12 and 5.13 plot the rise of temperature induced by the range of shockwave impulse added to each respective system. Figure 5.12: Rise of temperature of 1:1 POPC/CHL system at varying shock impulses. Figure 5.13: Rise in temperature in 9:1 POPC/CHL system at varying shock impulses.
As the plots show, the temperature increase is proportional to the strength of the shockwave impulse added. After all these validations are made, the analysis on the shockwave simulations and the structural effects can now be made.

Shockwave Simulation Analysis: Change in Bilayer Thickness
As the Koshiyama group discussed in their shockwave simulations performed on their DPPC (1,2-dipalmitoyl-sn-glycero-3-phoshocholine) models, one of the major structural effects due to shockwaves is the collapse and subsequent rebound of the lipid bilayer. This effect can be analyzed by measuring the change in the average bilayer thickness throughout the course of the shockwave simulations. In figure 5.14 and 5.15 respectively, the change in average bilayer thickness induced at varying shockwave impulses is plotted for the 1:1 POPC/CHL and the 9:1 POPC/CHL model.
An initial collapse phase, where the average membrane thickness decreased, is followed by a rebound phase where the membrane rebounds to its original thickness value. These two distinct phases are observed for both models at all shock impulses applied.  In more detail, after the initial collapse phase, the shockwave has propagated through the lipid bilayer and therefore the rebound phase begins and occurs at a range of 500-1000 fs for the respective models.
In addition, the higher the shock impulse applied, the greater the reduction in membrane thickness. Also at higher shockwave impulses, the rebound phase begins earlier because the shockwave has interacted and passed through the lipid bilayer faster than at lower shockwave impulses as expected. Finally, another observation is that the membrane thickness returns close to its normal value faster at higher shock impulses (7.5, 10 mPa * s) then at the lower range of impulses tested (2.5, 5.0 mPa * s).
In addition to analyzing the evolution of bilayer thickness of both models due to shockwaves separately, a comparison of the bilayer thickness changes due to shockwave impulses was carried out as well. To allow for this comparison, the percent decrease in bilayer thickness due to shockwave at similar shock impulses was calculated using equation 4.3 Since the starting membrane thickness in both models is different, calculating percent decrease in bilayer thickness by incorporating the starting membrane thickness values allows for a comparison between models to be made. Once the percent decrease in bilayer thickness was calculated for both systems at the range of shockwave impulses tested, plots were constructed comparing the percent decrease in bilayer thickness of the respective bilayer models at the range of shockwave impulses tested (figures 5.16, 5.17, 5.18 and 5.19) .    When comparing the percent decrease in bilayer thickness at varying shock impulses between models, there is a small difference observed. In general, the 1:1 POPC/CHL model shows a 1-2% less decrease in bilayer thickness corresponding to approximately a 1 Ǻ less reduction in average bilayer thickness when compared to the 9:1 POPC/CHL model. This small difference is not enough to draw any conclusions about the different models and the effect increased cholesterol concentration has on mitigating this effect.

Instantaneous Average Deuterium Order Parameter of Upper and Lower Lipid Monolayers During Shockwave
One of the primary effects of shockwaves is the collapse and rebound of the lipid bilayer caused by a change in the membrane thickness. Membrane thickness is known to be highly reliant on the length of acyl chains contained in the bilayer.   Overall, by calculating the temporal change of the instantaneous averaged order parameter S CDi, the time at which the collapse and rebound phase of the lipid membrane occur can be identified. In addition, it can be said that the decrease and increase of the membrane thickness due to shock impulse was due to an initial shortening of the chain lengths ( greater disorder) followed by an increase in chain lengths (decreased disorder) observed during the rebound phase.

Lateral Diffusion of Lipids During Shockwave
Generally, lipids undergo two types of movements in the lipid bilayer. The first type, tranverse diffusion, involved lipids flip-flopping from monolayer to monolayer and is a generally slow process. The second type of movement, the lateral diffusion, consists of lipids moving laterally (side to side) in the lipid monolayer they reside in.
In contrast to transverse diffusion, the lateral diffusion of lipids occurs rapidly and can  As can be seen in the figures above, the lateral fluidity of lipids during shockwave impulse increases rapidly and is proportional to the strength of the impulse applied.
On average, the 1:1 model shows a slightly higher lateral fluidity of lipids during shockwave impact at all shockwave impulses tested when compared to the 9:1 model.
Generally, the diffusion coefficient is a good measurement of the fluidity of the lipid bilayer. This value is usually derived from the lateral displacement values obtained over a long-time average. Since these shockwave simulations are in the picosecond time length, the calculation of the diffusion coefficient in these simulations would be irrelevant (Almedia & Vaz 1983).

Section 6 -Shock Slab Distance from Bilayer
The results presented above for both models consisted of applying an added velocity to water molecules in a water slab of thickness 5 Ǻ placed adjacent to the lipid bilayer where the far right end of the slab was placed approximately 1 Ǻ away.
One variable not tested by the Koshiyama group and other researchers during in-silico shockwave simulations on lipid bilayers is whether the distance placed between the lipid bilayer and the shockwave water slab is significant. In reality, the implementation of shockwaves using the methods presented in this research rely heavily on the efficiency of momentum transfer between water molecules, specifically water molecules in the water slab and in close proximity to it. Therefore, a similar set of simulations was carried out where the water slab was placed 10 Ǻ away from the bilayer. The velocities applied at both distances for each respective model are very similar. The only variable is the distance of the water slab from the lipid bilayer.
As with the systems tested above where water slabs were placed 1Ǻ away from the lipid bilayer, the evolution of bilayer thickness at varying shockwave impulses is   where the water slab is placed 10 Ǻ away from the lipid bilayer as opposed to 1 Ǻ.
This increased effect on the membrane thickness size can be observed in both models.
To take a closer look, the change in bilayer thickness of the 1:1 POPC/CHL model at 5 mPa * s and 10 mPa *s at both distances of water slab is plotted (figures 5.28, 5.29)  The same analysis for the different water slab distances tested on the 9:1 POPC/CHL system were done and are depicted in figures 5.30 (5 mPa * s) and 5.31 (10 mPa *s).  There are two trends observed in both models when comparing the two different water slab distances tested. First, unexpectedly, there is a greater reduction in average membrane thickness in the systems where the water slab was placed 10 Ǻ away from the lipid bilayer. Expectedly, the maximum decrease in bilayer thickness occurs later in time than the systems where the shock slab is placed 1 Ǻ away. This is expected because a shock placed further away from the bilayer will take longer to reach it.
In an effort to clarify why this increased bilayer thickness reduction is seen in the 10 Ǻ distance simulations, the rise and dissipation of kinetic energy is compared between both systems. The rise in kinetic energy is compared for the 1:1 model at 7.5 and 10 mPa * s impulses and different slab distances (figures 6.7, 6.8). The same comparison was carried out for the 9:1 POPC/CHL model at different water slab distances (figures 6.9, 6.10). Figure 6.7: Rise in kinetic energy of the 1:1 POPC/CHL system due to a 7.5 mPa shock impulse applied at different water slab distances (1 Ǻ vs 10 Ǻ) Figure 6.8: Rise in kinetic energy of the 1:1 POPC/CHL system due to a 10 mPa shock impulse applied at varying water slab distances (1 Ǻ vs 10 Ǻ). Figure 6.9: Rise in kinetic energy of the 9:1 POPC/CHL system due to a 7.5 mPa shock impulse applied at different water slab distances (1 Ǻ vs 10 Ǻ). Figure 6.10: Rise in kinetic energy of the 9:1 POPC/CHL system due to a 10.0 mPa shock impulse applied at different water slab distances (1 Ǻ vs 10 Ǻ).
As the kinetic energy plots comparing distance of water slab show, at a 10 Ǻ distance, the kinetic energy of the system is slightly higher than at a 1 Ǻ distance in both models at both shock impulses compared. This is representative of an increased momentum transfer efficiency that occurs when there is a greater distance between the lipid bilayer and the water slab that contains water molecules with the added average velocities.

Rise in Pressure Due to Shockwave Compared
As noted earlier, the implementation of a shockwave using the methods described, leads to a rise in the pressure of the system. Therefore, a comparison of pressures in both the 1:1 POPC/CHL and 9:1 POPC/CHL models at different shock slab distances is carried out. Figures 5.36-5.39 compare the rise in pressure due to shockwave in both models at 1 Å and 10 Å shock slab distance at shock impulses of 7.5 and 10 mPa * s. Figure 6.11: Rise in pressure of the 1:1 POPC/CHL system due to a 7.5 mPa shock impulse applied at different shock slab distances (1 Ǻ vs 10 Ǻ). Figure 6.12: Rise in pressure of the 1:1 POPC/CHL system due to a 10 mPa * s shock impulse applied at different shock slab distances (1 Ǻ vs 10 Ǻ). Figure 6.13: Rise in pressure of the 9:1 POPC/CHL system due to a 7.5 mPa * s shock impulse applied at different shock slab distances (1 Ǻ vs 10 Ǻ). Figure 6.14: Rise in pressure of the 9:1 POPC/CHL system due to a 10 mPa * s shock impulse applied at different shock slab distances (1 Ǻ vs 10 Ǻ).
In both models, at the 10 Ǻ distance of shock slab from the lipid bilayer, the system pressure is higher throughout simulations at both shock impulses plotted, as can be seen in figures 6.
Therefore, the placement of the shock slab has an impact on the momentum transfer efficiency of water molecules in the slab as well as in between the slab and the lipid bilayer. It can be seen in the results above that, applying similar shock impulses, but increasing the distance between the shock slab and the lipid bilayer has implications on system properties such as the kinetic energy and pressure as well as structural properties of the lipid bilayer such as the change in bilayer thickness.
Regardless, the goal of application of shockwaves to cells and more specifically the lipid bilayer relies heavily on alteration of the lipid membrane structure. By conducting this comparison of the distance of the shock slab relative to the lipid bilayer, this gives more insight to researchers on how distance of shock from target, in this case the lipid bilayer, effects the strength of the shock impulse applied.

Section 7 -Conclusions and Future Work
The goal of this research was to first analyze the effects of shockwaves of varying impulses on the structural effects on in-silico lipid bilayer models containing Upon applying a range of shock impulses (2.5 -10 mPa * s) to the respective bilayer models, the average bilayer thickness change was measured for both models. In addition, the percent decrease in bilayer thickness was calculated to allow a comparison among both models tested. On average, the 1:1 POPC/CHL model showed a lower percent decrease (1-2 %) at all shock impulses than the 9:1 POPC/CHL model. This was not enough of a difference to make any conclusions of how the level of cholesterol in a membrane might mitigate the impact of shock. In addition to bilayer thickness change, the temporal change in the instantaneous averaged deuterium order parameter S CDi due to shock impulse was calculated and plotted for the upper and lower monolayers separately. By doing this, it was shown that monitoring these values during these simulations can be indicative of when the collapse and rebound phase start and terminate. It was also shown that these values go hand in hand with the change in membrane thickness values and that the change in membrane thickness due to shockwave can be linked to the increase and decrease of the instantaneous order parameter of the upper and lower monolayers. Finally, the lateral displacement of both models at varying impulses was measured in order to give a representation of the fluidity of the lipids in these models during various shock impulses.
In addition to the effects of shockwaves, the distance of shock slab placement relative to the lipid bilayer was compared for a 10 Ǻ distance and for a 1 Ǻ distance. On average, at the same shock impulses, the 10 Ǻ water slab system demonstrated greater membrane thickness decrease, a slight increase in kinetic energy of the systems tested, and a slight pressure increase. This can be attributed to the increased momentum transfer efficiency that occurs between water molecules in the slab and near it at the 10 Ǻ distance when compared to the 1 Ǻ. At the 1 In future work, the effect of shock impulses on cancerous-like lipid membranes will be probed. Cancerous cell membranes have a different lipid bilayer composition than mammalian cells. Ideally, by doing this type of research, this is directly applicable to available cancer cell drug therapy methods that rely on permeabilizing the cancerous cell membrane in order to introduce therapeutic drugs into it.