Active filaments on curved surfaces: from single filaments to dilute suspensions

Abstract

Curvature plays a central organizational role in active polymer dynamics. Using large-scale Langevin-dynamics simulations, we study active semiflexible filaments confined to smooth curved surfaces and map how curvature, bending rigidity, and activity interact. We find geodesic alignment, curvature lensing, and curvature-induced trapping. In particular, regions of negative Gaussian curvature localize filaments and hinder global surface exploration. These results show how surface geometry can be used to control the organization and transport of active matter on curved substrates.

---

1 Introduction

Curvature is integral to the architecture and function of living systems, paradigmatic examples of active matter, across scales. During morphogenesis, epithelial sheets fold into tubes and furrows, generating structures such as the vascular system, brain and spinal cord.^1–4 Surface geometry contributes to long-range coordination in the spiral drift of cells across the hemispherical cornea.^5,6 At the subcellular level, organelles such as mitochondria remodel their membranes during apoptosis.^7–10 Understanding the role of curvature on active dynamics is, therefore, an important challenge.^11–25

To investigate these effects in controlled conditions, synthetic systems have been developed using microtubules and ATP-powered kinesin motors confined to curved oil–water interfaces. In both spherical²⁶ and toroidal²⁷ geometries, these systems exhibit rich curvature–activity coupling, such as oscillatory defect dynamics²⁶ and defect unbinding transitions.²⁷ Notably, such phenomena can be captured by continuum models^19,28,29 that do not explicitly resolve the filamentous structure of the constituents.

However, little is known about how the interplay of activity and curvature affects long, semi-flexible chains, when the extended nature of the agents cannot be neglected and leads to new motion patterns. Only recently, experiments have demonstrated that confining cytoskeletal filaments to the inner surface of spherical lipid vesicles can induce robust collective motion patterns, such as streams, polar vortices, equatorial bands, and globally arrested states.³⁰ While this work established a link between curvature and pattern selection, it was limited to spherical surfaces (i.e., surfaces of constant curvature), leaving open how spatially heterogeneous or sign-changing curvature affects the dynamics, and how filament flexibility and activity interact with the geometry.

Since experimental control of the vesicle shape and curvature distribution is not simple, computational models can be used to map out key properties of active filaments confined to move on curved surfaces. On flat substrates, numerical simulations of semiflexible filaments revealed interesting effects such as coiling and clustering, as well as nontrivial unwinding transitions.^31–37 Confinement to curved surfaces introduces additional elastic penalties, since passive filaments are generically misaligned with local geodesics.^38–40 All these combine to form intricate physics that remains poorly understood.

To address how spatially varying curvature influences the dynamics of active semiflexible filaments, we perform large-scale Langevin dynamics simulations on closed curved surfaces with fixed geometry. Our aim is to go beyond spherical confinement and investigate how sign-changing curvature interacts with filament flexibility and activity to produce new dynamical regimes. By systematically varying surface shape from spheres to Gaussian bumps and necked Cassini ovals and exploring densities from isolated filaments to dilute suspensions, we disentangle the effects of curvature, stiffness, and collisions. We find that curvature can significantly guide filament dynamics, with activity and bending rigidity competing to determine alignment with geodesics. At the collective level, we observe curvature-induced confinement and the emergence of self-organized rotating bands. Notably, on Cassini surfaces, regions of negative curvature promote filament trapping and hinder global exploration. These results highlight surface geometry as a powerful mechanism to direct active filament organization.

The paper is organized as follows. Section 2 introduces the theoretical framework for the bending energy of passive filaments on curved surfaces. Section 3 describes the simulation model for active semiflexible chains. Section 4 presents our main findings, focusing on both single-filament behavior and collective dynamics on curved surfaces. Finally, Section 5 summarizes our conclusions and outlines future directions.

2 Overview of the behavior of filaments on curved surfaces

Point-like active particles confined to a smooth, curved two-dimensional surface embedded in Euclidean space exhibit nontrivial dynamics due to the effects of the surface curvature.^18,19,23,25 For example, the curvature induces a geometric torque on the particles, directing their motion along the geodesics.⁴¹ Geodesics are the curved space equivalent of a straight line, i.e. they correspond to the path which a body constrained to a surface would follow under force-free conditions (beyond forces enforcing the constraint).

In many experimental setups, however, active agents are better represented as semi-flexible or stiff filaments rather than point particles. The influence of activity on filament behavior remains only partially understood, even in the flat case.^{33,35,36,42–45} For example, in contrast to passive filaments, where generic features of the system do not depend on the specifics of the discrete bond potentials, in the presence of activity, one needs to carefully distinguish between different types of discrete description of filaments.³⁷

For a filament embedded in a curved surface, in the absence of activity, its configuration is governed by intrinsic and extrinsic curvature; describing these effects presents a significant theoretical challenge.^{38–40,46–48} Theoretical approaches typically consider an inextensible filament modeled as a smooth curve γ of length L embedded in the smooth, orientable surface [scr S, script letter S]. The curve is parametrized by arclength s, while the surface is parametrized by curvilinear coordinates (u, v). Hence, γ(s) = r(u(s), v(s)), where r is the position vector in Euclidean space. At each point of the curve, one can assign the tangent vector (|t| = 1 since s is the arclength parametrization), and define the surface unit normal n. In general, n is not aligned with the normal to the curve, . However, since t lies in the tangent plane of the surface, n ⊥ t everywhere. One then defines the unit tangent normal l = t × n. The set of vectors {t, n, l} forms the Darboux frame,⁴⁹ an embedded curve analogue of the Frenet-Serret frame for curves in .

The scalar curvature of the curve is , i.e., it measures how the tangent vector changes its direction as one moves along the curve. The bending energy of the filament is given by ref. 50

where α is the bending rigidity. This expression is typically accompanied by the fixed-length requirement, e.g., by setting |t| = 1 everywhere, which is imposed using Lagrange multipliers. For an embedded curve, the curvature can be decomposed into the geodesic curvature, k_g, and the normal curvature, k_n, such that k² = k_g² + k_n². The normal curvature measures the curvature of the surface in the direction of t, while the geodesic curvature measures the intrinsic curvature of the curve, i.e., how much the curve bends along l.

Using the Darboux frame, the Euler–Lagrange equation that determines the equilibrium configurations of the filament is ref. 46,

where τ_g denotes the geodesic torsion, i.e., how the surface normal rotates around the tangent to the curve as one moves along the curve. Primes denote derivatives with respect to s, and the boundary conditions set the constant c.^38,40,46,51 Assuming free ends,⁴⁶ the boundary conditions at the filament endpoints (s = 0 and s = L) are

Eqn (2) is very hard to solve even on simple surfaces. It, however, highlights the coupling between the geometry of the surface and the elasticity of the filament, making even the passive case far from trivial. A key result from ref. 46 shows that geodesics (i.e., curves for which k_g = 0) are solutions to the Euler–Lagrange equation only when the condition (k_n²τ_g)′ = 0 is satisfied. Furthermore, a related theorem^38,46 states that on the plane or the sphere, an open-ended filament adopts an energy-minimizing configuration if and only if it follows a geodesic. This result follows from the fact that on both the plane (where for all curves k_n = 0) and the sphere (where k_n is constant), the geodesic torsion τ_g vanishes identically, so that (k_n²τ_g)′ = 0. However, under confinement, a closed loop on a sphere will deform away from the circular (geodesic) path, adopting a noncircular shape with increased bending energy.³⁹

In general, the condition (k_n²τ_g)′ = 0 is not satisfied, and the energy-minimizing configuration of the filament typically deviates from a geodesic.^{38,40,46,47,51} Since our goal is to study active filaments on surfaces with nonuniform curvature, where bending energy, curvature, and self-propulsion compete, the analytical treatment is not feasible and we instead resort to numerical simulations. We briefly note that recent work has shown that the mapping between the continuous and discrete behavior of active polymers is not necessarily straightforward.³⁷

3 Simulation model

We consider a collection of active filaments confined to a closed surface, with each filament composed of N_b beads. We ignore long-range hydrodynamic interactions, i.e., assume the dry limit.⁵² The motion of bead i is governed by the Langevin equation:where m_i is the bead mass, r_i is the position in Euclidean space of bead i at time t, and γ is the friction coefficient with the embedding surface. The active propulsion force f^act_i = f_pt_i propels bead i with magnitude f_p in the direction t_i set by the normalized sum of the unit-length vectors r̂_i+1,i = (r_i+1 − r_i)/|r_i+1 − r_i| and r̂_i,i−1 = (r_i − r_i−1)/|r_i − r_i−1|.^32,53 The interaction between beads i and j is given as f_ij = −∇_riϕ(r_ij), where r_ij = |r_i − r_j| and ϕ is the pair potential that describes both bonded and non-bonded interactions. R_i(t) represents thermal fluctuations and follows a Gaussian distribution with zero mean and variance 〈R_i(t)·R_j(t′)〉 = 4γk_BTδ_ijδ(t − t′) where k_B is the Boltzmann constant, T is the temperature, and 〈⋯〉 denotes the thermal average.

Bonded interactions, ϕ_B, account for chain stretching and bending within the filament, with stretching described by a tether bond potential⁵⁴ and bending by a harmonic angle potential.⁵³ Nonbonded interactions, ϕ_NB, capture steric repulsion and are modeled using the Weeks–Chandler–Anderson (WCA) potential,⁵⁵ for r_ij < 2^1/6σ, where ε sets the interaction energy scale and σ is the monomer diameter. Bonded interactions are applied only between beads i and j belonging to the same filament, whereas nonbonded WCA interactions act between any pair of beads that are within the interaction cutoff 2^1/6σ regardless of which filament they belong to.

Surface confinement is enforced by projecting bead positions onto the surface and restricting velocities and forces to the local tangent plane. For a general smooth surface [scr S, script letter S] parametrized as r(u, v), a point can be projected on [scr S, script letter S] by minimizing , e.g., by using Newton minimization with gradient . Projecting vector a onto the tangent plane is simply, where n_i is the unit-length surface normal at r_i.

In this study, we simulate polymers with a degree of polymerization N_b = 80, and work in the set of units where the monomer size σ = 1.0, interaction energy scale ε = 1.0, thermal energy k_BT/ε = 0.1, and monomer mass m = 1.0. Although our simulations are conducted using Langevin dynamics, we do not expect the qualitative behaviors discussed below to differ significantly under Brownian dynamics. On flat surfaces, inertial effects have primarily been shown to influence filament compactness.^44,45 Therefore, in the overdamped limit, we expect the same qualitative behaviour as that observed using Langevin dynamics with γ = 1. Each polymer has a length L ≈ b_bond(N_b − 1)σ, where b_bond = 0.86. For simulations of filament suspensions, the systems were first equilibrated for 10⁵ time units using a time step dt = 10⁻⁴ before data collection. All simulations were performed using the GPU-accelerated SAMoS molecular dynamics package⁵⁶ with the BAOAB Langevin integrator.⁵⁷ Data were analyzed using custom Python scripts and visualized with ParaView.⁵⁸

We introduce a geometry-dependent dimensionless Péclet number defined as where l_s represents the characteristic length scale along the surface over which the effects of curvature become significant to a filament's motion. Specifically, we define where |K|_max is the maximum of the absolute value of the Gaussian curvature. κ is the discrete equivalent of the continuum bending rigidity α. This Péclet number quantifies the impact of surface curvature on filament dynamics. This approach is compatible with findings on spherical surfaces, where curvature introduces an additional timescale, , with R being the sphere's radius and v₀ the self-propulsion speed of the active particles. This timescale is in addition to the conventional rotational diffusion time and affects the dynamics of active Brownian particles on curved surfaces.²⁵

4 Results and discussion

4.1 Single active filaments on a Gaussian bump

We begin by studying the deterministic (i.e., noiseless) active dynamics of a single active filament confined to a Gaussian bump. The Gaussian bump has non-uniform curvature, with positive curvature at the apex and a narrow band of negative curvature at its base (Appendix B). The curvature also quickly decays to zero as one moves away radially, making it a good model to study the effects of non-uniform curvature in a flat background. In addition, a generic curve has both nonconstant normal curvature and nonzero geodesic torsion, so the filament configurations that minimize bending energy are typically not geodesics.

The surface is axisymmetric about the vertical axis through its peak and can naturally be parameterized by (r, θ), where r measures the radial distance from the axis of symmetry and θ is the azimuthal angle. The position vector r(r, θ) = (r cos θ, r sin θ, h(r)), describes points on the surface, with a height function where A is the height at the peak, and σ_g sets the width of the bump. The height function h(r) describes the elevation above the xy-plane in Euclidean space . This sets l_s = σ_g²/A.

We initialize filaments sufficiently far from the bump that the Gaussian curvature can be neglected by placing them along geodesics (i.e., straight lines) along x at different y-positions and track their trajectories under varying activity. If geodesics minimize the bending energy, the filaments remain on them, with geodesic curvature k_g remaining zero as the filament evolves. This can be understood as follows. Active force is tangent to the filament and the surface at every point, and as such, forces beads to follow geodesics.⁴¹ Bending energy is minimized by remaining on a geodesic, and any deviation from it results in a restoring force that pushes each bead back to geodesics. Therefore, the entire filament follows a geodesic. For fully floppy filament, i.e., κ = 0, there is no bending penalty and the filament trajectory is fully determined by active force, i.e., it follows the initial geodesic it was placed on ref. 41 with k_g ≈ 0 (Fig. 1a). Note that due to the spatial discretization Δs = 1 (Appendix C), the computed k_g fluctuates between 10⁻⁸ and 10⁻⁴, but decreasing Δs reduces these fluctuations.

Fig. 1 Deviation of active filaments from geodesics at activity f_p = 10¹ as a function of bending rigidity κ on a Gaussian bump with amplitude A = 2 and width σ_g = 1. Each panel shows trajectories of the head monomer of eight filaments initially placed along geodesics at different y-coordinates in the x–y plane. Trajectories are color-coded by the instantaneous geodesic curvature k_g averaged over filament length, and the grayscale background indicates the Gaussian bump height profile h(r). The geodesic curvature is computed numerically, and details of the numerical procedure are provided in Appendix A. (a) For zero bending rigidity (κ = 0), filaments closely follow their initial geodesics, with k_g ranging from 10⁻⁸ to 10⁻⁴. (b) At κ = 5 × 10¹ (Pe = 10⁻¹), only the filament initialized at y = 0 remains on its geodesic, which has zero geodesic torsion (τ_g = 0) and minimizes bending energy. The others deviate from their initial geodesics and do not transition to alternative geodesics, as shown by the marked increase in k_g when crossing the bump. (c) At high bending rigidity (κ = 10², Pe = 5 × 10⁻²), all filaments except the one at y = 0 exhibit substantial deviation from geodesics. The higher k_g indicates that bending energy dominates filament dynamics, driving them away from geodesic paths.

As κ increases, filaments deviate from their initial geodesics and k_g increases (Fig. 1b). Finally, for very stiff filaments, these deviations and the associated k_g become more pronounced (Fig. 1c). This leads to the wavelike trajectories observed in Fig. 1c, where the filament pays a small penalty by increasing its geodesic curvature k_g in order to reduce its normal curvature k_n, causing deviations away from the geodesic path. Only the filament initialized at y = 0 continues to follow its geodesic because this geodesic has τ_g = 0, corresponding to a bending energy minimum (Section 2). In this regime, bending energy dominates, and filaments only follow geodesics that also minimize this energy, as with the central geodesic. These observations highlight that tangential active forces drive filaments along geodesics. When activity is strong enough to overcome bending rigidity, filament dynamics are dictated by surface geometry, leading to alignment with geodesic paths.

We can understand this competition between activity and bending rigidity through a simple scaling argument. We move to the overdamped limit of the filament dynamics and compare the time taken for the filament to pass over the bump, driven by its active self-propulsion t_a, with the time taken for the filament to deform due to its elastic bending forces t_b. When t_a ≪ t_b, we expect the filament to have passed smoothly over the bump on a geodesic trajectory before the elastic force has time to alter its dynamics. Whereas when t_b ≪ t_a, we expect the filament to be significantly pushed away from the geodesic trajectory by the bending forces, faster than it can move past the bump, leading to the trajectories seen in Fig. 1b and c. In Appendix D, we estimate the ratio . This means at large Pe we expect t_a ≪ t_b, so the filament is expected to remain on a geodesic, consistent with our numerical results.

4.2 Motion on surfaces with spherical topology

Motivated by experiments,³⁰ we focus on surfaces with spherical topology. A key distinction compared to the Gaussian bump is that the sphere is compact, and geodesics are closed loops, i.e., great circles. The geodesic torsion τ_g = 0, and an open active filament will follow a great circle (Fig. 2a). It, however, remains unclear how filament dynamics are affected if the spherical topology is retained but the surface is deformed to have non-uniform curvature. To address this question, we gradually deform the sphere to introduce spatial variations in curvature while preserving its closed topology. Specifically, we study the behavior of single active filaments on Cassini ovals revolved around the x-axis. These surfaces are implicitly defined by

f(x, y, z) = (x² + y² + z²)² − 2a²(x² − y² − z²) + a⁴ − b⁴ = 0,

where the parameters a and b control the shape. In the limit a = 0, this equation reduces to a sphere of radius b. For , the surface resembles a prolate spheroid, while for , it transitions to a peanut-shaped geometry with two lobes connected by a neck of negative curvature (Appendix B). For a sphere (a = 0), the natural length scale is simply the radius (i.e., l_s = b). For other Cassini ovals (a > 0), the maximum value of the Gaussian curvature, |K|_max, can be computed numerically. As shown in Fig. 2e, for peanut-shaped surfaces this maximum typically occurs in the neck region. The corresponding curvature values in this area are indicated in the color bar of Fig. 2e.

Fig. 2 Trajectories of the head monomers of eight non-interacting active filaments, each initially positioned along a different great circle on a sphere. Each filament is shown in a distinct color: the blue trajectory corresponds to the filament initialized on the equatorial great circle (i.e., the one passing through the center of the sphere), while the others lie on parallel great circles below the equator. As the surface gradually deforms, the filaments continue to follow geodesic paths. (a) Filaments on an undeformed spherical surface (a = 0). (b) Filaments on a prolate spheroid (a = 0.5b). (c) Filaments on a peanut-shaped surface (a = 0.9b). (d) Filaments on a peanut-shaped surface with a narrow channel (a = 0.99b). (e) Gaussian curvature K on Cassini ovals for increasing values of a. Surfaces are color-coded by Gaussian curvature, with the formula for K given in Appendix B.

To explore the stability of geodesic trajectory, we gradually tune the surface shape by tuning values of a. We begin by placing a filament along a great circle on the sphere (a = 0) as a reference case. After allowing it to complete a full revolution, we incrementally deform the surface by increasing a up to a = 0.99b, repeating the simulation and considering geodesic curvatures k_g ≤ 10⁻¹² effectively zero. Starting from a sphere, increasing a causes the filament to continue following a geodesic. At a = 0.5b (prolate spheroid, Fig. 2b), k_g and τ_g remain near zero, indicating it follows a geodesic that also minimizes bending energy. This remains true at a = 0.9b and a = 0.99b (peanut-shaped surfaces, Fig. 2c and d). In the most deformed case (Fig. 2d), k_g and τ_g increase slightly but stay very small, suggesting the filament remains close to a geodesic. This slight increase is likely due to numerical errors.

Having established how curvature variations, activity, and bending rigidity govern the dynamics of individual filaments, we now turn to the collective behavior of many interacting active filaments confined to curved surfaces. This shift is motivated by experimental observations of coordinated motion on closed geometries, where interactions become central to the emergent dynamics.³⁰ While our earlier analysis focused on trajectories initialized along specific geodesics, it is important to note that filament motion on a Cassini oval is also sensitive to initial conditions. We revisit this point at the end of the section, demonstrating how different initializations can lead to distinct patterns of spatial exploration. To connect with experimental relevance and gain insight into the collective behaviors, we begin by examining suspensions of active filaments.

Previous studies of dry active particles with polar or nematic alignment and short-range repulsion on spherical surfaces have shown that these systems can spontaneously form rotating polar bands around the equator.^18,19,41 This collective behavior arises from a balance between the particles’ tendency to follow geodesic paths (i.e., great circles) and their mutual steric repulsion, with bandwidth narrowing as activity increases. Building on our analysis of single-filament dynamics on curved spherical surfaces, we now explore the collective behavior of stiff active filaments at intermediate densities. In this regime, the filaments self-organize into a rotating band near the equator (Fig. 3a), closely mirroring the dynamics reported in prior simulations of active particles on curved geometries.^18,19,41 Notably, recent experiments³⁰ have observed analogous phenomena: actin filaments gliding along the inner surface of spherical vesicles exhibit surface–density–dependent transitions. As the filament concentration increases, these systems develop polar bands and display two additional dynamical regimes: polar vortices at off-equatorial latitudes at intermediate coverage, and disordered filament clusters at high coverage. Our simulations, performed under comparable conditions, reproduce the emergence of equatorial bands at low densities, in agreement with these experimental observations. We do not not observe the additional dynamical regimes, as this modeling approach is not designed to explore the high coverage states where hydrodynamics are likely to be an important factor.

Fig. 3 Steady-state configurations of stiff active filaments κ = 100 at intermediate activity for four values of the Cassini oval parameter a: (a) sphere (a = 0), (b) egg-like shape (a = 0.5b), (c) peanut-like shape (a = 0.9b), and (d) peanut-like shape with a narrow channel (a = 0.99b). Snapshots correspond to an intermediate density of 0.16, with filaments individually color-coded. On spherical-like surfaces (panels (a)–(c)), a rotating band state emerges. In contrast, for the narrow-channel peanut (panel (d)), where the curvature gradient between lobes is pronounced, filaments become trapped on one side. (e)–(h) Single-filament dynamics initialized from the corresponding configurations in panels (a)–(d). Filament colors match those in the first row for visual consistency.

To understand how non-uniform curvature affects collective dynamics, we further explore filament organization on Cassini ovals by varying the parameter a. For a = 0.5b, the surface has the shape of a prolate spheroid, and similar to their behavior on a sphere, filaments form a rotating band (Fig. 3b). A comparable pattern emerges at a = 0.9b, where the surface becomes peanut-shaped (Fig. 3c). However, at a = 0.99b, the neck connecting the lobes narrows substantially, creating a pronounced curvature gradient. In this regime, filaments tend to become trapped on one side of the surface (Fig. 3d). As a result, the equatorial band no longer forms, and motion becomes localized to a single lobe. Note that due to noise in the simulations, the filament configurations shown in Fig. 3a–d do not exactly follow geodesics (Appendix E).

As discussed in Section 4.1, some geodesics on Cassini ovals with a > 0 resemble those on the sphere, so the similar filament organization observed for 0 < a ≤ 0.9b is not surprising. In contrast, the qualitative change at a = 0.99b, where filaments become confined to one lobe, suggests a transition in the dynamics. This raises the question of whether the behavior in Fig. 3d is driven solely by surface curvature or also influenced by crowding.

To address this, we analyzed the dynamics of single filaments initialized with the same configurations as in Fig. 3a–d. As shown in Fig. 3e–h, single filaments explore the entire surface without preferring specific regions. We verified that this single-filament behavior persists even in the absence of noise (see Appendix E). This indicates that, in the absence of interactions, individual filaments do not favor any particular area. These results suggest that the trapping observed at a = 0.99b in Fig. 3d does not arise from the intrinsic dynamics of single filaments, but likely results from crowding effects.

Overall, we conclude that in the dilute regime, filaments on the Cassini oval move collectively across a wide range of the parameter a. Notably, the flocking behavior observed on curved spherical surfaces³⁰ is absent in the planar case at the same density,⁴⁵ suggesting that curvature promotes flocking in this dilute regime, consistent with findings for polar and nematic active systems.^18,19,41

4.3 Trapping effect as a function of density

To further investigate the trapping effect observed for a = 0.99b, we compute the probability that filaments cross the neck between the two lobes (Fig. 4). In the “peanut-shape” geometry, the left lobe corresponds to negative x values, the right lobe to positive x, and the neck is located at x = 0. To quantify the crossing probability, we first let the system equilibrate and then monitor it over a long observation window of t = 10⁵ time units. This duration is chosen to ensure that filaments have sufficient time to explore the surface, so that those which appear confined are truly trapped. During this period, we count how many filaments change the sign of their x-coordinate at least once. The crossing probability is defined as the fraction of filaments that exhibit such a sign change. Since we compare systems with different densities, for ρ < 0.16 we average over multiple initial conditions to match the total number of filaments analyzed in the highest-density case (ρ = 0.16).

Fig. 4 Trapping–escaping transition. (a) Probability of crossing P_c as a function of inverse density ρ⁻¹. At high densities (ρ⁻¹ < 4.5 × 10¹), filaments are trapped on one side of the peanut-shaped surface. At lower densities, all filaments escape and explore both lobes. (b) Probability of crossing P_c as a function of Péclet number Pe at fixed density ρ = 0.16. The behavior is similar to panel (a): at low Pe, filaments are trapped, while increasing Pe enables exploration of both lobes. (c) Band-like behavior observed at ρ⁻¹ = 2.2 × 10¹; similar dynamics occur at higher densities. This configuration matches that seen for a = 0.9b (see Fig. 3c). For ρ⁻¹ < 4.5 × 10¹, behavior resembles Fig. 3d. (d) Snapshot at ρ = 0.16 and Pe = 8 × 10⁻³. Unlike lower densities, filaments explore the full surface but do not form a coherent rotating band.

Measurements were performed at fixed Péclet number (Pe = 10⁻³) across densities from ρ = 0.16 down to ρ = 10⁻³, i.e. single filament (Fig. 4a). We observe that P_c, and thus surface exploration, increases as density decreases. In other words, lower densities promote broader exploration. At low densities, filaments form band-like trajectories (Fig. 4c) similar to those on the peanut-shaped surface (Fig. 3c), collectively crossing the neck. These findings suggest that at this Péclet number, filaments tend to move cohesively. When density is high and the neck is too narrow to accommodate all filaments, they remain confined to a single lobe. At lower densities, collective neck crossing becomes more likely, allowing exploration of the full surface. Moreover, we measured the time required for a filament to cross to the opposite lobe. When crossings occur, this time remains approximately constant across densities, with an average value of t_c = 2.5 × 10³ time units. Moreover, we verified that the trapping effect persists even in the overdamped limit (see Appendix F).

At fixed density (ρ = 0.16), as Pe increases, the probability of channel crossing P_c increases, allowing filaments to explore the surface uniformly (Fig. 4b). Unlike the low-density regime, where filaments form a single rotating band around the “equator” of the peanut-shaped surface, at higher Pe they move more uniformly without organizing into a coherent structure (Fig. 4d). They are no longer constrained to move collectively. These results explain the distinct behavior of the Cassini oval with a = 0.99b compared to those with a < 0.99b. Overall, these findings highlight how surface geometry can strongly influence filament dynamics and suggest that curvature and topology could be used to control or restrict the motion of active filaments or particles.⁴¹

5 Conclusions

In this paper, we investigated the dynamics of both individual and collective active filaments confined to curved surfaces using numerical simulations. Our goal was to understand how surface curvature, filament bending energy, and activity interact to shape filament behavior.

For single filaments on axisymmetric surfaces such as Gaussian bumps, we found that filaments follow geodesic paths when activity dominates. In contrast, when bending rigidity is strong, filaments deviate from geodesics, as these no longer minimize their elastic energy. On spherical surfaces, where geodesics coincide with energy-minimizing configurations, filaments consistently follow great circles regardless of activity. When deforming the sphere into Cassini ovals, filaments continue to follow geodesic paths that also minimize bending energy, highlighting a robustness of geodesic guidance in this regime.

We introduced a geometry-dependent Péclet number, which captures how the timescale of filament motion depends on surface curvature. In dilute suspensions on Cassini ovals, individual filaments explore the full surface, while collectively they self-organize into rotating bands, typically aligned along the equator. For highly deformed geometries (e.g., a Cassini oval with a = 0.99b), we observed a transition: at low density, rotating bands still form, but at higher densities, filaments become trapped in one lobe and fail to explore the entire surface due to spatial crowding. A similar transition occurs as a function of activity: at low Péclet number, filaments localize, while higher activity enables full surface exploration. However, in this case, filaments exhibit more disordered motion, and coherent bands do not emerge. Finally, we show that on both spherical and Cassini oval surfaces, the filament dynamics remain consistent with observations in flat space^32,45,53 (see Appendix G).

While our model captures the essential physics of active semiflexible filaments constrained to curved geometries, it abstracts away many biochemical and structural complexities present in filaments found in biological systems. As such, it can be informative for understanding the key physical mechanisms underlying experimental observations of cytoskeletal filaments confined to move on the surface of lipid vesicles.³⁰ More broadly, it may provide conceptual insight into cytoskeletal in vivo organization during processes such as cell division or polarity establishment. For example, the formation of a coherent band around one of the lobes of the Cassini oval suggests a curvature-guided localization mechanism that could partially underlie the accumulation of actomyosin at the cleavage furrow during cytokinesis.^59,60 Conversely, in asymmetrically dividing cells, similar curvature or confinement effects might promote cortical actomyosin segregation toward one pole,^61,62 contributing to polarity establishment. In both cases, the same underlying principle—activity-driven filament alignment modulated by surface geometry may help bias cytoskeletal organization in space, even in the absence of explicit biochemical cues.

However, while this simplified model is useful for understanding the key underlying physical principles, it is important to remember that real biological active filaments are more complex than discussed here. For example, cytoskeltal filaments such as actin and microtubules are inherently dynamic and polar due to their well-known treadmilling behavior, which leads to directed growth of the filament in a preferred direction. Additionally, the motion of the cytoskeleton is largely driven by molecular motors, which act to produce relative motion of the filaments, rather than an absolute self-propulsion as considered in this work.⁶³ Both of these effects, as well as many others, are not captured by our model, and so we expect, in many cases, a more complex model will be required to understand specific biological processes.

Overall, our results demonstrate that surface curvature and topology can be leveraged to guide or constrain the behavior of active filaments.⁴¹ This highlights the potential of geometric design as a mechanism to control active matter in soft materials, microfluidic devices, and bio-inspired systems.⁶⁴

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are available from the corresponding author upon reasonable request.

Appendices

A Geodesic and normal curvature

To understand the behavior of filaments moving on a Gaussian bump, we compute the normal and geodesic curvatures, as well as the geodesic torsion. When working with curves embedded in a surface, it is convenient to use the Darboux frame. In this framework, we define the surface unit normal n, the tangent vector to the curve t, and the unit tangent normal l = t × n.

The Darboux frame is characterized by structure equations analogous to the Frenet–Serret equations.⁴⁹ These are,

From these relations, the geodesic curvature, normal curvature, and geodesic torsion can be defined aswhere the curvature vector is defined by .

In simulations, each filament is represented as a sequence of discrete points, and all spatial derivatives are approximated using finite-difference schemes.⁶⁵ Specifically, the arc-length increments are computed as

Δs_i = |r_i+1 − r_i|, (A3)

where r_i denotes the position of the ith point along the filament. To compute derivatives with respect to arc length s, we use NumPy's⁶⁶ np.gradient function, which applies second-order accurate central differences.

B Gaussian curvature: Cassini ovals and Gaussian bump

The Cassini surface is defined implicitly, f(x, y, z) = 0. For such implicitly defined surfaces, the Gaussian curvature K can be expressed in terms of the gradient ∇f and the Hessian matrix H(f) as⁶⁷

We compute the curvature using sympy⁶⁸ to get

where r² = x² + y² + z². In the special case a = 0, corresponding to a sphere of radius r, this expression reduces to the expected result for a sphere, i.e., .

For a Gaussian-bump surface z = h(r), the Gaussian curvature is

where h′(r) and h″(r) denote the first and second derivatives of the height function h(r).

C Comparing continuum and discrete bending energy

Fig. 5 shows the relative error in the discrete bending energy E_d compared to the continuum limit (E_c), for a filament positioned on top of the Gaussian bump. The continuum bending energy is computed by evaluating eqn (1) numerically, using Simpson's rule to approximate the integral.⁶⁹ As the monomer spacing Δs increases, the error grows accordingly.

Fig. 5 Error in the discrete bending energy E_d compared to the continuum bending energy E_c, plotted as a function of monomer spacing Δs, evaluated when the polymer head is located at the peak of the Gaussian bump (x = 0). The dashed line is the case of a discrete polymer in a plane. while colored curves correspond to different bump amplitudes, A.

D Scaling argument for filament on Gaussian bump

In Section 4.1, we saw through numerical results that depending on the value of Pe, when an active filament is propelled towards a Gaussian bump by its active propulsion, it may pass directly over, maintaining a geodesic trajectory, or deviate around the bump. In this section, we seek to understand this behavior by examining the competition between activity and bending rigidity through a simple scaling argument.

We move to the overdamped limit of the filament dynamics and compare the time taken for the filament to pass over the bump, driven by its active self-propulsion t_a, with the time taken for the filament to deform due to its elastic bending forces t_b.

In the overdamped limit, where internal forces are typically balanced by external friction γv ≈ f. We can estimate the time to move over the entire filament over the bump when propelled by a constant force of magnitude f_p as . (Note this is valid as long as the size of the filament is much larger than that of the bump.) To compare this to t_b we need to estimate the force caused by bending when passing over the bump. If we assume that the curvature angle induced by the bump is approximated by , then the bending energy penalty the filament pays for passing over the bump will be . The force on the filament is then given , where l_s is the length scale discussed in Section 3, over which the filament senses the curvature. The filament is likely to experience this bending force until its path has deviated by a distance l_s, since this is approximately how far it must move to avoid sensing the bump. Therefore, the time taken minimize the bending penalty by avoiding the bump will be , since .

This means whether the filament avoids or passes over the bump is quantified by the dimensionless number , where Pe is the earlier defined geometry-dependent Péclet number. Using the expression for l_s for Gaussian bump, we have . This means for large Pe we expect t_a ≪ t_b, so the filament is expected to remain on a geodesic, consistent with our numerical results.

E Single active filaments on Cassini ovals

Fig. 6 shows the trajectories of single filaments without noise, initialized from the configurations shown in Fig. 3a–d. In the spherical case (Fig. 6a), filaments remain confined to the equator. In contrast, for nonzero values of the parameter a, the filaments progressively explore the entire surface and are no longer restricted to a specific region (Fig. 6b–d).

Fig. 6 Steady-state configurations of single filaments without noise at intermediate activity for four different values of the Cassini oval parameter aa: (a) sphere (a = 0); (b) prolate spheroid (a = 0.5b); (c) peanut-shaped surface (a = 0.9b); and (d) peanut-shaped surface with a narrow channel (a = 0.99b). Each filament is initialized from the corresponding configuration shown in Fig. 3a–d.

Fig. 7 presents the probability distributions of geodesic curvature corresponding to the configurations shown in Fig. 3. For the spherical case (Fig. 3a and e), the distributions for a single filament (Fig. 3e) and for the filament suspension (Fig. 3a) are nearly identical, both peaking around k_g ≈ 0. Due to the presence of noise, the values of k_g in the single-filament case are slightly higher than in the noiseless case. In contrast, for Cassini ovals with a > 0, the probability of observing zero geodesic curvature is noticeably lower in the filament suspension than in the single-filament case (Fig. 3b–d).

Fig. 7 Probability distributions of the geodesic curvature corresponding to the snapshots shown in Fig. 3. Green squares represent the case of filament suspensions (Fig. 3a–d), while blue circles correspond to single filaments with noise (Fig. 3e–h). (a) Sphere, (b) prolate spheroid, (c) peanut-shaped surface, and (d) peanut-shaped surface with a narrow channel.

F Trapping effect in the overdamped limit

Fig. 8 shows the probability of crossing, P_c, as a function of the friction coefficient for ρ = 0.16, using the same Pe = 10⁻³ as in Fig. 4(a). The results indicate that, in the overdamped limit (γ = 10²), the probability of crossing is slightly smaller than in the γ = 10⁰ case. This suggests that trapping is more pronounced in the overdamped regime. Nevertheless, consistent with previous observations on the plane,^44,45 the qualitative behavior (i.e., the presence of trapping) remains unchanged.

Fig. 8 Probability of crossing P_c as a function of the friction coefficient γ for filaments moving on a narrow peanut-shaped surface.

G Behavior of dilute suspensions as a function of bending rigidity

Fig. 9 shows the probability distributions of the geodesic curvature for two values of the bending rigidity, κ = 10⁰ and κ = 10², on the spherical and Cassini oval surfaces. The probability of motion along geodesics, P(k_g), is significantly lower for flexible filaments (κ = 10⁰) than for stiffer ones (κ = 10²). At low bending rigidity, consistent with observations in flat space,^32,45,53 the filaments tend to form spiral structures.

Fig. 9 Probability distributions of the geodesic curvature for bending rigidities κ = 10⁰ and 10², shown as blue dots and purple triangles, respectively. Panels: (a) sphere; (b) peanut-shaped surface; (c) peanut-shaped surface with a narrow channel. Snapshots of the systems at low bending rigidity are also shown for reference.

Acknowledgements

E. D. M. acknowledges funding by the endowed E. N.&M. N. Lindsay PhD studentship. R. S. acknowledges support from the UK Engineering and Physical Sciences Research Council (Award EP/W023946/1). D. M. F. and G. J. acknowledge support from the Comunidad de Madrid and the Complutense University of Madrid (Spain) through the Atracción de Talento program (grant no. 2022-T1/TIC-24007). D. M. F. also acknowledges support from MINECO (grant no. PID2023-148991NA-I00).

Notes and references

1. M. C. Simon, Dev. Cell, 2004, 6, 479–482.
2. A. Fatima, A. Culver, F. Culver, T. Liu, W. H. Dietz, B. R. Thomson, A.-K. Hadjantonakis, S. E. Quaggin and T. Kume, Dev. Dyn., 2014, 243, 957–964.
3. J. Whisler, S. Shahreza, K. Schlegelmilch, N. Ege, Y. Javanmardi, A. Malandrino, A. Agrawal, A. Fantin, B. Serwinski, H. Azizgolshani, C. Park, V. Shone, O. O. Demuren, A. D. Rosario, V. L. Butty, N. Holroyd, M.-C. Domart, S. Hooper, N. Szita, L. A. Boyer, S. Walker-Samuel, B. Djordjevic, G. K. Sheridan, L. Collinson, F. Calvo, C. Ruhrberg, E. Sahai, R. Kamm and E. Moeendarbary, Sci. Adv., 2023, 9, eadg9781.
4. D. S. Vijayraghavan and L. A. Davidson, Birth Defects Res., 2017, 109, 153–168.
5. K. Kostanjevec, R. Sknepnek, J. M. Collinson and S. Henkes, bioRxiv, 2024, preprint, 2024-08.
6. B. Schamberger, R. Ziege, K. Anselme, M. Ben Amar, M. Bykowski, A. P. Castro, A. Cipitria, R. A. Coles, R. Dimova and M. Eder, et al. , Adv. Mater., 2023, 35, 2206110.
7. M. Chuai, W. Zeng, X. Yang, V. Boychenko, J. A. Glazier and C. J. Weijer, Dev. Biol., 2006, 296, 137–149.
8. P. Buske, J. Galle, N. Barker, G. Aust, H. Clevers and M. Loeffler, PLoS Comput. Biol., 2011, 7, e1001045.
9. S. J. Callens, R. J. Uyttendaele, L. E. Fratila-Apachitei and A. A. Zadpoor, Biomaterials, 2020, 232, 119739.
10. G. Cereghetti and L. Scorrano, Oncogene, 2006, 25, 4717–4724.
11. L. Metselaar, J. M. Yeomans and A. Doostmohammadi, Phys. Rev. Lett., 2019, 123, 208001.
12. D. Matoz-Fernandez, F. A. Davidson, N. R. Stanley-Wall and R. Sknepnek, Phys. Rev. Res., 2020, 2, 013165.
13. A. Mietke, F. Jülicher and I. F. Sbalzarini, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 29–34.
14. D. Khoromskaia and G. Salbreux, eLife, 2023, 12, e75878.
15. L. A. Hoffmann, L. N. Carenza, J. Eckert and L. Giomi, Sci. Adv., 2022, 8, eabk2712.
16. F. Vafa and L. Mahadevan, Phys. Rev. Lett., 2022, 129, 098102.
17. Z. Wang, M. C. Marchetti and F. Brauns, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2220167120.
18. R. Sknepnek and S. Henkes, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 022306.
19. S. Henkes, M. C. Marchetti and R. Sknepnek, Phys. Rev. E, 2018, 97, 042605.
20. S. Ehrig, J. Ferracci, R. Weinkamer and J. W. Dunlop, Phys. Rev. E, 2017, 95, 062609.
21. I. R. Bruss and S. C. Glotzer, Soft Matter, 2017, 13, 5117–5121.
22. P. Castro-Villarreal and F. J. Sevilla, Phys. Rev. E, 2018, 97, 052605.
23. P. W. A. Schönhöfer and S. C. Glotzer, Soft Matter, 2022, 18, 8561–8571.
24. L. M. Janssen, A. Kaiser and H. Löwen, Sci. Rep., 2017, 7, 5667.
25. P. Iyer, R. G. Winkler, D. Fedosov and G. Gompper, Phys. Rev. Res., 2023, 5, 033054.
26. F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic and A. R. Bausch, Science, 2014, 345, 1135–1139.
27. P. W. Ellis, D. J. Pearce, Y.-W. Chang, G. Goldsztein, L. Giomi and A. Fernandez-Nieves, Nat. Phys., 2018, 14, 85–90.
28. F. Alaimo, C. Köhler and A. Voigt, Sci. Rep., 2017, 7, 5211.
29. D. Pearce, P. W. Ellis, A. Fernandez-Nieves and L. Giomi, Phys. Rev. Lett., 2019, 122, 168002.
30. C.-P. Hsu, A. Sciortino, Y. A. de la Trobe and A. R. Bausch, Nat. Commun., 2022, 13, 1–8.
31. H. Jiang and Z. Hou, Soft Matter, 2014, 10, 1012–1017.
32. R. E. Isele-Holder, J. Elgeti and G. Gompper, Soft Matter, 2015, 11, 7181–7190.
33. Ö. Duman, R. E. Isele-Holder, J. Elgeti and G. Gompper, Soft Matter, 2018, 14, 4483–4494.
34. V. Bianco, E. Locatelli and P. Malgaretti, Phys. Rev. Lett., 2018, 121, 217802.
35. R. G. Winkler and G. Gompper, J. Chem. Phys., 2020, 153, 040901.
36. G. Janzen, J. P. Miranda, J. Martín-Roca, P. Malgaretti, E. Locatelli, C. Valeriani and D. A. M. Fernandez, J. Chem. Phys., 2025, 162, 114905.
37. R. G. Winkler, J. Chem. Phys., 2025, 162, 154903.
38. G. S. Manning, Q. Appl. Math., 1987, 45, 515–527.
39. J. Guven and P. Vázquez-Montejo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 026603.
40. J. Guven, D. M. Valencia and P. Vázquez-Montejo, J. Phys. A: Math. Theor., 2014, 47, 355201.
41. E. D. Mackay, G. Janzen, D. Matoz-Fernandez and R. Sknepnek, arXiv, 2025, preprint, arXiv:2505.24730 DOI:10.48550/arXiv.2505.24730.
42. R. E. Isele-Holder, J. Elgeti and G. Gompper, Soft Matter, 2015, 11, 7181–7190.
43. K. Prathyusha, S. Henkes and R. Sknepnek, Phys. Rev. E, 2018, 97, 022606.
44. M. Fazelzadeh, E. Irani, Z. Mokhtari and S. Jabbari-Farouji, Phys. Rev. E, 2023, 108, 024606.
45. G. Janzen and D. A. Matoz-Fernandez, Soft Matter, 2024, 20, 6618–6626.
46. H. Nickerson and G. S. Manning, Geometriae Dedicata, 1988, 27, 127–136.
47. R. Capovilla, C. Chryssomalakos and J. Guven, J. Phys. A: Math. Gen., 2002, 35, 6571.
48. D. María Valencia, G. Torres-Vargas and R. Castañeda Priego, J. Phys. A: Math. Theor., 2019, 52, 395203.
49. M. P. Do Carmo, Differential Geometry of Curves and Surfaces, Prentice Hall, Englewood Cliffs, NJ, 1976.
50. J. F. Marko and E. D. Siggia, Macromolecules, 1995, 28, 8759–8770.
51. P. Castro-Villarreal and J. E. Ramírez, Phys. Rev. E, 2019, 100, 012503.
52. M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143.
53. K. R. Prathyusha, S. Henkes and R. Sknepnek, Phys. Rev. E, 2018, 97, 022606.
54. H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 011901.
55. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
56. R. Sknepnek, SAMoS: Self-propelled Agent-based Models with SAMoS, https://github.com/sknepneklab/SAMoS, 2024, Accessed: 2024-04-29.
57. B. Leimkuhler and C. Matthews, Molecular Dynamics: With Deterministic and Stochastic Numerical Methods, Springer International Publishing, 2015.
58. N. Fabian, K. Moreland, D. Thompson, A. C. Bauer, P. Marion, B. Gevecik, M. Rasquin and K. E. Jansen, 2011 IEEE symposium on large data analysis and visualization, 2011, pp. 89–96.
59. P. Chugh and E. K. Paluch, J. Cell Sci., 2018, 131, jcs186254.
60. T. D. Pollard and B. O'Shaughnessy, Annu. Rev. Biochem., 2019, 88, 661–689.
61. T. T. Pham, A. Monnard, J. Helenius, E. Lund, N. Lee, D. J. Müller and C. Cabernard, iScience, 2019, 13, 9–19.
62. C. H. Oon and K. E. Prehoda, eLife, 2019, 8, e45815.
63. N. Chaffey, B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts and P. Walter, Molecular biology of the cell, 4th edn, 2003.
64. E. Cereceda-López, M. De Corato, I. Pagonabarraga, F. Meng, P. Tierno and A. Ortiz-Ambriz, Soft Matter, 2024, 20, 7615–7622.
65. B. Fornberg, Math. Comput., 1988, 51, 699–706.
66. C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke and T. E. Oliphant, Nature, 2020, 585, 357–362.
67. R. Goldman, Comput.-Aided Geometric Des., 2005, 22, 632–658.
68. A. Meurer, C. P. Smith, M. Paprocki, O. Čertík, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh, T. Rathnayake, S. Vig, B. E. Granger, R. P. Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johansson, F. Pedregosa, M. J. Curry, A. R. Terrel, V. Roučka, A. Saboo, I. Fernando, S. Kulal, R. Cimrman and A. Scopatz, PeerJ Comput. Sci., 2017, 3, e103.
69. W. H. Press, Numerical recipes: The art of scientific computing, Cambridge University Press, 3rd edn, 2007.

---

Footnote

† G. J. and E. D. M. contributed equally to this work.

---