membranes Article In Silico Assessment of the Lipid Fingerprint Signature of ATP2, the Essential P4-ATPase of Malaria Parasites Mario López-Martín 1,2 , Pedro Renault 2,3,4,5 , Jesus Giraldo 2,3,4,5 , José Luis Vázquez-Ibar 6, * and Alex Perálvarez-Marín 1,2, * 1 Biophysics Unit, Department of Biochemistry and Molecular Biology, School of Medicine, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain; mario.lopez@uab.cat 2 Institut de Neurociències, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain; pedro.renault@uab.cat (P.R.); jesus.giraldo@uab.cat (J.G.) 3 Laboratory of Molecular Neuropharmacology and Bioinformatics, Unitat de Bioestadística, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain 4 Unitat de Neurociència Traslacional, Parc Taulí Hospital Universitari, Institut d’Investigació i Innovació Parc Taulí (I3PT), Institut de Neurociències, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain 5 Instituto de Salud Carlos III, Centro de Investigación Biomédica en Red de Salud Mental, CIBERSAM, 08193 Bellaterra, Spain 6 Institute for Integrative Biology of the Cell (I2BC), Université Paris-Saclay, CEA, CNRS, 91198 Gif-sur-Yvette, France * Correspondence: jose-luis.vazquez-ibar@i2bc.paris-saclay.fr (J.L.V.-I.); alex.peralvarez@uab.cat (A.P.-M.); Tel.: +33-1-69-08-89-44 (J.L.V.-I.); +34-93-581-4504 (A.P.-M.) Abstract: ATP2, a putative type 4 P-type ATPase, is a phosphatidylinositol-4-phosphate (PI4P)- regulated phospholipid transporter with an interesting potential as an antimalarial drug target due to its conservation across Plasmodium species and its essential role in the life cycle of Plasmodium falciparum. Citation: López-Martín, M.; Renault, Despite its importance, the exact mechanism of its action and regulation is still not fully understood. P.; Giraldo, J.; Vázquez-Ibar, J.L.; In this study we used coarse-grained molecular dynamics (CG-MD) to elucidate the lipid–protein Perálvarez-Marín, A. In Silico interactions between a heterogeneous lipid membrane containing phosphatidylinositol and Plas- Assessment of the Lipid Fingerprint Signature of ATP2, the Essential modium chabaudi ATP2 (PcATP2), an ortholog of P. falciparum ATP2. Our study reveals structural P4-ATPase of Malaria Parasites. information of the lipid fingerprint of ATP2, and provides structural information on the potential Membranes 2022, 12, 702. phosphatidylinositol allosteric binding site. Moreover, we identified a set of evolutionary conserved https://doi.org/10.3390/ residues that may play a key role in the binding and stabilization of lipids in the binding pocket. membranes12070702 Keywords: P-type ATPase; Plasmodium; flippases; coarse-grain molecular dynamics; lipid–protein Academic Editor: Matthew A. Coleman fingerprinting; modelling Received: 15 June 2022 Accepted: 8 July 2022 Published: 12 July 2022 1. Introduction Publisher’s Note: MDPI stays neutral In addition to the building blocks of membranes, phospholipids are also signaling with regard to jurisdictional claims in molecules that control the development and differentiation processes during the intraery- published maps and institutional affil- throcytic cycle of Plasmodium species, the parasite that causes malaria. For instance, iations. erythrocyte-derived lysophosphatidylcholine controls P. falciparum’s intraerythrocytic cycle by repressing the parasite’s sexual differentiation [1]. Moreover, a tight control of phos- phatidic acid homeostasis is necessary for the parasite’s egress and subsequent invasion of a new red blood cell [2]. In addition, most of the genes encoding membrane transport Copyright: © 2022 by the authors. proteins involved in lipid transport and in the maintenance of membrane lipid asymmetry Licensee MDPI, Basel, Switzerland. are essential for the parasite [3]. One of these essential transporters is ATP2, a putative This article is an open access article distributed under the terms and type 4 P-type ATPase (also called aminophospholipid-transporting P-type ATPase or lipid conditions of the Creative Commons flippase). In eukaryotes, P4-ATPases maintain the asymmetric distribution of lipids in Attribution (CC BY) license (https:// the membranes, necessary for many important physiological processes, including vesicle creativecommons.org/licenses/by/ budding and protein trafficking [4]. While the pivotal role of ATP2 is still unknown, a recent 4.0/). study has linked the activity of ATP2 with the uptake of the erythrocyte’s hemoglobin into Membranes 2022, 12, 702. https://doi.org/10.3390/membranes12070702 https://www.mdpi.com/journal/membranes Membranes 2022, 12, 702 2 of 13 the parasitic digestive vacuole, a fundamental process for parasite progression [5]. ATP2 has also recently emerged as a potential candidate for antimalarial drug targeting since in vitro evolution followed by whole genome analysis has situated the P. falciparum ATP2 as the potential target of two antimalarial drugs [6]. Recently, by using the heterologously- produced ATP2 ortholog encoded by the mouse malaria model P. chabaudi, PcATP2 (~60% sequence identity to the P. falciparum ATP2), it was possible to identify ATP20 s potential Cdc50 β-subunits [7], mandatory proteins for trafficking the majority of lipid flippases to the final membrane destination and, sometimes, for functional activity. Moreover, two phospholipid substrates of PcATP2, phosphatidylserine and phosphatidylethanolamine, were also identified by functional assessment of recombinant PcATP2, as well as the func- tional upregulatory effect of phosphatidylinositol 4-phosphate [7], a lipid known to be indispensable during the intraerythrocytic development of P. falciparum. As ATP2 is highly conserved among Plasmodium species, PcATP2 represents a fair functional and structural paradigm of Plasmodium ATP2. In view of the relevance of lipid–protein interactions for flippases’ membrane traffick- ing, stabilization, and activity, knowledge on the nature of the interacting lipids and the specific residues involved in the interaction is key towards the potential design of efficient antiparasitic strategies. Here, we use coarse-grain molecular dynamic simulations to assess and identify lipid–protein interactions in PcATP2. Using model membranes of asymmetric composition, we aim to fingerprint the lipids and the lipid-binding cavities in PcATP2, which may provide key information for the function and regulation of this transporter, and its therapeutic exploitation. From the perspective of the membrane protein’s structural biology, this approach allows insight to be gained on lipid–protein interactions that can help develop understanding of the structure–function relationships in a potential antimalarial target. 2. Materials and Methods PcATP2 homology modelling. The PcATP2 model is described in a previous study [7]. Briefly, the membrane topology was predicted using the TOPCONS algorithm [8]. The two large intracellular domains (ICD) were independently modeled in MODELLER [9] using 6KG7, 6ROH, 1IWO (ICD1), and 6KG7, 6ROH, 1IWO, 4HQJ, 5MPM (ICD2) as models. Transmembrane regions were modeled using 6K7G and 6ROH as templates in USCF Chimera [10] by alpha carbon molecular replacement. PcATP2 sequence was realigned to the ICD1, ICD2, and TMs models using MAFFT [11] and manually refined. The resulting alignment was used to build the final PcATP2 model in MODELLER [9] using the individual modelled regions as structural templates. Coarse-Grained Martini simulations. The PcATP2 model was converted to the Mar- tini 2.2 coarse-grained representation; in addition, an elastic network (ElNeDyn22) was used to preserve the protein stability and its secondary structures [12,13]. The protein model was inserted into a membrane using the CHARMM-GUI Martini Maker web server [14,15]. The inner leaflet composition consisted of phosphatidylcholine (POPC), phosphatidylethanolamine (POPE), phosphatidylserine (POPS), sphingomyelin (POSM), and phosphatidylinositol (POPI), while the outer leaflet contained the same lipids, except for POPI. Both leaflets contained approximately 300 lipids each (Table 1). The systems were solvated with Martini water and the charge was neutralized with NaCl to a final concentration of 150 mM. A control system was also built with the same membrane compo- sition but excluding PcATP2. In both cases, four replicates of each system with identical characteristics were run individually. Membranes 2022, 12, 702 3 of 13 Table 1. System initial set-up information, simulation conditions, and membrane composition. PcATP2–Membrane Membrane Only Number of 4 4 replicas Simulation time 25 µs 25 µs per replica Time step 20 fs 20 fs Number of 32,181 12,591 atoms Initial box size 14.68 nm × 14.68 nm × 17.99 nm 14.05 nm × 14.05 nm × 8.57 nm (xyz) Outer leaflet Inner leaflet Outer leaflet Inner leaflet POPC 175 100 175 100 POPS 10 36 10 36 POPE 26 103 26 103 POSM 86 21 86 21 POPI 0 45 0 45 Each replicate was energy minimized and equilibrated using the default CHARMM- GUI protocol for Martini 2.2 [16]. Briefly, a two-step minimization of 5000 steps each was carried out using the steepest descent algorithm. The first step was performed with soft-core potentials, and for the second step standard hard-core potentials were used. A 5-step equilibration was performed using an NPT ensemble with an increasing integration step (1, 5, 10, 15, and 20 fs, respectively) and decreasing positional restraints in the lipid heads and backbone beads of the protein for a total equilibration time of 4.75 ns. This 5-step equilibration protocol has been extensively tested for membrane and protein–membrane systems [15,17] and allows gradual relaxation of the solvent and the lipids around the protein, as well as relaxation of the lipid tails, prior to the production run by slowly decreasing the restraints in each step. During this, equilibration temperature was kept at 310 K using a v-rescale thermostat, while the Berendsen barostat was used to maintain the pressure at 1 bar. The reaction field and potential shift Verlet methods were used to calculate electrostatics and Van der Waals interactions, respectively. To ensure the system was in equilibrium during the production run, we plotted the temperature, density, and potential energy of the system for all four replicates to monitor the stability of the system (Figure S1). A 25 µs production for each replica was generated for both the PcATP2–membrane and membrane only replicates using the same parameters as in the equilibration steps but using the Parrinello–Rahman barostat instead. No restraints were applied during the production run. One frame of every 50,000 steps (1 ns time step) of the trajectory was saved for the analysis. Analysis. Trajectories were PcATP2-centered using the GROMACS built-in gmx trjconv tool [18]. The trajectories of the four replicates of each system were concatenated using gmx trjcat for the analysis when it was required [18]. Lipid density maps were generated using gmx densmap [18]. The relative depletion/enrichment (DE) of lipids around PcATP2 was calculated using in-house Python scripts with the MDAnalysis Python package [19,20] applying a contact cut-off distance of 0.55 nm between any bead of PcATP2 and any bead of each lipid species using a similar approach to a previously reported method [21]. The choice of a cut-off distance is dependent on the minimum distance between two beads, which in the case of the Martini 2.2 force field is ~0.50 nm for standard beads (excluding the small beads used for ring structures) [22]. Therefore, the definition of a contact between two beads must be made considering a distance higher but close to 0.5 nm. In our case, Membranes 2022, 12, 702 4 of 13 the choice of a 0.55 nm cut-off follows this criterion and has also been widely used for the contact analysis between Martini CG beads, including many lipid–protein interaction studies [23–27]. In order to calculate the DE indices, the fraction of each lipid species around 0.55 nm PcATP was calculated for each frame of the trajectory. The ratio between this fraction and the fraction of the corresponding lipid species in the leaflet composition was calculated to obtain the DE index for each species in each frame L X (0.55)/L T (0.55) DEX = L X (membrane)/L T (membrane) where DE is the depletion/enrichment index, X is the lipid species analyzed, T refers to all lipid species, and L is the number of lipids. Between brackets the area under consideration, 0.55 nm around PcATP2 or the whole membrane. For each species, DE mean and standard deviation of the whole trajectory were calculated for each individual replicate. The analysis was performed separately for the outer and inner leaflet of the membrane, which were assigned using the MDAnalysis LeafletFinder algorithm [19,20] on the first frame of the trajectory. As POPI was shown to accumulate around PcATP2 and phosphatidylinositol is a known regulator of PcATP2 activity [7], it was decided to further explore POPI–PcATP2 interactions in our simulations by analyzing occupancies, residence times, and putative binding pockets. POPI heads occupancy (O), defined as the fraction of frames in which a contact is present, were calculated with MDAnalysis [19,20] in-house Python scripts applying a cut-off distance of 0.55 nm between any POPI head atoms (NC1, NC2, or NC3) and PcATP2 NA (POPI) OA = NT If A is the residue under analysis, NA (POPI) the number of frames the residue A is 0.55 nm or less from any POPI head bead, and NT is the total number of frames. POPI residue residence times and binding sites analysis were performed with PyLipID [28] using a dual cut-off of 0.5 and 0.65 nm, analyzing 2500 frames (10 ns time step) of each PcATP2 trajectory. Residence times measure the dynamics of bound lipids and provide an insight on the decay periods of bound lipids to residues/binding pockets, revealing specific (high residence times/long decay periods) or transient (low residence times/short decay periods) interactions. Briefly, PyLipID calculates a dissociation constant (koff ) from a survival time correlation function σ (t) Nj T − t 1 1 σ(t) = Nj T − t ∑ ∑ ñj (v, v + t) j =1 v =0 If T is the length of the trajectory, Nj is the number of lipid contacts between any bead −t of PcATP2 and POPI, and ∑vT= 0 ñ j ( v, v + t ) is a binary function that adopts a value of 1 or 0 depending on whether the lipid j contact lasts from time v to time v + t or not, respectively. The σ (t) values are calculated for every time step of the trajectory, from t = 0 to t = T, and divided by σ (0) to normalize the survival time correlation function so σ (t) = 1 at t = 0. The normalized function is fitted to a biexponential model σ (t) ∼ Ae−k1 t + Be−k2 t (k1 ≤ k2 ) where the first term (Ae−k1 t ) models long decays of POPI relaxation, and the second term (Be−k2 t ) models the short decays. A, B, k1 , and k2 are constants. From this model, k1 is taken as the dissociation constant koff . PyLipID calculates the residence time (τ) as τ = 1/ko f f Membranes 2022, 12, 702 5 of 13 To calculate binding sites for a given lipid species, PyLipID detects clusters of residues interacting at the same time with a common lipid molecule. Briefly, for each PcATP2 residue a distance vector is calculated containing the distances to each POPI molecule as a function of time. A Pearson correlation matrix is calculated for those residues binding the same lipid molecule. Using this matrix, a weighted network is constructed in which the protein residues are the nodes and the correlation coefficients of pairs of residues are the weights. The Louvain algorithm [29] is then used to decompose the network into communities (i.e., subsets of nodes/residues more densely connected between them than with the rest of the network). Essentially, Louvain finds the best partitions for the network by iterative optimization of the modularity (Q), which is a measure of the quality of the network partitions by comparing the density of the edges within a community and the rest of the network ki k j   1 ∑  Q= A ij − ∂ ci , c j 2m i,j 2m where m is the sum of all edges in the network, Aij is the edge weight between nodes i and j, ki and kj are the sum of the weights of the edges attached to the nodes i and j, respectively; ci and cj are communities of nodes i and j, respectively, and ∂ is a Kronecker delta function   0 i f i 6= j that ∂ ci , c j = . 1 if i = j Each optimization iteration is separated in two phases. After assigning each node to its own community, modularity is optimized by placing each node to the community that increases the modularity the most. This is repeated until no modularity increase is possible. In the second phase, all nodes within the same community are collapsed into a new node and a new network is created. This network is then used to start a new iteration. Iterations are repeated until no changes in the modularity can be observed. The result of this community analysis is a set of nodes (binding sites). Occupancies and residence times are then calculated for the binding sites, allowing identification of those with more stable interactions with POPI. Residue evolutionary conservation of PcATP2 was performed with ConSurf [30,31] using the default settings. ConSurf conservation scores, PcATP2 residue occupancy by POPI heads and POPI residence time, and POPI-predicted binding pockets were mapped to the all-atom PcATP2 structure b-factor column. Data were plotted using in-house Python scripts and Seaborn [32,33]. Molecular graphics of PcATP2 were generated using UCSF ChimeraX [34]. 3. Results and Discussion In this study, we computationally analyzed the lipid–protein interactions of PcATP2, a close Plasmodium ortholog of the essential P. falciparum P4-ATPase, PfATP2. Live fluorescent imaging of the P. berghei ATP2 ortholog localizes ATP2 at the parasite’s plasma membrane (PPM) [35], situated a few nanometers away from the parasitophorous membrane at the interface parasite/host during the intraerythrocytic stage. In our simulations, we embedded PcATP2 in an asymmetric membrane composed of five different lipid species using coarse- grained molecular dynamic simulations, where the outer leaflet mimics the luminal leaflet of the PPM and the inner leaflet of the cytoplasmic-faced leaflet. All these lipids have been previously identified in a lipidomic analysis of isolated P. falciparum parasites from infected erythrocytes [36,37]. PC and PE are the majority species accounting for nearly 50% of the total lipids of the parasite’s membranes. Sphingomyelin (SM) is the third most abundant lipid, representing ~15% of the total membrane’s lipids. This lipid is important for building and maintaining the tubulovesicular network of membranes, a system used by the parasite to traffic its own proteins and then remodel the parasitized red blood cell membrane [38]. Finally, phosphatidylserine (PS) and phosphatidylinositol (PI) represent, respectively, ~6% and ~5% of all phospholipids present in the parasite’s membranes [36,37]. Interestingly, the membrane of the apicoplast, an internal organelle present in apicomplexan parasites such as Plasmodium is enriched in PI, reaching up to ~15% of total phospholipid content [37]. Membranes 2022, 12, x FOR PEER REVIEW 6 of 13 present in apicomplexan parasites such as Plasmodium is enriched in PI, reaching up to Membranes 2022, 12, 702 6 of 13 ~15% of total phospholipid content [37]. Figure 1a and Table 1 show the system set-up for the simulations. In this asymmetric membrane system, the outer leaflet contained mostly phosphatidylcholine (POPC) and Figure 1a sphingomyelin and Table (POSM), with1 show minorthe system of amounts set-up for the simulations. In this phosphatidylethanolamine asymmetric (POPE) and membrane system, phosphatidylserine the outer (POPS). Theleaflet innercontained leaflet, onmostly phosphatidylcholine the other hand, contained (POPC) higherand sphingomyelin amounts (POSM), of POPE and POPS,with minor as well asamounts of phosphatidylethanolamine phosphatidylinositol (POPE) aand (POPI), which generates phosphatidylserine (POPS). The inner leaflet, on the other hand, contained negative charge at the cytoplasmic side and is also a species that has been involved in higher amounts the of POPE and POPS, as well as phosphatidylinositol (POPI), which modulation of PcATP2 activity [7]. In the middle of the membrane we embedded the generates a negative charge at thePcATP2 coarse-grained cytoplasmic sidewhich model, and iswas alsoobtained a speciesfrom that has been involved a previous in the study [7]. modulation After 25 μs of PcATP2 of MD activity simulation, we[7]. In the middle analyzed of the the whole membrane molecular we embedded dynamic trajectorythe coarse-grained to characterize thePcATP2 model, lipid shell which formed was obtained around from aand the protein previous study establish the[7]. After 25 unique µs of lipid MD simulation, fingerprint of we PcATP2.analyzed the whole molecular dynamic trajectory to characterize the lipid shell formed around the protein and establish the unique lipid fingerprint of PcATP2. Figure Figure 1. (a) 1. (a) System System set-up. set-up. Lateral Lateral view.view. EachEach lipidlipid species species is colored is colored distinctively distinctively (POPC: (POPC: pink, pink, POPE: POPE: orange, orange, POPS: POPS: grey,grey, POSM: POSM: cyan, cyan, andand POPI: POPI: green). green). Only Only lipidlipid phosphate phosphate groups groups are shown. are shown. (b) (b) All-atom All-atom cartoon cartoon representation representation of of PcATP2.Protein PcATP2. Proteinorientation orientationininthe thelipid lipidbilayer bilayerdetermined determined by by the the PPM PPMserver server[39] [39]isisindicated indicatedby byred red(corresponding (correspondingtoto the theouter/luminal outer/luminal leaflet) andand leaflet) blue dotsdots blue (inner/cytoplasmic leaflet). Transmembrane regions of helices 7, 8, and 10 are (inner/cytoplasmic leaflet). Transmembrane regions of helices 7, 8, and 10 are colored in green.colored in green. TheThe remaining transmembrane regions are colored in purple. (c) Average depletion/enrichment remaining transmembrane regions are colored in purple. (c) Average depletion/enrichment (DE) (DE) of of lipids as the ratio between the proportion of each lipid species around 0.55 nm of PcATP2 and the lipids as the ratio between the proportion of each lipid species around 0.55 nm of PcATP2 and the ratio ratio of the same species in the membrane composition. Each individual replicate is represented for of the same species in the membrane composition. Each individual replicate is represented for each each lipid with a different color. Error bars represent standard deviations across the whole simula- tionlipid time.with Theablack different color. dashed Error line (DEbars = 1)represent indicates standard deviations a homogeneous across theofwhole distribution simulation time. the depletion/en- The black richment dashed index. line Values (DE =1 1) under indicates indicate a homogeneous lipid depletion around distribution PcATP2, of the depletion/enrichment while values over 1 indi- cateindex. Values under 1 indicate lipid depletion around PcATP2, while values over 1 indicate an an enrichment. enrichment. Figure 1b shows the average relative lipid depletion/enrichment around PcATP2 Figure with respect to 1b theshows the average membrane relative bulk across lipid depletion/enrichment the whole simulation duration using around a 0.55PcATP2 nm with respect distance cut-off to the membrane between the lipidbulk andacross the protein wholeInsimulation beads. duration the outer leaflet, weusing a 0.55anm observed distance cut-off between the lipid and protein beads. In the outer leaflet, we observed a net net enrichment of POPS, the only negatively charged lipid present in the initial composi- enrichment of POPS, the only negatively charged lipid present in the initial composition of tion of this leaflet. Interestingly, after the simulation we observed a higher variability in this leaflet. Interestingly, after the simulation we observed a higher variability in the lipid the lipid composition of the inner leaflet surrounding PcATP2. This leaflet is consistently composition of the inner leaflet surrounding PcATP2. This leaflet is consistently enriched enriched with POPI while the presence of the other lipid species are depleted, with the with POPI while the presence of the other lipid species are depleted, with the exception of POPS, which remains approximately at the same proportion as in the membrane bulk. On the other hand, in the inner leaflet POPC displayed the highest degree of depletion around PcATP2. In terms of lipid enrichment through the simulation time, the enhancement Membranes 2022, 12, x FOR PEER REVIEW 7 of 13 Membranes 2022, 12, 702 7 of 13 exception of POPS, which remains approximately at the same proportion as in the mem- brane bulk. On the other hand, in the inner leaflet POPC displayed the highest degree of depletion around PcATP2. In terms of lipid enrichment through the simulation time, the of POPS’s presence enhancement around of POPS’s the protein presence aroundwasthehighly protein variable (Figure was highly S2), while variable POPI’s (Figure S2), enrichment while POPI’s at enrichment the inner layer of PcATP2 at the reached inner layer its maximum of PcATP2 reachedvery early in thevery its maximum simulation early in and the remained simulationstable during thestable and remained whole simulation during for allsimulation the whole four replicates (Figure for all S3). The four replicates density plots in Figure 2 that illustrate the redistribution of the different lipid (Figure S3). The density plots in Figure 2 that illustrate the redistribution of the different species around PcATP2 compared to empty membranes clearly demonstrate lipid species around PcATP2 compared to empty membranes clearly demonstrate thatthat this redistribution (enrichment or depletion) this redistribution is driven (enrichment or by PcATP2–lipid depletion) interactions, is driven rather than by PcATP2–lipid spontaneous interactions, ra- lipid–lipid interactions. Similarly, ther than spontaneous lipid–lipidradial distribution interactions. functions Similarly, of the radial lipids around distribution PcATP2 functions of also revealaround the lipids a similar enrichment PcATP2 pattern also reveal (Figure enrichment a similar S4). pattern (Figure S4). Figure2.2.Density Figure Densitymaps maps(in (inlipids nm−−33 )) of lipidsnm of each each lipid lipid species species of of both both leaflets leaflets across acrossthe themembrane membrane planeof plane ofthe thesimulation simulationbox boxfor forboth boththe the control control (without (without PcATP2) PcATP2)membranes membranes(top) (top)and andthe thePcATP2 PcATP2 systems (bottom). systems (bottom). Theenrichment The enrichment in inthe theanionic anioniclipidslipids(especially (especiallyPIs)PIs)around aroundthe theintegral integralmembrane membrane proteinshas proteins hasbeen beenextensively extensivelyreported reportedbefore beforebybymolecular moleculardynamic dynamicsimulations simulationsin inaawide wide varietyof variety ofmembrane membraneproteins, proteins,including including GPCRs, GPCRs, ion ion channels, channels, and and membrane membrane receptors, receptors, among others among others[21,40–42]. [21,40–42]. Indeed, Indeed, CG-MD CG-MDsimulations simulations usingusinghighly highlycomplex complexmembrane membrane compositions with more than 60 different species have shown compositions with more than 60 different species have shown that, despite revealing that, despite revealing unique unique lipid lipid fingerprints fingerprints for eachfor each protein, protein, there isthere is a general a general trend across trend across most proteins most proteins to shellto shell themselves themselves in PIs,with in PIs, along along otherwith other negatively negatively charged charged lipids, lipids,PS generally generally or PA [21]. PS Inorour PA [21]. In our computational computational approach, weapproach, reached similarwe reached similarusing conclusions conclusions a simpleusing a simple membrane mem- system brane system containing onlycontaining only lipid five different five different species. lipid This species. suggestsThis thatsuggests that this this simplistic simplistic membrane system membrane is able to provide system is ableato proper provide lipid environment a proper to the PcATP2 lipid environment in order to the PcATP2 to generate in order accurate to generatepredictions, thus reducing accurate predictions, thus computational requirements reducing computational without compromising requirements without com- the biophysical promising properties of the biophysical the system. properties of the system. The Thelipid lipidhead head group group is is the the main main determinant determinant of of substrate substrate specificity in P4-ATPases. Phosphatidylserine, Phosphatidylserine, one one ofof the the most most prevalent prevalent substrates substrates of of P4-ATPases P4-ATPases acrossacrossdifferent different species species [43], is also one of the putative physiological substrates of is also one of the putative physiological substrates of PcATP2 [7]. Phosphati-PcATP2 [7]. Phos- phatidylinositol, dylinositol, on the other hand, has been extensively characterized in simulation andand on the other hand, has been extensively characterized in simulation ex- experimental studiesacross perimental studies acrossdifferent differentmembrane membraneproteins proteins [21,44,45]. [21,44,45]. PIs PIs are are involved involved in in aa wide widevariety varietyof ofprocesses processes in inthethecellcellmembrane membrane including including protein protein trafficking trafficking and and anchor- anchor- ing. ing. Moreover, Moreover,they theyarearesignaling signalingmolecules moleculeswithwiththethecapacity capacityto tomodulate modulatethe theactivity activityof of integral and membrane-associated proteins. The inositol head integral and membrane-associated proteins. The inositol head group can either be non- group can either be non- phosphorylated phosphorylated or or mono-, mono-, bi-,bi-, or or tri-phosphorylated tri-phosphorylated (PIPs)(PIPs) at at positions positions 3, 3, 4, 4, or or 55 of of the the ring by different kinases, leading to seven different PIPs with ring by different kinases, leading to seven different PIPs with unique subcellular mem-unique subcellular mem- brane brane distribution distribution [46]. [46]. Interestingly, Interestingly,phosphatidylinositol phosphatidylinositol 4-phosphate 4-phosphate (PI4P) (PI4P) is is able able to to upregulate upregulate the ATPase activity of PcATP2 in the presence of PS or PE lipid substrates [7], the ATPase activity of PcATP2 in the presence of PS or PE lipid substrates [7], an observation previously reported in other PcATP2 homologs, such as the yeast flippase Drs2p [47,48]. PI4P is an essential lipid for the malaria parasite found at the plasma mem- brane and at the Golgi apparatus in all stages of the erythrocytic cycle [37]. The inhibition Membranes 2022, 12, x FOR PEER REVIEW 8 of 13 Membranes 2022, 12, 702 8 of 13 an observation previously reported in other PcATP2 homologs, such as the yeast flippase of PI4P Drs2p synthesis [47,48]. PI4P iscompletely an essential blocks lipid forthe thedevelopment malaria parasite and progression found of the at the plasma mem-parasite inside the erythrocyte brane [49].apparatus in all stages of the erythrocytic cycle [37]. The inhibition and at the Golgi of PI4PThe synthesis cryo-EMcompletely structure blocks of the Drs2pdevelopment in its activeand form progression revealedof the a parasite putativein-binding site side for the PI4Perythrocyte between[49]. the TMs 7, 8, and 10 [43,50]. Interestingly, most of the Drs2p residues thatThe cryo-EM establish structure contact of Drs2p with PI4P in areitswell-conserved active form revealed a putative in PcATP2 binding [7]. site forevaluate the To further PI4P between the TMs 7, 8, and 10 [43,50]. Interestingly, most of the Drs2p residues that robustness of our method and gain more structural insights about the PI binding sites establish contact with PI4P are well-conserved in PcATP2 [7]. To further evaluate the ro- of PcATP2, we used PyLipID to analyze the protein–lipid interactions and characterize bustness of our method and gain more structural insights about the PI binding sites of PcATP2, wePOPI putative bindingtosites used PyLipID [28].theBy analyze using thisinteractions protein–lipid more in-depth analysis,pu- and characterize we observed that, POPI tative despite POPIsites binding molecules [28]. Bybeing using well-distributed this more in-depthacross analysis,thewetransmembrane observed that, residues of PcATP2 despite exposed POPI to thebeing molecules innerwell-distributed leaflet of the lipid across membrane (Figure 3a), the transmembrane the residence residues of times of bound PcATP2 POPItolipids exposed are relatively the inner leaflet of the low (<1 lipid µs) in most membrane of the (Figure 3a),positions (Dataset the residence times S2). Indeed, offurther bound analysis POPI lipids are relatively of each replicate low (<1 μs) in most comparing of theoccupancies PyLipID positions (Dataset and the S2).residence In- times deed, further of each analysis residue of eachthat, revealed replicate comparing despite PyLipID observing occupancies a weak positive and the residence correlation between both times of eachmost variables, residue revealed with residues that, despite observing high (>0.5) a weak positive occupancies correlation displayed betweenlow (<1 µs) relatively both variables, most residues with high (>0.5) occupancies displayed relatively low (<1 residence times (Figure S6). Residues with the highest (>2 µs) residence times, on the μs) residence times (Figure S6). Residues with the highest (>2 μs) residence times, on the other hand, tended to also possess high occupancies, but there was also variability across other hand, tended to also possess high occupancies, but there was also variability across residues residues andand replicates. replicates. Interestingly, Interestingly, the highest the highest residenceresidence times in all times fourin all four simulations simulations of replicates were located in in the predicted PI4P binding of replicates were located in in the predicted PI4P binding pocket of PcATP2 (Figure pocket of PcATP2 3b (Figure 3b andDataset and Dataset S2)S2) [7].[7]. Specifically, Specifically, some someof the of the highest highest residenceresidence times of POPItimeswithin of POPIthe within the predicted predicted PI4P PI4P binding binding pocket, pocket, and consistent and consistent in the different in the different replicates,replicates, were for residues were for residues Trp1288, Trp1288,His1295, His1295,Lys1340, Lys1340,Ala1404, and Arg1407, Ala1404, conserved and Arg1407, in Drs2p as conserved in well Drs2p as in as all well as in all Plasmodium ATP2 orthologs, along with residues located around those, Plasmodium ATP2 orthologs, along with residues located around those, such as Tyr1208, such as Tyr1208, Tyr1279, Tyr1279,Phe1280, Phe1280,Asn1292, Asn1292,Phe1294, Ser1359, Phe1294, and Asp1408 Ser1359, (Dataset (Dataset and Asp1408 S2). S2). Figure Figure (a)(a) 3. 3. Total Total occupancy occupancy of theofPcATP2 the PcATP2 residues residues by POPI byhead POPI head represented groups, groups, represented as percent as percent of oftime timethe the residue residue isisinin contact contact (≤0.55 (≤0.55 nm) the nm) with with the inositol inositol group, group, for for the time the combined combined time of all four of all four replicates. replicates.(b)(b) Residence Residence timestimes (in μs)(in of µs) PcATP2 residuesresidues of PcATP2 with POPI, for each with POPI,of for the each four replicates. of the four replicates. The selected plane for visualization contains the residues with the highest residence times, located The selected plane for visualization contains the residues with the highest residence times, located in in the putative PI4P/POPI binding site. the putative PI4P/POPI binding site. Following these observations, we decided to further analyze the POPI–PcATP2 stable interaction sites by calculating putative binding pockets using PyLipID. The four replicates showed similar binding sites, with the binding pocket with the highest residence time being consistent across replicates and also located between TM7–8–10 (Figure 4c), consistent with our previous analysis and experimental data [43]. To further investigate the conservation of key residues of the POPI binding pocket, we performed a ConSurf conservation analysis of PcATP2. The results revealed some hot spots of conserved residues across PcATP2, interaction sites by calculating putative binding pockets using PyLipID. The four repli- cates showed similar binding sites, with the binding pocket with the highest residence time being consistent across replicates and also located between TM7–8–10 (Figure 4c), consistent with our previous analysis and experimental data [43]. To further investigate the conservation of key residues of the POPI binding pocket, we performed a ConSurf Membranes 2022, 12, 702 9 of 13 conservation analysis of PcATP2. The results revealed some hot spots of conserved resi- dues across PcATP2, including many of the POPI binding site (Figure 4a). Indeed, a cor- relation plot between occupancies (>0.1) and conservation scores revealed that, despite including there many not being anyofcorrelation the POPI binding betweensite (Figure 4a).and conservation Indeed, a correlation occupancy, the keyplot between residues of occupancies the pocket were (>0.1) and conservation well-conserved, scores revealed independently of theirthat, despite there occupancies not 4b). (Figure being any Con- correlation sistent with allbetween conservation these analyses, whenand occupancy, looking the key residues at the relationship betweenof residence the pockettimeswere well-conserved, and conservation,independently most of the highof their occupancies residence time (>1(Figure 4b). Consistent μs) residues belonged withto theallPOPI these analyses, binding whenand pocket looking at the relationship were well-conserved. betweenthe Moreover, residence residuestimes of theand conservation, binding pocket most of the high residence time (>1 µs) residues belonged to the POPI with low residence times (<1 μs) also possessed a good degree of conservation (Figurebinding pocketS5).and were well-conserved. Altogether, these resultsMoreover, suggest that thePcATP2 residuescould of the binding have pocket a similar PI4Pwith low residence allosteric activa- times tion site(<1 thatalso to µs) possessed a good experimentally degree ofinconservation characterized Drs2p [43], (Figure and thisS5). Altogether, could these be well-con- results suggest that PcATP2 could have a similar PI4P allosteric activation served across other PcATP2/Drs2p orthologs with a similar allosteric modulation mecha- site to that experimentally nism. characterized in Drs2p [43], and this could be well-conserved across other PcATP2/Drs2p orthologs with a similar allosteric modulation mechanism. (a)ConSurf Figure4.4.(a) Figure ConSurfconservation conservationscores scoresfor forPcATP2 PcATP2mapped mappedonto ontothetheprotein proteinstructure. structure.Golden Golden residuesrepresent residues representresidues residueswith withinsufficient insufficientdata datatotoassign assigna aconservation conservationscore. score.(b) (b)Correlation Correlationplot plot between betweenoccupancies occupanciesand andconservation conservationscores. scores.Orange Orangedotsdotsrepresent representresidues residuesofofthe thepredicted predictedPOPI POPI binding bindingsite. site.Regression Regressionline line with with confidence intervalsisisrepresented confidence intervals representedininblack. black. The The Pearson Pearson coeffi- coefficient cient valuevalue (r) is(r) is shown shown in theintop theright. top right. (c) Predicted (c) Predicted POPI binding POPI binding site in all-atom site in all-atom representation. representation. Residues Residues commoncommon to all fourtoreplicates’ all four replicates’ (core BP) predictions (core BP) predictions are coloredare coloredResidues in green. in green.predicted Residuesinpre- only dicted in only some of the replicates (extended BP) are represented in blue. In purple are Tyr1423 some of the replicates (extended BP) are represented in blue. In purple are Tyr1423 and His1424, the and His1424, the conserved key residues that confer binding specificity for PI4P in Drs2p. conserved key residues that confer binding specificity for PI4P in Drs2p. Interestingly, structural and experimental data in Drs2p revealed that residues Tyr1235 and His1236 confer the specificity for PI4P binding by interacting with the inositol-4- phosphate group of PI4P, and the subsequent activation of the ATPase activity of the transporter [43]. These two residues are highly conserved in Plasmodium orthologues of Drs2p and are present in PcATP2 as Tyr1423 and His1424. Consistently, ConSurf conser- vation scores showed that both these residues were highly conserved across all analyzed sequences (Dataset S1). However, in our PcATP2 model the secondary structure in the corresponding sequence region is not a helix but is instead an extended, disordered confor- mation. Indeed, in Drs2p the helix switch domain is proposed to change its conformation Membranes 2022, 12, 702 10 of 13 upon PI4P binding from an extended conformation in the inactive state of Drs2p to a fully helical structure, which also induces a movement of this region closer to the binding site [43]. The interaction between the negatively charged PI4P phosphate group with Tyr1235 and His1236 stabilizes the binding of the lipid and the helix structure, triggering the release of the autoinhibitory C-terminal domain [43]. Our PcATP2 model likely represents an autoinhibited conformation where Tyr1423 and His1424 are located far from the predicted PI4P binding pocket, and therefore no interactions between the lipid headgroup of POPI and residues Tyr1423 and His1424 are observed in our simulations. Moreover, the elastic network in the coarse-grained model (ElNeDyn) does not allow changes in secondary structures, as they are fixed by defining additional harmonic bonds [51]. Therefore, the set of conformational changes proposed to occur upon PI binding cannot be observed in our current approach. Nevertheless, as with Drs2p, in this autoinhibited conformation we observe the stabilization of the glycerol backbone of POPI through the interactions with positively charged residues, mainly Arg1407. Indeed, density maps of the Cryo-EM structure of the autoinhibited conformation of Drs2p reveals that other lipids can bind to this pocket with less selectivity than PI4P [43]. In our study we only used non-phosphorylated phosphatidylinositol (POPI in Martini 2.2) as a representative of this lipid species. Despite the lack of phosphorylation in the inositol group, POPI shows stable interactions in specific regions of PcATP2, specifically the predicted binding pocket. Nonetheless, a more extensive study using different PIs species mixed together could help an understanding to be gained of which phosphorylation state, if any, is preferred to interact with the putative binding site of PcATP2. Moreover, the use of all-atom MD or a more flexible coarse-grained model that allows changes in secondary structure, such as SIRAH [52], would be particularly useful to determine whether some of the other residues that play a key role in the PI4P-mediated activation of Drs2p, but do not show a relevant role in POPI binding to the PcATP2-predicted PI4P binding pocket (Tyr1411, Lys1412, Lys1415, Lys1416, Tyr1423, and His1424), also play a role in stabilizing the phosphate group of the inositol group in the Plasmodium flippases. In this study we used a computational approach combining homology modeling of PcATP2 with coarse-grained molecular dynamics using the Martini 2.2 force field to study lipid–protein interactions of this flippase with a complex asymmetric membrane containing five different lipid species. We observed that PcATP2 creates a local environment rich in anionic lipids, specifically PI. Moreover, we found a conserved PI binding site situated between transmembrane helices TM7, TM8, and TM10 corresponding to the experimentally- proposed PI4P binding site of Drs2p, the Saccharomyces cerevisiae ortholog of PcATP2. Our results can help in better understanding the mechanism used by PI4P to upregulate the ATPase activity of PcATP2 in the presence of putative substrates. Moreover, our approach also suggests that the use of coarse-grained systems, as well as simplified membranes with a reduced but relevant selection of different lipid species, can give useful insights into protein– lipid interactions that may be functionally significant, as seen, for example, in the study of lipid-modulated GPCR dimerization [53]. Though coarse graining inevitably leads to a loss of detail in comparison to all-atom models, it allows an important gain in sampling, and therefore it represents a good compromise in the investigation of protein–lipid interactions in PcATP2. Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/membranes12070702/s1, Figure S1: Plots of temperature, pressure, density, and potential energy of the systems over the trajectory time. Figure S2: Outer leaflet average depletion–enrichment index; Figure S3: Inner leaflet average depletion–enrichment index; Figure S4: Radial distribution function of each lipid species; Figure S5: Regression plots of the ConSurf conservation score of each PcATP2 residue; Figure S6: Regression plots of the overall contact time of each PcATP2 residue and its residence time with POPI; Dataset S1: PcATP2 ConSurf scores; Dataset S2: PcATP2 POPI PyLipid_Analysis. Figure S7: Structural comparison between PcATP2 model and Drs2p (6ROH). Sequence alignment of the transmembrane regions of PcATP2 and Drs2p [54]. Membranes 2022, 12, 702 11 of 13 Author Contributions: Conceptualization, A.P.-M. and J.L.V.-I.; methodology and formal analysis, M.L.-M. and P.R.; investigation, all authors; resources, J.G. and A.P.-M.; writing—original draft preparation, M.L.-M.; writing—review and editing, all authors. All authors have read and agreed to the published version of the manuscript. Funding: Authors acknowledge financial support by the Spanish Government MCIN/AEI/ 10.13039/ 501100011033 (Project PID2020-120222GB-I00 to A.P.-M.), by the French Agence Nationale de la Recherche, ANR-18-C811-0009-02 to JLVI, and a Margarita Salas Award (MGSD2021-10 to M.L.-M.). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. References 1. Brancucci, N.M.B.; Gerdt, J.P.; Wang, C.; De Niz, M.; Philip, N.; Adapa, S.R.; Zhang, M.; Hitz, E.; Niederwieser, I.; Boltryk, S.D.; et al. Lysophosphatidylcholine Regulates Sexual Stage Differentiation in the Human Malaria Parasite Plasmodium falciparum. Cell 2017, 171, 1532–1544.e15. [CrossRef] 2. Sah, R.K.; Garg, S.; Dangi, P.; Ponnusamy, K.; Singh, S. Phosphatidic acid homeostasis regulated by a type-2 phosphatidic acid phosphatase represents a novel druggable target in malaria intervention. Cell Death Discov. 2019, 5, 107. [CrossRef] 3. Martin, R.E. The transportome of the malaria parasite. Biol. Rev. Camb. Philos. Soc. 2020, 95, 305–332. [CrossRef] 4. Andersen, J.P.; Vestergaard, A.L.; Mikkelsen, S.A.; Mogensen, L.S.; Chalat, M.; Molday, R.S. P4-ATPases as Phospholipid Flippases—Structure, Function, and Enigmas. Front. Physiol. 2016, 7, 275. [CrossRef] 5. Patel, A.; Nofal, S.D.; Blackman, M.J.; Baker, D.A. CDC50 orthologues in Plasmodium falciparum have distinct roles in merozoite egress and trophozoite maturation. bioRxiv 2022. [CrossRef] 6. Cowell, A.N.; Istvan, E.S.; Lukens, A.K.; Gomez-Lorenzo, M.G.; Vanaerschot, M.; Sakata-Kato, T.; Flannery, E.L.; Magistrado, P.; Owen, E.; Abraham, M.; et al. Mapping the malaria parasite druggable genome by using in vitro evolution and chemogenomics. Science 2018, 359, 191–199. [CrossRef] 7. Lamy, A.; Macarini-Bruzaferro, E.; Dieudonné, T.; Perálvarez-Marín, A.; Lenoir, G.; Montigny, C.; le Maire, M.; Vázquez-Ibar, J.L. ATP2, The essential P4-ATPase of malaria parasites, catalyzes lipid-stimulated ATP hydrolysis in complex with a Cdc50 β-subunit. Emerg. Microbes Infect. 2021, 10, 132–147. [CrossRef] 8. Tsirigos, K.D.; Peters, C.; Shu, N.; Käll, L.; Elofsson, A. The TOPCONS web server for consensus prediction of membrane protein topology and signal peptides. Nucleic Acids Res. 2015, 43, W401–W407. [CrossRef] 9. Fiser, A.; Šali, A. Modeller: Generation and Refinement of Homology-Based Protein Structure Models. Methods Enzymol. 2003, 374, 461–491. [CrossRef] 10. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; Ferrin, T.E. UCSF Chimera-a visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605–1612. [CrossRef] 11. Katoh, K.; Misawa, K.; Kuma, K.I.; Miyata, T. MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002, 30, 3059–3066. [CrossRef] 12. de Jong, D.H.; Singh, G.; Bennett, W.F.D.; Arnarez, C.; Wassenaar, T.A.; Schäfer, L.V.; Periole, X.; Tieleman, D.P.; Marrink, S.J. Improved Parameters for the Martini Coarse-Grained Protein Force Field. J. Chem. Theory Comput. 2013, 9, 687–697. [CrossRef] 13. Periole, X.; Cavalli, M.; Marrink, S.-J.; Ceruso, M.A. Combining an Elastic Network With a Coarse-Grained Molecular Force Field: Structure, Dynamics, and Intermolecular Recognition. J. Chem. Theory Comput. 2009, 5, 2531–2543. [CrossRef] 14. Jo, S.; Kim, T.; Iyer, V.G.; Im, W. CHARMM-GUI: A web-based graphical user interface for CHARMM. J. Comput. Chem. 2008, 29, 1859–1865. [CrossRef] 15. Qi, Y.; Ingólfsson, H.I.; Cheng, X.; Lee, J.; Marrink, S.J.; Im, W. CHARMM-GUI Martini Maker for Coarse-Grained Simulations with the Martini Force Field. J. Chem. Theory Comput. 2015, 11, 4486–4494. [CrossRef] 16. Lee, J.; Cheng, X.; Swails, J.M.; Yeom, M.S.; Eastman, P.K.; Lemkul, J.; 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. [CrossRef] 17. Hsu, P.C.; Bruininks, B.; Jefferies, D.; de Souza, P.C.T.; Lee, J.; Patel, D.S.; Marrink, S.J.; Qi, Y.; Khalid, S.; Im, W. CHARMM-GUI Martini Maker for modeling and simulation of complex bacterial membranes with lipopolysaccharides. J. Comput. Chem. 2017, 38, 2354–2363. [CrossRef] 18. Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A.E.; Berendsen, H.J.C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701–1718. [CrossRef] Membranes 2022, 12, 702 12 of 13 19. Michaud-Agrawal, N.; Denning, E.J.; Woolf, T.B.; Beckstein, O. MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. J. Comput. Chem. 2011, 32, 2319–2327. [CrossRef] 20. Gowers, R.J.; Linke, M.; Barnoud, J.; Reddy, T.J.E.; Melo, M.N.; Seyler, S.L.; Domanski, J.; Dotson, D.L.; Buchoux, S.; Kenney, I.M.; et al. MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations. In Proceedings of the 15th Python in Science Conference, Austin, TX, USA, 11–17 July 2016; pp. 98–105. 21. Corradi, V.; Mendez-Villuendas, E.; Ingólfsson, H.I.; Gu, R.-X.; Siuda, I.; Melo, M.N.; Moussatova, A.; DeGagné, L.J.; Sejdiu, B.I.; Singh, G.; et al. Lipid–Protein Interactions Are Unique Fingerprints for Membrane Proteins. ACS Cent. Sci. 2018, 4, 709–717. [CrossRef] 22. Marrink, S.J.; Risselada, H.J.; Yefimov, S.; Tieleman, D.P.; de Vries, A.H. The MARTINI force field: Coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812–7824. [CrossRef] [PubMed] 23. López, C.A.; Agarwal, A.; Van, Q.N.; Stephen, A.G.; Gnanakaran, S. Unveiling the Dynamics of KRAS4b on Lipid Model Membranes. J. Membr. Biol. 2021, 254, 201–216. [CrossRef] [PubMed] 24. Prakaash, D.; Cook, G.P.; Acuto, O.; Kalli, A.C. Multi-scale simulations of the T cell receptor reveal its lipid interactions, dynamics and the arrangement of its cytoplasmic region. PLoS Comput. Biol. 2021, 17, e1009232. [CrossRef] [PubMed] 25. Prasanna, X.; Jafurulla, M.; Sengupta, D.; Chattopadhyay, A. The ganglioside GM1 interacts with the serotonin 1A receptor via the sphingolipid binding domain. Biochim. Biophys. Acta (BBA)-Biomembr. 2016, 1858, 2818–2826. [CrossRef] 26. De Vecchis, D.; Reithmeier, R.A.F.; Kalli, A.C. Molecular Simulations of Intact Anion Exchanger 1 Reveal Specific Domain and Lipid Interactions. Biophys. J. 2019, 117, 1364–1379. [CrossRef] 27. Chakraborty, S.; Wagh, K.; Gnanakaran, S.; López, C.A. Development of Martini 2.2 parameters for N-glycans: A case study of the HIV-1 Env glycoprotein dynamics. Glycobiology 2021, 31, 787–799. [CrossRef] 28. Song, W.; Corey, R.A.; Ansell, T.B.; Cassidy, C.K.; Horrell, M.R.; Duncan, A.L.; Stansfeld, P.J.; Sansom, M.S.P. PyLipID: A Python Package for Analysis of Protein–Lipid Interactions from Molecular Dynamics Simulations. J. Chem. Theory Comput. 2022, 18, 1188–1201. [CrossRef] 29. Blondel, V.D.; Guillaume, J.-L.; Lambiotte, R.; Lefebvre, E. Fast unfolding of communities in large networks. J. Stat. Mech. Theory Exp. 2008, 2008, P10008. [CrossRef] 30. Ashkenazy, H.; Abadi, S.; Martz, E.; Chay, O.; Mayrose, I.; Pupko, T.; Ben-Tal, N. ConSurf 2016: An improved methodology to estimate and visualize evolutionary conservation in macromolecules. Nucleic Acids Res. 2016, 44, W344–W350. [CrossRef] 31. Landau, M.; Mayrose, I.; Rosenberg, Y.; Glaser, F.; Martz, E.; Pupko, T.; Ben-Tal, N. ConSurf 2005: The projection of evolutionary conservation scores of residues on protein structures. Nucleic Acids Res. 2005, 33, W299–W302. [CrossRef] 32. Hunter, J.D. Matplotlib: A 2D graphics environment. Comput. Sci. Eng. 2007, 9, 90–95. [CrossRef] 33. Waskom, M.L. Seaborn: Statistical data visualization. J. Open Source Softw. 2021, 6, 3021. [CrossRef] 34. Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Meng, E.C.; Couch, G.S.; Croll, T.I.; Morris, J.H.; Ferrin, T.E. UCSF ChimeraX: Structure visualization for researchers, educators, and developers. Protein Sci. 2020, 30, 70–82. [CrossRef] [PubMed] 35. Kenthirapalan, S.; Waters, A.P.; Matuschewski, K.; Kooij, T.W.A. Functional profiles of orphan membrane transporters in the life cycle of the malaria parasite. Nat. Commun. 2016, 7, 10519. [CrossRef] [PubMed] 36. Gulati, S.; Ekland, E.H.; Ruggles, K.V.; Chan, R.B.; Jayabalasingham, B.; Zhou, B.; Mantel, P.-Y.; Lee, M.C.; Spottiswoode, N.; Coburn-Flynn, O.; et al. Profiling the Essential Nature of Lipid Metabolism in Asexual Blood and Gametocyte Stages of Plasmodium falciparum. Cell Host Microbe 2015, 18, 371–381. [CrossRef] 37. Ebrahimzadeh, Z.; Mukherjee, A.; Richard, D. A map of the subcellular distribution of phosphoinositides in the erythrocytic cycle of the malaria parasite Plasmodium falciparum. Int. J. Parasitol. 2018, 48, 13–25. [CrossRef] [PubMed] 38. Elmendorf, H.G.; Haldar, K. Plasmodium falciparum exports the Golgi marker sphingomyelin synthase into a tubovesicular network in the cytoplasm of mature erythrocytes. J. Cell Biol. 1994, 124, 449–462. [CrossRef] 39. Lomize, M.A.; Pogozheva, I.D.; Joo, H.; Mosberg, H.I.; Lomize, A.L. OPM database and PPM web server: Resources for positioning of proteins in membranes. Nucleic Acids Res. 2012, 40, D370–D376. [CrossRef] 40. Buyan, A.; Cox, C.D.; Barnoud, J.; Li, J.; Chan, H.S.; Martinac, B.; Marrink, S.J.; Corry, B. Piezo1 Forms Specific, Functionally Important Interactions with Phosphoinositides and Cholesterol. Biophys. J. 2020, 119, 1683–1697. [CrossRef] [PubMed] 41. Sejdiu, B.I.; Tieleman, D.P. Lipid-Protein Interactions Are a Unique Property and Defining Feature of G Protein-Coupled Receptors. Biophys. J. 2020, 118, 1887–1900. [CrossRef] 42. Cabezudo, A.C.; Khoualdi, A.F.; D’Avanzo, N. Computational Prediction of Phosphoinositide Binding to Hyperpolarization- Activated Cyclic-Nucleotide Gated Channels. Front. Physiol. 2022, 13, 485. [CrossRef] 43. Timcenko, M.; Lyons, J.A.; Januliene, D.; Ulstrup, J.J.; Dieudonné, T.; Montigny, C.; Ash, M.-R.; Karlsen, J.L.; Boesen, T.; Kühlbrandt, W.; et al. Structure and autoregulation of a P4-ATPase lipid flippase. Nature 2019, 571, 366–370. [CrossRef] 44. Suh, B.-C.; Hille, B. Regulation of ion channels by phosphatidylinositol 4,5-bisphosphate. Curr. Opin. Neurobiol. 2005, 15, 370–378. [CrossRef] [PubMed] 45. De Matteis, M.A.; Wilson, C.; D’Angelo, G. Phosphatidylinositol-4-phosphate: The Golgi and beyond. BioEssays 2013, 35, 612–622. [CrossRef] [PubMed] 46. Dickson, E.J.; Hille, B. Understanding phosphoinositides: Rare, dynamic, and essential membrane phospholipids. Biochem. J. 2019, 476, 1–23. [CrossRef] Membranes 2022, 12, 702 13 of 13 47. Natarajan, P.; Liu, K.; Patil, D.V.; Sciorra, V.A.; Jackson, C.; Graham, T.R. Regulation of a Golgi flippase by phosphoinositides and an ArfGEF. Nat. Cell Biol. 2009, 11, 1421–1426. [CrossRef] 48. Jacquot, A.; Montigny, C.; Hennrich, H.; Barry, R.; le Maire, M.; Jaxel, C.; Holthuis, J.; Champeil, P.; Lenoir, G. Phosphatidylserine Stimulation of Drs2p·Cdc50p Lipid Translocase Dephosphorylation Is Controlled by Phosphatidylinositol-4-phosphate. J. Biol. Chem. 2012, 287, 13249–13261. [CrossRef] 49. McNamara, C.W.; Lee, M.C.S.; Lim, C.S.; Lim, S.H.; Roland, J.; Nagle, A.; Simon, O.; Yeung, B.K.S.; Chatterjee, A.K.; McCormack, S.L.; et al. Targeting Plasmodium PI(4)K to eliminate malaria. Nature 2013, 504, 248–253. [CrossRef] 50. Bai, L.; Kovach, A.; You, Q.; Hsu, H.-C.; Zhao, G.; Li, H. Autoinhibition and activation mechanisms of the eukaryotic lipid flippase Drs2p-Cdc50p. Nat. Commun. 2019, 10, 4142. [CrossRef] 51. Monticelli, L.; Kandasamy, S.K.; Periole, X.; Larson, R.G.; Tieleman, D.P.; Marrink, S. The MARTINI Coarse-Grained Force Field: Extension to Proteins. J. Chem. Theory Comput. 2008, 4, 819–834. [CrossRef] 52. Machado, M.R.; Barrera, E.E.; Klein, F.; Sóñora, M.; Silva, S.; Pantano, S. The SIRAH 2.0 Force Field: Altius, Fortius, Citius. J. Chem. Theory Comput. 2019, 15, 2719–2733. [CrossRef] [PubMed] 53. Gahbauer, S.; Böckmann, R.A. Comprehensive Characterization of Lipid-Guided G Protein-Coupled Receptor Dimerization. J. Phys. Chem. B 2020, 124, 2823–2834. [CrossRef] [PubMed] 54. Xavier, R.; Patrice, G. Deciphering key features in protein structures with the new ENDscript server. Nucleic Acids Res. 2014, 42, W320–W324. [CrossRef]