The GROMACS molecular package version 2021, with the CHARMM36 all-atom force field, was used to conduct simulations of various polymers interacting with a lipid bilayer membrane (
17). The coordinate and topological information related to the starting structures of the lipid bilayer and polymers were obtained from the CHARMM-GUI server (
https://www.charmm-gui.org/).
The lipid content and their constituent percentages were prepared to mimic those of the human cell membrane, including cholesterol, dipalmitoylphosphatidylcholine (DPPC), 1-palmitoyl-2-oleoylphosphatidylcholine (POPC), 1-stearoyl-2-oleoylphosphatidylcholine (SOPC), dipalmitoylphosphatidylethanolamine (DPPE), 1-palmitoyl-2-oleoylphosphatidylethanolamine (POPE), 1-stearoyl-2-oleoyl-phosphatidylethanolamine (SOPE), stearoyloleoylphosphatidylserine (SOPS), and 1-stearoyl-2-oleoyl phosphatidylglycerol (SOPG) (
Table 1). The percentage and number of each lipid in the model membrane were selected based on a literature review (
18-
20).
For polymer selection, two positively charged polymers, chitosan and PEI, with high and low molecular weights, two negatively charged polymers, alginate and PLGA, with high and low molecular weights, and cellulose as a neutral polymer were used to investigate the interaction of polymers with the bilayer membrane regarding their charges and molecular weights of monomers (
Figure 1).
A certain number of each polymer was randomly placed on the outer layer of the membrane so that the final systems were equal, based on their molecular weight. The TIP3P water model (
21) was used for solvation of the system, and to completely neutralize the charge, an appropriate number of counter ions (Na
+ and Cl
-) were added to the simulation boxes. Thus, Na
+ ions were added to the boxes containing negatively charged polymers, and Cl
- ions were added to the boxes containing positively charged polymers. The water model is recommended by the software developer as a more suitable model for working with the CHARMM force field.
Energy minimization of the systems was performed using the steepest descent algorithm until the maximum force on the atoms reached a final value of < 10 kJ/mol.nm (
22). Periodic boundary conditions (PBC) were applied in all X, Y, and Z directions to eliminate energy drifts resulting from edge effects in the small simulation boxes. To simulate real body conditions, the temperature and pressure were adjusted to 310 K and 1 bar using NVT and NPT ensembles for 200 ps. In this regard, the temperature and pressure were equilibrated using the Hoover thermostat, Nose thermostat and Parrinello and Rahman barostat, respectively. Both methods are precise algorithms for coupling the conditions of molecular simulation (
23-
25). The linear constraint solver (LINCS) algorithm was employed to keep all bonds around their equilibrium value (
26). Short-range van der Waals and Coulombic interactions were calculated for all particles within a cut-off distance of 1 nm. The particle-mesh Ewald (PME) method (
27) was also applied for the calculation of long-range electrostatic interactions. This method helps manage the atomic forces between particles outside the central box when PBC is applied in the simulation.
Eventually, all simulations were conducted for 100 nanoseconds under the default leap-frog algorithm in the GROMACS package. The output structure of the simulated systems was prepared with the visual molecular dynamics (VMD) package (
28). The interaction energies were calculated using the MM/PBSA method (
29), and the GridMAT-MD tool (
30) was used for thickness parameter measurements.