Unraveling the Phase Behavior, Mechanical Stability, and Protein Reconstitution Properties of Polymer–Lipid Hybrid Vesicles

Hybrid vesicles consisting of natural phospholipids and synthetic amphiphilic copolymers have shown remarkable material properties and potential for biotechnology, combining the robustness of polymers with the biocompatibility of phospholipid membranes. To predict and optimize the mixing behavior of lipids and copolymers, as well as understand the interaction between the hybrid membrane and macromolecules like membrane proteins, a comprehensive understanding at the molecular level is essential. This can be achieved by a combination of molecular dynamics simulations and experiments. Here, simulations of POPC and PBD22-b-PEO14 hybrid membranes are shown, uncovering different copolymer configurations depending on the polymer-to-lipid ratio. High polymer concentrations created thicker membranes with an extended polymer conformation, while high lipid content led to the collapse of the polymer chain. High concentrations of polymer were further correlated with a decreased area compression modulus and altered lateral pressure profiles, hypothesized to result in the experimentally observed improvement in membrane protein reconstitution and resistance toward destabilization by detergents. Finally, simulations of a WALP peptide embedded in the bilayer showed that only membranes with up to 50% polymer content favored a transmembrane configuration. These simulations correlate with previous and new experimental results and provide a deeper understanding of the properties of lipid-copolymer hybrid membranes.

This standard representation for polymers in CG scale was initially employed, and bonds, angles, and dihedrals for this model were also obtained (as discussed in Section B).This method resulted in a satisfactory description of the PBD-b-PEO all-atom (AA) structure, but when polymers and phospholipids were mixed to form hybrid bilayers, we observed a phase separation (Figure S2.A).This outcome conflicted with the expectation of homogeneous dispersion, as shown by Seneviratne et al. 8 , indicating that the interactions between the molecules were not accurately computed.To determine which part of the molecule needed to be better modeled, we conducted tests by altering the Lennard-Jones parameters of the C4 and SN3a beads.We found that the C4 beads exerted a considerable repulsion on the acyl chains of the phospholipids, which caused the observed separation.To address this issue, a two-bead parameterization using Tiny Martini beads was proposed for PBD.The main chain, which consists of two bonded carbon atoms, was represented by the TC2 bead while the side chain, composed of two double bonded carbons, was represented by TC4 (Polybutadiene -Model 2 in Figure S1).T beads were used for this polymer since each particle contains only two heavy atoms.Despite the higher computational cost, this approach produced well-mixed hybrid bilayers (Figure S2.B).Martini beads can be further divided into attractive, normal and repulsive beads, and all different conformations were tested, and electron density profiles were compared with experimentally obtained ones by Ref. 9 .The normal version of the SN3a, along with the repulsive versions of the TC2 and TC4 beads (TC2r and TC4r) were used, which resulted in electron density profiles very similar to the observed in experiments, as discussed in the main manuscript.The developed model also introduced a small amount of hydrophilic PEO beads into the hydrophobic core.This contradicts experimental findings 10 , and a similar issue had been previously reported in another study that used the Martini force field to model a linear conformation of the PBD-b-PEO polymer 1 .To address this issue, Grillo and co-workers 1 reduced the ϵ parameter describing the interactions between PBD and PEO beads by 11.4 %.Following a similar approach, we decreased this parameter by ≈ 9 %, resulting in ϵ TC2r−SN3a changing from 1.6 to 1.45 kJ/mol and ϵ TC4r−SN3a changing from 1.75 to 1.6 kJ/mol.This small adjustment was sufficient to achieve proper separation of hydrophobic and hydrophilic polymer regions.

B. Bonded parameters
The bonded parameters were derived by comparing the distributions of the CG beads with those from all-atom (AA) simulations and further refined to represent experimental features visualized in experiments.PBD and PEO chains with 10 monomers were analyzed to derive parameters for the pure polymers, while PBD 5 -b-PEO 5 was used to derive parameters for the copolymer hydrophobic-hydrophilic interface.The AA simulations were carried out in the NPT ensemble (T= 310 K, P = 1 atm) with OPLS-AA FF force field, using the Nosé-Hoover barostat and thermostat to control pressure and temperature, respectively.Each polymer chain was solvated with 8000 TIP3P water molecules, using 1.4 nm cut-offs for both nonbonded and electrostatic potentials, and PPPM for long-range electrostatic interactions.The simulations included a 50 ns equilibration and a 100 ns data collection run.The CG representations were generated from AA structures using the CG builder tool https://jbarnoud.github.io/cgbuilder/.For PEO, bonds, angles, and dihedrals were parametrized, while only bonds and angles were used for PBD due to the lack of a well-defined dihedral distribution for this monomer.PEO models have been previously developed for the Martini force field, and initial parameters for this molecule were adopted from literature sources 1,11,6 and refined based on our simulation results.To match the AA bonded parameters with the CG representation, a combination of optimization algorithms and trial and error were employed until the CG distributions were in close agreement with those of the AA simulations.
The harmonic potential was used to model the bonds following the Martini force-field formulation.Bond distributions of the simulated polymers are shown in Figure S3, indicating good agreement between the CG and AA simulations.However, the bonds between TC2r-TC2r and TC2r-SN3a beads showed a bimodal behavior in the AA representation which was not captured by the CG model.As a result, the two configurations were combined into one due to the simplified phase-space sampling typical of mesoscale models.The use of non-charged Martini 3 beads in modeling can lead to underestimation of structural properties when using bond distances directly from atomistic simulations, as their attractive interactions are slightly underestimated 12 .To address this issue, bond lengths can be adjusted for a more accurate representation.One proposed solution is to slightly increase (≈ 15 %) the atomically adjusted bond length of main chains 12 .Following this suggestion, we applied this correction for SN3a-SN3a bonds and TC2r-TC2r bonds.TC2r-TC2r TC2r-TC4r TC2r-SN3a SN3a-SN3a Fig. S3.AA (blue) and CG (red) bond length distributions.
The angles were also modeled using a harmonic potential, and the results are depicted in Figure S4, which shows close agreement with the atomistic reference data.Five different types of angle were considered, including main chain angles (TC2r-TC2r-TC2r, TC2r-TC2r-SN3a, TC2r-SN3a-SN3a, SN3a-SN3a-SN3a) and the main and side chain of a PBD bead with its neighboring particle (TC4r i -TC2r i -TC2r i+1 ).Again, angular distribution that exhibited a bimodal behavior were combined into a single distribution in the CG scale.Simulations became unstable when the SN3a-SN3a-SN3a angles approach 180°if only an harmonic potential is used, as also reported elsewhere ( 6, 7 ).To correct this, the constrained bending potential developed by Bulacu et al. 13 was applied along with the harmonic potential for this angle.With this correction, a timestep of 2 femtoseconds could be used, whereas without it, only a step of 1 femtosecond could be applied.TC2r-TC2r-TC2r TC4r-TC2r-TC2r SN3a-SN3a-SN3a TC2r-TC2r-SN3a TC2r-SN3a-SN3a The dihedral distribution of PBD did not exhibit a well-defined behavior, and therefore dihedral terms were not included for these beads.Regarding PEO, there are some peculiarities in modeling this term, such as bimodalities, and for simplicity we used the CG dihedral model implemented in the Poliply Python suite 11 .A summary of all the bond parameters is provided in Table S1.

C. Validation of the proposed model
To validate the efficacy of the proposed parameterization, we conducted calculations of the radius of gyration (R g ) and surface area of available solvent (SASA) for both the all-atom and CG representations.The SASA was calculated based on the system used for parameterization, as it represents the polymers fully immersed in water.Meanwhile, the R g was calculated through independent simulations of 30 polymer chains in a simulation box to assess the ability of the proposed model to capture the degree of packing of the bulk polymer.It is important to note that both bonded and non-bonded parameters can impact the R g and SASA, and the presented results are based on the final set of parameters utilized in the simulations.The radius of gyration is a measure of the overall size of the molecule and is defined by Eq.S1: where r(i) and r com are the coordinates of atom i and the center of mass of the molecule, respectively, and N is the number of atoms.Figure S5 illustrates the radius of gyration of the three polymer chains systems studied during parameterization.The excellent agreement between the distributions of the AA and CG trajectories suggests that the Martini parameterization proposed in this work accurately represents the overall size of the molecule.SASA measures the area of a molecule available to establish contact with its surrounding solvent or media.As polar and nonpolar residues are being parameterized, SASA is a valuable metric for evaluating their behavior in a solvated environment.According to Lee and co-workers 14 , SASA is calculated based on the atomic van der Waals radius by: where R is the radius, L i is the length of the arc drawn in a given section i, Z i is the perpendicular distance of the section i to the center of the sphere, ∆Z is the space between sections and ∆ ′ Z is (∆Z/2) or (R − Z j ) (whichever is smaller).The SASA distributions of the three polymers were evaluated for both their AA and CG representations (Fig. S6).In general, the CG model underestimated the SASA values by about 5 %, which is in good agreement with the trends reported by the Martini developers 3 .As this is a small deviation, our analysis of the radius of gyration (R g ) and SASA provides strong evidence that the CG model is a good representation of the AA structure.

INCORPORATION OF PEPTIDE A. Umbrella Sampling
Umbrella sampling (US) simulations were conducted to analyze the interactions of HVs with a peptide model as discussed in the main text.Figure S7 shows the system representation with the peptide embedded in the center of the membrane and histograms of the US windows.The good overlap between adjacent windows suggests that the simulations were properly spaced and converged.The peptide can be pulled out of the membrane from both its carboxyl and amine termini.While the main manuscript reported results for the N terminus, Figure S8 illustrates the results for both termini.The curves obtained for each terminus differ slightly, primarily due to differences in the peptide's behavior within the membrane.The data indicate that the energy in the center of the bilayer and in the bulk water are equal for both termini, within the error margin.The findings suggest that the charged termini primarily influence the behavior of the peptide within the bilayer.These results are in line with the findings of Yue et al. 15 , who proposed that the N terminus is thermodynamically preferred for peptide reconstitution.In line with these results, our research revealed an additional energy barrier for the C terminus.

B. Equilibrium simulations
We further investigated the conformations of the transmembrane peptide that showed a favorable transmembrane configuration using equilibrium simulations.To examine local properties, we discretized the membrane into a grid with a resolution of 0.5 nm along the x and y axes.We applied the nonlocal weighted averaging method proposed by 16 to ensure that no empty cells were present.This approach involves determining the local value of a variable by averaging the values of that property for all molecules within a radius of σ from the center of the cell, known as the horizon.The local value of a generic variable f (x i ) is then determined by the following relation: where w(r x j ) is the weight assigned to a particle at radius r x j away from the center of the element xi , f (x j ) is the value of the measured property for this particle and N c is the number of particles that lie within the horizon.The weights assigned to each molecule are given by S4: It can be observed that the average weight of molecules decreases exponentially with distance from the center at a rate defined by α.The value of α can be adjusted to achieve a desired level of resolution and, in this study, α = 8 and σ = 0.5 nm, which have been found to provide an adequate sampling by preventing empty cells.The profiles for the total, POPC, and polymer-normalized density, membrane thickness, and POPC order parameter obtained at all analyzed temperatures and bilayer compositions are shown in Figures S9-S13.

FIRST MOMENT OF THE PRESSURE PROFILE
Table S2.First moment of pressure profiles for all temperatures and membrane compositions.Distance (nm)

ADDITIONAL ANALYSIS REGARDING THE STRUCTURE OF THE MEMBRANE
To analyze the degree of membrane interdigitation, we utilized the MembPlugin developed by Ref. 17 .This tool employs a correlation-based approach to evaluate the overlap of the mass distribution of the two leaflets along the membrane normal.
The fraction of overlap, denoted by the symbol I ρ , is computed using the equation , where ρ a and ρ b represent the density profiles of the two leaflets.The resulting interdigitation values at a temperature of 310 K are presented in Figure S17 for the complete membrane, as well as for the POPC and PBD 22 -b-PEO 14 molecules.At low polymer concentrations, the degree of interdigitation increases due to an increase in the fraction of polymers that interdigitate.Simultaneously, the lateral area of the membrane expands, which has been discussed in detail in the main manuscript.The combined effect of these phenomena is responsible for the observed thin membrane configuration, which is further elaborated on in the main document.It is noteworthy that both the 25:75 and 50:50 cases exhibit an increase in membrane interdigitation of approximately 25% compared to membranes with higher polymer concentration.Interestingly, the average lateral area of the membrane also increases by a similar amount in these compositions.These findings suggest that interdigitation and collapsed structures may have comparable influence on the 'thin' configuration of the membrane Furthermore, the indicate that as the polymer fraction increases, the interdigitation of POPC decreases, indicating that more lipids are concentrated in the hydrophilic interface of the thicker conformation.

Fig. S2 .
Fig. S2.Hybrid membranes formed by PBD modeled by (A) one and (B) two beads.
Fig. S7.(A) Representation of a system used in the US simulations and (B) Histogram showing the overlap between adjacent windows.

Fig. S8 .
Fig. S8.Potential of mean force for different vesicles (lipid:polymer %) for the peptide being pulled from the N (red) and C (blue) terminus.
Fig. S17.Fraction of overlap between leaflets for different HVs at 330 K.

FigureFig. S18 .
Figure S18 displays the phosphate-phosphate distance of a 50:50 HV.Thin regions (with a of approximately 40 Å) are present throughout the membrane, with greater concentration observed in the upper part.It is noteworthy that this configuration constitutes only about 20% of the total area and is embedded in a much thicker membrane (approximately 70 Å), providing further evidence of the existence of two distinct configurations in this condition.

Table S3 .
4PL parameters and regression coefficient for varied polymer/lipid ratios.