Describing Membranes As Meshes

Overview

Questions

  • What is the Triangulated Mesh Model?

  • How can I create a mesh modeling a flat surface?

  • How can I create a mesh modeling a sphere?

Objectives

  • Introducing the Triangulated Mesh Model.

  • Create a triangulation of the xy-plane.

  • Create a triangulation of a sphere surface. ## Boilerplate Code

[1]:
import numpy
from scipy.spatial import ConvexHull

Triangulated Mesh Model

Modeling membranes and other deformable interfaces as triangulated meshes is a common coarse-grained approach in molecular simulations, particularly for studying lipid bilayers and biological membranes at macroscopic scales. The triangulated mesh model represents the interface surface as a set of vertices connected by edges forming a mesh of triangles. Each vertex represents a coarse-grained membrane patch (e.g., a group of lipids). This approach allows users to integrate membranes with specific elasticities and interactions with other particles.

Creating a triangulation of the xy-plane.

To create a triangulation of the xy-plane the positions of the vertices have to defined first.

[2]:
# box size
L = 10

# Number of vertices per row
row = 10

# Number of vertices
N = row * row

# Create meshgrid for a rectangular pattern
x = numpy.arange(-L / 2, L / 2, L / row)
y = numpy.arange(-L / 2, L / 2, L / row) * numpy.sqrt(3) / 2
X, Y = numpy.meshgrid(x, y)

# Offset every second row to create a triangular pattern
X[1::2, :] += 0.5

# Combine X and Y coordinates to define vertex positions
vertex_positions = numpy.column_stack([X.flatten(), Y.flatten(), numpy.zeros(N)])

print(vertex_positions[:10, :])
print("...")
[[-5.         -4.33012702  0.        ]
 [-4.         -4.33012702  0.        ]
 [-3.         -4.33012702  0.        ]
 [-2.         -4.33012702  0.        ]
 [-1.         -4.33012702  0.        ]
 [ 0.         -4.33012702  0.        ]
 [ 1.         -4.33012702  0.        ]
 [ 2.         -4.33012702  0.        ]
 [ 3.         -4.33012702  0.        ]
 [ 4.         -4.33012702  0.        ]]
...

In addition to the vertex positions also the triangulation data is needed. The mesh information is given by an array of triplets which define each triangle in the mesh and consist of the vertex indices.

[3]:
# Create array of vertex indices
v_index = numpy.arange(N)

# Find index of vertex to the right in same row
vr_index = v_index + 1

# Consider periodic boundaries
vr_index[row - 1 :: row] -= row

# Find index of vertex above and below the edge between idx1 and idx2
e_index = numpy.repeat(v_index, 2).reshape(-1, 2)

# Find index of vertex above the edge between v_index and vr_index
e_index[:, 0] += row
# Find index of vertex below the edge between v_index and vr_index
e_index[:, 1] -= row

# Consider periodic boundaries
e_index = e_index.reshape(row, -1, 2)
# periodic boundaries in x direction
e_index[1::2] += 1
# periodic boundaries in y direction
e_index[1::2, -1] -= row
e_index = e_index.reshape(-1) % N

# Create duplicates of the vertices to create triangle above and below each edge
v_index = numpy.repeat(v_index, 2)
vr_index = numpy.repeat(vr_index, 2)

# Switch vr_index and e_index for triangle below each edge to make mesh orientable
vr_index[1::2] = e_index[1::2]
e_index[1::2] = vr_index[0::2]

# Create all triangles making up the mesh
triangles = numpy.column_stack([v_index, vr_index, e_index])

print(triangles[:5, :])
print("...")
[[ 0  1 10]
 [ 0 90  1]
 [ 1  2 11]
 [ 1 91  2]
 [ 2  3 12]]
...
Plane-Triangulation

Creating a triangulation of a sphere surface.

First, we distribute vertices on a sphere using the Fibonacci sphere algorithm.

[4]:
# Number of points on the sphere
num_pts = 1000

# Radius of the sphere
R = 10

# Fibonacci algorithm
indices = numpy.arange(0, num_pts, dtype=float) + 0.5

phi = numpy.arccos(1 - 2 * indices / num_pts)
theta = numpy.pi * (1 + 5**0.5) * indices

x, y, z = (
    R * numpy.cos(theta) * numpy.sin(phi),
    R * numpy.sin(theta) * numpy.sin(phi),
    R * numpy.cos(phi),
)

vertex_positions = numpy.column_stack((x, y, z))

Decoration-Sphere We determine the triangles of the mesh by computing the convex hull on the vertex positions

[5]:
hull = ConvexHull(vertex_positions)

triangles = hull.simplices

# calculate normal of each triangle and orient mesh so all triangle normals
# point away from the center
for i in range(len(triangles)):
    normal = numpy.cross(
        vertex_positions[triangles[i, 0]] - vertex_positions[triangles[i, 2]],
        vertex_positions[triangles[i, 0]] - vertex_positions[triangles[i, 1]],
    )
    if numpy.dot(vertex_positions[triangles[i, 0]], normal) < 0:
        triangles[i, 1], triangles[i, 2] = triangles[i, 2], triangles[i, 1]

print(triangles)
[[459 438 493]
 [404 438 459]
 [178 199 144]
 ...
 [703 669 648]
 [703 648 682]
 [703 682 737]]
Triangulation-Sphere