Feb 23, 2026
Membrane permeation by molecular dynamics simulation
- Tamir Dingjan1
- 1Weizmann Institute of Science
Protocol Citation: Tamir Dingjan 2026. Membrane permeation by molecular dynamics simulation. protocols.io https://dx.doi.org/10.17504/protocols.io.kxygx8qy4v8j/v1
License: This is an open access protocol distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Protocol status: Working
We use this protocol and it's working
Created: February 23, 2026
Last Modified: February 23, 2026
Protocol Integer ID: 243832
Keywords: molecular dynamics, membrane permeation, membrane permeability of aldopentose stereoisomer, membrane permeation by molecular dynamics simulation, membrane permeability, membrane permeation, aldopentose stereoisomer, permeability coefficient, molecular dynamics simulation, acid bilayer
Abstract
This method is an approach to the simulation of membrane permeability of aldopentose stereoisomers. It yields potential of mean force measurements across a fatty acid bilayer, which can be used as described to calculate permeability coefficients.
Troubleshooting
General description of simulated systems
The study simulated the permeation of four stereoisomers of ribose (β-d-ribopyranose, β-d-arabinopyranose, β-l-xylopyranose, and β-l-ribopyranose) through six different bilayer compositions. Four of the bilayers simulated contain a single lipid species: oleic acid (18:1 Δ9), myristoleic acid (14:1 Δ9), lauric acid (12:0), and POPC. The remaining two bilayers were composed of mixtures of short chain fatty acids: an equimolar mixture of myristoleic acid, lauric acid, and myristic acid (14:0), and a mixture containing excess lauric acid (8 : 1 : 1 lauric : myristoleic : myristic acids). Each bilayer contained a total of 240 lipid molecules. Each bilayer was positioned in the x, y plane in a box of dimensions 6.0 nm × 6.0 nm × 12.0 nm, with a water height of 4.0 nm above and below the bilayer. All fatty acid species were represented as a mixture of protonated and deprotonated in 1 : 1 ratio. Sodium and chloride ions were added as counter ions to neutralize the electric charge of the system, to a final concentration of 200 mM.
The stereoisomers of ribose were all simulated as β-anomeric forms, owing to the greater abundance of the β-pyranose anomer in solution for ribose and xylose. While α-arabinopyranose is more abundant in solution than β-arabinopyranose, we have simulated the permeation of β-d-arabinopyranose for consistency with previous computational work on aldopentose permeation by Wei and Pohorille. Additionally, in agreement with previous spectroscopic and computational experiments determining saccharide conformations in solvent, β-d-ribopyranose was simulated in the 4C1 conformation, β-l-ribopyranose and β-l-xylopyranose in the 1C4 conformation, and β-d-arabinopyranose in the 1C4 conformation for consistency with the previous computation permeation study by Wei and Pohorille.
Molecular dynamics simulation approach
All MD simulations in this study were performed using GROMACS 2022.1. The Accelerated Weight Histogram (AWH) method was used to bias sampling of the solute-bilayer systems in order to calculate bilayer permeability for each solute-bilayer pair. In all simulations, van der Waals interactions were evaluated pairwise within a cutoff of 1.2 nm, with a force-switch between 1.0 and 1.2 nm. Coulomb interactions were calculated within a cutoff of 1.2 nm using PME. Bonds to hydrogen atoms were constrained using the P-LINCS algorithm. Water molecules were parameterized as TIP3P. Lipid and solute molecules were parameterized using the CHARMM36 July 2022 force field. Simulation was performed using a 2 fs integration time step, except where noted otherwise below.
Preparation of simulated bilayer and solute systems
Bilayer preparation
Bilayers were prepared for simulation using CHARMM-GUI. The exact composition of each bilayer is as follows. Oleic acid bilayers contained 60 deprotonated (OLE) and 60 protonated (OLEP) fatty acid molecules in both the upper and lower leaflets. Myristoleic acid bilayers contained 60 deprotonated (MYRO) and 60 protonated (MYROP) fatty acid molecules in both the upper and lower leaflets. Lauric acid bilayers contained 60 deprotonated (LAU) and 60 protonated (LAUP) fatty acid molecules in both leaflets. The equimolar mixture of short chain fatty acids contained 20 molecules of MYRO, MYROP, LAU, LAUP, deprotonated myristic acid (MYR) and protonated myristic acid (MYRP) in both leaflets. The mixture of short chain fatty acids with excess lauric acid contained 48 molecules of LAU and 48 molecules of LAUP in both leaflets, with an additional 6 molecules each of MYRO, MYROP, MYR, and MYRP in both leaflets. Bilayer systems were equilibrated prior to the addition of solute molecules (see below “Bilayer equilibration simulations”).
Bilayer equilibration simulations
Bilayer systems were minimized using the steepest descent integrator for 5000 steps to within a force tolerance of 1000 kJ mol-1 nm-1. Equilibration was performed in the NVT ensemble at a temperature of 296.15 K (coupling time constant of 1.0 ps) in two stages with progressively weaker position and dihedral restraints on all lipid molecules. Both stages used a 1 fs time step for a total duration of 125 ps, with the temperature maintained by the Berendsen thermostat. The first stage used restraints with force constants of 1000 kJ mol-1 nm-2, the second stage used 400 kJ mol-1 nm-2. Following this, bilayer systems were equilibrated in the NPT ensemble at a pressure of 1 atm (coupling time constant 5.0 ps) using the Berendsen barostat applied in a semiisotropic fashion (isotropic coupling in the x and y direction, different in the z direction) with a compressibility of 4.5e-5 bar-1 in both directions. Equilibration proceeded in multiple stages with decreasing force constants: 125 ps with 1 fs time step using 400 kJ mol-1 nm-2 position restraints and 200 kJ mol-1 nm-2 dihedral restraints, followed by 500 ps with 2 fs time step using 200 kJ mol-1 nm-2 for all restraints, 500 ps using 40 kJ mol-1 nm-2 position restraints and 100 kJ mol-1 nm-2 dihedral restraints, 500 ps with no restraints. Finally the system was equilibrated using the Nose-Hoover thermostat and Parrinello-Rahman barostat for 10,000 ps.
Solute system preparation
Ribose stereoisomers were parameterized using the CHARMM36 July 2022 force field version. Initial coordinates were sourced from CHARMM-GUI, using the following residue identifiers: β-d-ribopyranose, BRIBP; β-d-arabinopyranose, BDARBP; β-l-xylopyranose, BLXYL; β-l-ribopyranose, BLRBIP. Solute coordinates were minimized in vacuum (see below “Solute minimization in solvent”), and minimized coordinates inserted into equilibrated bilayer systems using the GROMACS 2022.149 utility gmx insert-molecules. For each combination of bilayer composition and solute, 4 replicate simulations were prepared each containing the bilayer with a single molecule of solute. Replicate bilayer-solute systems were minimized and equilibrated independently (see below “Bilayer-solute system equilibration simulations”), and subsequently simulated together as 4 communicating walker simulations.
Solute minimization in solvent
Solutes coordinates were placed in the center of a box of dimensions 2.5 nm × 2.5 nm × 2.5 nm and minimized in vacuum using the steepest descent integrator for 5000 steps to within a force tolerance of 1000 kJ mol-1 nm-1.
Bilayer-solute system equilibration simulations
A single molecule of minimized solute was inserted into pre-equilibrated bilayer systems using gmx insert-molecules with a van der Waals scale factor of 0.4. A total of 4 independent insertions were performed to yield 4 independent replicate bilayer-solute systems. Each bilayer-solute system was minimized using the steepest descent integrator for 5000 steps to within a force tolerance of 1000 kJ mol-1 nm-1. Equilibration was performed in the NPT ensemble for 1000 ps. Temperature was equilibrated to 310 K using the Nose-Hoover thermostat (coupling time constant 1.0 ps), with separate coupling groups for the bilayer and for the solvent with solute. Pressure was equilibrated to 1 atm using the C-rescale barostat (coupling time constant 5.0 ps), with no compressibility in the z-axis to ensure each independent replicate system maintained the same z-axis length for subsequent AWH multi-walker simulations.
AWH sampling simulations
AWH simulations were used to sample the position of the solute along the z-axis to calculate the potential of mean force for the permeation of the solute through the bilayer. Simulation was performed using a convolved biasing potential with shared biases across the 4 walker simulations. A single bias coordinate was applied, with an estimated initial average error of 40 kJ mol-1. The bias histogram was equilibrated before entering the initial stage. The target distribution was defined as a uniform coordinate distribution across the entire z-axis of the bilayer-solute system. The estimated diffusion constant for the z-axis coordinate dimension was set to 5e-4 nm2 ps-1, with a cover diameter of 8 nm. The reaction coordinate was provided by the pull module, which was configured to use the external potential from the AWH module. The pull was applied in the positive z axis direction, with a force constant of 25,000 kJ mol-1 nm-2. AWH simulations were continued until the walkers had exited the initial stage, and at least one walker simulation had sampled a complete transition of the solute across the entire bilayer, leading to total simulation times between 400 ns to 2,800 ns depending on the solute-bilayer combination.
Permeability coefficient calculation
The free energy of solute permeation was extracted from the AWH simulation energy files using the gmx awh tool. The permeability coefficient was calculated from the free energy of permeation according to the relation from Marrink and Berendsen:
\frac{1}{P_{m}}=\int_{z1}^{z2}\frac{e^{A(z)/k_{B}T)}}{D(z)}dz
Where:
• Pm is the permeability coefficient
• z1 and z2 are the limits of the solute’s position in bulk water on either side of the bilayer
• A(z) is the change in the free energy of the solutes as a function of their position along the z-axis perpendicular to the plane of the bilayer
• kB is the Boltzmann constant (here used as the molar thermal energy value 8.314463 J K-1)
• T is temperature in K
• D(z) is the diffusion coefficient in the direction perpendicular to the bilayer surface.
In the present study, the diffusion coefficient D(z) value used is that determined for the β-d-ribopyranose 4C1 chair form by Wei and Pohorille, shown graphically in Figure 2A. To correct the baseline of the free energy profile for each permeation free energy profile, the average of the force reported in the bulk water region (the 0.5 nm region furthest from the bilayer) was subtracted from all force measurements to shift the reported force profile to 0 in the bulk water regions. Spearman ρ rank correlation measured using the scipy library.