Conformation of ultra-long-chain fatty acid in lipid bilayer: Molecular dynamics study
Kazutomo Kawaguchi, Koh M. Nakagawa, Satoshi Nakagawa, Hideo Shindou, Hidemi Nagao, Hiroshi Noguchi
aa r X i v : . [ phy s i c s . b i o - ph ] S e p Conformation of ultra-long-chain fatty acid in lipid bilayer:Molecular dynamics study
Kazutomo Kawaguchi, ∗ Koh M. Nakagawa, Satoshi Nakagawa, Hideo Shindou,
Hidemi Nagao, and Hiroshi Noguchi Institute of Science and Engineering, Kanazawa University, Kanazawa, Ishikawa 920-1192, Japan Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan Department of Lipid Signaling, National Center for Global Health and Medicine, Shinjuku-ku, Tokyo 162-8655, Japan Department of Lipid Science, Graduate School of Medicine,University of Tokyo, Bunkyo-ku, Tokyo, 113-0033, Japan (Dated: September 28, 2020)Ultra-long-chain fatty acids (ULCFAs) are biosynthesized in the restricted tissues such as retina, testis, andskin. The conformation of a single ULCFA, in which the sn-1 unsaturated chain has 32 carbons, in three typesof tensionless phospholipid bilayers is studied by molecular dynamics simulations. It is found that the ultra-longtail of the ULCFA flips between two leaflets and fluctuates among an elongation into the opposite leaflet, lyingbetween two leaflets, and turning back. As the number ratio of lipids in the opposite leaflet increases, the ratioof the elongated shape linearly decreases in all three cases. Thus, ULCFAs can sense the density differencesbetween the two leaflets and respond to these changes.
I. INTRODUCTION
Cellular membranes consist of various types of lipids,cholesterol, and proteins . The most abundant lipids are phos-pholipids, which have a polar head group and two hydro-carbon tails (fatty acids). From the combinations, diversephospholipids (>1,000 molecular species) are biosynthesizedand have numerous structural and functional roles in cells .Each tail typically contains between 14 and 22 carbon atoms.Around C22 fatty acids such as arachidonic acid (C20:4) anddocosahexaenoic acid (C22:6, DHA) are called very long-chain fatty acids (VLCFAs). Moreover, much longer chainswith 32 to 36 carbons with 6 double bonds were found atthe sn-1 position of phosphatidylcholine (PC) in photorecep-tors, spermatocytes, fibroblasts, and keratinocytes . Theseare called ultra-long-chain fatty acids (ULCFAs). Until now,the biosynthestic mechanisms and biological roles of ULCFAscontaining phospholipids are still unclear.In cells, fatty acids are activated to acyl-CoAs, which areesterified to lysophospholipids to form phospholipids. Us-ing acyl-CoAs as substrates, phospholipids are biosynthesizedby the Kennedy pathway (de novo pathway) and maturedby the Lands’ cycle (remodeling pathway) to generate mem-brane asymmetry and diversity . ULCFA-CoA is elongatedfrom DHA-CoA by elongation of very long-chain fatty acid4 (ELOVL4) and used for the biosynthesis of phospholipid-containing ULCFA as a substrate. ELOVL4 mutations havebeen implicated in Stargardt disease, a type of juvenile macu-lar degeneration ? . Recently, it was reported that ULCFA isstored as a precursor of bioactive lipid mediators. Derivativesof C32:6 and C34:6 are neuroprotective in the retina ? . De-pite these findings, further studies are needed to understandthe importance of these molecules.Computer simulations have been widely used to study lipidmembranes . Among these approaches, all-atom moleculardynamics (MD) simulations are a suitable tool to investigatethe detailed interactions between lipids. They can reproducethe membrane properties well . VLCFAs have been simu-lated by a few groups . They have shown that the long tail is interdigitated into the opposite leaflet . However, themaximum tail length in these studies is C24. Thus, ULCFAshave not yet been simulated. Since ULCFAs have a longer hy-drophobic tail, they are expected to interact strongly with theopposite leaflet. Such an interaction may be relevant to theULCFA function in living cells.In the present study, we examine the conformation of UL-CFA in a fluid bilayer using all-atom MD simulation. AsULCFA, dotriacontahexaenoic acid (C32:6) containing phos-phatidylcholine (dTSPC, C32:6-C18:0) is employed. Themembrane containing the dTSPC is analyzed with othermembrane phospholipids, i.e., distearoyl PC (DSPC, C18:0-C18:0), stearoyl-DHA PC (SDPC, C18:0-C22:6), or stearoyl-oleoyl PC (SOPC, C18:0-C18:1). Moreover, we describethe effects of the difference in lipid density between the twoleaflets. In vesicles, such differences in the lipid density oftwo leaflets can induce membrane bending, as described byan area-difference-elasticity model . We show how thisdensity difference modifies dTSPC conformation. II. MATERIALS AND METHODA. Membrane systems
The molecular structures considered in this study are shownin Fig. 1. DSPC (C18:0-C18:0) contains two saturatedstearoyl chains. SDPC (C18:0-C22:6) contains a steroyl chainand a docosahexaenoyl chain with six double bonds at the sn-1 and sn-2 positions, respectively. SOPC (C18:0-C18:1) con-tains a steroyl chain and an oleoyl chain with a double bondat the sn-1 and sn-2 positions, respectively. dTSPC (C32:6-C18:0) has an ulta-long sn-1 chain of 32 carbons with sixdouble bonds and a steroyl sn-2 chain. In this study, we con-sider a dTSPC molecule inserted into pure DSPC, SDPC, andSOPC membranes to investigate the conformation of the longsn-1 chain of dTSPC. A membrane bilayer consisting of 100lipid molecules per a leaflet was prepared for each membranesystem. The membrane was connected by its periodic images
FIG. 1. Molecular structures of (a) DSPC, (b) SDPC, (c) SOPC, and(d) dTSPC. in the xy plane under the periodic boundary conditions. Onelipid molecule in the upper leaflet was replaced by one dT-SPC molecule for each pure membrane bilayer. In order tostudy the influence of the difference between two leaflets onthe long sn-1 chain of dTSPC, 10 lipid molecules in the up-per (lower) leaflet were removed in the ’189u’ (’189l’) system.Four lipid molecules in the upper (lower) leaflet were removedin the ’195u’ (’195l’) system. Two lipid molecules in the up-per (lower) leaflet were removed in the ’197u’ (’197l’) system.Labels 189, 195, and 197 represent the total number of lipidmolecules except for dTSPC. The number of water moleculesper lipid was fixed at 50 in all cases. Details of the mixedbilayer systems are shown in Table I. The lipids do not flip–flop to the opposite leaflets on a simulation time scale (the TABLE I. Model systems used in the present MD simulations.Model Upper leaflet Lower leaflet Water T / KDS189u 1 dTSPC, 89 DSPCs 100 DSPCs 9,500 343DS189l 1 dTSPC, 99 DSPCs 90 DSPCs 9,500 343DS197u 1 dTSPC, 97 DSPCs 100 DSPCs 9,900 343DS197l 1 dTSPC, 99 DSPCs 98 DSPCs 9,900 343SD189u 1 dTSPC, 89 SDPCs 100 SDPCs 9,500 343SD189l 1 dTSPC, 99 SDPCs 90 SDPCs 9,500 343SD197u 1 dTSPC, 97 SDPCs 100 SDPCs 9,900 343SD197l 1 dTSPC, 99 SDPCs 98 SDPCs 9,900 343SO189u 1 dTSPC, 89 SOPCs 100 SOPCs 9,500 303SO189l 1 dTSPC, 99 SOPCs 90 SOPCs 9,500 303SO195u 1 dTSPC, 95 SOPCs 100 SOPCs 9,800 303SO195l 1 dTSPC, 99 SOPCs 96 SOPCs 9,800 303SO197u 1 dTSPC, 97 SOPCs 100 SOPCs 9,900 303SO197l 1 dTSPC, 99 SOPCs 98 SOPCs 9,900 303 flip–flop time is typically hours or days ). This differencebetween the lipid number of two leaflets results in the devia-tion of the lipid density from a stable value even in a tension-less membrane, as described by an area-difference-elasticitymodel . Thus, the lipids in the leaflet with higher densityare more compressed, although a flat membrane connected bythe periodic boundary does not bend because of the symmetry. B. Molecular dynamics simulations
The CHARMM 36 force field and TIP3P water model were adopted for lipid and water molecules, respectively.The system pressure was controlled at 0.101 MPa by us-ing a semi-isotropic Parrinello-Rahman . The membranewas maintained in a tensionless state. The system tempera-ture in each system is shown in Table I and was controlledusing a Nose-Hoover thermostat . The van der Waals in-teractions were truncated within a radius range of 1 to 1.2nm by using a switching scheme. Electrostatic interactionswere calculated by using the particle mesh Ewald method .The DSPC/dTSPC and SDPC/dTSPC mixtures were equili-brated for 200 ns, which was followed by 800-ns productionruns at 343 K. Since dTSPC is moved more slowly in theSOPC/dTSPC mixtures at 303 K, twice longer time periodswere employed for the SOPC/dTSPC mixtures: 400 ns for theequilibration and 1.6 µ s for production runs. All MD simula-tions were performed by using GROMACS version 2016.4 .Initial configurations were generated using CHARMM-GUIMembrane Builder . Images were visualized using VisualMolecular Dynamics (VMD) software . III. RESULTS AND DISCUSSIONA. Conformation of the sn-1 chain of dTSPC
Figure 2 shows snapshot structures obtained from theSD189u trajectory. Large conformational changes are ob-served in the sn-1 chain of dTSPC. In Fig. 2(a), the sn-1 chain
TABLE II. Transit time of C in dTSPC through the lipid bilayer.Standard deviations are described in parentheses.lipid nsDS197u 5.3 (3.4)DS197l 4.7 (3.6)SD197u 6.4 (4.9)SD197l 5.5 (4.3)SO197u 13.7 (9.0)SO197l 12.0 (7.8) is folded, and the C -C -C angle in sn-1 chain of dTSPCis approximately 0 ◦ ; The terminal carbon atom (C ) is lo-cated in the upper leaflet, and the z coordinate of C is 0.70nm. Here, the z -coordinate is defined as the unit vector paral-lel to the bilayer normal, and the origin is taken at the middleof the lipid bilayer. An L-shaped conformation of the sn-1chain is shown in Fig. 2(b), where the C -C -C angle isapproximately 90 ◦ , and the z coordinate of C is − .
10 nm.The C atom is located between the upper and lower leaflets.In Fig. 2(c), the sn-1 chain exhibits a stretched conformationand C is located in the lower leaflet: the C -C -C angleis approximately 180 ◦ , and the z coordinate of C is − . z -coordinate ofthe C of dTSPC in SD189l and SO189l. The magnitudes offluctuation of C in SD189l and SO189l are similar to eachother. It fluctuates between z = − z = moves across from the membranesurface to another surface. We estimated the transit time ofC to move across the lipid bilayer. Table II shows the aver-age transit time of C to move from an upper boundary to alower boundary and vice versa. The upper and lower bound-aries are defined as z = . − . ∼
10 ns in all threetypes of the lipid bilayers. It is quite fast in comparison withthe flip–flop of the phospholipids between the two leaflets.
B. Distribution of atoms
The distribution of atoms in the sn-1 chain of dTSPC wascalculated along the z -axis. Figure 4 shows the probabilitydistribution of C , C , and C atoms of the sn-1 chain of dT-SPC, phosphate (P) of DSPC, SDPC, and SOPC (see Fig. 1).The distribution of P in DS189u is similar to that in DS189l asshown in Fig. 4(a). Similar distributions are observed in the FIG. 2. Sequential snapshots of the SD189u membrane. The C -C -C angle of dTSPC is around (a) 0 ◦ , (b) 90 ◦ , and (c) 180 ◦ .SDPC molecules are shown in gray, and the gray spheres representthe phosphate atoms of SDPC. The colored spheres represent the dT-SPC molecule. The terminal carbon C32 of the sn-1 chain is depictedin orange. The other carbon atoms and oxygen atoms are depicted incyan and red, respectively. Water molecules are not shown for clarity. cases of SD189u/l (Fig. 4(b)) and SO189u/l (Fig. 4(c)). Theseresults indicate that the membrane thickness is not affectedby the lipid-density difference between two leaflets so that thehost lipid bilayer structure is not significantly modified in theexamined range of the density difference.In contrast, the conformation of the sn-1 chain of dTSPCis strongly changed by the lipid-density difference betweenthe two leaflets. The terminal carbon C of the sn-1 chain -3.0-2.0-1.0 0.0 1.0 2.0 3.0 0.2 0.4 0.6 0.8 1.0 z / n m t / µ s SD189lSO189l
FIG. 3. Time evolution of the z -coordinate of C of dTSPC dur-ing MD simulations. The red and green lines represent SD189l andSO189l, respectively. is widely located from z of − largely shifts toward the lower leaflet, onaverage 0.36, 0.28, and 0.33 nm in DSPC, SDPC, and SOPC,respectively, from 189u to 189l. The middle carbon C islocated almost at the center of the lipid bilayer. The peakposition of C shifts toward the lower leaflet in all cases, asobserved in C , but the shift magnitude is smaller (average of0.22, 0.16, and 0.27 nm in DSPC, SDPC, and SOPC, respec-tively). The distribution of C in the sn-1 chain only slightlyshifts toward the lower leaflet in all cases, as the lipid den-sity relatively decreases in the lower leaflet. Thus, the longerregion (C ∼ C ) of the sn-1 chain exhibits larger changesowing to the lipid-density difference.To quantitatively investigate the effects of the lipid-densitydifferences on the dTSPC conformation, we calculated thenormalized z position as a function of the lipid-density differ-ence as shown in Fig. 5. Positive linear correlations are foundin all cases. From the least-squares fitting, it is found that theterminal of the sn-1 chain moves to the upper leaflet ≃ . N low / N lip . At N low / N lip = . z C / z P is almost 0 for C in all cases. Thus, C is located in themiddle of the lipid bilayer without lipid-density differencesand moves the upper or lower leaflets as the lipid density ofthe upper or lower leaflets relatively decreases, respectively.These conformation changes in dTSPC reduce the lipid den-sity differences. C. Order parameters
Lipid order parameters were calculated to investigate theorientation of the acyl chains. The order parameters S CD aredefined as S CD = (cid:28) α − (cid:29) , (1)where α is the angle between the C-H bond vector and the bi-layer normal. The bracket represents the average over time p ( z ) z / nm(a) dTSPC-sn1-C3 dTSPC-sn1-C18dTSPC-sn1-C32DSPC-P p ( z ) z / nm(b) dTSPC-sn1-C3 dTSPC-sn1-C18dTSPC-sn1-C32SDPC-P p ( z ) z / nm(c) dTSPC-sn1-C3 dTSPC-sn1-C18dTSPC-sn1-C32SOPC-P FIG. 4. Probability distribution of atoms in (a) DS189u/l, (b)SD189u/l, and (c) SO189u/l systems. The solid and dotted lines rep-resent the probability distributions in ’u’ and ’l’, respectively. and lipid molecules. S CD is experimentally measurable bydeuterium nuclear magnetic resonance ? . Figure 6 shows − S CD calculated for SD189u and SO197u.First, we describe the order of the host lipids. The overallorder profile of the sn-1 chain of SDPC in the upper leaflet(blue filled triangles in Fig. 6(a)) is lower than that in thelower leaflet (green open triangles in Fig. 6(a)). For SD189l, ahigher order is obtained in the upper leaflet (data not shown).Hence, the sn-1 chain is more disordered in the leaflets withthe lower lipid density. The order of the sn-2 chain of SDPCis lower than that of the sn-1 chain in both upper and lowerleaflets, since the sn-2 chain has six double bonds. The orderof the sn-2 chain of SDPC in the upper leaflet (orange filledinverse triangles in Fig. 6(a)) is similar to that in the lower -0.3-0.2-0.1 0.0 0.1 0.2 0.3 0.47 0.48 0.49 0.50 0.51 0.52 0.53 DS189l DS197l DS197u DS189u z C / z P N low / N lip (a) C C y = 2.0 x -0.9 y = 3.2 x -1.6 -0.3-0.2-0.1 0.0 0.1 0.2 0.3 0.47 0.48 0.49 0.50 0.51 0.52 0.53 SD189l SD197l SD197u SD189u z C / z P N low / N lip (b) C C y = 1.6 x -0.7 y = 2.8 x -1.4 -0.3-0.2-0.1 0.0 0.1 0.2 0.3 0.47 0.48 0.49 0.50 0.51 0.52 0.53 SO189l SO195lSO197l SO197uSO195u SO189u z C / z P N low / N lip (c) C C y = 2.5 x -1.2 y = 3.4 x -1.7 FIG. 5. Correlation between the position of C and C atomsand the number ratio of lipid molecules in the lower leaflet for (a)DPSC/dTSPC, (b) SDPC/dTSPC, and (c) SOPC/dTSPC mixtures.The horizontal axis represents the number N low of lipid moleculesin the lower leaflet normalized by the total number N lip of lipidmolecules. The vertical axis represents the vertical positions z C ofthe carbons C and C of dTSPC are normalized by the position z P of the phosphate atoms of the host lipids. Error bars are calcu-lated from the standard deviations of eight time-windows in singlesimulation runs. The solid and dashed lines are obtained by the least-squares fitting. leaflet (red open inverse triangles in Fig. 6(a)). As shown inFig. 6(c), for SO197u, the orders of the sn-1 and sn-2 chainsof SOPC are quite similar in the upper and lower leaflets, be-cause the lipid-density difference is small. For SO189u/l, ahigher order is obtained for the higher lipid-density leaflet, asobserved in SD189u (data not shown). The order of the sn-2 chain takes lower value at C and C owing to the double - S C D Carbon(a)
SDPC-sn1 upperSDPC-sn2 upperSDPC-sn1 lowerSDPC-sn2 lower - S C D Carbon(b) dTSPC-sn1dTSPC-sn2 - S C D Carbon(c)
SOPC-sn1 upperSOPC-sn2 upperSOPC-sn1 lowerSOPC-sn2 lower - S C D Carbon(d) dTSPC-sn1dTSPC-sn2
FIG. 6. Order parameter profile − S CD in each leaflet of (a) SDPCin SD189u, (b) dTSPC in SD189u, (c) SOPC in SO197u, and (d)dTSPC in SO197u. bond between C and C .The longer region of the sn-1 chain of dTSPC exhibits loworder in both the SDPC and SOPC membranes, as shown inFigs. 6(b) and (d) (blue filled triangles). A similar dependenceis also observed in the DSPC membranes. The order of the sn-1 chain rapidly decreases to null at C and is maintained untilthe terminal. Thus, this region of the sn-1 chain is randomlyoriented. It corresponds with the conformation of dTSPC inthe snapshots and the distribution of C described above. Theoverall order of the sn-2 chain of dTSPC (orange filled inversetriangles in Figs. 6(b) and (d)) is similar to those of the sn-1chains of SDPC and SOPC. This is reasonable because the sn-2 chain of dTSPC is a stearoyl chain, which is the same as thesn-1 chains of SDPC and SOPC. Thus, the conformation ofthe sn-2 chain of dTSPC is not modified by that of the longersn-1 chain.The main difference among the three types of host lipids isthe number of the double bonds, and clear effects appear in theorder profiles. Nevertheless, the conformation of dTSPC ex-hibits no qualitative differences. A minor influence is found inthe distribution of C of dTSPC, as shown in Fig. 4. Differentshapes are obtained between SDPC and the others, whereasthose of C and C are not. In the SDPC membrane, C hasa rounded triangular distribution. In contrast, a small peak orshoulder shape appears around z ≃ ± IV. CONCLUSION
In this study, we have performed MD simulations for dT-SPC/DSPC, dTSPC/SDPC, and dTSPC/SOPC mixtures to in-vestigate the conformation of dTSPC. The conformation ofthe ultra-long sn-1 chain of dTSPC largely fluctuates in bothleaflets and forms a straight shape deeply interdigitated intothe opposite leaflet, an L-shape bending at the interface be-tween the two leaflets, and a turned shape where the wholechain remains in one leaflet. The ratio of these three states de-pends on the lipid-density difference between the two leaflets.The sn-1 chain is located at the opposite leaflet more fre-quently, as the lipid density of the opposite leaflet relativelydecreases. We have clarified the linear relationships betweenthe position of the sn-1 terminal of dTSPC and the lipid-density difference in all three types of membranes. The timescale of the conformational change between the elongated andturned shapes is ∼
10 ns. Thus, ULCFA can rapidly respondto the lipid-density differences to reduce such differences; thisresponse may be essential for the functions of ULCFAs in liv-ing cells.We investigated a single ULCFA molecule embedded ina fluid membrane consisting of a single type of phospho-lipid. Biomembranes consist of many types of lipids, andthe two leaflets exhibit different lipid compositions. In thepresent study, we found differences in the vertical positions of the sn-1 terminal between the SDPC membrane and theothers. This difference might indicate that ULCFAs prefercontact with a specific type of lipid in multi-component mem-branes. Additionally, ULCFAs might exhibit specific interac-tions with membrane proteins via the conformational changesin the ultra-long chain. Accordingly, further studies are nec-essary to investigate the conformations of ULCFAs in variousmembrane systems.
ACKNOWLEDGMENTS
This research was supported by JSPS KAKENHIJP17K05607 (H.N.), AMED-CREST 20gm0910011 (H.S.),AMED-P-CREATE 20cm0106116 (H.S.), AMED Programfor Basic and Clinical Research on Hepatitis JP20fk0210041(H.S.). MD simulations were carried out by using the facilitiesof the Supercomputer Center, Institute for Solid State Physics,University of Tokyo. We would like to thank CHARMM-GUIdevelopers to add dTSPC into the lipid list.
DATA AVAILABILITY
The data that support the findings of this study are availablefrom the corresponding author upon reasonable request. ∗ [email protected] G. van Meer, D. R. Voelker, and G. W. Feigenson, Nature Rev.Mol. Cell Biol. , 112 (2008). H. Shindou and T. Shimizu, J. Biol. Chem. , 1 (2009). B. Anthony, S. Vanni, H. Shindou, and T. Ferreira, Trends CellBiol. , 427 (2015). J. P. SanGiovanni and E. Y. Chew, Prog. Retin. Eye Res. , 87(2005). A. McMahon and W. Kedzierski, Br. J. Ophthalmol. , 1127(2009). P. Barabas, A. Liu, W. Xing, C.-K. Chen, Z. Tong, C. B. Watt,B. W. Jones, P. S. Bernstein, and D. Križaj, Proc. Natl. Acad. Sci.USA , 5181 (2013). M. Müller, K. Katsov, and M. Schick, Phys. Rep. , 113(2006). M. Venturoli, M. M. Sperotto, M. Kranenburg, and B. Smit, Phys.Rep. , 1 (2006). H. Noguchi, J. Phys. Soc. Jpn. , 041007 (2009). R. M. Venable, F. L. H. Brown, and R. W.Pastor, Chem. Phys.Lipids , 60 (2015). S. J. Marrink, V. Corradi, P. C. T. Souza, H. I. Ingólfsson, D. P.Tieleman, and M. S. Sansom, Chem. Rev. , 6184 (2019). A. P. Ramos, P. Lagë, G. Lamoureux, and M. Lafleur, J. Phys.Chem. B , 6951 (2016). R. Gupta, B. S. Dwadasi, and B. Rai, J. Phys. Chem. B ,12536 (2016). T. Róg, A. Orłowski, A. Llorente, T. Skotland, T. Sylvänne,D. Kauhanen, K. Ekroos, K. Sandvig, and I. Vattulainen,Biochim. Biophys. Acta , 281 (2016). M. Manna, M. Javanainen, H. M.-S. Monne, H.-J. Gabius, T. Rog,and I. Vattulainen, Biochim. Biophys. Acta , 870 (2017). E. Wang and J. B. Klauda, J. Phys. Chem. B , 2757 (2018). U. Seifert, Adv. Phys. , 13 (1997). S. Svetina and B. Žekš, Adv. Colloid Interface Sci. , 189(2014). R. D. Kornberg and H. M. McConnell, Biochemistry , 1111(1971). J. B. Klauda, I. M. Venable, J. A. Freites, J. W. O’Connor, D. J.Tobias, C. Mondragon-Ramirez, I. Vorobyov, A. D. M. Jr., andR. W. Pa, J. Phys. Chem. B , 7830 (2010). A. D. M. Jr. and R. W. Pa, J. Phys. Chem. B , 7830 (2010). M. Parrinello and A. Rahman, J. Appl. Phys , 7182 (1981). W. G. Hoover and B. L. Holian, Phys. Lett. A , 253 (1996). T. Darden, D. York, and L. Pedersen, J. Chem. Phys , 10089(1993). M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith,B. Hess, and E. Lindahl, SoftwareX , 19 (2015). S. Jo, T. Kim, V. G. Iyer, and W. Im, J. Comput. Chem. , 1859(2008). J. Lee, X. Cheng, J. M. Swails, M. S. Yeom, P. K. Eastman, J. A.Lemkul, S. Wei, J. Buckner, J. C. Jeong, Y. Qi, et al. , J. Chem.Theory Comput. , 405 (2015). W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics14