LiPyphilic: A Python toolkit for the analysis of lipid membrane simulations

Molecular dynamics simulations are now widely used to study emergent phenomena in lipid membranes with complex compositions. Here, we present LiPyphilic - a fast, fully tested, and easy to install Python package for analysing such simulations. Analysis tools in LiPyphilic include the identification of cholesterol flip-flop events, the classification of local lipid environments, and the degree of interleaflet registration. LiPyphilic is both force field and resolution agnostic, and thanks to the powerful atom selection language of MDAnalysis it can handle membranes with highly complex compositions. LiPyphilic also offers two on-the-fly trajectory transformations to i) fix membranes split across periodic boundaries and ii) perform nojump coordinate unwrapping. Our implementation of nojump unwrapping accounts for fluctuations in box volume under the NPT ensemble — an issue that most current implementations have overlooked. The full documentation of LiPyphilic, including installation instructions, is available at https://lipyphilic.readthedocs.io/en/latest. Graphical TOC Entry


Introduction
The plasma membrane was once thought to be a passive divider between a cell and its external environment. We now understand that it is in fact a dynamic interface upon which many cellular processes, from cell-signaling to membrane transport, depend. These processes are emergent phenomena that arise from a complex interplay between the molecular species that comprise the plasma membrane. As such, there is a great interest in understanding how the lipids, proteins, and carbohydrates of the plasma membrane interact with one another.
Molecular dynamics (MD) simulations are routinely used to study lipid-lipid and lipidprotein interactions at a molecular level, and there exists many excellent tools for analysing the trajectories of such simulations. FATSLiM 1 and MemSurfer, 2 for example, both specialise in the analysis of non-planar membranes such as buckled bilayers or vesicles. PyLipID 3 and ProLint 4 are designed for the easy and efficient analysis of lipid-protein interactions.
ML-LPA is a recently developed Python package that employs various machine learning algorithms to identify the phase -L o or L d -of lipids in a bilayer. 5 LOOS, on the other hand, is a C++ library with a Python interface for analysing MD simulations. 6,7 Unlike the above packages, LOOS handles the trajectory reading internally whilst also offering a large set of analysis tools, some of which are for lipid membranes. Between them, these software packages provide an extensive analysis suite for MD simulations of lipid membranes.

LiPyphilic
LiPyphilic is an object-oriented Python package for analysing MD simulations of lipid membranes. It is built directly on top of MDAnalysis, and makes use of NumPy 43 and SciPy 44 for efficient computation. It is both resolution and force field agnostic, working with any file format that MDAnalysis can load so long as the topology contains residue names. All analysis classes in LiPyphilic share the same application programming interface (API) as those in MDAnalysis. This shared API makes it simple for users of MDAnalysis to learn how to use LiPyphilic.
At its core, LiPyphilic is designed to easily integrate with the wider scientific Python stack. Results are typically stored in a two dimensional NumPy array of shape N lipids by N f rames (Figure 1), making it simple to post-process the results for further analysis. Some analysis tools also take a two dimensional NumPy array of the same shape as input. This input array may contain information about each lipid, such as which leaflet they belong to or their phase (L o or L d ). Alternatively, the input array may be a boolean mask -an array of True and False values specifying which lipids to include in the analysis. As these inputs are generic NumPy arrays instead of types specific to LiPyphilic, it is possible to use the output from other membrane analysis tools as input to LiPyphilic. For example, you may assign lipids to leaflets using FATSLiM, 1 determine their phase state using ML-LPA, 5 or calculate local membrane normals using MemSurfer, 2 and extract the results to perform further analysis with LiPyphilic. The workflow for using LiPyphilic generally involves the following steps: 1. import MDAnalysis along with the required LiPyphilic analysis modules 2. load a topology and trajectory as an MDAnalysis universe 3. create an analysis object using the MDAnalysis universe and specifying the relevant input options 4. use the run() method to perform the analysis 5. store the results, either by serialising the analysis object itself or saving the results data as a NumPy array Below we will discuss the implementation and usage of some of the analysis tools currently available in LiPyphilic. We will then discuss the on-the-fly transformations that LiPyphilic can perform on MDAnalysis trajectories. We then provide benchmarks of the analysis tools and transformations. Finally, we will briefly discuss the software engineering best practices used in developing LiPyphilic. If you are more interested in learning how to use LiPyphilic, rather than how LiPyphilic works per se, we recommend working through the examples in the documentation at https://lipyphilic.readthedocs.io/en/latest/reference/ analyses.html.

Assign leaflets
For many analyses, such as calculating the area per lipid, it is necessary to know the leaflet within which a lipid is found. LiPyphilic has two tools for assigning lipids to lealfets. The class lipyphilic.lib.assign_leaflets.AssignLeaflets assigns each lipid to a leaflet based on the distance in z to its local membrane midpoint. This is suitable only for planar bilayers. On the other hand, the class lipyphilic.lib.assign_leaflets.AssignCurvedLeaflets can be used to identify leaflets in a buckled bilayer or a micelle. This uses the MDAnalysis leaflet finder 45,46 to assign non-translocating lipids to leaflets, then at each frame assigns the remaining lipids based on their minimum distance to each leaflet. AssignLeaflets remains useful for planar bilayers, especially if the rate of cholesterol translocation is of interest, which is typically measured by assigning lipids to leaflets based on their z-coordinate.
LiPyphilic can assign molecules not just to the upper or lower leaflet, but also to the midplane. This is useful for studying, for example, the local lipid environment of midplane cholesterol 24 or its role in the registration of nanodomains. 47 Assigning cholesterol to the midplane also creates a buffer zone for determining whether a flip-flop event was successful (i.e. crossed the buffer zone) or not. 25 AssignLeaflets and AssignCurvedLeaflets have opposing approaches to assigning molecules to the midplane. The former considers a molecule's distance to the midplane whereas the latter considers its distance to each leaflet. There is naturally an inverse relationship between these two measures -the further a molecule is from the midplane the closer it is to a leaflet. However, the distance to the midplane is typically employed when studying flip-flop in planar bilayers, 17,[21][22][23][24][25]31 whereas the distance to each leaflet is used for studying flip-flop in undulating bilayers. 47 As with all analysis tools in LiPyphilic, the assigning of lipids to leaflets is both resolution and force field agnostic. Instead of reading atom selections hard coded into the package, the analysis tools rely on the powerful selection language of MDAnalysis. Figure 2A shows how AssignLeaflets may be used to determine the leaflet membership of all lipids in the 58 component neuronal plasma membrane studied by Ingólfsson et alii. 31 First, an MDAnalysis Universe must be created. Lipids are then assigned to leaflets by passing this Universe to AssignLeaflets along with an atom selection of lipids in the bilayer, using the universe and lipid_sel arguments respectively. Optionally, to allow molecules to be in the midplane, we can use the midplane_sel and midplane_cutoff arguments. In the example shown in Figure   2A, cholesterol will be assigned to the midplane if its ROH bead is within 8 Å of its local midpoints. Local midpoints are computed by first splitting the membrane into an n by n grid in xy, where n is specified using the n_bins argument. The local midpoint of a grid cell is then given by the center of mass of all atoms selected by lipid_sel that are in the grid cell. Through calculating local membrane midpoints, this algorithm can account for small undulations in a bilayer. However, for bilayers with large undulations or for non-bilayer membranes, AssignCurvedLeaflets should be used. After creating leaflets as described above, the analysis is performed by calling the run method. Here the start, stop and step arguments are used to specify which frames of the trajectories to use, and a progress bar can be displayed on the screen by setting verbose=True.
Leaflet data are then stored in the leaflets.leaflets attribute as a two-dimensional NumPy array. Each row in the results array corresponds to an individual lipid and each column to an individual frame. For example, leaflets.leaflets [i, j] contains the leaflet membership of lipid i at frame j. leaflets.leaflets [i, j] is equal to 1 if the lipid is in the upper leaflet, −1 if the lipid is in the lower leaflet, or 0 if the lipid is in the midplane.

Flip flop
Cholesterol is unevenly distributed across the plasma membrane, although the precise distribution is still under debate. 48 Figure 3A illustrates how to do so using the output from AssignLeaflets.
The same MDAnalysis Universe that is used for assigning leaflets is passed to FlipFlop. An atom selection that specifies which molecules to consider when identifying flip-flop events is passed to the lipid_sel argument. The leaflet membership of each lipid selected by lipid_sel is passed to the leaflets argument. In this example, this is achieved by filtering the results array of AssignLeaflets to include only the leaflet membership of cholesterol molecules. The leaflet membership must be a NumPy array, of shape ((N lipids , N f rames )), in which each element is equal to:  In this example, the frame_cutoff argument is used to specify that a molecule must remain in its new leaflet for at least two consecutive frames in order for the flip-flop to be considered successful.
We again call the run method to perform the analysis. For each molecule, LiPyphilic will then identify the frames at which it leaves one leaflet and enters another for at least frame_cutoff frames. If the new leaflet is different to the original leaflet, the flip-flop was successful.
The success or failure of each flip-flop event is stored in a one-dimensional NumPy array, accessible via the flip_flop.flip_flop_success attribute. Elements in this array are strings, equal to either "Success" or "Failure". From this results array, it is easy to calculate the flip-flop rate based on the number of observed events, the number of cholesterol molecules in the membrane, and the total simulation time ( Figure 3B).
The example in Figure 3A shows how to use the results from AssignLeaflets to identify flip-flop events. However, leaflets need not be identified using LiPyphilic in order to use the FlipFlop analysis tool. FlipFlop expects a NumPy array of leaflet membership as described above. Flip-flop events can thus be found even if leaflets were assigned using, for example,

Registration
The translocation of cholesterol across leaflets is thought to be important in several cellular process, including the modulation of the lateral heterogeneity of the membrane. 47 Recently, transient nanodomains of L O phase lipids were observed in live mammalian cell plasma membranes. 50 These nanodomains are thought to be enriched in sphingomyelin and cholesterol, and to act as functional platforms for cell signaling. However, their nature, formation, and roles in cellular processes are still not fully understood.
There is particular interest in understanding under what conditions nanodomains in apposing leaflets are spatially aligned. Such alignment is known as interleaflet registration.
This has been the subject of several MD simulations, which have revealed registration to be a complex process. 11,[34][35][36][37][38][39][40][41][42] Registration is modulated by many factors, including the length and saturation of lipid tails as well as the relative affinity of cholesterol for the different lipid species in a domain-forming mixture.
The class lipyphilic.lib.registration.Registration can be used to quantify the degree of interleaflet registration in a planar bilayer. Registration is an implementation of the registration analysis described by Thallmair et alii. 40 The degree of registration is calculated as the Pearson correlation coefficient of molecular densities in the upper and lower leaflets.
First, the two-dimensional density of each leaflet is calculated: where the (x, y) positions of lipid atoms in leaflet L are binned into two-dimensional histograms with bin lengths of 1 Å. L is either the upper (u) or lower (l) leaflet. The twodimensional density is then convolved with a circular Gaussian density of standard deviation σ. The registration between the two leaflets, r u/l , is then calculated as the Pearson correlation coefficient between ρ(x, y) u and ρ(x, y) l . Values of r u/l = 1 correspond to perfectly registered domains and values of r u/l = −1 correspond to perfectly anti-registered domains. is useful to know the degree of registration of L o domains, rather than the registration of a specific molecular species. Registration can be used to calculate this, providing the phase of each lipid at each frame has been determined using for example, ML-LPA, 5 a hidden Markov Model, 27,33 or the deuterium order parameter of lipid tails. If the lipid phase data is stored in a two-dimensional NumPy array of shape (N lipids , N f rames ), it can be used to create a boolean mask that will tell Registration which lipids to include in the analysis. For example, if our array is named lipid_phase_data, and its elements are strings of either "Lo" or "Ld", then we can select the L o lipids for analysis by passing the boolean mask lipid_pahse_data == "Lo" to the filter_by argument of Registration.

Neighbour matrix
The plasma membrane is comprised of hundreds of different lipid species. In this complex mixture, lateral heterogeneities and aggregates of specific lipid species arise spontaneously.
Over the past decade, this compositional complexity has begun to feature in MD simulations of membranes. 17,25,31,[51][52][53][54][55][56] In these simulations, the lateral organization of the membrane is typically quantified via a consideration of local lipid environments. Specifically, the lipid enrichment index of species B around species A, E AB , may be defined as: 17 This is because the lil_matrix is efficient at adding new neighbours on-the-fly, whereas the csc_matrix is efficient for slicing columns of the matrix. The latter is used frequently in the largest_cluster and count_neighbours helper methods of Neighbours.

Largest cluster
Some glycolipids in the plasma membrane are known to aggregate, forming platforms for cell-signalling. [57][58][59] The size of the largest cluster of glycolipids in the neuronal plasma membrane 31 can be calculated using the largest_cluster method of Neighbours ( Figure 5B). For this, an atom selection must be provided to the cluster_sel argument. In the example in Figure 5B, it is specified that only lipids in the upper leaflet should be included in the cal-culation by passing a boolean mask to the filter_by keyword. The return_indices argument is used to specify that the residue indices of the lipid molecules in the largest cluster at each frame are also to be returned. There is no need to call a run method nor to specify which frames of the analysis to use -the same frames specified in 5A will be used for the cluster analysis.
The To find the largest cluster at a given frame, the neighbours.neighbours sparse adjacency matrix is first sliced to give a matrix for the current frame only, A f rame . The connected_components function of SciPy is then used to find all connected components at the current frame. NumPy's unique function, with return_counts set to True, is then used to identify the largest connected component and thus the largest cluster size.

Enrichment index
After constructing the adjacency matrix, the count_neighbours method can be used to determine the local environment of each lipid molecule ( Figure 5C). For a single lipid, its local lipid environment is defined as the number of neighbours of each species. In the example in Figure 5C, return_enrichment=True is set to specify that the lipid enrichment index is also to be returned.
As with neighbours, the results are not stored in an attribute in our neighbours object.
The neighbour counts and the enrichment index are each returned as a Pandas DataFrame.

User-defined counts
By default, the count_neighbours method will calculate the number of neighbouring species around each individual lipid. This is done using the residue name of each lipid. However, it is also possible to use any ordinal or string data for counting lipid neighbours. For example, the enrichment index of lipids in the neuronal plasma membrane can be calculated based on their tail saturation ( Figure 6). For this, first a two-dimensional NumPy array of shape (N lipids , N f rames ) that contains the saturation of each lipid needs to be created ( Figure S1).
Then this array is passed to the count_by argument of count_neighbours, and the local lipid environment and enrichment index will be determined based on the information in this array.

On-the-fly transformations
MDAnalysis has a powerful set of on-the-fly trajectory transformations. These transformations can do away with the need to create multiple instances of the same trajectory using, for example, the GROMACS trjconv tool. Instead, the transformations are applied each time a frame is loaded into memory by MDAnalysis. LiPyphilic extends the set of transformations available in MDAnalysis to include the ability to repair membranes split across periodic boundaries and to perform nojump trajectory unwrapping. The latter prevents an atom from being wrapped into the primary unit cell when it crosses a periodic boundary.

Center membranes
The callable class lipyphilic.transformations.center_membrane can be used to center a membrane -or any supramolecular structure -in a box, providing it is not self-interacting across periodic boundaries. For each frame, all atoms in the system are iteratively shifted along a specified set of dimensions until the membrane is no longer split across the periodic boundaries ( Figure 7A). After each translation, all atoms are wrapped back into the primary unit cell. For example, to check if a bilayer is split across the periodic boundary in the z-dimension, its extent in z, H, could be compared to the box length in z, L z . If H is within a user-specified cutoff value of L z , the bilayer is split across z and is thus translated in this dimension. Once the bilayer is no longer split across boundaries, it is then moved to the center of the box in z.
This transformation can be applied to an MDAnalysis universe, u, using the u.trajectory .add_transformation method ( Figure 7A). This method takes as input a transformation, or a list of transformations. The center_membrane callable class is passed to add_transformation along with the arguments for center_membrane. The atoms that comprise the membrane are specified using the ag argument. The atom selection should include all atoms in the membrane, not a subset of atoms -otherwise the extent of the membrane cannot be calculated accurately. In the example in Figure 7A, the bilayer is centered in only the z-dimension; however, the membrane can be made whole in each dimension independently. This is controlled using the center_x, center_y and center_z arguments. The shift argument is used to specify the distance in Å that the membrane will be translated at each iteration. Too small a value of shift would require many iterations to make a membrane whole. Too large a value, on the other hand, may result in a membrane being translated nearly the length of the unit cell and thus remaining broken. We have found a tranlsation of 10 Å to be suitable for bilayers, but the optimal value will depend on the membrane structure and the size of the system.

NoJump trajectory unwrapping
lipyphilic.transformations.nojump can be used to prevent atoms from jumping across periodic boundaries. It is analogous, but not equivalent, to using the GROMACS command trjconv with the flag -pbc nojump. This transformation can be applied to an MDAnalysis universe in much the same way as lipyphilic.transformations.center_membrane. We must pass an atom selection to the ag argument of nojump, and specify the dimensions to which the transformation should be applied (Figure 8). Upon adding this transformation to your trajectory, nojump will perform an initial pass over the trajectory. It will determine the frames at which each atom crosses a boundary, keeping a record of the net movement of each atom at each boundary. This net movement across each boundary is used to determine the distance an atom must be translated in order to be moved from its wrapped position to its unwrapped position. Subsequently, every time a new frame is loaded into memory by MDAnalysis, such as when iterating over the trajectory, the relevant translation is applied to each atom to move it to its unwrapped coordinates.
Below we describe the nojump unwrapping algorithm implemented in LiPyphilic. We also explain how this algorithm avoids the artefacts introduced by the standard unwrapping scheme that, to our knowledge, is employed by all molecular dynamics simulation related software. Specifically, the standard unwrapping scheme fails to account for fluctuating systems sizes caused by barostats in NPT ensemble simulations.

Unwrapping scheme
The unwrapped position of a particle at frame N , denoted x u N , is given by: L n c n where x w N is the wrapped position of the particle at frame N , and N n=0 L n c n accounts for the displacement that results from all jumps across periodic boundaries from frame 0 to frame N . L n is the box length at frame n, with the box centered at L/2. This box length is multiplied by a factor c n , which depends on whether an atom crossed a periodic boundary from frame n − 1 to frame n: where x w n−1 is the particle's wrapped position at frame n − 1. At frame n = 0, for which there is no previous position, we take x w −1 to be the particles raw atomic coordinate at frame n = 0. That is, if the particle is not in the primary unit cell at frame n = 0, we calculate the displacement required to move from x w 0 to x u 0 . This unwrapping scheme is equivalent to that recently described by von Bulow et al., 60 although it was derived independently. Both our unwrapping algorithm and that of von Bulow et al. avoid the problems that the standard unwrapping scheme suffers. To calculate x u N , the standard unwrapping scheme iteratively adds the box length at frame N to the wrapped coordinated at frame N until |x w n − x u n−1 | < L n /2. However, in the example in Figure 9, this would result in x u 2 = x w 2 − L 2 = −0.5 instead of the correct value x u 2 = −1.5. von Bulow et al. demonstrated clearly the effect that this inaccurate unwrapping of atomic coordinates has on the calculated diffusion coefficient. 60 The unwrapping scheme described above, and previously by von Bulow et al., 2 correctly accounts for the fluctuating box size in the NPT ensemble. However, it is only accurate in the case where coordinates are stored every timestep. In fact, it is impossible to correctly unwrap coordinates unless we store them at every timestep. This logically follows from the same argument made above -to correctly unwrap coordinates we must know the length of the box at the timestep at which the jump occurred. See the Supporting Information for further details.

Benchmarking
For benchmarking the LiPyphilic analysis tools, we have performed a simulation of an equimolar mixture of DPPC/DOPC/Cholesterol using the MARTINI 2 force field. 62, 63 We created a symmetric bilayer of 12,000 lipids in total using the CHARMM-GUI MARTINI Maker. 64 The production run was performed for 8.0 µs and coordinates were stored every 5.0 ns, giving a trajectory of 1,600 frames. Where analogous analysis tools are available in either FATSLiM 1 0.2.2 or GROMACS 65 2020.4, we compare their performance with that of LiPyphilic. We benchmark against FATSLiM as this is generally the fastest membrane analysis tool available. 1 The FATSLiM benchmarks were performed using eight OpenMP threads. All other benchmarks were performed in serial. All of the benchmarks can be seen in Table 1.
LiPyphilic is generally very fast, with most analysis tools and trajectory transformations taking on the order of 10 ms per frame for the 12,000 lipid membrane. For example, both methods of assigning leaflets are faster than the corresponding implementation in FATSLiM.
There are two further benefits of assigning leaflets with LiPyphilic: i) molecules may reside in the midplane and ii) the results are stored in a single NumPy array whereas FATSLiM creates a GROMACS index file for each frame of the trajectory. LiPyphilic is significantly faster in calculating the bilayer thickness, although the imeplmetation in FATSLiM is more sophisticated. The FATSLiM thickness tool is able to handle both bilayers and non-bilayer membranes, whereas the one in LiPyphilic is only appropriate for planar bilayers. LiPyphilic is also faster than GROMACS when it comes to calculating the coarse-grained order parameter, S CC . Using LiPyphilic to calculate S CC is also simpler than using GROMACS, which requires the creation of a separate index group for each unique atom along a tail of each lipid species.
The calculation of interleaflet registration, construction of a neighbour matrix, and calculation of the lipid enrichment index are slower, on the order of 100 ms per frame. This is still relatively fast -although they are different analyses, these times are comparable to The slowest analysis in LiPyphilic is the calculation of the area per lipid. This takes over 1.5 s per frame, and is 7.6 times slower than FATSLiM. We therefore recommend using FATSLiM for the area per lipid calculation if you have a large membrane (>1,000 lipids). It should be remembered, however, that the FATSLiM analyses were run in parallel over eight cores whereas all other benchmarks were performed in serial. Parallelising the analysis in LiPyphilic could result in a similar performance boost. Whilst this is out of the scope of the package in its current state, we do have plans to parallelise the slower tools in LiPyphilic in due course. In this respect, we are particularly interested in the development of Parallel MDAnalysis (PMDA). 66 PMDA is based on MDAnalysis and uses Dask to parallelise the analysis modules. 67 The analysis classes inherit a modified abstract base class that is specifically designed to make the parallelisation straightforward. We will wait until PMDA is out of the alpha stage of development before assessing which modules would benefit from being parallelised.
Finally, the on-the-fly transformations are fast, with center_membrane and nojump taking 23.37 ms and 23.89 ms per frame, respectively. Upon applying the nojump transformation, a first pass over the trajectory is performed to calculate the translations that need to be applied at each frame. This first pass takes 13.98 ms per frame for the 12,000 atoms selected in the benchmark. After this first pass, the time to load a frame into memory and apply the translations is 9.91 ms. This performance is put into perspective by considering that iterating over the trajectory with MDAnalysis takes 8.54 ms per frame itself, with no transformations applied. The total time per frame of nojump (23.89 ms) is approximately twice that required by the GROMACS trjconv tool to do the same transformation. However, using nojump has two benefits over GROMACS trjconv in that it i) prevents the need to create duplicate trajectories and ii) accounts for box size fluctuations caused by barostats.

Software Engineering
LiPyphilic is free, open-source software licensed under the GNU General Public License v2 or later. In developing LiPyphilic, we have followed the best practices in modern scientific software engineering. 68 We use version control, unit testing and continuous integration, and We encourage users to submit feature requests and bug reports via GitHub, and we welcome any question about usage on our discussion page at https://github.com/p-j-smith/ lipyphilic/discussions.
Unit testing in LiPyphilic is performed using Pytest. 70 We have constructed a set of toy systems for testing each analysis tool. These systems are typically composed of two sets of atoms each arranged on a hexagonal lattice in xy, with the two lattices separated vertically in z. This setup approximates the topology of the headgroups of a lipid bilayer. Using these toy systems, we know what the results of each analysis tool should be a priori. We can thus test each analysis tool with full confidence in the results if the tests pass, without relying on regression tests that involve highly complex systems. Further, using Pytest-cov, 71 we have ensured that all analysis tools and trajectory transformations have 100% test coverage.
Finally, LiPyphilic is simple to install. We have packaged LiPyphilic to make it available for installation using widely-used package managers. The easiest way to install LiPyphilic along with all of its dependencies is through Anaconda. 72 Alternatively, it can be installed via the Python Package Index using Pip. 73

Conclusion
We have developed a new Python package for the analysis of lipid membrane simulations. We have focused on providing functionality not available in other membrane analysis tools, such as calculating the lipid enrichment/depletion index, the degree of interleaflet registration, and the flip-flop rate of molecules between leaflets. LiPyphilic is a modular object-oriented package that makes extensive use of NumPy, 43 SciPy, 44 and MDAnalysis 45,46 for efficient computation. For analysing a 12,000 lipid MARTINI membrane, the analysis classes typically take on the order of 10-100 ms per frame. This is comparable to the performance of analysis tools in GROMACS 65 and FATSLiM 1 -a very fast package for membrane analysis. All analysis tools in LiPyphilic share the same API as those of MDAnalysis. This shared API makes LiPyphilic simple to learn for current users of MDAnalysis.
The modularity of LiPyphilic, along with its focus on integrating with the wider scientific Python stack, means the output of other analysis tools such as FATSLiM 1 or ML-LPA 5 can be used as input for further analysis in LiPyphilic. Further, the output of LiPyphilic is in the form of NumPy arrays, Scipy sparse matrices, or Pandas Dataframes. This means the results can readily be plotted or further analysed using the standard libraries of the scientific Python stack.
LiPyphiic is built upon sound software engineering principles. It uses version control, is fully unit-tested, employs continuous integration, and has extensive documentation.
LiPyphilic is also trivial to install -it can be installed using either Anaconda or Pip. 72,73 We encourage users to submit feature requests and bug reports via GitHub, and are always open to new contributors to the project.

Supporting Information Available
Python script for calculating the lipid enrichment index based on tail saturation. Discussion of problems with nojump trajectory unwrapping.