High Genus Membranes - weria-pezeshkian/FreeDTS GitHub Wiki
High Genus Membranes
Created by Beatrice Geiger (July 2025)
Introduction
On this page we provide detailed tutorials for simulating High Genus Membranes with FreeDTS, specifically to reproduce results to understand the bimodal mechanical response of membrane necks (Geiger and Pezeshkian, 2025). We will guide you through setting up a high-genus system as whole stomatocytic vesicle and in periodic boundary conditions (PBC) in Part1. Based on these simulations you will then be able to simulate systems under internal pressure, using a core bead in Part2 as well as systems in PBC under lateral tension or over a range of constant frame areas in Part3. Lastly, we provide you with initial files to start all kinds of simulations or skip steps here.
Part1
Preparing High Genus Systems
Part1.1: High Genus Stomatocytes
To simulate a whole vesicle system we need to chose a genus>17 so that it will adopt a stomatocytic shape naturally without a need for volume control. We use "stomatocyte" as a generalized term here: not only referring to a genus=0 shape as in red blood cells (typical use) but also to shapes like the one of the nuclear envelope, which has a high-genus surface that can be described as two quasi-spherical membranes connected by multiple membrane necks. To setup such system we use the GEN property of FreeDTS to create a closed cuboid surface with 21 necks (genus=20) with the command:
$path/GEN -box 100 100 100 -type high_gen -g 21 -N 30 -o gen20.q
We now need to relax and equilibrate this membrane mesh. To do so we need an input file defining the paramemeters for the simulation. It should contain:
Integrator_Type = MC_Simulation
Min_Max_Lenghts = 1 3
MinfaceAngle = -0.5
Temperature = 1 0
Box_Centering_F = 0
Set_Steps = 1 20000000
EnergyMethod = FreeDTS1.0_FF
Kappa = 10 0 0
Edge_Parameters = 5 0 0
VertexArea = 0 0.7 0 0
TimeSeriesData_Period = 100
VertexPositionIntegrator = MetropolisAlgorithmOpenMP 1 1 0.05
VolumeCoupling = SecondOrder 0 0.0 0.0
AlexanderMove = MetropolisAlgorithmOpenMP 1
VisualizationFormat = VTUFileFormat VTU_F 2000
NonbinaryTrajectory = TSI TrajTSI 2000
Restart_Period = 1000
You can also find this file here as "input_relaxG20.dts". As you can see we do not use a volume constraint here (we have set the coupling constant to 0). You can now start the simulation, after making a .top file:
echo "gen20.q 2" > top.top
$path/DTS -in input_relaxG20.dts -top top.top
Please check your energy output files (dts-en.xvg) afterwards to see whether your system has equilibrated within the 20.000.000 steps or whether a longer simulation is needed. The last simulation output provides you with a genus-20-stomatocyte with which you can proceed in Part2. You should see be able to see the following development of your system, using a visualization software such as paraview:
Fig.1 Progression from initial configuration to high genus stomatocyte.
Part1.2: A Neck in PBC
To uncover more details one may want to isolate one neck from the whole stomatocyte. To do so, we can simulate a neck embedded in two quasi-flat membrane segments in PBC. To setup such system we use the GEN property of FreeDTS to create an opensided cuboid surface with 1 neck in PBC with the command:
$path/GEN -box 30 30 30 -type high_genpbc -g 1 -N 20 -o gen1_PBC.q
We now need to relax and equilibrate this membrane mesh. To do so we need an input file defining the paramemeters for the simulation. It should contain:
Integrator_Type = MC_Simulation
Min_Max_Lenghts = 1 3
MinfaceAngle = -0.5
Temperature = 1 0
Box_Centering_F = 0
Set_Steps = 1 20000000
EnergyMethod = FreeDTS1.0_FF
Kappa = 10 0 0
Edge_Parameters = 5 0 0
VertexArea = 0 0.7 0 0
TimeSeriesData_Period = 100
VertexPositionIntegrator = MetropolisAlgorithmOpenMP 1 1 0.05
AlexanderMove = MetropolisAlgorithmOpenMP 1
Dynamic_Box = IsotropicFrameTensionOpenMP 1 0 XY
VisualizationFormat = VTUFileFormat VTU_F 2000
NonbinaryTrajectory = TSI TrajTSI 2000
Restart_Period = 1000
You can also find this file here as "input_relaxG1PBC.dts". Here we use a frame tension of $\tau=0 [kT/l_{DTS}^2]$ to ensure that the membrane neck can relax in a tension-free setup. You can now start the simulation after making a .top file:
echo "gen1_PBC.q 2" > top.top:
$path/DTS -in input_relaxG1PBC.dts -top top.top
Please check your energy output files (dts-en.xvg) afterwards to see whether your system has equilibrated within the 20.000.000 steps (this is likely the case, it may even require fewer than 5.000.000 steps). The last simulation output provides you with a relaxed neck with which you can proceed in Part3. You should see be able to see the following development of your system, using a visualization software such as paraview:
Fig.2 Progression from initial configuration to a relaxed neck in PBC.
Part2
Core Bead Simulations: Neck Response in Stomatocytes under Internal Pressure
We now want to investigate the response of the membrane necks in a high-genus stomatocye to an increase in internal pressure, which manifests as an increase in lateral tension on the membrane. Achieving this is in our simulations is challenging. Unlike the well-defined volume within a triangulated surface, the internal compartment of the stomatocyte remains connected to the environment through the necks, hence there is no unambiguous definition of its volume. We can however control this internal volume in a different way: we simulate the stomatocyte containing a large bead that acts as an infinite potential barrier for the vertices, meaning any shape update that places a vertex inside the bead is rejected. This effectively creates a fixed volume for the stomatocyte’s internal compartment. During the simulation, this core bead is gradually expanded, representing a slow and continuous increase in the volume of the inner compartment, analogous to the effects of internal pressure increase.
Part2.1: Simple Membrane System
We start with the output of Part 1.1. You can copy the last .tsi file and rename it to "G20Stomatocyte_NoPro.tsi" (or get this file here), indicating that this is a simple membrane system without proteins. We will consider protein addition in Part 2.2. below.
The simulations are performed at near constant surface area by coupling to a harmonic stretching potential ($K_A=1000 kT$). We use a membrane bending rigidity of $\kappa=10 kT$, but you may change this at will to lower or higher values. Using the following input file, you can start the expanding core bead simulation. Please note that this is a non-equilibrium simulation, as we actively expand the core bead.
Integrator_Type = MC_Simulation
Min_Max_Lenghts = 1 3
MinfaceAngle = -0.5
Temperature = 1 0
Box_Centering_F = 0
Set_Steps = 1 5000000
EnergyMethod = FreeDTS1.0_FF
Kappa = 10 0 0
Edge_Parameters = 5 0 0
VertexArea = 0 0.7 0 0
TimeSeriesData_Period = 100
VertexPositionIntegrator = MetropolisAlgorithmOpenMP 1 1 0.05
VolumeCoupling = SecondOrder 0 0.0 0.0
TotalAreaCoupling = HarmonicPotential 1000 0.37
AlexanderMove = MetropolisAlgorithmOpenMP 1
Boundary = EllipsoidalCore 5 1 1 1
NonequilibriumCommands ExpandEllipsoidalCoreWall 100 0.005
VisualizationFormat = VTUFileFormat VTU_F 2000
NonbinaryTrajectory = TSI TrajTSI 2000
Restart_Period = 1000
You can also find this file here as "input_expCore.dts". You can now start the simulation:
$path/DTS -in input_expCore.dts -top G20Stomatocyte_NoPro.tsi
Before you draw final conclusions about the membrane necks' behaviour, you should actually equilibrate your system at certain core bead sizes, i.e. at certain points of the expanding core simulation. This is to avoid seeing only results of transient dynamics. Hence, we advise you to visually inspect your results, then chose files at certain increasing core bead radii. Here is an example of the progression you may see at $\kappa=20 kT$:
Fig.3 Stomatocyte with expanding core bead at $\kappa=20 kT$ and near-constant surface area (compressibility $K_A=1000 kT$)
If you used the exact parameters in the example input file above, radii of $5, 6, 7, 8, 8.5, 8.7, 9.0$ and $9.16 l_{DTS}$ are helpful. They then correspond to the .tsi files dts5.tsi, dts112.tsi, dts240.tsi, dts390.tsi, dts472.tsi, dts506.tsi, dts560.tsi and dts590.tsi, which you can also find here. With each of these .tsi files you should now start a simulation to equilibrate at their core bead radius using the following input:
Integrator_Type = MC_Simulation
Min_Max_Lenghts = 1 3
MinfaceAngle = -0.5
Temperature = 1 0
Box_Centering_F = 0
Set_Steps = 1 3000000
EnergyMethod = FreeDTS1.0_FF
Kappa = 10 0 0
Edge_Parameters = 5 0 0
VertexArea = 0 0.7 0 0
TimeSeriesData_Period = 100
VertexPositionIntegrator = MetropolisAlgorithmOpenMP 1 1 0.05
VolumeCoupling = SecondOrder 0 0.0 0.0
TotalAreaCoupling = HarmonicPotential 1000 0.37
AlexanderMove = MetropolisAlgorithmOpenMP 1
Boundary = EllipsoidalCore {} 1 1 1
VisualizationFormat = VTUFileFormat VTU_F 2000
NonbinaryTrajectory = TSI TrajTSI 2000
Restart_Period = 1000
And at {} fill in each of the exact values of bead radius:
5.049, 6.016, 6.999, 7.999, 8.496, 8.694, 8.999, 9.164
Frames from each of the now equilibrated systems can then be used to analyse neck diameters or the distance of inner and outer membrane. We provide our analysis routine here as a jupyter notebook "Analysis_NecksStomatocyte.ipynb".
Part2.2: Adding Proteins to the Stomatocyte
To see the effect of protein aggregations in the membrane necks, one can first add curvature preferring proteins to the system and equilibrate with these, before performing core bead simulations as above. To do so, follow these steps:
- Convert "G20Stomatocyte_NoPro.tsi" to a .q-file, to use the inbuilt features of FreeDTS for easy addition of proteins to the system:
$path/CNV -in G20Stomatocyte_NoPro.tsi -o G20Stomatocyte_NoPro.q
- Run a simulation with 10% curvature preferring proteins with $C_{0,P}=\frac{1}{l_{DTS}},C_{0,L}=\frac{-1}{l_{DTS}}$ and local impact on the bending rigidity of $\kappa_P=10 kT, \kappa_L=5 kT$, to let them assemble in the necks. The input file should contain:
Integrator_Type = MC_Simulation
Min_Max_Lenghts = 1 3
MinfaceAngle = -0.5
Temperature = 1 0
Box_Centering_F = 0
Set_Steps = 1 3000000
EnergyMethod = FreeDTS1.0_FF
Kappa = 10 0 0
Edge_Parameters = 5 0 0
VertexArea = 0 0.7 0 0
TimeSeriesData_Period = 100
VertexPositionIntegrator = MetropolisAlgorithmOpenMP 1 1 0.05
VolumeCoupling = SecondOrder 0 0.0 0.0
TotalAreaCoupling = HarmonicPotential 1000 0.37
AlexanderMove = MetropolisAlgorithmOpenMP 1
VisualizationFormat = VTUFileFormat VTU_F 1000
NonbinaryTrajectory = TSI TrajTSI 1000
Restart_Period = 1000
INCLUSION
GenerateInclusions
Selection_Type Random
TypeID 1 2
Density 0.1 0.0
Define 2 Inclusions
SRotation Type K KG KP KL C0 C0P C0L lambda lkg lkn cn0
0 Pro1 10 0 10 5 0.0 1 -1
Inclusion-Inclusion-Int
1 1 1 2 0.0 0.0
You can also find this file here as "input_G20_wPro.dts". You can now start the simulation after making a .top file:
echo "G20Stomatocyte_NoPro.q 2" > top.top:
$path/DTS -in input_G20_wPro.dts -top top.top
After checking that your system has been well equilibrated, you can rename the last .tsi output to "G20Stomatocyte_wPro.tsi" (or get it here) and use it to follow Part 2.1 - and see how proteins impact the behaviour of membrane necks under internal pressure in the stomatocyte. Remember to add the information about the proteins (inclusions) to the input file(s) given in Part 2.1., i.e. add the following lines at the end of the file(s):
INCLUSION
Define 2 Inclusions
SRotation Type K KG KP KL C0 C0P C0L lambda lkg lkn cn0
0 Pro1 10 0 10 5 0.0 1 -1
Inclusion-Inclusion-Int
1 1 1 2 0.0 0.0
Part3
Lateral Tension in PBC: Two-State Response of One Neck to Stretching
We now want to investigate the behaviour of one membrane neck in periodic boundary conditions under lateral tension to isolate the neck response. While in Part2 (whole stomatocyte) this tension was indirectly achieved through an internal pressure increase, we can now directly apply a frame tension on the simulation box in XY directions (Part 3.1 for a simple membrane, Part 3.2 with a protein complex). Moreover, to further characterise the response, one can perform constant frame area ensemble simulations $(N,A_p,T)$ for a wide range of frame areas. Because (positive) tension favours increasing $A_p$ this allows to track its effect on neck size in a more controlled manner, while sampling configurations that might be inaccessible in constant-tension simulations, where the area can expand too rapidly. Instructions are provided in Part 3.3.
Part 3.1: Frame Tension Simulations, Simple Membrane System
We continue with our output of Part1.2. You can copy the last .tsi file and rename it to "G1_PBC_NoPro.tsi" (or get this file here), indicating that this is a simple membrane system without proteins. We will consider protein addition in Part 3.2. below.
To induce a frame tension on the membrane, for example $\tau=1.5 \frac{kT}{l_{DTS}}$, your input file needs to be set up as below (you can ofcourse vary the value of $\tau$, and we encourage to explore the neck behaviour for a range of values).
Integrator_Type = MC_Simulation
Min_Max_Lenghts = 1 3
MinfaceAngle = -0.5
Temperature = 1 0
Box_Centering_F = 0
Set_Steps = 1 4000000
EnergyMethod = FreeDTS1.0_FF
Kappa = 10 0 0
Edge_Parameters = 5 0 0
VertexArea = 0 0.7 0 0
TimeSeriesData_Period = 100
VertexPositionIntegrator = MetropolisAlgorithmOpenMP 1 1 0.05
AlexanderMove = MetropolisAlgorithmOpenMP 1
Dynamic_Box = IsotropicFrameTensionOpenMP 1 1.5 XY
VisualizationFormat = VTUFileFormat VTU_F 2000
NonbinaryTrajectory = TSI TrajTSI 2000
Restart_Period = 1000
You can also find this file here as "input_tensionPBC_noPro.dts". You can now start the simulation(s):
$path/DTS -in input_tensionPBC_noPro.dts -top G1_PBC_NoPro.tsi
You will find that the neck size behaves differently below and above a certain threshold tension (constricting vs. dilating). If you are close to that tension, your simulation may get stuck in local minima and you might see none of your simulations dilating, even though you are above the threshold. Hence, for these systems it is helpful to set up a large amount of replica for each tension, e.g. 20 or more, to make sure you actually do or don't see dilation at that tension.
Fig.4 Tension dependent neck constriction and dilation for membranes at $\kappa=20 kT$ and a starting neck diameter of $D_{start}\approx 8 l_{DTS}$
Part 3.2: Frame Tension Simulations, Protein Addition
To simulate a protein complex and its influence on the necks' tension response, we proceed similarly as in the whole stomatocyte. Before applying tension to the membrane, we first add proteins to the relaxed membrane neck. Again we use the output of Part1.2. You can copy the last .tsi file and rename it to "G1_PBC_NoPro.tsi" (or get this file here). We recommend 2.5-5% of inclusions to cover the neck.You may use $C_{0,P}=\frac{-0.8}{l_{DTS}}$ and local impact on the bending rigidity of $\kappa_P=10 kT, \kappa_L=5 kT$, to let them assemble in the neck. To form a protein complex of varying rigidity you may use attractive protein-protein interactions with $A=1 kT$ and $B=0,-1,-2,-4 kT$ (more negative $B$ value leads to a more rigid complex). To do so, follow these steps:
- Convert "G1_PBC_NoPro.tsi" to a .q-file, to use the inbuilt features of FreeDTS for easy addition of proteins to the system:
$path/CNV -in G1_PBC_NoPro.tsi -o G1_PBC_NoPro.q
- Run a simulation to let the protein complex assemble in the neck. Your input file for one of these simulations needs to look like this (remember to adapt the protein percentage and $B$, which is the very last value -bottom right- in the file):
Integrator_Type = MC_Simulation
Min_Max_Lenghts = 1 3
MinfaceAngle = -0.5
Temperature = 1 0
Box_Centering_F = 0
Set_Steps = 1 4000000
EnergyMethod = FreeDTS1.0_FF
Kappa = 10 0 0
Edge_Parameters = 5 0 0
VertexArea = 0 0.7 0 0
TimeSeriesData_Period = 100
VertexPositionIntegrator = MetropolisAlgorithmOpenMP 1 1 0.05
AlexanderMove = MetropolisAlgorithmOpenMP 1
Dynamic_Box = IsotropicFrameTensionOpenMP 1 0 XY
VisualizationFormat = VTUFileFormat VTU_F 2000
NonbinaryTrajectory = TSI TrajTSI 2000
Restart_Period = 1000
INCLUSION
GenerateInclusions
Selection_Type Random
TypeID 1 2
Density 0.05 0.0
Define 2 Inclusions
SRotation Type K KG KP KL C0 C0P C0L lambda lkg lkn cn0
0 Pro1 10 0 10 5 0.0 1 0
Inclusion-Inclusion-Int
1 1 1 2 1.0 -2.0
You can also find this file here as "input_G1PBC_wPro.dts". You can now start the simulation(s) after making a .top file:
echo "G1_PBC_NoPro.q 2" > top.top
$path/DTS -in input_G1PBC_wPro.dts -top top.top
After checking that your system has been well equilibrated, you can rename the last .tsi output to "G1_PBC_wPro.tsi" (or get it here) and use it to follow Part 3.1 - and see how proteins impact the behaviour of a membrane neck under tension. Remember to add the information about the proteins (inclusions) to the input file(s) given in Part 2.1., i.e. add the following lines:
INCLUSION
Define 2 Inclusions
SRotation Type K KG KP KL C0 C0P C0L lambda lkg lkn cn0
0 Pro1 10 0 10 5 0.0 1 0
Inclusion-Inclusion-Int
1 1 1 2 1.0 -2.0
Fig.5 Curvature preferring proteins (5% membrane coverage, $C_{0,P}=-0.8 l_{DTS}^{-1}$, $\kappa_P=10 kT$,$\kappa_L=5 kT$) on a membrane neck.
Part 3.3: Constant Frame Area Simulations
Using states from a constant tension simulation in which the neck dilates, one can set up an array of simulations each at a fixed frame area corresponding to a certain state. (Positive) tension favours increasing frame area $A_p$. Therefore, performing constant $A_p$ simulations for a wide range of frame areas allows you to track its effect on neck size in a more controlled manner, and see configurations that might be inaccessible in constant-tension simulations, where the area can expand too rapidly. The workflow can be visualized as follows:
Fig.6 Principle of constant frame area ($A_p$) simulations.
To simulate with constant frame area instead of tension, just delete this line of the input file ("input_tensionPBC_noPro.dts"):
Dynamic_Box = IsotropicFrameTensionOpenMP 1 0 XY
A simulation of a given topology file (e.g. G1_PBC_NoPro.tsi) will then run with the fixed box side lengths of that topology file.
We provide an exemplary set of states as .tsi files here, which cover a range of scaled frame area $0.8<2A_p/A_0<1.3$. To uncover the details of the tension response with this approach, it is useful to set up 10 (or more) replica for each state, simulate it at fixed frame area, then analyse all replica for bending energy, neck diameters and membrane distance at all frame areas/states to obtain results as shown below. You may vary the bending rigidity (e.g. $\kappa=5,10,20 kT$) and also area compressibility (e.g. $K_A=0,100,1000 kT$) to see their effect. Exemplary results at $K_A=1000 kT$, i.e. almost constant surface area $A_0$:
Fig.7 Bending energy $E_B$ over frame area $A_p$, scaled by constant surface area (compressibility $K_A=1000 l_{DTS}$).
We provide an exemplary python script to set up such an array of simulations and analyse the results here).