RESEARCH PROPOSAL– BIOMEMBRANES
NAME: VIVEK GOVIND KUMAR (010846928)
TITLE: Elucidation of the role played by sphingomyelin in the oligomerization of lysenin, a ?-pore forming toxin
Specific Aims:
Recent X-ray crystallography and cryo-EM studies have shown that the lysenin pore is actually a nonamer (1, 2). While it is widely accepted that lysenin pore formation is dependent on the specificity of binding to sphingomyelin, a detailed mechanism of pore formation has yet to be published. Other recent studies (3) speculate on the potential roles of sphingomyelin in the pore-formation process, but there is still a lack of clarity. For instance, the sphingomyelin-binding regions of monomeric lysenin have not been thoroughly characterized either experimentally (extensive mutagenesis studies) or computationally (mutation free energy and binding free energy calculations). The primary aim of the proposed research project would be to formulate a well-rounded and definitive hypothesis for the mechanism of sphingomyelin-dependent pore formation in lysenin.

A secondary aim of the proposed project would be to look at the potential usage of lysenin as a tool to visualize the distribution and dynamics of sphingomyelin in cell membranes, under various pathophysiological conditions (4). An extensive characterization of the biochemical properties of lysenin would lead to an improvement in our understanding of pore-forming toxin assembly. It would also provide general insights into the reorganization of lipid bilayers by pore-forming toxins as well as other membrane proteins with similar structure-function relationships.
Research Strategy:
Significance
Pore-forming toxins (PFTs) generally interact with a protein or lipid binding partner to facilitate the recognition of specific target membranes (5, 6). This is followed by an oligomerization process that occurs on the surface of the lipid bilayer and insertion of the oligomeric complex into the target membrane (5, 6). This conversion from a soluble state to a membrane-inserted state requires complex structural and conformational rearrangement of the protein (5, 6). For instance, it may involve the conversion of a soluble monomer with an ?-helical structure to a ?-barrel pore (7, 8, 9). Several studies have demonstrated that the transmembrane region has an amphipathic character, which facilitates simultaneous interactions with the hydrophobic lipid tails of the bilayer and the aqueous pore (8, 10). Our understanding of specific protein-lipid interactions is still somewhat limited, as only a small number of membrane-protein structures have been resolved successfully (3).

Depending on the secondary structural elements that are present in the transmembrane region, PFTs can be classified as either ?- or ?-PFTs (11). The first ?-PFT to have its X-ray crystal structure determined was the soluble form of aerolysin (Aeromonas) (12). Studies have shown that proteins from the aerolysin family tend to form heptameric, octameric or nonameric oligomers, as reported for aerolysin, Dln1 (zebrafish) and monalysin (Pseudomonas entomophila), respectively (13, 14, 15).
Lysenin (from the earthworm Eisenia fetida) is a ?-PFT that specifically targets membranes with sphingomyelin (SM) (1, 2, 3, 4). It is thought to aid the innate immune response by attacking parasitic membranes pore-forming protein that specifically interacts with sphingomyelin (SM) and may confer innate immunity against parasites by attacking their membranes (16). Due to the specificity of its interactions, lysenin is an important biochemical marker for SM (17). Studies that investigate the structure of lysenin bound to SM might potentially throw light on the molecular mechanism of the protein-SM interaction and how it affects pore assembly .The crystal structures of monomeric lysenin were published in 2012 (PDB: 3ZX7, 3ZXD, 3ZXD), with the protein in its apo form and bound to sphingomyelin and phosphocholine (3). At the time, these crystal structures were the first representations of a direct and specific protein-sphingomyelin interaction (3).
Recent x-ray crystallography and cryo-electron microscopy studies have established that the lysenin pore is in fact a nonamer (1, 2) (PDB: 5EC5, 5GAQ).   Monomeric lysenin is a 297 amino acid protein that has two distinct domains – the elongated N-terminal domain (pore-forming module, PFM) which is conserved in all aerolysin ?-PFTs, and the C-terminal ?-trefoil lectin type domain (3). When exposed to membranes with SM, lysenin forms highly stable oligomers that are resistant to high SDS concentrations and temperatures in excess of 95?°C (18). The lysenin nonamer forms a mushroom-like transmembrane pore, with a central stem consisting of an elongated ?-barrel (1, 2).
During structure determination, every protomer was completely defined with the exception of 9 residues at the N-terminus and the ?-hairpin connector loops at the trans side of the barrel (1, 2). The upper half (collar region) of the mushroom cap consists of ?-strands that originate from the N-terminal domains, while the bottom half consists of the C-terminal domains of each protomer (1, 2). No interactions have been observed between ?-barrel transmembrane region and the C-terminal domains (1, 2). The ?-barrel is situated in the center of the oligomeric complex and consists of nine ?-hairpins that are derived from the unwinding of each protomer’s N-terminal domain (1, 2). The ?-barrel spans the entire length of the pore, which is a feature unique to this structure among those currently available for ?-PFTs (1, 2).
A prominent structural feature of aerolysin ?-PFTs is the presence of serine and threonine residues that are potentially involved in membrane binding (19), oligomerization (20) or positioning the amphipathic loops during pore formation (13). The lysenin nonameric structures have revealed that the serine and threonine residues are present along ?-hairpins of the entire tube and not just the transmembrane part (1, 2). Most side chains face the pore lumen, with the exception of some patches on the outside of the ?-barrel where they point into the gap between the barrel and the cap domain (1, 2). As these side chains are polar and relatively short, it is possible that they contribute to the hydrogen bonding network and thus the overall stability of the barrel (1, 2).

The inner surface of the lysenin channel is known to be almost completely negatively charged which explains the pore’s preference for cation transport (21). The shape of the pore in combination with the localization of charged and hydrophobic residues provides an indicator of its orientation in membranes (1, 2). The bottom portion of the cap domain contains positively and negatively charged pockets that could fit the charged phosphocholine (PC) heads of SM (1, 2). Superposition of the soluble monomer (PDB: 3ZX7) with bound POC molecules has revealed that two POC groups can fit into the pockets, with a third one to the side of the cap domain relatively close to the membrane surface (1, 2, 3). This multivalent binding of membrane lipid receptors is quite uncommon in PFTs and potentially explains the affinity of lysenin for SM clusters and the capacity to form oligomers in SM-rich lipid domains (22, 23).

A comparison of the monomeric structure and pore protomers reveals that complex conformational changes occur during pore formation (1, 2, 3). The monomeric N-terminal domain reassembles into the twisted ?-hairpin of the ?-barrel by extending strands ?3 and ?7 through reorganization of a segment comprised of three ?-strands (?4, ?5, ?6) and a single 310 ?-helix (roughly 25% of residues in the monomer) (1, 2). A stack of tilted and twisted ?-hairpins thus forms the final ?-barrel (1, 2). The other portion of the N-terminal domain, a twisted ?-sheet including strands ?2, ?8, ?10 and ?11, remains relatively intact in comparison to the monomer (1, 2, 3). However, the angle between ?11 and the C-terminal domain changes drastically during pore formation, resulting in this portion of the N-terminal domain being tilted by 45° (1, 2). As a consequence of these conformational changes, these ?-sheets form a collar region that connects the transmembrane region to the C-terminal domains in the mushroom cap (1, 2).
The C-terminal domain fold of the pore is completely preserved in comparison to the monomer (1, 2). It is known that lysenin binds specifically to SM through interactions with the C-terminal domain (3, 24). Each C-terminal domain also uses roughly 15% of its overall surface area to interact with C-terminal domains of neighboring protomers (1, 2). Therefore, this clearly indicates that the C-terminal domain plays two key roles. It is the first point of contact with the lipid bilayer and is also an important component of the oligomerization process. (1, 2, 3, 4).

It has been proposed that one of the first post-oligomerization steps is the displacement of the insertion loop (tongue region – Met44 to Gly67) from the rest of the molecule, which would then facilitate the breakage of hydrogen bonds between strands of the insertion loop and the linker region (pairs ?3-?9 and ?5-?11) (1, 2, 4). An important difference between the prepore and functional pore states is that the N-terminal twisted ?-sheet needs to tilt by 45° in order to obtain its final conformational position (1, 2). Further analysis of the lysenin pore structure has also shown that the ?-hairpin precursors are correctly aligned prior to the interaction of all the subunits. This effectively sets up a template for ?-barrel formation (1). While the larger events of the pore formation process are now understood to some extent, there is still a distinct paucity of information on the sphingomyelin binding-initialization step.

Innovation
Dislocation of the tongue region and adjacent parts of the pore-forming module after prepore formation is thought to be important for barrel formation in lysenin (1, 2, 4). The conformational changes seem to involve almost half of the lysenin monomer and involve complete refolding of one-third of the molecule with subsequent refolding of a functional ?–barrel (1, 2). On the other hand, a more typical ?-PFT like ? -hemolysin (Staphylococcus aureus) is thought to oligomerize without any significant structural disturbances to the monomer (25). A potential explanation for this could be lysenin’s unique multivalent SM-binding site on the C-terminal domain. This would control lysenin binding and pore formation solely through interaction with clustered sphingo-myelin in SM-rich domains (22, 23).

Previous studies have resulted in the determination of both monomeric and nonameric lysenin structures (1, 2, 3). Monomeric lysenin has been crystallized bound to the phosphatidylcholine (PC) head groups of SM molecules (3). The PC groups were localized in the C-terminal ?-trefoil domain, in a pocket lined by the residues Lys185, Ser227, Gln229, Tyr233 and Tyr282 (3). In the functional lysenin pore, this pocket would be ideally located for interaction with the SM head groups as it lies on the tip of the receptor-binding domain (2, 3). While, the SM head groups were also bound to residues in the monomeric N-terminal PFM (Lys21, Tyr24, Tyr26, Gln117 and Glu128), these residues would be buried within the cap of the fully-formed lysenin pore (1, 2). Therefore, it is possible that SM recognition at these residues actually occurs prior to oligomerization and membrane insertion (2). A major focus of the proposed project would be functional characterization of these SM binding residues. A review of the available literature has indicated that no studies focusing on the SM binding residues have been conducted. A combination of mutagenesis and pore-formation experiments would help determine the exact role played by each of these residues. For instance, it would be very interesting to study the effects of a K185E or Y233D mutation on sphingomyelin binding and pore formation. What effects would be observed if two or more residues are mutated? The proposed project seeks to address these issues in detail.

Another avenue of investigation would be the dependence of lysenin on the biophysical state of SM, as determined by the presence of other lipids in the membrane. It has been shown that the presence of glycolipids in a membrane causes a drastic reduction in the binding of lysenin to the membrane (18, 22). For example, in polarized MDCK epithelial cells, apical membranes have a high concentration of glycolipids and only basolateral membranes are accessible to lysenin (4). It is possible that glycolipid sugar chains sterically hinder the SM-lysenin interaction (4). Another possibility could be the change in distribution of SM due to preferential localization of both SM and glycolipids in lipid rafts (4).
Studies have also shown that lysenin is capable of targeting SM/DOPC liposomes, while it is unable to bind SM/DPPC liposomes (18). In the membrane, SM mixes well with DPPC but not DOPC, because the gel-to-liquid crystalline phase transition temperature of SM is very similar to that of DPPC and much higher than that of DOPC (4, 18). As a result, SM molecules form clusters in DOPC liposomes and tend to be more dispersed in DPPC liposomes (4, 18). This indicates that lysenin recognizes clustered SM but not dispersed SM (4). Isothermal titration calorimetry experiments have shown that one lysenin typically binds to five SM molecules, thus supporting the idea that lysenin binds to clustered SM (18). Cholesterol is also known to play a role in the oligomerization of lysenin (4). High lateral diffusion along the SM-cholesterol membrane facilitates hexagonal close packing of lysenin oligomers (4).
This project would use high-speed atomic force microscopy (AFM) to investigate the effects of varying membrane lipid composition on the lysenin-SM interaction. Molecular dynamics simulations could also be used to investigate the behavior of wildtype and mutant lysenin (monomers and nonamers) in different lipid environments. Therefore, a combination of experimental and computational techniques could be utilized to elucidate the role played by sphingomyelin in the oligomerization of lysenin. A long term goal would be to use the results of this research to improve our understanding of the properties and mechanisms of assembly of other PFTs.

Approach and methodology
Mutagenesis studies
Every known sphingomyelin binding residue would initially undergo an individual point mutation. The choice of mutant residue would be determined by the nature of the wildtype residue (charged, polar etc.). Two or more residues could also be mutated to investigate the effect of multiple mutations on sphingomyelin binding and pore formation. Genes encoding each mutant could be produced through polymerase chain reaction based mutagenesis (26). All constructs would be assembled in a pT7 expression vector and verified by DNA sequencing of the entire gene insert (27). A potential construct would be an N-terminal fusion with DNA encoding STrEP (II) ligand followed by DNA encoding thioredoxin (1, 3, 27). A tobacco etch virus protease cleavage site would be present between the thioredoxin and lysenin (1, 3, 27).

A culture of E. coli cells harboring the respective plasmids would then be grown in an appropriate broth, and protein expression would be induced following the addition of isopropyl 1-thio-?-d-galactopyranoside to proceed overnight at 18?°C (1, 27). The buffered cell lysate would be loaded on a resin column and the protein would be eluted. The N-terminal tag could be removed overnight at 4?°C using tobacco etch virus protease (1, 27). The lysenin monomers would then be purified further using gel filtration, followed by the removal of the non-cleaved protein as well as cleaved N-terminal tags (1, 27). The monomers could be stored at 4?°C (1, 27).

Lysenin oligomerization and pore formation
The lysenin pore could be assembled from the monomers using a four-step purification process (1). Initially, the monomer would be oligomerized by incubating with liposomes containing SM and cholesterol (1). The resulting precipitate would be collected by centrifugation and solubilized (1). After removing (particulates by centrifugation, the supernatant would be loaded onto a column and fractions containing the lysenin pore could be identified using SDS–PAGE (1). If needed, further purification could be achieved using another column (1). This whole procedure could then be repeated with liposomes containing other SM and other lipids like DOPC, DPPC, POPC etc.

Visualization of oligomerization
Lysenin oligomerization could be visualized using AFM (atomic force microscopy) (1, 4, 28). The effects of each mutation could potentially be characterized by studying the time-lapse AFM images (4, 28). The next logical step would be to characterize the effect of varying the lipid environment. This could be followed by an investigation of the compound effect of mutations and the lipid environment. Would a specific mutation cause a partial loss of sphingomyelin binding activity? If so, would this partial loss of function be amplified by adding glycolipids to the liposome or membrane? High-speed AFM would enable us to view the assembly of lysenin oligomers into a hexagonal close packed structure (28). The effect of phase coexistence on lysenin assembly could also be investigated using AFM (28). It has been proposed that pore formation in SM-rich domains causes phase mixing (4, 28). Without high-speed scanning, it would be exceedingly difficult to determine how lysenin pore formation causes a reorganization of the target membrane (4, 28). As the distribution of lysenin in lipid rafts may vary based on the lipid composition, performing pore assembly assays in live cells would be a potential future project. Lysenin has primarily been characterized in liposomes and model membranes (4, 28). Therefore, cell-based experiments with the correct lipid environment should be a priority for any future studies.

Molecular dynamics studies
Molecular dynamics (MD) simulations should also be performed to test the validity of any hypotheses based on experimental results. 4 different model types would be simulated – wildtype monomer, mutant monomer, wildtype nonamer, mutant nonamer. The lipid environment could also be changed for each simulation, in order to mimic the experimental conditions as closely as possible.
MD simulation trajectories can be analyzed using specialized software to provide an insight into the flexibility and stability of mutant and wildtype systems. Root mean square deviation and fluctuation can be calculated to quantify stability and flexibility respectively. The hydrogen-bonding and salt-bridge interactions within a system can also be analyzed both quantitatively and qualitatively. The number of sphingomyelin (and other lipid) molecules bound to each system could be counted and compared to the experimental setup.

The simulations would be based on the published cryo-EM and crystal structures (PDB: 5EC5. 5GAQ, 3ZX7, 3ZXD, 3ZXD) (1, 2, 3). An example simulation setup would work as follows: The CHARMM-GUI web interface (29,30) would be used to generate all the simulation models required. The simulations would then be performed using the NAMD 2.11 simulation package (31) and the CHARMM36 force field (32, 33). The input files for energy minimization, equilibration steps and production run would be generated using CHARMM-GUI’s Membrane Builder plugin (30). A heterogeneous lipid bilayer could be used, with 800 lipids in the upper leaflet and 800 in the lower leaflet. For instance, the lipid composition could be POPC (17%), POPE (16.5%), POPI (16.5%), cholesterol (25%) and sphingomyelin (25%). The models would then be solvated in a rectangular water box and 150 mM of NaCl ions would be inserted into each system using the Monte-Carlo ion placing method (30).
Initially, energy-minimization would occur for 10,000 steps using conjugate gradient algorithm (34). Then, the systems would be relaxed using restrained MD simulations in a stepwise manner (for a total of ?1 ns) (29, 30). Production runs would be carried out in an NPT ensemble (with constant number of particles, pressure, and temperature) at 310 K using a Langevin integrator with a time step of 2 fs and damping coefficient of 0.5 ps?1. The pressure would be maintained at 1 atm using the Nose??Hoover Langevin piston method (35, 36). The smoothed cutoff distance for nonbonded interactions would be set to 10?12 Å, and long-range electrostatic interactions would be computed with the particle mesh Ewald (PME) method (37). The production run for each model would last roughly 10 microseconds. Analysis of the simulation trajectories would carried out using Visual Molecular Dynamics (VMD) software (38).

Mutation and binding free energy calculations
Alchemical free energy calculations model the physically impossible but computationally realizable process of gradually mutating a subset of atoms of a system from one state to another, through a series of intermediate steps (39). Two alternative methods for alchemical calculation of free energies from molecular dynamics simulation are available in NAMD: Free energy perturbation (FEP) and thermodynamic integration (TI) (39).
Free energy perturbation (FEP) is a statistical mechanics based method that can be used in computational chemistry to calculate free energy differences from molecular dynamics simulations (39). A normal MD simulation can be run for state A, but each time a new configuration is accepted, the energy for state B is also computed (39). The difference between both states could potentially be in the atom types involved, in which case the ? ?G obtained is for “mutating” one molecule onto another, or it may be a geometrical difference where a free energy map can be obtained along the reaction coordinates (39).
If the difference between the two states is not small enough, the free energy perturbation calculations would fail to converge properly (40). Therefore, any given perturbation must be divided into a series of smaller compartments which are computed independently (40). As no communication is required to occur between simulation compartments, the process can be simplified by running each simulation compartment on a different CPU (40).

Mutation free energy calculations could be used to identify those mutations that are energetically unfavorable, thus providing an indication of the functional importance of that wildtype residue. It may even be possible to identify a mutation that could increase sphingomyelin binding affinity. Binding free energy calculations would quantify the direct interactions between wildtype/mutant residues and sphingomyelin.

The proposed project would thus use a mixture of experimental and computational techniques to investigate the role played by sphingomyelin in the oligomerization of lysenin, a ?-pore forming toxin.References:
Podobnik, M., Savory, P., Rojko, N., Kisovec, M., Wood, N., Hambley, R., Anderluh, G. Crystal structure of an invertebrate cytolysin pore reveals unique properties and mechanism of assembly. Nature Communications 2016; 7, 11598.

Bokori-Brown, M., Martin, T. G., Naylor, C. E., Basak, A. K., Titball, R. W., & Savva, C. G. Cryo-EM structure of lysenin pore elucidates membrane insertion by an aerolysin family protein. Nature Communications 2016; 7, 11293.

De Colibus, L., Sonnen, A. F. P., Morris, K. J., Siebert, C. A., Abrusci, P., Plitzko, J., Gilbert, R. J. C. Structures of Lysenin Reveal a Shared Evolutionary Origin for Pore-Forming Proteins And Its Mode of Sphingomyelin Recognition. Structure 2012; 20(9), 1498-1507.
Yilmaz, N., Yamaji-Hasegawa, A., Hullin-Matsuda, F., & Kobayashi, T. Molecular mechanisms of action of sphingomyelin-specific pore-forming toxin, lysenin. Seminars in Cell & Developmental Biology 2018; 73, 188-198.
Anderluh G., Lakey J.H. Disparate proteins use similar architectures to damage membranes. Trends Biochem. Sci. 2008; 33:482–490.

Gilbert R.J. Cholesterol-dependent cytolysins. Adv. Exp. Med. Biol. 2010; 677:56–66.

Gilbert R.J. Inactivation and activity of cholesterol-dependent cytolysins: what structural studies tell us. Structure. 2005; 13:1097–1106.

Shatursky O., Bayles R., Rogers M., Jost B.H., Songer J.G., Tweten R.K. Clostridium perfringens beta-toxin forms potential-dependent, cation-selective channels in lipid bilayers. Infect. Immun. 2000; 68:5546–5551.

Tilley S.J., Orlova E.V., Gilbert R.J., Andrew P.W., Saibil H.R. Structural basis of pore formation by the bacterial toxin pneumolysin. Cell. 2005; 121:247–256.

Song L., Hobaugh M.R., Shustak C., Cheley S., Bayley H., Gouaux J.E. Structure of staphylococcal alpha-hemolysin, a heptameric transmembrane pore. Science. 1996; 274:1859–1866.

Gilbert R. J., Dalla Serra M., Froelich C. J., Wallace M. I. & Anderluh G. Membrane pore formation at protein-lipid interfaces. Trends Biochem. Sci. 2014; 39, 510–516.

Parker M. W. et al. . Structure of the Aeromonas toxin proaerolysin in its water-soluble and membrane-channel states. Nature 1994; 367, 292–295.

Degiacomi M. T. et al. . Molecular assembly of the aerolysin pore reveals a swirling membrane-insertion mechanism. Nat. Chem. Biol. 2013; 9, 623–629.

Jia N. et al. . Structural basis for receptor recognition and pore formation of a zebrafish aerolysin-like protein. EMBO Rep. 2016; 17, 235–248.

Leone P. et al. . X-ray and cryo-electron microscopy structures of monalysin pore-forming toxin reveal multimerization of the pro-form. J. Biol. Chem. 2015; 290, 13191–13201.

Bruhn H., Winkelmann J., Andersen C., Andrä J., Leippe M. Dissection of the mechanisms of cytolytic and antibacterial activity of lysenin, a defence protein of the annelid Eisenia fetida. Dev. Comp. Immunol. 2006; 30:597–606.

Hullin-Matsuda F., Ishitsuka R., Takahashi M., Kobayashi T. Imaging lipid membrane domains with lipid-specific probes. Methods Mol. Biol. 2009; 580:203–220.

Yamaji-Hasegawa A. et al. . Oligomerization and pore formation of a sphingomyelin-specific toxin, lysenin. J. Biol. Chem. 2003; 278, 22762–22770.

Akiba T. et al. . Crystal structure of the parasporin-2 Bacillus thuringiensis toxin that recognizes cancer cells. J. Mol. Biol. 2009; 386, 121–133.

Cole A. R. et al. . Clostridium perfringens epsilon-toxin shows structural similarity to the pore-forming toxin aerolysin. Nature Struct. Mol. Biol. 2004; 11, 797–798.

Kwiatkowska K. et al. . Lysenin-His, a sphingomyelin-recognizing toxin, requires tryptophan 20 for cation-selective channel assembly but not for membrane binding. Mol. Membr. Biol. 2007; 24, 121–134.

Makino A. et al. . Visualization of the heterogeneous membrane distribution of sphingomyelin associated with cytokinesis, cell polarity, and sphingolipidosis. FASEB J 2015; 29, 477–493.

Kulma M. et al. . Sphingomyelin-rich domains are sites of lysenin oligomerization: implications for raft studies. Biochim. Biophys. Acta 2010; 1798, 471–481.

Yamaji A. et al. . Lysenin, a novel sphingomyelin-specific binding protein. J. Biol. Chem. 1998; 273, 5300–5306.

Sugawara T. et al. . Structural basis for pore-forming mechanism of staphylococcal alpha-hemolysin. Toxicon 2015; 108, 226–231.

Jones D. H. in: PCR Primer: A Laboratory Manual eds Dieffenbach C. W., Dveksler G. S. 1995; 591–601, Cold Spring Harbor Laboratory.

Cheley S. et al. . Spontaneous oligomerization of a staphylococcal alpha-hemolysin conformationally constrained by removal of residues that form the transmembrane beta-barrel. Protein Eng. 1997; 10, 1433–1443.

Yilmaz N. & Kobayashi T. Visualization of lipid membrane reorganization induced by a pore-forming toxin using high-speed atomic force microscopy. ACS Nano 2015; 9, 7960–7967.

Jo, S., Kim, T., Im, W. Automated builder and database of protein/membrane complexes for molecular dynamics simulations. PLoS One 2007; 2, e880.

Lee, J., Cheng, X., Swails, J. M., Yeom, M. S., Eastman, P. K., Lemkul, J. A., Wei, S., Buckner, J., Jeong, J. C., Qi, Y., et al. CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. J. Chem. Theory Comput. 2016; 12, 405?413.

Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., Chipot, C., Skeel, R. D., Kale, L., Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005; 26, 1781?1802.

Klauda, J. B., Venable, R. M., Freites, J. A., O’Connor, J. W., Tobias, D. J., Mondragon-Ramirez, C., Vorobyov, I., MacKerell, A. D., Jr., Pastor, R. W. Update of the CHARMM all-atom additive force field for lipids: Validation on six lipid types. J. Phys. Chem. B 2010; 114, 7830?7843.

Best, R. B., Zhu, X., Shim, J., Lopes, P. E. M., Mittal, J.; Feig, M., MacKerell, A. D. Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone ?, ? and side-chain ?1 and ?2 dihedral dngles. J. Chem. Theory Comput. 2012; 8, 3257?3273.

Darden, T., York, D., Pedersen, L. G. Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993; 98, 10089?10092.

Reid, J. K. In Large Sparse Sets of Linear Equations; Reid, J. K., Ed.; Academic Press: London 1971; 231?254.

Martyna, G. J., Tobias, D. J., Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994; 101, 4177?4189.

Feller, S. E., Zhang, Y., Pastor, R. W., Brooks, B. R. Constant pressure molecular dynamics simulation: The Langevin piston method. J. Chem. Phys. 1995; 103, 4613?4621.

Humphrey, W., Dalke, A., Schulten, K. VMD ? visual molecular dynamics. J. Mol. Graphics 1996; 14, 33?38.

Chipot, C.; Pearlman, D. A., Free energy calculations. The long and winding gilded road, Mol. Sim. 2002; 28, 1–12.

Chipot, C.; Pohorille, A., Eds., Free energy calculations. Theory and applications in chemistry and biology, Springer Verlag, 2007.