Coarse Grained (CG) simulations were carried out using the MARTINI 2 force field (Marrink, de Vries, and Mark, 2004; Marrink et al., 2007). Membranes were built with the insane.py script available in the Martini web page (Wassenaar et al., 2015). The script was adapted to add RLs in the membranes.

All CG simulations were performed with a time step of 14 fs under periodic boundary conditions applied in all directions for a total simulation time length of 15 µs. For long range interactions, the reaction-field option was used. For the short-range electrostatic and van der Waals interactions, the Verlet cutoffs (Páll and Hess, 2013) were fixed to 1.1 nm. All systems were studied at 299 K with the V-rescale thermostat (Bussi, Donadio, and Parrinello, 2007) in combination with a time constant of 1.0 ps. The Parrinello-Rahman barostat (Parrinello and Rahman, 1981) was used to maintain a mean average pressure of one bar in a semi-isotropic manner with a compressibility of 3.0 × 10^–4 bar^-1 and a coupling constant of 12 ps.

All simulations were minimised using the steepest descent algorithm and then equilibrated for a total time of about 4 ns during which the time step was gradually increased (1 fs, 5 fs, 10 fs, and 14 fs) while decreasing the constraints on polar head movements (from 200, 100, 50, 20 to 10 kJ mol^-1 nm^-2).

Lipidic moieties were described by common Martini parameters. RLs CG topologies (not available in the Martini FF) were built starting from those of monogalactosyldiacylglycerol (MGDG) and digalactosyldiacylglycerol (DGDG) in their anionic form (López et al., 2013; van Eerden et al., 2015). Out of many different models available, differing in the mapping, Martini bead types, rhamnose/s bead bonds, dihedral angles as well as constraint and/or exclusion definitions, we selected the one shown in Figure 2 on the basis of the comparison of results obtained by all-atom (AA) and CG simulations. In particular, a RL (mono- or di-RL) was simulated in both water (10 ns) and POPC membrane (128 total lipids, 500 ns) using both C36 (2 fs time steps) and Martini 2 (20 fs time step) FFs.

First, CG distributions of bonds, angles and dihedral angles were compared to the atomistic ones (for the comparison we used the centre of mass of atoms corresponding to each CG bead). The choice of the best model required several tests until a good reproduction of CG bonds, angles and dihedral angles was achieved (see Supplementary Figures S1, S2 for the optimised mono- and di-RLs topologies, shown in Supplementary Tables S1, S2).

To ensure that the chosen CG topologies could also reproduce membrane properties, we monitored both area per lipid and membrane thickness of (a) RL-containing POPC membrane as well as for (b) a more complex model containing POPC/POPG/ergosterol (53/23/25) and 1 RL/25 lipids (40/60, mono- and di-RL) (see Supplementary Figure S3). Both models allowed us to compare the results with previous works from our group (Monnier et al., 2019). These latter models were run with 256 lipids, 4 mono-RLs and 6 di-RLs for 500 ns (AA) and with 864 lipids, 9 mono-RLs and 27 di-RLs for 10 µs (CG).

For model (a), membrane properties were nicely reproduced with stable trajectories (20 fs time step). Good results were also achieved for model (b) (see Supplementary Figure S3) when reducing the time step (14 fs), as previously reported for other complex glycolipids (López et al., 2013).

The final parameters used for mono- and di-RLs topologies are listed in Supplementary Tables S1, S2. They correspond to those of the model which gave the best results in terms of system stability, bond, angle and dihedral angle distributions (Supplementary Figure S1), and structural membrane properties (Supplementary Figure S3). Mappings and Martini bead assignments are shown in Figure 2A.

In the present work, the following models have been studied by means of CG simulations (Figures 2B, C) (described further in the text in the Results and Discussion Section 3.1).

- Model 1: POPC (1-palmitoyl-2-oleoyl-glycero-3-phosphocholine), POPE(1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine), POPS (1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine), PIP1 (phosphatidylinositol phosphate) and ergosterol (30/18/6/6/40).

CG order parameters were calculated by the python script do-order-gmx5. py available on the MARTINI website (Ingolfsson, 2022) using the equation: P 2 = 1 2 3   cos ⁡ ²   < Θ > − 1

P₂ considers the angle θ formed by the bonds of two consecutive CG beads and the normal of the membrane (axis z). Values of 1, -0.5 and 0 correspond to a perfect alignment, anti-alignment and random orientation with respect to the membrane normal, respectively.

Data were analysed by GROMACS tools (gmx command): density, rdf (for radial distribution function), and (for lateral diffusion). More details about how the calculations are performed are described in detail in GROMACS manual (Lindahl and Spoel. 2022).

Ergosterol flip-flops were calculated using LiPyphilic toolkit (Smith and Lorenz, 2021). A successful flip-flop event was defined as a translocation from one leaflet to another followed by permanence in the destination leaflet for at least 10 ns, as previously done for other works (Gu, Baoukina, and Tieleman, 2019).

AA MD simulations were carried out with Charmm36 (C36) force field (Huang and MacKerell, 2013), using a time integration step of 2 fs under periodic boundary conditions applied in all directions and with an overall time length of 500 ns. The Particle Mesh Ewald (PME) method (Essmann et al., 1995) was used for long-range interactions. For short-range electrostatic and van der Waals interactions, the Verlet cutoff radii (Páll and Hess, 2013) was set at 1.2 nm. The systems were studied at 299 K. The Nosé-Hoover temperature coupling (Hoover, 1985; Nosé, 2002) was used with a time constant of 0.5 ps. The average pressure was maintained at one bar by the Parrinello-Rahman barostat (Parrinello and Rahman, 1981) in a semi-isotropic fashion (pressure in the XY plane independent of the pressure along Z) with a compressibility of 4.5 × 10⁻⁵ bar^-1 and a coupling constant of 5 ps. All bonds were constrained with the LINCS algorithm (Hess et al., 1997). Water molecules were described by the TIP3P model (Jorgensen et al., 1983). All simulations were minimised using the steepest descent algorithm. Equilibrations were carried out using the 6-step protocol provided by CHARMM-GUI for 375 ps.

A water layer of 45 Å thickness was added above and below the lipid bilayer of 128 lipids (64 per leaflet) which resulted in about 12,000 water molecules being incorporated, the exact number depending on the nature of the membrane. Systems were neutralised with K⁺ counterions. Binary systems phospholipid/sterol and the ternary model POPE/POPI/ERGO were mixed in ratio 70/30 and 35/35/30, respectively. The whole process (minimization, equilibration and production run) was repeated once.

“Ligand Reader and Modeller” from CHARMM-GUI was used to generate CHARMM-compatible topologies and parameters files for the anionic RLs form. Penalties reported by CgenFF were between the recommended values for well parametrized molecules (Kim et al., 2017). RLs were later incorporated to the membrane using the usual CHARMM-GUI workflow (Jo et al., 2009) and placed either over the upper leaflet at a non-interacting distance (>10 Å) or placed inside the bilayer at the beginning of the simulation.

Multiple lipid systems have been evaluated during the work.

- Single lipid models containing POPC, POPE, POPS, POPI or DLiPE (1,2-dilinoleoyl-sn-glycero-3-phosphoethanolamine).

Data were analysed by GROMACS tools (gmx command): density (for electron density profiles), potential (for dipole potential calculations), order (for order parameter calculations), rdf (for radial distribution function), and (for lateral diffusion). More details about how the calculations are performed are described in detail in GROMACS manual (Lindahl and Spoel, 2022). Polar and apolar contacts were analysed using radial distribution function by an in-house developed script.

The effect of rhamnolipids was investigated by introducing a mixture of mono and di-RLs in ratio 40/60 (1 RL for 25 lipids) to the model membranes (Monnier et al., 2019; 2020; Botcazon et al., 2022).

- Model 1: Our simulation of model 1 closely reproduced previous data (Ermakova and Zuev, 2017) in terms of both structural (membrane thickness, area per lipid) and dynamical properties (lipid diffusion and CG order parameters) (Supplementary Table S3). Such mono:di-RL and lipid:RL ratios reproduce previous studies where a fluidification effect was observed by 2H-NMR in some sterol-containing membrane models (Monnier et al., 2019) (see Figures 3A–D). For the model 1, which contains only one unsaturation per lipid, we observe a slight decrease in membrane thickness, consistent with an increase in membrane fluidity (Lemkul, 2019; Moradi, Nowroozi, and Shahlaei, 2019) as monitored by CG order parameters (Figure 3). Lipid diffusion appears unaltered (Supplementary Table S3).

An increased fluidity is in very good agreement with experimental observations reported to date, that also point to the modification of phase behaviour due to the RLs action (Sánchez et al., 2006; Oliva et al., 2020). Indeed, this was observed before for RLs in DMPC or DMPS models for which the presence of RLs broadens and shifts the transition to lower temperatures. Similar effect was detected for more complex models composed of DPPC/DPPG/DOPC/DOPG/cholesterol in molar ratio 45/5/20/5/25 (Herzog et al., 2020). Other studies using Atomic force microscopy (AFM) showed that the insertion of the RLs in the latter anionic model caused a reversal of the normal proportion between the liquid ordered (68.5%) and liquid disordered (31.5%) phases (Herzog et al., 2020). Finally, this is in coherence with fluorescence microscopy studies that revealed that the presence of mono-RLs in giant unilamellar vesicles (GUVs) causes the disappearance of pre-existing liquid ordered phase (Herzog et al., 2020).

We took advantage of the higher resolution of AA models to get insight on how RLs are able to get internalised into their target membranes. To this end, placing rhamnolipids (mono-RLs and di-RLs) well outside the surface, we simulated different systems for a long time (up to 3 µs).

RLs penetrate PE-containing membranes much more often and easily than PC- ones, including those containing sterols. This finding might be explained by the significant difference in the steric hindrance of the choline moiety of PC headgroup with respect to the amino group of PE. Moreover, we observe a quicker penetration when sterols are added to the membrane models as discussed further in the text.

As we observe some degree of internalisation in all PE-based systems, independently of the internalisation time, we decided to focus on systems where all RLs are internalised from the start.

To take into account the role of intrinsic fluidity of the membrane and its possible influence on the action of RLs (see Section 3.2.5), each phospholipid was studied with two different types of acyl chains, one with only one insaturation (oleic acyl chain) and one with 6 insaturations (3 insaturations per acyl chain). The effect of sterols was also taken into account as all simulations were run in the absence and in the presence of ergosterol. Besides the commonly studied cholesterol, β-sitosterol and stigmasterol were studied too, allowing to ascertain the differences each sterol could produce in the RLs action (Supplementary Figure S5). The effect of RLs was studied with both mono and di-RLs.

Figure 4 shows snapshots of RLs interacting with a variety of membranes. In most cases, RLs are found well inserted into the bilayers with the rhamnose head in close proximity and just beneath the phosphate groups. RLs acyl chains flank those of phospholipids, especially in the presence of PE or PS headgroups, both playing an important role in several fungal organisms (Cassilly and Reynolds, 2018).

In order to get insight into the interactions favouring the described localization, we analysed polar and apolar contacts (Figures 5–7, Supplementary Figures S6–S10) in the simulations with RLs initially placed inside (to mimic final RLs distribution in the membrane after their full internalisation). MD simulations have previously reported the formation of hydrogen bonds between water and the carbonyl group of DMPS (Oliva et al., 2020). Figure 5 clearly show that RLs only rarely form polar contacts with PC, accounting for RLs selectivity, that would be important to avoid toxicity effects on cells exposing more PC in their outer leaflets (Phoenix et al., 2015; Cassilly and Reynolds, 2018; Ramos-Martín and D'Amelio, 2022). On the contrary, a frequent and dominant interaction (Figure 5) is the formation of a salt bridge between the amino group of PE (or PS) and the carboxylate of RL (atoms O5 and O6 in mono-RLs and O9 and O10 in di-RLs, Figure 1).

PI behaves similarly to PC but in this case the rhamnose oxygen atoms form sporadic H-bonds with its phosphate or the hydroxyl groups of the inositol. Dominant interactions are more evident when we analyse a mixture of phospholipids representing a simplified fungal membrane (PE/PI/ERGO, Figure 5). In this case, the salt bridge with PE largely dominates over all other polar contacts. It is worth noticing that the localization of RLs in the membrane is ideal for the formation of H-bonds between hydroxyl groups of the rhamnose and that of ergosterol.

We believe that the formation of the salt bridge between RLs and the amino group of PE (or PS) plays a role in keeping RL at the observed position. When we analyse the apolar contacts (Figure 6) we find that the acyl chain of RLs make significant van der Waals interactions with those of phospholipids and such contacts become more frequent as we move away from the rhamnose. Interestingly, in the case of PC and PI, RL chains are capable of making contacts with atoms of phospholipid chains localised close to the core of the bilayer (in Figure 6 atom colours become darker as we move away from the headgroups), indicating that they are not well anchored to the surface and confirming the role of the salt bridge with the PE (or PS) headgroup. In all cases, contacts with unsaturated chains (oleoyl) seem slightly favoured over those with saturated ones (palmitic).

Our observations on RL localization in the bilayer are coherent with previous FTIR studies showing that the presence of di-RLs affects the C=O bonds of the polar heads of a DPPC model and, to a lesser extent, the bonds of the hydrophobic chain (Sánchez et al., 2009). As for the importance of the salt bridge, isothermal calorimetric titration (ITC) studies have shown that, compared to liposomes formed only by POPC, di-RLs have a higher affinity for those composed of POPC/POPE (1/1) but less affinity for liposomes containing POPC/cholesterol (1/1) (Aranda et al., 2007).

FTIR combined with molecular modelling of an implicit membrane (hypermatrix) showed that the presence of sterols tends to increase the interaction of RLs with membrane lipids (Nasir et al., 2017).

These data seem to indicate that sterols may play an important role in the interaction of RLs with biomembranes. We therefore investigated apolar contacts between RLs and several sterols (Figure 7, Supplementary Figures S9, S10). Our data show that, while RL acyl chains interact frequently with phospholipids with their terminal ends, they seem to adhere completely to sterols when they are present. Such interactions (in violet in Figures 7, Supplementary Figure S9) are important and tend to dominate over all other apolar contacts both with the more rigid (containing palmitoyl and oleoyl acyl chains) and more fluid membranes (containing dilinoleyl-acyl chains, shown in Figure 7). It should be noted that similar effects are found also in POPC/CHOL systems (Supplementary Figure S10), when forcing RLs to be inserted in the membrane at the start of our simulations. Conversely, while free floating RLs spontaneously insert in PE-containing bilayers, they are unable to integrate mammalian models made of POPC/CHOL at least in the time range explored (2 μs), thus providing a possible explanation for the selectivity.

Indeed, our data are consistent with previous literature and provide a molecular picture of the interaction (Figure 8). RLs approach the membrane surface by its polar head. In the case of POPE the driving force is the electrostatic interaction between its carboxylate and the amino group of POPE. Although RLs tend to aggregate before they get inserted, they usually penetrate the membrane one-by-one. After reaching the surface, RLs insert one acyl chain in the hydrophobic part of the bilayer (Figure 8A). In the case of POPE, the bridge with the amino group is maintained as the molecule further descends into the bilayer, facilitating the entrance of the second RL chain. Once inserted, the acyl chains of RLs recruit two (or more) sterol molecules (Figure 8B). The remarkably similar length of RL acyl chains and ring structure of sterols maximises van der Waals interactions. In this localization the rhamnose rings can form H-bonds with polar atoms of phospholipid head groups while the salt bridge between the RL carboxylate and the amino group persists. Occasionally the carboxylate can also form H-bonds with the hydroxyl of sterols, especially in phospholipids lacking the amino groups. Rarely it can also make intramolecular H-bonds with its rhamnose and descend deeper in the bilayer.

Our coarse grained simulations involving complex fungal membrane models (Figures 2B, C) have highlighted a fluidification of biomembranes caused by the addition of RLs (Figure 3), an effect observed also in other models by Laurdan-based fluorescence spectroscopy (Herzog et al., 2020), 2H-NMR and MD simulations (Monnier et al., 2019).

In general, mono-RLs and di-RLs effects are similar in trends and slightly more intense for di-RLs. Therefore, for simplicity we will focus on DiRLs.

Fluidification can be monitored by changes in the order parameter of phospholipid acyl chains (Figure 9 and Supplementary Figure S11). As it was observed in our CG simulations, also in AA simulations we observe a larger fluidification for membranes with only one unsaturation (POP-phospholipids) rather than those bearing 6 unsaturations (DLi-phospholipids).

Analysing electron density profiles, we can evaluate lipid interdigitation processes or an increase in the amount of water molecules near the polar head groups due to loose packaging. The POPE model displays stronger changes than POPC model (Supplementary Figure S13), showing an increase in the density in the central part of the membrane due to lipid interdigitation and a reduction of the density in the headgroups zone, probably due to the intercalation of additional water molecules. The presence of sterols significantly reduces the effect in mixed POPE models (especially sitosterol and stigmasterol). Similar phenomena are observed for other phospholipids often found in fungi such PS and PI (Supplementary Figure S14). On the other hand, no big changes are observed between single and mixed POPC-based models (Supplementary Figure S15).

The dipole potential ΨD originates from the alignment of dipolar parts of the lipids and water molecules (Wang, 2012; Clarke, 2001; Pearlstein, Dickson, and Hornak, 2017; Brockman, 1994; Ramos-Martín and D’Amelio, 2022). The magnitude of the dipole potential is dependent on the structure of the lipid, the degree of unsaturation and the nature of the linkage between headgroup and hydrocarbon chain (ester or ether). Intercalation of dipolar molecules such as some antimicrobial peptides can greatly modify the potential, leading to significant membrane damage (Dreyer et al., 2013; Pearlstein, Dickson, and Hornak, 2017). Changes in the dipole potential are in fact able to affect the translocation rates of ions across the bilayer and the partition and translocation of macromolecules. The gating kinetics of voltage-gated ion channels can also be altered by the modulation of ΨD. This mechanism of action is employed by some antifungal agents (Zakharova et al., 2019) due to the fact that ΨD is important for fungal survival (Zhao and Tombola, 2021), as in eukaryotes for the regulation of cardiac or neuronal processes (Brockman, 1994; Pearlstein, Dickson, and Hornak, 2017).

Given the importance of ΨD for fungal survival, we examined changes in dipole potential in single-lipid models. The most important reduction is observed for pure POPE (Supplementary Figure S16). Interestingly, the presence of sterols reduces the effect of RLs in all different types of membranes examined (Supplementary Figures S16–S19).

Regarding other parameters, such as the area per lipid, we can observe a slight increment due to RLs incorporation to the membranes, specially in POP-models (Supplementary Figure S20). Membrane thickness (Supplementary Figure S21) slightly decreases in all the membrane models, coherently with the reduction of the hydrophobic thickness that can be monitored by electron density profiles (Supplementary Figure S13).