ABSTRACT

Title of Dissertation: SIMULATING MEMBRANE-BOUND
CYTOSKELETAL DYNAMICS

Haoran Ni
Doctor of Philosophy, 2023

Dissertation Directed by: Professor Garegin A. Papoian
Department of Chemistry and Biochemistry

The cell membrane defines the shape of the cell and plays an indispensable role in bridging

intra- and extra-cellular environments. The membrane, consisting of a lipid bilayer and various

attaching proteins, mechanochemically interacts with the active cytoskeletal network that dynam-

ically self-organizes, playing a vital role in cellular biomechanics and mechanosensing. Compre-

hensive simulations of membrane-cytoskeleton dynamics can bring insight in understanding how

the cell mechanochemically responds to external signals, but a computational model that cap-

tures the complex cytoskeleton-membrane with both refined details and computational efficiency

is lacking. To address this, we introduce in this thesis a triangulated membrane model and in-

corporate it with the active biological matter simulation platform MEDYAN (“Mechanochemical

Dynamics of Active Networks”). This model accurately captures the membrane physical prop-

erties, showing how the membrane rigidity, the structure of actin networks and local chemical

environments regulate the membrane deformations. Then, we present a new method for simulat-

ing membrane proteins, using stochastic reaction-diffusion sampling on unstructured membrane



meshes. By incorporating a surface potential energy field into the reaction-diffusion sampling,

we demonstrate rich membrane protein collective behaviors such as confined diffusion, liquid-

liquid phase separation and membrane curvature sensing. Finally, in order to capture stretching,

bending, shearing and twisting of actin filaments which are not all available with traditional acto-

myosin simulations, we introduce new finite-radius filament models based off Cosserat theory of

elastic rods, with efficient implementation using finite-dimensional configurational spaces. Us-

ing the new filament models, we show that the filaments’ torsional compliance can induce chiral

symmetry breaking in a cross-linked actin bundle. All the new models are implemented in the

MEDYAN platform, shedding light on whole cell simulations, paving way for a better under-

standing of the membrane-cytoskeleton system and its role in cellular dynamics.



SIMULATING MEMBRANE-BOUND CYTOSKELETAL DYNAMICS

by

Haoran Ni

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy

2023

Advisory Committee:
Professor Garegin A. Papoian, Chair/Advisor
Professor Arpita Upadhyaya
Professor Pratyush Tiwary
Professor Sergei Sukharev
Professor Jeffery B. Klauda



© Copyright by
Haoran Ni

2023



Acknowledgments

I owe my gratitude to all the people who have made this thesis possible and because of

whom my graduate experience in the past 6.5 years has been one that I will cherish forever. First

and foremost I’d like to thank my advisor, Garegin Papoian. I have gained a lot by discussing with

you about scientific, research and coding topics, and more importantly the way of prioritizing

things as a researcher. I am grateful to be mentored by such an inspiring scientist. I would also

like to thank my committee members Arpita Upadhyaya, Pratyush Tiwary, Sergei Sukharev and

Jeffery Klauda. I appreciate all your great help in guiding this thesis and the opportunity to learn

from you.

I would extend my thanks to the Papoian lab members, who have been colleges and friends

along the way. Haiqing Zhao, Hao Wu, Mary Pitman, Ryan Nival, as well as the MEDYAN de-

velopers James Komianos, Aravind Chandrasekaran, Carlos Floyd, Qin Ni, Nathan Zimmerberg

and Mengxin Gu have made the my lab experience joyful. It has been a great pleasure to talk

about scientific problems, debug erroneous codes or chat about life in general. I have learned a

lot from you and am grateful for all the helps I received, thank you.

I would also like to acknowledge and thank the University of Maryland organizations I

have been a member of, including the Biophysics program, the Institute for Physical Science

and Technology and the Physics of Living Systems network. The faculty and students in these

organizations have broaden my views with research in a wide range of disciplines and I am

ii



grateful for inspiring discussions with people outside my lab. I would also like to thank the staff

members in these organizations for delivering excellent services.

Finally, I owe my deepest thanks to my family. To my parents, Yongping Ni and Xunjuan

Qiu, my grandparents and my cousins, thank you for being my constant source of strength, even

you are overseas thousands of miles away. To all my friends along the way, thank you all.

iii



Table of Contents

Acknowledgements ii

Table of Contents iv

List of Tables vii

List of Figures viii

List of Abbreviations x

Chapter 1: Introduction 1
1.1 Biological membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Cytoskeleton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Cytoskeleton-membrane interactions . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Outline of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Chapter 2: Simulating deformable vesicles containing complex cytoskeletal networks 12
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.2.1 MEDYAN Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.2 Membrane mechanics and discretization . . . . . . . . . . . . . . . . . . 15
2.2.3 Adaptive surface remeshing . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2.4 Interactions between the membrane and the cytoskeleton . . . . . . . . . 26
2.2.5 Computational efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.6 Code availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Chapter 3: Simulating membrane proteins and their dynamics 39
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.2.1 Thermodynamics of interacting proteins on the surface . . . . . . . . . . 41
3.2.2 Surface discretization and stochastic diffusion sampling . . . . . . . . . . 43
3.2.3 Membrane-related reactions . . . . . . . . . . . . . . . . . . . . . . . . 46
3.2.4 Free energy of the particle system . . . . . . . . . . . . . . . . . . . . . 48
3.2.5 Membrane curvature sensing and generation model . . . . . . . . . . . . 49

3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

iv



3.3.1 Energy barriers acting as protein confinements decrease diffusivity . . . . 50
3.3.2 Pairwise protein interactions promote LLPS . . . . . . . . . . . . . . . . 53
3.3.3 Concentration sorting by local curvature . . . . . . . . . . . . . . . . . . 54

3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Chapter 4: Cosserat rod model of actin filaments 61
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.2.1 Cosserat theory of elastic rods . . . . . . . . . . . . . . . . . . . . . . . 66
4.2.2 Variational approach to rod mechanics . . . . . . . . . . . . . . . . . . . 69
4.2.3 Geodesic models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.2.4 MEDYAN model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.2.5 Dynamical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.2.6 Binding to surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.3.1 Model comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.3.2 Bundle study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.3.3 MEDYAN Implementation . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Chapter 5: Summary and future directions 93

Appendix A: Supporting information for Chapter 2 96
A.1 Simulation parameters and initial setups . . . . . . . . . . . . . . . . . . . . . . 96

A.1.1 Simulations for interactions between actin and membrane vesicle . . . . . 96
A.1.2 Perpendicular bundle protrusion against the membrane . . . . . . . . . . 97
A.1.3 Non-perpendicular bundle protrusion against the membrane . . . . . . . . 97

A.2 Representation of membrane surface using a triangular mesh . . . . . . . . . . . 98
A.3 Benchmark on membrane adaptive remeshing . . . . . . . . . . . . . . . . . . . 100
A.4 Membrane bending energy on discretized surface . . . . . . . . . . . . . . . . . 100

A.4.1 Edge based curvature and bending energy estimation . . . . . . . . . . . 101
A.4.2 Vertex based curvature and bending energy estimation . . . . . . . . . . 103
A.4.3 Comparing curvature estimation algorithms . . . . . . . . . . . . . . . . 104

A.5 Derivation of triangle-bead volume exclusion interaction . . . . . . . . . . . . . 108
A.5.1 Evaluating the integral . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
A.5.2 Geometric interpretation of Eq. A.17 . . . . . . . . . . . . . . . . . . . . 111
A.5.3 Resolving numerical issues for coplanar cases . . . . . . . . . . . . . . . 115

A.6 Performance of bundle protrusion simulations . . . . . . . . . . . . . . . . . . . 119

Appendix B: Supporting information for Chapter 3 122
B.1 Estimating the relationship between mesh size and the maximum gradient of energy122
B.2 Membrane protein cluster nucleation . . . . . . . . . . . . . . . . . . . . . . . . 125
B.3 Membrane bending energy with protein inclusion . . . . . . . . . . . . . . . . . 127

Appendix C: Supporting information for Chapter 4 129

v



C.1 Euler buckling study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
C.2 Description of test cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
C.3 Strain profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
C.4 Parameterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
C.5 Q(ŝ) and r(ŝ) in the geodesic models . . . . . . . . . . . . . . . . . . . . . . . . 136
C.6 Energies in the geodesic models . . . . . . . . . . . . . . . . . . . . . . . . . . 138
C.7 List of symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

Bibliography 145

vi



List of Tables

3.1 Membrane diffusion barrier simulation parameters. . . . . . . . . . . . . . . . . 52
3.2 Protein LLPS simulation parameters. . . . . . . . . . . . . . . . . . . . . . . . . 54
3.3 Curvature-induced protein sorting simulation parameters. . . . . . . . . . . . . . 57

4.1 The difference metrics Cr and Cd between each model and the dynamical model. . 80
4.2 The equilibrated filament energies E, external energies from bound linkers Eext,

and total energy Etot for the three test cases. . . . . . . . . . . . . . . . . . . . . 80

A.1 The parameters used in simulation. . . . . . . . . . . . . . . . . . . . . . . . . . 96
A.2 Parameters for vesicle bundle deformation simulation. . . . . . . . . . . . . . . . 97
A.3 Parameters for filament bundle polymerization with different contact angles with

the membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
A.4 Time consumption of 0.2 s simulation time window when N = 7. Simulating

16.1 s in total took 1 d 6 h 22 min. . . . . . . . . . . . . . . . . . . . . . . . . . . 120
A.5 Time consumption of 0.2 s simulation time window when N = 10. Simulating

10.4 s in total took 3 d. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
A.6 Time consumption of 0.2 s simulation time window when N = 14. Simulating

4.01 s in total took 3 d. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

C.1 The deviation of the fitted prefactor a from 1 is shown for the five variational
models tested. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

vii



List of Figures

1.1 Existing models that involve membrane curvature generation in various cellular
processes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2 Schematic representation of the cell with various cytoskeleton structures. . . . . . 6
1.3 Membrane regulates actin network dynamics. . . . . . . . . . . . . . . . . . . . 8

2.1 Membrane manifold representation and deformation modes. . . . . . . . . . . . 17
2.2 An illustration of membrane as a chemical regulator in 2D space. . . . . . . . . . 28
2.3 Evolution of the vesicle shape is shown under increasingly hyperosmotic conditions. 31
2.4 Vesicle deformation by actin filaments crucially depends on the spatial organiza-

tion of actin filaments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.5 Filament bundle deforming the membrane. . . . . . . . . . . . . . . . . . . . . . 34
2.6 Membrane protrusion behaviors under various membrane rigidities. . . . . . . . 36

3.1 Geometric elements used in calculating a jk. . . . . . . . . . . . . . . . . . . . . 45
3.2 Membrane species diffusion with a barrier. . . . . . . . . . . . . . . . . . . . . . 52
3.3 Membrane protein phase separations. . . . . . . . . . . . . . . . . . . . . . . . . 55
3.4 Surface protein concentration sorting for membrane vesicles with various radii. . 58

4.1 The quantities used to specify a thin rod’s configuration in the Cosserat theory. . . 65
4.2 The cross-section of a filament with a bound cross-linker. . . . . . . . . . . . . . 76
4.3 A comparison of the variational and dynamical models. . . . . . . . . . . . . . . 81
4.4 The structural details of the chiral bundle simulation. . . . . . . . . . . . . . . . 83
4.5 Results from the chiral bundle simulation. . . . . . . . . . . . . . . . . . . . . . 86
4.6 A visualization of a contractile actomyosin network from a MEDYAN simulation. 89

A.1 Examples of fans around a vertex. . . . . . . . . . . . . . . . . . . . . . . . . . 99
A.2 The remeshing algorithm generates anisotropic mesh suitable for simulation. . . . 101
A.3 Ways of estimating the area at a vertex. . . . . . . . . . . . . . . . . . . . . . . . 105
A.4 Result of curvature estimate benchmark on a red blood cell. . . . . . . . . . . . . 107
A.5 The area element on a triangle facet. . . . . . . . . . . . . . . . . . . . . . . . . 109

B.1 Two equilateral triangles under a uniform energy gradient. . . . . . . . . . . . . 123
B.2 Large membrane mesh size reduces the accuracy. . . . . . . . . . . . . . . . . . 124
B.3 Snapshot of membrane protein LLPS with too large mesh size. . . . . . . . . . . 125
B.4 The area concentration profile among membrane cells with the phase-separated

initial state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

viii



C.1 The Euler critical buckling force Fc for various filament lengths L̂ for five varia-
tional models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

C.2 A comparison of strain profiles along an equilibrated filament in different models. 134

ix



List of Abbreviations

ALPS Amphipathic lipid packing sensor
Arp2/3 Actin-related protein 2/3 complex
ATP Adenosine triphosphate
BCR B cell receptor
GC Geodesic Cosserat
GEK Geodesic extensible Kirchoff
GTP Guanosine-5’-triphosphate
MEDYAN Mechanochemical dynamics of active networks
MRR Membrane-related reactions
NMII Non-muscle myosin II
PLC Phosphorylated phospholipase C
VASP Vasodilator-stimulated phosphoprotein
WASP Wiskott-Aldrich syndrome protein
WAVE WASP-family verprolin-homologous protein

x



Chapter 1: Introduction

The cell is the fundamental unit of life, and its shape is defined by the plasma membrane.

The membrane guards the frontier of the cell, allowing the cell to communicate with its sur-

rounding environment [1]. Inside the cell there is a ubiquitous polymeric network known as

the cytoskeleton. Crosslinked filament networks serve as a scaffold stiffening the cell [2]. With

chemical energy supplies, the polymeric networks actively self-reorganize into various structures,

converting chemical energy to mechanical work [3,4]. Through active remodeling, the cytoskele-

ton can exert forces on the membrane, initiating the cell’s mechanical response. On the other

hand, mechanical or chemical cues transmitted from the membrane can affect this cytoskeleton

reorganization processes, complementing the cell’s mechanosensing process [5, 6]. Understand-

ing how cytoskeleton works with the membrane to help the cell adapt to external environment is

crucial in understanding how cell works.

Computer simulations provide a powerful way to understanding such complex systems

[7]. Knowledge of individual components in the system can be obtained from past theoretical,

computational and experimental studies, and is used as the inputs of the simulation programs.

The simulation programs then evolve the system based on fundamental physics and chemical

principles, reproducing mesoscopic emergent behaviors that are difficult to obtain from analytical

studies. Moreover, intermediate information that is generally inaccessible in experiments could

1



be easily extracted from computer simulations, providing more insights into the complex systems

[4,8]. Over last two decades, efforts on simulation software on both the cytoskeleton systems and

the membranes greatly enhanced understanding of these systems. With past advances, we aim to

build a computation model that captures the essential mechanochemical properties of the entire

cell, allowing cell dynamics predictions.

In this thesis, we bring new simulation capabilities to the active matter simulation platform

MEDYAN (Mechanochemical Dynamics of Active Networks) [9], including a new cytoskeleton-

interacting membrane model with membrane protein dynamics. We also introduce a more ac-

curate Cosserat-rod-based filament model, and use it to investigate symmetry breaking in actin

bundles. In the remainder of the introduction, we first discuss the cytoskeleton and the membrane

systems in greater detail, including their biological functions and existing models to describe

them. Next, we discuss the cytoskeleton-membrane interactions and the efforts to simulate these

processes. Finally, we will briefly outline the subsequent chapters in the thesis.

1.1 Biological membrane

The membrane defines the boundary of a cell as well as many important organelles such

as the mitochondria and the nucleus [10]. The membrane consists of a thin lipid bilayer as the

main structural component as well as proteins attaching on the bilayer [11]. Mechanically, the

membrane resembles a thin fluid film, resisting in-plane and out-of-plane deformations, which,

in turn, influences its cellular functions. For example, membrane tension has been shown to

influence exocytosis, endocytosis and cell motility [12, 13]. The membrane proteins can diffuse

on the bilayer due to its in-plane fluidity [14]. The membrane proteins can be categorized into

2



a variety of classes, carrying out important cellular functions. For example, ion channels are

transmembrane proteins that control flows of specific ions across the membrane and, in turn,

control cell volume, action potential, exocytosis and many others [15]. Membrane proteins can

also undergo clustering via 2D liquid-liquid phase separation (LLPS), a thermodynamic effect

arising from intermolecular attraction [16–19]. Computation models for biological membranes

spread across multiple orders of magnitude in terms of length scale and resolution, ranging from

all-atom simulations to micrometer-scale meshless or manifold representations [20].

Membrane curvature generation processes have profound impact on cellular processes such

as membrane trafficking and endocytosis [21, 22]. Many different membrane proteins influ-

ence membrane curvature with various pathways [23]. For example, the BIN-Amphiphysin-

Rvs (BAR) domain proteins locally curve the membrane by both adapting the lipid bilayer to

its intrinsic curvature [24], and by inserting hydrophobic helices into the lipid layer [23, 25].

Conversely, many curvature generating proteins are known to sense the membrane curvature as

well. In the BAR domain protein example, the helix insertion was also identified to be the key

factor for curvature sensing processes [25]. Advances in computer modeling have enabled bet-

ter understanding of mechanisms underlying membrane curvature sensing and generation (see

Fig. 1.1) [23, 26–29]. However, how the curvature sensing and generating effects synergistically

reinforce each other and modulate the membrane shape is still not well understood [23].

1.2 Cytoskeleton

Cytoskeleton exists ubiquitously in eukaryotic cells, consisting of interconnected poly-

mer networks with regulatory proteins [1]. The cytoskeletal networks can actively undergo non-

3



Figure 1.1: Existing models that involve membrane curvature generation in various cellular pro-
cesses. (A) Membrane budding process during endocytosis [26, 30]. (B) Tubular membrane
structure formation and stabilization [28, 31]. (C) Nuclear envelope topology change [32]. (D)
Membrane deformation by caveolae [33, 34]. (E) Membrane protrusion by actin force [27, 35].
(F) Mitochondrial fission [36, 37]. Image reprinted from Ref. [23] from Molecular Diversity
Preservation International under Creative Common Attribution License.

4



equilibrium reorganization, powered by chemical energies [38]. For example, with ATP supplies,

an actin filament — a major type of polymer in the cytoskeleton — can have net growth on the

barbed and and net destruction on the pointed end at steady state [39,40]. The energy-consuming

polymerization of the filaments can also be converted to useful work under external loads [41].

Another example is the non-muscle myosin II (NMII), which a crosslinking protein with heads

that can bind to actin filaments. NMII heads can consume ATP to move along actin filaments [42],

which may create stress in cytoskeleton and lead to network contraction [9, 43]. With highly dy-

namic components, the cytoskeletal networks can self-organize into various structures (Fig. 1.2),

and control important cellular activities such as cell stiffening, migration, transport and mechan-

otransduction [38, 44–48].

Computer simulation packages have been developed for mechanochemical models of cy-

toskeletal networks. Examples include MEDYAN [9], Cytosim [50], AFiNeS (Active Filament

Network Simulation) [51], etc. In this work, we base our work on the MEDYAN platform. In

MEDYAN, diffusing chemical species in each compartment are considered “well-mixed”, but

the locations for force bearing species (such as actin filaments and bound myosin motors) are

explicitly modeled and tracked [9]. MEDYAN simulations are based on iteratively interleaving

stochastic draws of chemical events, using the Gillespie algorithm, followed by a mechanical

equilibration step. Importantly, chemistry influences mechanics and vice versa, creating intricate

mechanochemical feedbacks.

In MEDYAN, polymeric filaments are modeled as straight, compressible cylinders con-

nected at hinge points [9]. Since these filaments have high aspect ratio (an actin filament can

reach several micrometers in length and has only ∼ 7nm in diameter [52]), this simplified ap-

proach captures salient filament properties such as stretching and bending rigidities, and is com-

5



Figure 1.2: Schematic representation of the cell with various cytoskeleton structures. (i) The
cell cortex. (ii) A contractile stress fiber. (iii) The lamellipodium. (iv) The filopodium. Image
reprinted from Ref. [49] with permission from the American Physiological Society.

6



putationally cheap enough to allow for micron-scale actomyosin simulations to reach a timescale

of ∼1 hour. However, this model neglects finite width of the filaments and omits the chiral

asymmetry of helical actin filaments [53], so simulations using it cannot reproduce nor explain

emergent chiral symmetry broken behavior observed in experiments [54]. To address this is-

sue, in this work, we present new filament models in MEDYAN based on Cosserat rod theory,

capturing all allowed deformation modes, while minimizing the impact on computation times.

1.3 Cytoskeleton-membrane interactions

A complete cellular mechanochemical sensing process often involves the membrane, a

mostly inactive flexible surface, cooperating with the cytoskeleton, a highly active dynamic poly-

mer network [55–57]. Membranes can interact with the cytoskeleton mechanically via repulsion

and adhesion [58]. Cells with the plasma membrane attached to a dense cross-linked actin cor-

tex have their boundaries significantly rigidified [12]. On the contrary, disruption of filament-

membrane attachments leads to cellular blebbing [59,60]. Membrane proteins also participate in

numerous signaling processes, regulating, in particular, many cytoskeletal interactions (Fig. 1.3).

For instance, the membrane attached WASP-VASP complex regulates Arp2/3 activation, which

is an actin branching protein essential in lamellipodia formation [61].

From a different perspective, membrane receptor proteins and their spatial organizations

play an indispensable role in cellular signaling processes, which modulate, for example, the

formation of immunological synapse [63, 64]. Clustering of receptors may alter the signaling

cascade and affect cytoskeleton dynamics. For example, transmembrane complex nephrin-NCK-

N-WASP clustering regulates actin assembly [65, 66]. Reciprocally, the cytoskeleton dynamics

7



Figure 1.3: Membrane regulates actin network dynamics. (A) Membrane lipids affect profilin
activity. Profilin preferably binds to PI(4,5)P2, allowing for more free actin monomers locally,
while phosphorylated phospholipase C (PLC) hydrolyzes PI(4,5)P2 and allows for more profilin-
bound actin monomers locally. Free actin monomers preferably nucleate with Arp2/3 complex
to form branched networks, while profilin-bound actin monomers prefer formin-mediated fila-
ment polymerizations. (B) Small activated GTPases of the Rho superfamily proteins attach on
the membranes, activating WASP/WAVE factors, which modulates Arp2/3 activity. F-BAR pro-
teins interact with WASP, regulating actin polymerization. These interactions lead to processes
described in the text boxes. (C) Small activate GTPases of the Rho superfamily on the membrane
can bind to formins, activating actin polymerization. F-BAR proteins also regulate formin activ-
ities, which, in turn, modulates actin activities and drives membrane protrusion and cytokinesis.
Formins binding directly to the membrane can be found in eukaryotes that lack formins with
obvious Rho-binding domains, such as plants. Image reprinted from Ref. [62] with permission
from Rockefeller University Press.

8



in membrane’s vicinity can influence the distribution of membrane proteins as well, facilitating

or hindering their clustering. For example, the spatial distribution of the transmembrane recep-

tor CD44 is guided by the underlying cytoskeleton dynamics [67]. Moreover, these complex

interactions may further deform the membrane into highly curved shapes, and, in turn, affect

the dynamics of curvature-sensing proteins. Specifically, How membrane-actomyosin dynamics,

driven by far-from-equilibrium processes, and the passive intermolecular interactions synergisti-

cally regulate the aggregation and disassembly of surface receptors remains a central question at

the interface of physical chemistry and cell biology. A simulation framework that captures these

complex surface-receptor/membrane/cytoskeleton interactions may greatly improve the under-

standing of mechanism of surface protein phase behavior.

Efforts of using computational models to study cell periphery dynamics involving membrane-

cytoskeleton interactions have been made in the past decades. For example, both analytical and

computational methods have been used to study filopodia and lamellipodia systems [68–70]. A

recent work combined branched actin network simulation with hemisphere-shaped membrane

endocytic pit, and emphasized cytoskeleton’s role in initiating endocytosis [71]. However, it is

still desirable to develop a detailed, flexible, 3D membrane model which interacts with cytoskele-

tons. In this thesis, we will introduce our new membrane model in MEDYAN, capturing salient

properties such as membrane surface proteins and membrane-filament interactions. This work

enables cytoskeleton mechanochemical simulations enclosed in a movable membrane domain,

paving way for whole cell modeling.

9



1.4 Outline of Thesis

In Chapter 2, we present a new computational model for simulating membranes and their

interactions with the cytoskeleton. We introduce a mesh-based membrane model in MEDYAN,

capturing the membrane’s elastic properties, membrane-filament steric interactions and volumet-

ric effects such as the osmotic pressure. The membrane-MEDYAN simulations reveal how the

growing actin bundles generate tubular membrane protrusions. We discuss our investigations

on the protrusion initiation and dynamics under various filament bundle geometries, membrane

rigidities and local G-Actin concentrations. Using simulations, we found that the bundles’ protru-

sion propensities depend on the synergy between bundle thickness and inclination angle at which

the bundle approaches the membrane. The new membrane model paves the way for whole-cell

simulations using MEDYAN.

In Chapter 3, we present a novel way to simulate surface protein dynamics using a compartment-

based reaction-diffusion scheme on membrane meshes. We show how reaction rates and diffusion

coefficients for proteins in a potential field can be translated to reaction/hopping propensities on

unstructured 2D meshes. Then, we show membrane proteins undergoing thermodynamics driven

LLPS, which can be studied using this model by incorporating the protein inter-molecular in-

teractions. Moreover, utilizing the curvature-dependent membrane binding energy in the protein

potential field, we demonstrate curvature-induced protein curvature sorting on membrane vesi-

cles with various sizes and analyze the sorting strength under various membrane protein elastic

properties.

In Chapter 4, we introduce new variational models for cytoskeletal filaments. In order to

simulate large filament networks that incorporate stretching, shearing and twisting of filaments

10



with finite radius, we propose new filament models based on the Cosserat rod theory. We compare

the cylinder-bead filament model in the original MEDYAN, a polynomial spline-based Cosserat

filament model, as well as a Cosserat model where geodesic curves in the configurational space is

used, and find that the last model both captures all the relevant degrees of freedom and allows for

arbitrarily large deformations without systematic error. We implement new procedures necessary

for these models to work properly with MEDYAN, and show that the new models reveal chiral

symmetry breaking for filament bundles with torsional compliance.

In chapter 5, we provide the summary of the work in this thesis and discuss future directions

of research.

11



Chapter 2: Simulating deformable vesicles containing complex cytoskeletal net-

works

This chapter and Appendix A are adapted with permission from: Haoran Ni and Garegin

A. Papoian (2021) Membrane-MEDYAN: simulating deformable vesicles containing complex cy-

toskeletal networks. The Journal of Physical Chemistry B, 125(38), 10710-10719. Copyright

2021 American Chemical Society.

2.1 Introduction

Cytoskeletal fibers, along with regulatory molecules, such as α-actinin cross-linking nearby

actin filaments, and myosin walking on the actin filaments, work together to establish cell’s main

mechanochemical engine, continuously consuming energy. The plasma membrane, on the other

hand, is mostly inactive. However, it defines the boundary of a cell, playing an important role in

being a mechanical and chemical window to the outside world, regulating cytoskeletal functions

[72].

Mechanically, the membrane is a thin fluid film, resisting in-plane and out-of-plane defor-

mations, which, in turn, influences its cellular functions. For example, membrane tension has

been shown to influence exocytosis, endocytosis and cell motility [12, 13]. Membranes can in-

teract with the cytoskeleton mechanically via repulsion and adhesion [58]. Cells with the plasma

12



membrane attached to a dense cross-linked actin cortex have their boundaries significantly rigid-

ified [12]. On the contrary, disruption of filament-membrane attachments leads to cellular bleb-

bing [59, 60]. Membrane proteins also participate in numerous signaling processes, regulating,

in particular, many cytoskeletal interactions. For instance, the membrane attached WASp-VASP

complex regulates Arp2/3 activation, which is an actin branching protein essential in lamellipodia

formation [61].

Computer simulations can shed light on mechanisms of membrane-cytoskeleton interac-

tions, building on and complementing extensive experimental observations of these processes

[72]. A variety of such computational models have been proposed during the last two decades.

Both analytical and computational models were used to study the filopodia formation and prop-

erties [68, 69]. To study lamellipodial protrusions, the membrane leading edge was modeled as

a quasi 1D curve, interacting mechanochemically with a branched actin filament network [70].

Following up on these advances, it is desirable to develop a flexible, scalable, 3D model for mem-

branes and their interactions with cytoskeleton that has detailed mechanochemical components

essential for simulating whole-cell dynamics. With this goal in mind, we based our work on

MEDYAN (Mechanochemical Dynamics of Active Networks, a software designed for carrying

out detailed cytoskeletal simulations [9]), extending it such that cytoskeletal networks can be

enclosed within a movable membrane domain.

MEDYAN can simulate cytoskeletal behaviors at high spatial and temporal resolutions.

It couples cytoskeletal network mechanics with spatially resolved chemical dynamics, enabling

studies of cytoskeletal processes at the micrometer scale and the time scale of tens of minutes [9].

Simulations using MEDYAN have shed light on the turnover and treadmilling dynamics for actin

filaments [73], the stability and shape transformations of actin bundles [74] and the entropy pro-

13



duction and cytoskeletal avalanches (cytoquakes) of evolving actomyosin networks that are far

from equilibrium [4,75,76]. Currently, however, the enclosing domain for cytoskeletal evolution

in MEDYAN is rigid and immovable, imposing inflexible steric restriction on cytoskeletal fila-

ments. In this work, we incorporate into MEDYAN an explicit membrane model, which accounts

for membrane’s salient mechanical properties, volumetric effects such as osmotic pressure, as

well as filament-membrane interactions. The new membrane model paves the way for whole-cell

simulations using MEDYAN.

In this work, we applied membrane-MEDYAN to simulate cytoskeletal networks enclosed

in a vesicle, finding that actin network’s architecture plays an important role in generating mem-

brane tubular protrusions. Furthermore, our analysis indicates that bundle protrusion initiation

sensitively depends on local G-Actin concentration, bundle thickness, membrane rigidity and the

inclination angle of the bundle approaching the membrane.

The remainder of the chapter has the following organization. In Section 2.2, we introduce

the new membrane model, including its discrete representation, the mechanical energy model and

the membrane-cytoskeleton interactions. In Section 2.3, we discuss the application of membrane-

MEDYAN to several relatively simple model systems. In particular, we report on our simulations

of vesicles containing various actomyosin networks, exploring conditions that lead to formation

of successful protrusions. Finally, in Section 2.4, we summarize our work, and further elaborate

on the limitations and potential applications of the new model.

14



2.2 Methods

2.2.1 MEDYAN Simulations

The new membrane model was incorporated into MEDYAN. The latter can be used to

simulate active matter systems using a compartment-based reaction-diffusion scheme, where dif-

fusing chemical species are considered “well-mixed” in each compartment, but the locations for

force bearing species (such as actin filaments and bound myosin motors) are explicitly modeled

and tracked [9]. MEDYAN simulations are based on iteratively interleaving stochastic draws

of chemical events, using the Gillespie algorithm, followed by a mechanical equilibration step.

Importantly, chemistry influences mechanics and vice versa, creating intricate mechanochemical

feedbacks.

2.2.2 Membrane mechanics and discretization

Typically, the phospholipid bilayer has a thickness of 3− 7 nm [77]. At the scale of a

whole cell (1− 10 ţm), the thickness of the plasma membrane is much smaller than its lateral

dimensions. Therefore, we model the plasma membrane as a 2D Riemannian manifold, M ,

embedded in the 3D Euclidean space. Some of the mechanical properties of the lipid bilayer,

arising from in-plane and out-of-plane deformations and from interactions with the cytoskeleton,

are mapped onto this 2D manifold. We note here that from the perspective of a point on the

surface, we use the term “lateral directions” to denote the directions in its tangent plane.

Next, we discretize the 2D surface M into a triangulated meshwork M̂ (Fig. 2.1a). The

mesh is comprised from sets of vertices V , edges E and triangles T . Each vertex vi ∈ V is a

15



representative point on the original surface. The 3D coordinates of all vertices control the shape

of the surface. Each triangle t ∈ T is a set of 3 vertices {vi,v j,vk} ⊂ V . These piecewise linear

triangles tile the entire surface of M̂ , which, in turn, approximates the original manifold M .

Each edge e is a set of 2 vertices {vi,v j} ⊂ V , and the set of all edges E := {e = {vi,v j} : e ⊂

t for some t ∈ T}. We define the 1-ring neighbor vertices Nv,1(v) := {vi ∈ V : {v,vi} ∈ E}, and

the 1-ring neighbor triangles Nt,1(v) := {t ∈ T : v ∈ t}. Several other constraints are added to

the structure of the mesh to enforce a manifold-like structure. See the SI for more information.

Below, all membrane related potentials take a discretized form as well, using coordinates of the

vertices as corresponding arguments for the potentials and force functions elaborated below.

Excluding interactions with the cytoskeleton, we define the total mechanical energy for the

membrane as the sum of the following potentials (see Fig. 2.1b top-left, top-right, bottom-left),

Emembrane = Earea +Etension +Ebending +Evolume, (2.1)

where Earea is the area elastic energy, Etension is the surface tension energy, Ebending is the bending

energy, and Evolume is the free energy associated with the 3D volume enclosed by the membrane,

which can arise, for example, as a result of the intracellular osmotic pressure. The discretized

form of these free energies are functions of mesh vertex coordinates, hence, one can compute the

forces on vertices, with fi j = −∂E/∂x j(vi), where fi j indicates the j-th component of the force

on vertex vi, x j(vi) is the j-th coordinate of vertex vi, for vi ∈V and j ∈ {1,2,3}.

In the above outlined discretization scheme, mesh vertices do not have to represent the

material points: the time evolution of an individual vertex does not represent the trajectory of any

lipid molecule in the membrane. Therefore, one may also freely remesh the original surface at any

16



Figure 2.1: (a) The plasma membrane is represented as a 2D manifold M , and is implemented
as a triangulated mesh M̂ . (b) The potential terms in our membrane model are illustrated. The
shapes with a dashed stroke indicate allowed deformations that lead to: (top-left) surface tension;
(top-right) membrane bending rigidity; (bottom-left) volumetric potentials in case of closed sur-
face; (bottom-right) membrane-filament repulsion. The actin filaments are shown in red.

17



time as long as the resulting shape is a good approximation of the original shape. Alternatively,

we can enforce the vertices to be material points, where each vertex represents a patch of lipid

molecules on the membrane that always stay together. In this case, the local area elasticity can be

incorporated naturally in the membrane model, which will be discussed below. Our membrane

model supports both options.

In addition to closed surfaces, we can also represent membranes with fixed boundaries. The

boundaries can optionally connect the membrane surface to an implicit lipid reservoir, allowing

the total number of lipid molecules to change.

Membrane area elasticity A piece of free open membrane spontaneously relaxes to an equilib-

rium area that has zero surface tension [78]. Assuming homogeneity of lipid composition of the

membrane, in the continuum limit, the membrane is considered to have a constant equilibrium

mass density (mass per unit surface area) ρ0, so the free energy density (free energy per unit

surface area) earea can be written in a quadratic form

earea =
karea

2

(
ρ

ρ0
−1
)2

, (2.2)

where ρ is the local density on the surface, and karea is a positive constant with dimension MT−2,

where M and T are the dimensions for mass and time, respectively. Total elastic energy is then,

Earea =
∫
M

eareadS. For a membrane not connected to a lipid reservoir, the mass conservation

equation
∫
M

ρdS = m must be satisfied, where m is the constant total mass. In these integrals, dS

is the surface area element.

Assuming that vertices represent material points, there is a natural discretization of surface

18



elastic energy using Eq. 2.2. Each vertex vi carries a fixed amount of lipid molecules around it,

and, therefore, has a well-defined equilibrium area A0(vi) around it. The actual area around a

vertex A(vi) is defined as a third of the sum of areas of its 1-ring neighbor triangles, so A(vi) =

∑t∈Nt,1(vi)A(t)/3, where A(t) is the area of triangle t. Assuming that the surface area density in

the vicinity of a vertex is uniform, the area elastic energy of the discretized membrane can be

written as

Êarea = ∑
v∈V

karea

2A0(v)
(A(v)−A0(v))

2 , (2.3)

where the “hat” on E indicates that it is a discretized form.

On the other hand, if the vertices only parametrically represent the surface, not being con-

nected to material points, the area elasticity cannot be naturally discretized around each vertex.

However, if the membrane is not buffered by a lipid reservoir, we can use the mass conservation

equation, with the assumption that the surface density is uniform across the whole membrane, to

obtain the following relation,

Êarea =
karea

2A0(M̂ )

(
A(M̂ )−A0(M̂ )

)2
, (2.4)

where A(M̂ ) = ∑v∈V A(v) is the area of the discretized membrane, and A0(M̂ ) is the equilibrium

area of the membrane. One can check, using the Cauchy-Schwarz inequality, that the energy

in Eq. 2.4 is a minimum of the energy reachable in Eq. 2.3, if A0(M̂ ) = ∑v∈V A0(v). The

uniform density assumption is valid, because at the smallest time resolution that we are interested

in (∼1 ms), the inhomogeneity of density on the membrane is quickly relaxed throughout the

membrane due to its in-plane fluidity [12, 79].

19



In fact, the value of karea for biomembranes is so large, that often the membrane area in-

compressibility is assumed, and the Lagrange multipliers are used in calculations to fix the total

area [80]. However, since we employ energy minimization as the main mechanical simulation

method, we use Eq. 2.3 or Eq. 2.4 to enforce a softer constraint, by allowing the area of mem-

brane to change but with large energy penalty, instead of strictly fixing the area.

Membrane surface tension Many types of cells can actively maintain the membrane tension

within a certain range via various biological processes, such as membrane undulation and caveo-

lae formation and dissolution, which function as implicit lipid reservoirs [12]. The experimentally

measured tensions typically range from 0.01 mNm−1 to 0.04 mNm−1 [81].

If an implicit lipid reservoir is considered, the total mass of the membrane is no longer

conserved, and in this case, the vertices in our model will not be associated with material points.

Assuming a constant surface mass density and a constant chemical potential of the lipid reservoir,

the total area-dependent energy is reduced to a linear form,

Etension = γA(M ), (2.5)

being further discretized to Êtension = γA(M̂ ) where γ is the surface tension, a property of the

external reservoir.

Membrane bending elasticity The bending energy can be described using the Helfrich Hamil-

tonian [82, 83],

Ebending =
∫
M

dS
[
2κH(H −H0)

2 +κGK
]
, (2.6)

20



where κH and κG are the bending modulus and the saddle-splay modulus respectively, H is the

mean curvature, H0 is the spontaneous curvature and K is the Gaussian curvature.

According to the Gauss-Bonnet theorem, the contribution from the second term (the Gaus-

sian curvature term) is constant with fixed membrane topology and boundary conditions, so we

can ignore the second term for most purposes. We discretize the first term in Eq. 2.6 on a per-

vertex basis.

Êbending = 2κH ∑
v∈V

(Ĥ(v)−H0(v))2A(v), (2.7)

where Ĥ(v) and H0(v) are the mean curvature and the spontaneous mean curvature estimated

around vertex v, and are assumed to be uniform in the neighborhood of the vertex with area A(v).

Many different ways of estimating the surface curvature on a triangulated surface mesh

are discussed in [84, 85]. Based on our numerical experiments, we decided to allow for two

different ways of estimating the mean local curvature at a vertex in MEDYAN, which provide

good numerical stability for the systems that we studied. The first way is described in the Surface

Evolver [86, 87],

Ĥ(vi) =
〈∇iA,∇iV 〉
2〈∇iV,∇iV 〉

, (2.8)

where A is the total area of the membrane, V is the total volume enclosed by the membrane, ∇i

means “taking the gradient with respect to coordinate of vertex vi”, and 〈·, ·〉 denotes taking the

inner product of the vectors. This definition gives a continuous and signed curvature value, the

sign depending on the orientation of the surface, because the volume calculation is orientation-

dependent. Therefore, this curvature estimate is suitable for bending energy computation when

the spontaneous curvature H0 is not zero. In the continuum limit, ∇A and ∇V are both parallel

to the normal direction, but computed numerically on a discretized mesh. Hence, there can be

21



misalignments between these vectors. Under extreme conditions these two vectors can even

become almost perpendicular, leading to unphysical geometries.

When H0 = 0, the sign of the curvature no longer matters in Eq. 2.7, and we rely on a way

to directly estimate the unsigned curvature as,

Ĥ(vi) =
‖∇iA‖

2‖∇iV‖
. (2.9)

This equation provides better numerical stability under extreme conditions. For detailed discus-

sion on how the bending energy was implemented, see the SI.

The volume for a closed discretized surface mesh can be computed as V (M̂ ) = ∑t∈T V (t),

where V (t) is the signed volume of the tetrahedron formed by the triangle t and an arbitrary fixed

point p ∈ R3, and the sign depends on the orientation. It is often convenient to use the origin of

the coordinate system as the fixed point.

In-plane shearing of the membrane is ignored We do not explicitly include the in-plane shearing

stress of the membrane in our model. Because the membrane behaves as a fluid in the lateral

directions, the in-plane shearing stress cannot be sustained, while there can be a dynamic stress,

which depends on fluid’s properties. In MEDYAN, mechanical energy minimization steps imply

that the system is nearly equilibrated [9]. Therefore, we also ignore the dynamic shearing stress.

Volumetric energy A membrane vesicle is also resistant to a change of its enclosing volume.

The following quadratic expression captures the volume dependent free energy near the equilib-

rium volume,

Evolume =
kvolume

2V0(M )
(V (M )−V0(M ))2 , (2.10)

22



where V (M ) and V0(M ) are the volume and the equilibrium volume enclosed by the membrane

M respectively, and kvolume is a positive constant, with dimension ML−1T−2, where M, L and T

are the dimensions for mass, length and time, respectively. Next, we show how this formula can

be derived using the osmotic pressure difference across the membrane vesicle.

Assuming a constant osmotic pressure outside the vesicle, the free energy due to osmotic

pressure inside an enclosed volume is,

Eosmotic =−
∫ V

V ∗
∆Π(V )dV =−

∫ V

V ∗
(Πin(V )−Πout)dV, (2.11)

where Πin and Πout are the osmotic pressure inside and outside the vesicle respectively, ∆Π =

Πin −Πout, V is the volume enclosed by the membrane and V ∗ is an arbitrary reference volume.

Under an assumption that the membrane is perfectly permeable to water and impermeable

to solutes, the minimum osmotic energy occurs when V =V0 and Πin(V0) = Πout. Under further

assumptions that the total mole fraction of the solutes is far less than the mole fraction of the

solvent, and that the solutes are well mixed in the system, one can write Πin = NkBT/V , where

N is the total number of solute molecules in the vesicle. In this case, V0 = NkBT/Πout, and by

choosing V ∗ =V0, the energy can be evaluated as

Eosmotic = Πout

[
−V0 ln

V
V0

+(V −V0)

]
≈ Πout

2V0
(V −V0)

2,

(2.12)

where the approximation is the result of expanding the energy to the second order with respect to

V around V0.

23



Additional discussion of the membrane force field Our membrane model is based on Helfrich’s

scheme, which assumes that the lipid molecules are perpendicular to the membrane surface [82].

However, lipid tilting may also contribute to the total free energy, which was shown to be im-

portant at length scales of lipid molecules themselves [88]. In our current model, we are not

including the tilting effect because the resolution of the membrane mesh in our model (i.e., the

triangle element length) is typically 20− 70 nm, which is much larger than the length scale at

which tilting plays an important role.

Our current model also does not include a repulsive potential among the membrane ele-

ments. Although this may hypothetically allow unphysical geometries for folded membranes, in

our simulations reported below, we have not observed any membrane self-crossing. In part this is

explained by the steric interactions between the membrane and actomyosin networks that engen-

der a positive tension, preventing crossing of membrane’s segments. In future, a self-avoidance

potential can be naturally incorporated if self-crossing becomes a problem, however, at the ex-

pense of added computational effort.

For 2D manifold, all free energy terms mentioned above, except for the neglected in-plane

shearing energy, are parametrization invariant, which means that the energy of the surface de-

pends only on the shape of the surface. How the surface is triangulated or how the surface

coordinates are chosen only affect the discretization error. More specifically, in the absence of

membrane shearing energy, as long as the triangulations represent the same shape of the surface

up to some precision (the similarity between two surfaces can be measured using the Hausdorff

distance [89]), the total discretized free energies computed on those triangulations should be

approximately the same. This is, in fact, very useful in a simulation, because it allows us to

re-triangulate the membrane meshwork to improve the triangulation quality, without changing

24



membrane’s physical behaviors. This approach effectively addresses the problem of the quality

of triangulation becoming poor over time due to membrane’s large deformations, which, in turn,

leads to undesirable numerical instabilities.

2.2.3 Adaptive surface remeshing

The aforementioned membrane potentials (especially bending energy, that is based on a

discretized estimate for curvature) suffer numerical issues when the membrane deformation be-

comes large. To resolve this problem, we rely on the method of triangulated surface remeshing,

which essentially optimizes the membrane quality without changing its shape [90, 91].

The quality of the mesh is related to the following: the triangle shape qualities and the dihe-

dral qualities. For best numerical stability, the triangles used in the mesh should be as equilateral

as possible, while the dihedral angles between the triangles sharing an edge should be as close to

π as possible (or the angle between oriented normal vectors of the two triangles are as close to 0

as possible).

Our implementation of such an algorithm was inspired by [90, 91]. We apply local opera-

tions on the meshwork in an iterative manner to optimize both triangle and dihedral quality. With

this algorithm, we are able to get an isotropic surface triangular mesh with good quality, where

the vertex density adapts to the local curvature, making numerical simulations accurate for high

curvature regions, while being efficient for low curvature regions. See the figures in the SI for an

example of the remeshing procedure.

In MEDYAN, the adaptive mesh quality optimization for a membrane is carried out after

every energy minimization. We rely on the Marching Tetrahedra algorithm as the mesh gen-

25



erator [92], which is fast, robust and easy to implement. However, the quality of the mesh by

Marching Tetrahedra may lead to triangles that are far from being equilateral. Therefore, the

mesh optimization algorithm is essential for achieving numerical stability.

2.2.4 Interactions between the membrane and the cytoskeleton

2.2.4.1 Membrane-filament repulsion

Without tethering proteins connecting the cytoskeletal filaments to the plasma membrane,

the filaments and the membrane sterically repel each other (Fig. 2.1b bottom-right). A hypothet-

ical repulsive interaction between filament beads and membrane vertices would be insufficient

in this case, because our simulations showed that the filaments are able to occasionally pierce

through the triangle. Therefore, we developed an excluded volume interaction via repulsion be-

tween filament’s tip and a meshwork triangle, where essentially every point on the surface of

the membrane interacts with the tip. The interaction energy between a given filament tip and

membrane takes the following form,

Eexcl = kexcl

∫
M

dS
1
dn = kexcl ∑

t∈T

∫
t
dS

1
dn , (2.13)

where kexcl is the force constant to the interaction and n describes the steepness of the

potential, and d is the distance between the filament tip and a point on the membrane. This

form of the potential generalizes the potential for segmented interactions from [93] to higher

dimensional discrete elements. We found that choosing n= 4 produces a repulsive force field that

is steep enough while still amenable to analytical evaluation of forces. For detailed derivation of

this volume exclusion energy, see the SI.

26



The force acting on the filament tip from this repulsive interaction is f =−∇Eexcl, where ∇

means “taking the gradient with respect to the coordinate of the filament tip”. The reaction force

(repulsive forces on the membrane by the filament tip) is balanced by the load on the membrane

when the mechanical energy is minimized. Therefore, we directly use the value of f to be the

load force when a filament polymerizes against a membrane. The polymerization of a filament

near the membrane ratchets the membrane thermal fluctuation [94], and statistically work is done

against the membrane load, resulting in membrane protrusion or deformation [41]. In the absence

of explicit modeling of membrane fluctuation, we use the value of the load force f to compute a

scaled filament polymerization rate [9, 41],

kpoly,eff = kpolye
〈 f ,δ 〉
kBT , (2.14)

where 〈·, ·〉 indicates an inner product, and δ is a vector pointing in the direction of polymeriza-

tion. The average monomer length projected along an actin filament is, ‖δ‖ ∼2.7 nm. This cou-

pling introduces a non-linear feedback that allows instantaneously generated mechanical stresses

to locally modulate chemical dynamics. A schematic diagram of membrane-filament interaction

can be seen in Fig. 2.2.

2.2.4.2 Membrane and diffusing species

Diffusing large molecules, such as G-actin monomers, cannot move across the membrane

without assistance. Therefore, the membrane must confine the space where the protein diffu-

sion and chemical reactions take place. In a compartment based reaction-diffusion system, the

effect of membrane confinement of reaction volume is achieved by full or partial deactivation

27



Figure 2.2: An illustration of membrane as a chemical regulator in 2D space. The membrane
disallows big diffusing molecules to move across, such as G-Actin molecules or the α-actinin
linkers. Therefore, the reaction volume is defined by the region inside the membrane (approx-
imated by the straight dotted lines cutting through the compartments). The polymerization of
filaments against the membrane is, however, possible but with reduced rates due to the Brownian
ratchet effect. The actual implementation is carried out with 3D geometrical elements.

28



of reaction networks in affected compartments. As shown in Fig. 2.2, a completely deactivated

compartment does not allow chemical reactions or diffusion with neighboring compartments. A

partially deactivated compartment, however, works as a normal compartment, but the reaction

and diffusion rates are scaled by the diminished volume as well as the diminished contact area

between neighboring compartments.

To simplify computation, if a part of a membrane cuts through a certain compartment, we

fit the shape of that part of the membrane into a plane, and perform a planar cut through that

compartment to obtain the sliced volume and the sliced area on each side of the compartment.

This latter method introduces a small inaccuracy in sliced volume computation. To estimate

the corresponding error, we initialized a 0.8 ţm radius spherical vesicle at the center of a 6×6×6

compartment grid, where each compartment is a cube of size 500 nm. For each partially activated

compartment sliced by the membrane mesh, we first estimated the actual sliced volume using a

Monte-Carlo integration with 10000 samples. We found that the average relative error between

the planar approximation of the sliced volume and the more accurate Monte Carlo estimate is

below 4%. If needed, this error can be further reduced by choosing a smaller compartment size.

2.2.5 Computational efficiency

Every simulation in this work was run on a single CPU core. Detailed benchmarks on

the performance of perpendicular bundle protrusion simulations (discussed in the next section)

are reported in the SI. These benchmarks indicate that mechanical force relaxation is the most

time-consuming process in our simulations. Hence, parallelization of the mechanical energy

computation will enable accessing significantly larger time and length scales.

29



2.2.6 Code availability

The source code for this work is published with MEDYAN version 5 on http://medyan.

org, where example input files for simulating filament-vesicle interactions are provided in the

examples directory. Additionally, simulation input files and trajectories for results presented

below have been archived in Digital Repository at the University of Maryland [95].

2.3 Results and Discussion

We first applied the new force field to investigate how a spherical vesicle would deform

as its volume starts to shrink under hyperosmotic conditions. The latter causes a mismatch of

internal and external pressures, which affects our new volumetric term, leading to corresponding

reduction of system’s volume. We started simulations with an initially spherical vesicle having a

radius of 1 ţm, applying potentials given in Eqs. 2.4, 2.7 and 2.12. The equilibrium area, A0, was

fixed to be the initial spherical surface area. We continuously decreased the equilibrium volume

V0 over time, which effectively creates an increasingly hyperosmotic chemical environment, since

V0 = NkBT/Πout under dilute conditions.

Several snapshots of the simulation trajectory are shown in Fig. 2.3, indicating that un-

der a constant area constraint, the membrane is unable to maintain its original spherical shape,

crumpling as the vesicle volume shrinks. The local surface deformations, however, are smooth,

because the bending potential prevents formation of sharp features where the mean curvature

would be too large. We also notice that, over time, the regions with high convex curvature (as

shown in blue in Fig. 2.3) become more distinct and highly localized.

Next we investigated how two classes of actomyosin architectures interact with the vesi-

30

http://medyan.org
http://medyan.org


Figure 2.3: Evolution of the vesicle shape is shown under increasingly hyperosmotic conditions.
Scale bar: 1 ţm. The red-white-blue coloring scale marks the local mean curvature, with red
being the most concave and blue being the most convex. The membrane surface becomes more
crumpled over time, and the bending stress is balanced by the forces arising from volume and
area potentials. Using adaptive remeshing, the mesh is sampled more densely where the absolute
values of local curvatures are higher.

cle membrane boundary. Our simulations of vesicles containing disordered actomyosin networks,

which are not bundled, indicate that randomly initialized filaments growing toward the membrane

tend to bend instead of making a membrane protrusion, continuously inducing small membrane

ruffles (Fig. 2.4a). At a steady state, the filaments accumulate just beneath the membrane sur-

face. On the other hand, as displayed in Fig. 2.4b, the initially bundled filaments grow toward

two poles of the vesicle, forming strong tubular protrusions. Interestingly, in some cases, a few

filaments grow away from the axis of the bundle, become eventually buckled, hence, resembling

the behavior of filaments in disordered actomyosin networks. Finger-like membrane protrusions,

such as the filopodia, are ubiquitous structures in eukaryotic cells [96, 97]. Although the bio-

logical functions of the filopodia were widely explored [97–99], a clear understanding of how

filopodia form as a result of membrane-cytoskeleton interactions is still lacking. Overall, these

results demonstrate that distinct seeding patterns lead to different actomyosin network architec-

tures, where formation of bundles is crucial for generating tubular membrane protrusions.

To study what bundle properties most facilitate protrusion, we explored how the tubula-

31



Figure 2.4: Vesicle deformation by actin filaments crucially depends on the spatial organization
of actin filaments. The color codes are the same for two images: Cyan: Membrane, Red: Actin
filaments, Green: α-actinin cross linker, Blue: Myosin motor. Filaments are randomly initial-
ized in (a), and initialized as a bundle in (b). The boxed images on the bottom right show the
initial filament arrangements for each condition. Bundled filaments can generate tubular protru-
sion, while filaments in disordered actomyosin networks accumulate under the membrane, being
largely parallel to the surface.

32



tion process depends on the number of actins filaments and G-Actin concentrations. To sus-

tain stable tubular growth, the number of filaments, N, should satisfy kpoly[A]exp(−‖ f‖‖δ‖
NkBT ) ≥

kdepoly [68, 69], where ‖ f‖ = 2π
√

2κγ , which is the load force at the tip of a stable filopodium.

This imposes a minimum G-Actin concentration threshold which increases exponentially as the

number of actin filaments decreases. For example, a bundle with 7 actin filaments requires at

least 3.9× 10−4 mol/m3 of G-Actin at the tip to sustain protrusion against the membrane with

γ = 0.02pNnm−1 and κ = 100pNnm, but with 1 filament alone, the concentration has to be

at least 0.474 mol/m3. Furthermore, for a tube to form in the first place, the filaments need to

overcome the initial resistive force, which is greater in magnitude than subsequent steady state

load forces. Another salient requirement for successful protrusion formation is for the bundle to

avoid buckling [68].

To elaborate on these ideas, we simulated an actin bundle polymerization against a sheet

of planar membrane at a constant surface tension with fixed edge boundaries. The results are

summarized in Fig. 2.5. Our simulations show that a bundle with 4 filaments cannot effectively

protrude against the membrane because it is not strong enough to resist buckling (Fig. 2.5a/b).

With more actin filaments comprising the bundle, the protrusion can form at sufficiently high

G-Actin concentrations (Fig. 2.5c/d). When a finger-like protrusion is successfully formed,

its rate of growth increases with both the number of filaments in the bundle, and the G-Actin

concentration (Fig. 2.5e).

We assumed in this work that the pointed ends of the filaments are fixed, so there is no

overall retrograde flow [100] of actin filaments. We also assumed that the G-Actin concentration

is uniform along a membrane tube. Under these assumptions, the dependencies of the growth rate

on the number of filaments in the bundle and the G-Actin concentration are qualitatively in accord

33



Figure 2.5: Filament bundle deforming the membrane. Color codes for (a) and (c): Cyan: Mem-
brane, Red: Actin filaments. (a) An actin bundle with 4 filaments cannot make a successful pro-
trusion. (b) The maximum height of the membrane corresponding to (a) is shown. The growing
filaments slightly deform the membrane in the initiation phase, subsequently becoming buckled
and, hence, failing to form a tubular protrusion. (c) A bundle with 7 filaments form a finger-like
membrane protrusion. (d) The protrusion height on the membrane corresponding to (c) is shown.
Interestingly, bundles tend to buckle less when surrounded by the membrane, in agreement with
a previous work [103]. (e) The growth rate for a nascent tubular protrusion increases almost lin-
early as the G-Actin concentration increases, and also increases with the number of filaments in
the bundle.

with previous theoretical models [27, 68]. However, these assumptions allow the filopodia to

grow at a constant speed indefinitely, which is unphysical in real filopodia. The existence of actin

retrograde flow and the finite diffusion rate of G-Actin molecules lead to a steady state filopodial

structure having a finite length and a dynamically generated G-Actin concentration gradient along

the tube [69]. Furthermore, transportation of G-Actin by motor molecules, such as Myosin-X,

may lead to a non-monotonic G-Actin concentration profile along the filopodium [101, 102].

These additional aspects of growth of mature filopodia are potentially amenable to simulations

based on further elaboration of our new model.

In real cells, the polymerization directions of filament bundles are not necessarily perpen-

dicular to the membrane. At sufficiently high G-Actin concentrations, we anticipate that the

34



tubular protrusion is less likely to form when the filaments are inclined to the membrane surface,

in particular, because those filaments would be more susceptible to buckling.

To investigate how bundle’s orientation affects protrusion formation, we simulated a fila-

ment bundle polymerizing towards a planar membrane sheet, where we varied the incident angle,

α , which is the angle between the axis of the filament polymerization and the normal vector of

the membrane. We also varied the load force at the tip of a stable protrusion, by modifying both

the membrane tension and the bending rigidity simultaneously, while keeping their ratio con-

stant, i.e. γ = χγ0, κ = χκ0 where γ0 and κ0 are the reference membrane tension and bending

rigidity, and χ is a unitless parameter. In this way, the tubes formed would have the same radius

R =
√

κ0/2γ0 [68], but the load force at the tip for a stable protrusion changes as ‖ f‖ = χ2 f0.

The results are demonstrated in Fig. 2.6. When the filament bundle is perpendicular to the mem-

brane surface (α = 0), as the membrane rigidity increases (with increasing χ), the filaments first

bend, and then eventually buckle, resulting in a protrusion failure. Furthermore, our simulations

show that when filaments are not perpendicular to the membrane (α > 0), they easily buckle even

under low membrane rigidities.

The above-discussed results suggest that bundle’s orientation relative to the membrane is

important in controlling the process of protrusion initiation. These insights help to explain the

prior experimental observations showing that bundles growing within a vesicle can either be bent

by the membrane and form ring-like structures, or protrude the membrane to form finger-like

structures [104].

35



Figure 2.6: Membrane protrusion behaviors are shown under different relative membrane rigidi-
ties characterized by a unitless parameter χ (where a larger χ means larger membrane tension
and bending modulus) and the incident angle α between the filament bundle and the membrane
surface normal. All protrusions are made by a bundle comprised of 7 actin filaments. Color
codes: Cyan: Membrane, Red: Actin filament, Green: α-actinin cross-linker.

2.4 Conclusions

In this work we have addressed the problem of the interactions between the active matter

and the plasma membrane, which are important for cell functions and have been a growing field

of study. To study the rich behaviors arising from these mechanochemical interactions, we have

developed a new force field for simulating the plasma membrane and its interaction with the ac-

tive network. These algorithms were incorporated into MEYDAN, which is a broad simulation

platform for active matter. The new membrane model enables simulations of rich cytoskeletal

behaviors near the cell membrane, accounting for geometrical, mechanical and chemical details.

On a single CPU, membrane-MEDYAN simulations can treat systems on the length scale up to

several micrometers and time scales on the order of thousand seconds. Our software implementa-

36



tion allows to easily tune various model parameters, such as the elastic coefficients characterizing

a membrane, and its initial shape. Our new model enables simulations of various cellular systems

where membrane-cytoskeleton interactions play a salient role, such as lamellipodia, filopodia and

dendritic spines.

In the current work, we applied membrane-MEDYAN to simulate a membrane vesicle,

as well as the membrane-cytoskeleton systems under a variety of conditions. Simulations of

an empty membrane vesicle under hyperosmotic stress lead to membranes’ intricate crumpling.

Simulations of vesicles containing actomyosin networks show that bundling of actin filaments

is crucial for generating tubular membrane protrusions. We found that the initial growth rate

of successful protrusions is determined by the thickness of the protruding actin bundle and the

G-Actin concentrations at the growth tip. Furthermore, unsurprisingly, lowering membrane’s

surface tension or bending modulus facilitates the tubulation process. Finally we showed that

the inclination angle between the filaments and the membrane is also critically important for the

formation of finger-like protrusions. However, softer membranes can significantly mitigate this

constraint.

One limitation of the current membrane model is the absence of membrane-associated

proteins and membrane-cytoskeleton adhesions, which modulate both mechanical and chemical

properties of the membrane-cytoskeleton system. Membrane based proteins, such as the BAR-

domain proteins, can cause membrane remodeling by generating local membrane curvatures,

being present in many cellular processes such as membrane trafficking and endocytosis [21, 22].

Membrane-cytoskeleton adhesions could increase the apparent rigidity of membrane attached to

the actin cortex, which is important in cell blebbing [22, 58, 72]. They are also important in cre-

ating membrane invaginations, which are important in initiating endocytosis [105]. Moreover,

37



chemical signaling via membrane receptors are known to regulate cytoskeleton behaviors [62].

Our current model can be naturally extended in future to enable studying of such intricate pro-

cesses.

38



Chapter 3: Simulating membrane proteins and their dynamics

3.1 Introduction

The cell has the membrane as its boundary, and uses the membrane to communicate with

external environments. Cell membranes use the lipid bilayer as the “backbone”, while various

types of proteins can attach to the lipid bilayer, shaping the mechanical and chemical proper-

ties of the membranes. Membrane proteins are vital in cellular functions such as morphology

control and mechanochemical signaling. For example, cadherin is a transmembrane protein that

promotes and stabilizes cell-cell contact, while also regulating the contractility of cortical cy-

toskeleton [106]. Curvature generating proteins, such as BIN-Amphiphysin-Rvs (BAR) domain

proteins, can regulate membrane deformations, which are important, for example, in initiating

endocytosis or dendritic spines [23,105,107,108]. The spatial distribution of membrane proteins

can also regulate cell activities. As an example, B cell receptors (BCRs) form microclusters dur-

ing B cell activation, mediating intracellular signaling [109,110]. Understanding how membrane

protein dynamics affect the membrane-cytoskeleton dynamics helps unveiling their mechanisms

in cellular sensing processes.

The dynamics and spatial distribution of membrane proteins can be affected by protein in-

termolecular interactions, protein-lipid interactions, as well as interactions with the cytoskeleton.

Pairwise attractions among proteins could cluster the protein molecules, forming 2D liquid-liquid

39



phase separation (LLPS) on the membrane [16–19]. Protein-lipid interactions allow proteins

such as BAR domains and Amphipathic lipid packing sensor (ALPS) motifs to exhibit higher

affinity on specific local instantaneous curvature, by matching membrane geometry and sensing

lipid packing defects, respectively, effectively acting as membrane curvature sensor [23, 111].

Computer simulations allow us to access the details in the complex interactions among the lipid

bilayers, the membrane proteins and the cytoskeleton, thus making it possible to reveal the dy-

namics of membrane proteins during important cellular processes. To achieve this, we base off

this work of the cellular simulation platform, MEDYAN (Mechanochemical Dynamics of Active

Networks), capable of simulating the cytoskeleton and its interactions the membrane (discussed

in Chapter 2) [9, 112].

In this work, we first discuss the membrane protein dynamics with pairwise interactions and

under external potential fields. Then we introduce our new method of simulating protein reaction-

diffusion processes on the unstructured membrane surface mesh in the MEDYAN platform, while

also discussing how the potential energy field can be accounted for by adjusting the reaction rates.

The we discuss how this model can be applied to study the behavior of membrane curvature

sensing and generating proteins. Simulations using the new models show that membrane protein

diffusion can be hindered by energy barriers and confinement is stronger with higher energy

barrier and slower diffusion coefficients. Apart from confinements, the phase separation can also

occur in simulations if pairwise protein interactions are included. Furthermore, we will study how

the curvature sensing proteins preferably dock on membrane pieces with various curvature values,

and we show how curvature generation affect cytoskeleton’s ability to deform the membrane.

40



3.2 Methods

3.2.1 Thermodynamics of interacting proteins on the surface

Consider N interacting particles living in Rd . Particle i lives is at position x(i). Particles

have pairwise interactions that only depend on their distances: w(r) where r = |x(i)−x( j)|. Let X

be a vector in RdN with X = (x1
(1),x

2
(1), ...,x

d
(1), ...,x

d
(N)). The overall interaction energy is

W (X) = ∑
i< j

w(|x(i)− x( j)|). (3.1)

Let V (x) be a potential energy function that is the same for every particle. Then let U(X) =

W (X)+∑iV (x(i)) to be the total potential energy of the system.

Assume that the particles undergo overdamped dynamics, such that

Ẋ i =
1
γ

∂iU +ξ (t), (3.2)

where γ is the damping constant, ξ (t) is the thermal white noise satisfying Cξ (t) := 〈ξ (s)ξ (s+

t)〉 = 2Dδ (t), where D is the diffusion constant. The Einstein relation holds: βγD = 1. Let

p(X , t) be the probability density of the system at state X given the time t, and satisfies
∫

V N p(X , t)dV N =

1. The Fokker-Planck equation is

∂t p(X , t) = ∑
i

1
γ

∂i((∂iU)p)+D∂i∂i p. (3.3)

Eq. 3.3 is extremely hard to solve due to the high dimensionality of X . Therefore, we

41



assume that with a large enough number of particles, the particles form a density field n(x, t), such

that we can separate one particle which interacts with this density field. Naively, the interaction

energy as a result of adding this particle is V (x(i)) +
∫

V w(|x − x(i)|)n(x, t)dV . However, this

energy may be inaccurate or even blow up if w(r) rapidly varies or becomes singular when r

approaches 0, due to the inclusion of self-interaction. Instead, we can decompose w(r) into two

parts: A short-range strong repulsive interaction and a long-range weaker interaction. For the

strong repulsive interaction we can model it as an upper limit on local particle density n0. We

will use µ(r) to represent the long-range weak interaction.

At x= x(i), n(x, t) can be expanded to the second order as n(x, t)≈ n(x(i), t)+∑ j ∂ jn(x(i), t)(x j−

x j
(i))+

1
2 ∑ jk ∂ j∂kn(x(i), t)(x j − x j

(i))(x
k − xk

(i)). The function µ(|x− x(i)|) has spherical symmetry

centered at x(i), so the particle energy from interaction can be written as R(x, t) :=
∫

V µ(|x−

x(i)|)n(x, t)dV ≈ n(x(i), t)
∫

V µ(|x − x(i)|)dV + 1
2 ∑ j ∂ j∂ jn(x(i), t)

∫
V µ(|x − x(i)|)(x j − x j

(i))
2dV =

n(x(i), t)
∫

V µ(|x− x(i)|)dV + 1
2d ∆n(x(i), t)

∫
V µ(|x− x(i)|)|x− x(i)|2dV .

Under this formulation, let p(x(i), t) be the probability density for particle to be at position

x(i) at time t. The Fokker-Planck equation for this particle is

∂t p(x(i), t) = ∑
j

1
γ

∂ j
[
∂ j(R(x(i), t)+V (x(i)))p(x(i), t)

]
+D∂ j∂ j p(x(i), t), (3.4)

where

R(x, t) = n(x, t)µ0 +∆n(x, t)µ2,

µ0 =
∫

V
µ(|x|)dV,

µ2 =
1

2d

∫
V

µ(|x|)|x|2dV,

n(x, t)≤ n0,∀x, t,

(3.5)

42



The formulation above can be extended to particle systems on a differential manifold as

well. Instead of Rd , we have a d-dimensional manifold M . Locally one can use coordinates

x ∈ Rd with diffeomorphism φ : Rd → M . Assume that M can be embedded in Rm where

m ≥ d, with the inclusion map i : M → Rm. A Riemannian metric g can defined on M by

pulling back the metric δi j in Rm. Now, with two points p(1), p(2) ∈ M , we can define their

distance using either the geodesic distance or the Euclidean distance in the embedding space. We

assume that the pairwise interaction function w(r) decays fast enough such that the difference

between the two is negligible at relatively flat extrinsic curvature. Moving forward, we will use

d(p(1), p(2)) to represent the geodesic distance between p(1) and p(2). With similar arguments,

we may continue using v1 and v2 obtained in the flat space.

With the volume form dV =
√

det(g)dx1∧ ...∧dxd , we still require that
∫
M

p(p(i), t)dV =

1. Eq. 3.4 now becomes

∂t p(p(i), t) =
1
γ

div
[
grad(R(p(i), t)+V (p(i)))p(p(i), t)

]
+D∆p(p(i), t). (3.6)

where ∆ = div◦grad is the Laplace operator.

3.2.2 Surface discretization and stochastic diffusion sampling

Next, we will discretize the system for d = 2 and m = 3, that is, a 2D manifold embedded

in R3. Let {q(i)} be the set of vertices of a triangulation M̂ of M , where the edges E are

straight lines in R3 connecting neighboring vertices. Let y(i) be the coordinates of q(i) in a certain

coordinate system. Let N1,v(q(i)) be the set of all 1-ring neighbor vertices of q(i) and N1,t(q(i))

be the set of triangles that has vertex q(i). Let t(q(i),q( j),q(k)) be the triangle formed by 3 unique

43



vertices. Now, we can discretize a function g : M → R supported by basis functions f j. f j :

M̂ → R is defined on the triangulation, which evaluates to 1 at q( j) and 0 at all other vertices,

and is linearly interpolated in each triangle in R3. Therefore, g(p) ≈ ∑ j g j f j(p) where p ∈ M .

Let F be the vector space spanned by these basis functions.

Eq. (3.6) can now be evaluated in the weak form using the inner product with the test

function: 〈 f j,g〉 :=
∫
M̂

f j(p)g(p)dV ≈∑k gk
∫
M̂

f j(p) fk(p)dV . Define M̃ jk :=
∫
M̂

f j(p) fk(p)dV

as the elements of the mass matrix, which is non zero only if j = k or { j,k} ∈ E. To simplify

computation, mass lumping technique is applied to diagonalize the mass matrix [113]. After mass

lumping, the new mass matrix can have non zero terms only on the diagonal: M j j = |C j|, where

|C j| is the area around vertex q( j) and can be defined as a third of the 1-ring area ∑t∈N1,t(q( j))
A(t).

Consequently, the weak form with mass lumping becomes 〈 f j,g〉= g j|C j|.

With such finite element discretization, the first order derivatives (gradients) are well-

defined on triangles, but the higher order derivatives are zero on the triangles and undefined

on edges or vertices. To estimate higher order derivatives, we again use the weak formulation.

Consider two functions g,h : M →R, and a second order derivative of form i := div(ggradh). If

g is constant on M , then i = g∆h. Using the divergence theorem, we have

〈 f j, i〉=
∫
M̂

f j(p)div(ggradh)dV,

=
∫
M̂

[
div( f jggradh)−g〈grad f j,gradh〉

]
dV

=
∫

∂M̂

f jg〈gradh, n̂〉dl −
∫
M̂

g〈grad f j,gradh〉dV,

(3.7)

where n̂ is the outward normal at the boundary of M̂ . With vanishing boundary conditions, the

first term is zero. Therefore, 〈 f j, i〉=−
∫
M̂

g〈grad f j,gradh〉dV ≈−∑k,l gkhl
∫
M̂

fk〈grad f j,grad fl〉dV .

44



(a)

q( j) q(n)

αn

βn

(b)

q( j) q(k)

q(n1)

q(n2)

α

γ

β

δ

Figure 3.1: Geometric elements used in calculating a jk. (a) j = k. (b) j 6= k.

Define a jk :=−∑l hl
∫
M̂

fk〈grad f j,grad fl〉dV , and then 〈 f j, i〉 ≈ −∑k a jkgk. When j = k,

a j j =
1
6 ∑

q(n)∈N1,v(q( j))

(cotαn + cotβn)(h j −hn). (3.8)

When j 6= k, a jk is non-zero only if {q( j),q(k)} ∈ E. In this case,

a jk =
1
6
[
(h j −hk)(cotα + cotβ )+(h j −hn1)cotδ +(h j −hn2)cotγ

]
. (3.9)

The angles used in the equations above are illustrated in Fig. 3.1. With mass lumping, i j =

−∑k a jkgk/|C j|. As a special case where g j = 1,∀ j, one has i = ∆h, and obtains the famous

cotangent formula for the discrete Laplacian on surface mesh: (∆h) j =
1

2|C j| ∑q(n)∈N1,v(q( j))
(cotα+

cotβ )(h j −hn). Let S be a matrix such that ∑k S jkhk = |C j|(∆h) j,∀ j,{h j}, and S is also known

as the “stiffness matrix” [113].

With the finite element discretization, Eq. 3.6 becomes

|C j|∂t p j(t) = ∑
k

1
|Ck|

(
−

a jk

γ
+DS jk

)
|Ck|pk, (3.10)

45



where a jk is defined with respect to functions g = p and h = R+V .

For a small time period, the above equation can be applied to all particles in the system.

From Eq. 3.10, the transition rate q jk for a particle to move from C j to Ck can be found as

q jk =
1

|C j|

(
−

ak j

γ
+DSk j

)
. (3.11)

If the potential energy profile R+V is smooth, for a fine enough mesh, a jk scales with the

mesh size, while S jk is invariant to the mesh size. Therefore, the physical transition rate q jk being

positive puts a constraint on the maximum size of the mesh. In practice, we use a small enough

mesh to ensure that the transition rate is positive most of the time. Considering rare cases where

negative transition rate occur due to extraordinary local mesh shapes or a large energy gradient,

we manually set the transition rate to zero. On the other hand, the transition rate is also set to

zero if the transition causes the target cell to violate the nk ≤ n0 constraint.

3.2.3 Membrane-related reactions

Apart from membrane species diffusions, we also incorporate in our model the membrane-

related reactions (MRRs), which are the reactions happening at the membrane surface, but in-

cluding at least one 3D diffusing species. Specifically, we want to model the reactions of form

Mr1 +Mr2 + ...+Cr1 +Cr2 + ...
M−→ Mp1 +Mp2 + ...+Cp1 +Cp2 + ... (3.12)

where Mri, Mpi are 2D diffusing membrane species as reactants and products, respectively, and

Cri, Cpi are 3D diffusing species as reactants and products, respectively. Assuming that the law

46



of mass action applies to the membrane domain, and given a small patch of membrane with area

A, the membrane species and nearby 3D diffusing species are uniformly distributed, then the

propensity of such a reaction is related to the rate constant r as

p = rA∏
i
[Mri]∏

i
[Cri], (3.13)

where the operator [·] returns the area density for membrane species and the volume density for

3D diffusing species.

Some common cellular events can be described using MRRs:

Adsorption / membrane binding C M−→ M

Desorption / membrane unbinding M M−→C

Membrane species activator M+C1
M−→ M+C2

Membrane activator (implicit) C1
M−→C2

Membrane species dimerization M1 +M1
M−→ M2

Similar to membrane diffusion processes, the potential energy field on the membrane also

affects reaction rates. For example, let the reaction rates of a reversible adsorption/desorption

reaction be rad and rde, respectively, and their ratio encodes the free energy change between the

bound and the unbound states. If the protein is subject to extra potential fields on the membrane,

we will adjust the desorption rate to make up for that:

rde,eff = rdeeβV . (3.14)

47



3.2.4 Free energy of the particle system

Following the derivation in Ref. 114, the total free energy can be found by

E =
∫
M

(
− 1

β
φ(n(p))− 1

4

∫
M

µ(d(p, p′))
[
n(p)−n(p′)

]2 dS′+V (p)n
)

dS, (3.15)

where φ(n) = n ln
[n0−n

n

]
+ n− 1

2β µ0n2. Assume that in the neighborhood of p, one can use a

geodesic frame x̃ to represent surface positions, such that it is isometric to δi j in R3 at p, and

close to being isometric near p. Then one can expand the second term above, such that

E ≈
∫
M

− 1
β

φ(n)− 1
4

∫
M

µ(|x̃− x̃′|)

[
∑

i
∂x̃in(x̃)(x̃− x̃′)i

]2

dS′+V n

dS

=
∫
M

(
− 1

β
φ(n)− 1

4
〈gradn,gradn〉

∫
M

µ(|x̃− x̃′|)
[
(x̃− x̃′)‖

]2
dS′+V n

)
dS

=
∫
M

(
− 1

β
φ(n)− µ2

2
〈gradn,gradn〉+V n

)
dS,

(3.16)

where (x̃− x̃′)‖ is the unit vector along the gradn direction.

In MEDYAN, we select the enthalpic part of the energy above, which enters energy mini-

mization process. The entropic energy minimization is implicitly done via the stochastic diffusion-

reaction sampling mentioned above, so that the dynamics of protein distribution equilibration is

retained.

48



3.2.5 Membrane curvature sensing and generation model

A protein binding to the membrane locally with curvature H results in free energy change

∆E = Eb +VH(n,H), (3.17)

where Eb < 0 is the fixed portion of binding energy, and VH(n,H) is the energy dependent on

local protein concentrations and the curvature. Various computational models were developed

to provide an estimate to VH(n,H) [23]. Two simplest examples are the spontaneous curvature

model and the curvature mismatch model [115, 116]. In the spontaneous curvature model, the

bending energy density of the protein-bound membrane surface is still in the Helfrich form:

ebend = 2κ(H −H0)
2, where the bending rigidity κ is a constant, but the spontaneous curvature

H0 depends on the membrane protein’s mechanical characteristics and its local concentration. On

the other hand, the curvature mismatch model treats the lipid bilayer and the proteins separately,

with ebend = 2κl(H −H0)
2 + 2φpκp(H −Hp)

2, where κl and κp are the bending rigidities of the

lipid bilayer and the protein, respectively, Hp is the spontaneous curvature of the membrane

bound protein, and φp is the fraction of area occupied by the protein. If the area of individual

membrane protein is Ap, then φp = nAp. Hybrid models where both the bending rigidity and the

spontaneous curvature depend on protein concentration are used as well [29]. In this work, we

adopt the following bending energy density expression [117]

ebend = 2κ(H −H0φp)
2, (3.18)

49



where

κ =

(
1−φp

κl
+

φp

κl +κp

)−1

H0 = Hp
κp

κl +κp

(3.19)

When φp = 0, Eq. 3.18 reduces to the simple Helfrich form for lipid bilayer with zero equilibrium

curvature. For this reason, the constant used in Eq. 3.18 is scaled from 1/2 to 2 to match the

Helfrich form conventions in this work. The appendix provides more discussions on this energy

formula.

In the discretized mesh model, VH(n,H) at each vertex can be obtained from the difference

in bending energy before and after adding an additional molecule to the cell around this vertex,

ie VH(n,H) = |C |[ebend(n+ |C |−1)− ebend(n)]. If a cell already contains many such molecules

(where n|C | � 1), one can approximate VH(n,H) using VH(n,H)≈ ∂nebend(n).

3.3 Results

3.3.1 Energy barriers acting as protein confinements decrease diffusivity

With our new membrane-protein model, we simulate membrane protein diffusions with

diffusion barriers, where we also benchmark the appropriate mesh size under specific potential

energy gradients. In reality, the surface diffusion barriers can arise due to dense clustered mem-

brane proteins [67] or nearby actomyosin cortical networks [57].

We initialized a piece of membrane of 1500nm× 1500nm parallel to the xy plane. The

membrane species M can diffuse on the membrane under a potential energy field. The potential

energy field is a constant along the y direction and is a Gaussian function in the x direction with

center at x = 750nm and width 400 nm, which represents a linear diffusion barrier parallel to the

50



y axis. The maximum value of the potential energy is varied from 0 to 4 kBT (Fig. 3.2ab). The

membrane domain is triangulated, and the characteristic mesh size is controlled during both mesh

generation and adaptation processes, with a value in 20 – 300 nm. We added 22500 copies of M

uniformly at region x < 500nm at initialization (t = 0). Then, we start the simulation and let the

species diffuse in the surface. We track the number of molecules on the right half of the plane

(x > 750nm) at every frame, which starts at 0 and goes toward the equilibrium value over time.

At equilibrium, due to symmetry, the average molecule count on the right should be a half of the

initialized copy number. Unless otherwise specified, all parameters used are listed in Table. 3.1.

Fig. 3.2c show trajectories of molecule counts on the right under various barrier heights.

Notice that when the height is 0, the potential energy is zero everywhere on the surface, meaning

that the proteins undergo free diffusion. The equilibration rate is slower with higher barriers,

resembling a slower state transition rate under higher activation energy. With a constant barrier

height, the equilibration rate increases with increasing diffusing coefficient of M (Fig. 3.2d).

In Fig. 3.2e, we also show that the equilibrium can be reached much faster if the protein

can unbind from the membrane and diffuse in 3D in addition to in-membrane diffusion, circum-

venting the energy barrier in this pathway.

In a triangulated meshwork, the transition rate in Eq. 3.11 should be non-negative, which

sets an upper limit to the mesh size given potential energy gradients. As discussed in Appendix

B.1 and Fig. B.2, the upper limit of permissible mesh sizes decreases with the height of the

energy function.

51



Table 3.1: Membrane diffusion barrier simulation parameters.

Description Value
Default 2D diffusing coefficient (D2D) 5×106 nm2/s a

Default 3D diffusing coefficient (D3D) 1×108 nm2/s b

Protein binding rate (rad) 6000 nm/s c

Protein unbinding rate (rde) 10 s−1 d

aFrom Ref. 118, typical diffusion rates for membrane proteins with radius under 5 nm are between 3×106 nm2/s
and 6×106 nm2/s.

bEstimated using Stokes-Einstein equation D3D = kBT/6πηr, where η is the dynamic viscosity of water (∼
1×10−6 pNs/nm2) and r is the protein radius (∼ 2.5nm) [119].

cProtein binding and unbinding rates are chosen such that in equilibrium, ∼ 22500 molecules are bound on the
membrane with ∼ 40000 molecules in total, and that this equilibrium can be reached before 0.4 s.

dSee footnote c.

Figure 3.2: Membrane species diffusion with a barrier. (a) On a piece of square membrane, the
energy barrier is initialized as a Gaussian function of various heights with 400 nm width along
the x direction and is uniform along the y direction. The triangular mesh discretization is shown
in gray. (b) The energy function described in (a) in the x direction. The proteins initialized on the
left must overcome the energy barrier to reach the right. (c) Fraction of equilibrium protein count
on the right as a function of time, varying barrier heights. Equilibration is slower with higher
barrier energy. Diffusion coefficient: 5×106 nm2/s. (d) Fraction of equilibrium protein count on
the right as a function of time, varying diffusion coefficients. Equilibration is faster with higher
diffusion coefficients. Barrier height: kBT . (e) If the protein can unbind from the membrane
and has a higher 3D diffusion rate, it can overcome the barrier easily, and greatly accelerates the
equilibration process.

52



3.3.2 Pairwise protein interactions promote LLPS

To incorporate protein pairwise interactions in our simulations, we need a pairwise at-

tractive potential µ(r) that µ0 and µ2 defined in Eq. 3.5 are finite. If one normalizes the

pairwise interaction function resulting in a 2D distribution f (x) := µ(|x|)/µ0, its standard de-

viation σ characterizes the range of the interaction. σ can be represented by µ0 and µ2 as

σ2 =
∫

V f (x)|x|2dV = (1/µ0)
∫

V µ(|x|)|x|2dV = 4µ2/µ0. For simplicity, we model µ(r) as a

Gaussian function with µ(r) = −hexp(−r2/σ2), where h ≥ 0 is the height of this energy func-

tion. In this case, µ0 =−πhσ2 and µ2 =−πhσ4/4. We note that if the attractive potential takes

the form of a power function µ(r) =−hrp where p <−1, µ0 necessarily diverges when h > 0.

We initialized a square membrane of size 400nm×400nm, and uniformly initialized 1500

molecules bound to the membrane. Then, we ran chemistry simulations, allowing those proteins

to diffuse. The simulation parameters are displayed in Table 3.2. The simulation lasts 0.5 s

and snapshots are captured every 0.01 s starting at 0.11 s. When there is no protein attractive

interaction, the protein shows gaseous behavior. As the protein interactions become stronger,

protein condensates that span multiple membrane cell volumes are formed, showing a liquid

phase separation behavior. Fig. 3.3b and c shows two snapshots when the attractive interaction

is weak (h = 0.01kBT ) and strong (h = 0.21kBT ), respectively. The equilibrium area fraction

distribution is collected at a variety of energy heights and is displayed in Fig. 3.3a. The phase

separation starts at h ∼ 0.09kBT , where the bifurcation showing phase separation starts to appear

beyond that. In Appendix B.2, we show that the initial protein distribution has no observable

effect on the final phase separation results.

In LLPS-manifesting simulations, the range of mesh sizes faces stronger restrictions. It

53



Table 3.2: Protein LLPS simulation parameters.

Description Value
Diffusing coefficient 106 nm2/s a

Protein area (Ap) 19.6 nm2 b

Interaction range (σ ) 18 nm c

Target mesh size 12 nm

aThe same order of magnitude as described in footnote a in Table 3.1.
bEstimated using protein with a circular membrane projection of radius 2.5 nm.
cWe choose pairwise interaction ranges with the order of magnitude ∼ 10nm. This is backed by the experimental

results in Ref. 120, where the minimum pairwise interaction energy between ATP synthase c-rings happens at a
distance of 10.3 nm and drops below 1kBT for distances larger than 12.5 nm.

should not be greater than the protein interaction range σ , otherwise high potential energy gra-

dient (mentioned in the previous section) will emerge on the cluster boundaries, leading to er-

roneous clustering behaviors. An example is presented in Appendix B.1 and Fig. B.3. On the

other hand, the mesh size should still be greater than the protein size. Consider a hypothetical

membrane mesh where every cell has size 1.5Ap. Since each cell contains an integral number

of proteins, a cell can never contain 2 such proteins simultaneously and the remaining 0.5Ap can

never be utilized by the proteins, which effectively increases the protein area by 50%. When the

size is larger, such discrepancy (1−b|C |/ApcAp) can be reduced. Due to the range restrictions

of mesh sizes, this model is more useful for small proteins with long-range attractive pairwise

interactions.

3.3.3 Concentration sorting by local curvature

To simulate protein curvature sensing behavior, we use Eq. 3.18 to determine the surface

potential energy. We initialized membrane vesicles with various radii while keeping their shapes

fixed after initialization. Proteins were added to the system as 3D diffusing species, that can bind

to and unbind from the membrane surface (Fig. 3.4a). The effective unbinding rate depends on

54



Figure 3.3: Membrane protein phase separations. (a) A snapshot for the protein distributing on
the membrane in a gas phase. h = 0.01kBT . (b) A snapshot of the protein with liquid conden-
sations on the membrane. h = 0.21kBT . (c) The probability distribution of area fraction (φp) for
different interaction strengths (h). When h = 0.21kBT , the distribution becomes bimodal, dis-
playing phase separation behavior. (d) The area concentration profile among membrane cells, at
various protein attractive interaction strengths quantified by its height. The brighter color indi-
cates a higher probability in the distribution. When the pairwise interaction is weak, the diffusion
takes dominance and the system is in a gas phase. When the height of the energy function goes
above ∼ 0.09kBT in