◌

# Complex motion of steerable vesicular robots filled with active colloidal
rods

Sophie Y. Lee  Department of Materials Science and Engineering, University of
Michigan, Ann Arbor, Michigan 48109, USA  Philipp W. A. Schönhöfer  Department
of Chemical Engineering, University of Michigan, Ann Arbor, Michigan 48109,
USA  Sharon C. Glotzer<EMAIL_ADDRESS>Department of Materials Science and
Engineering, University of Michigan, Ann Arbor, Michigan 48109, USA
Department of Chemical Engineering, University of Michigan, Ann Arbor,
Michigan 48109, USA  Biointerfaces Institute, University of Michigan, Ann
Arbor, Michigan 48109, USA

(2062023)

###### Abstract

While the collective motion of active particles has been studied extensively,
effective strategies to navigate particle swarms without external guidance
remain elusive. We introduce a method to control the trajectories of two-
dimensional swarms of active rod-like particles by confining the particles to
rigid bounding membranes (vesicles) with non-uniform curvature. We show that
the propelling agents spontaneously form clusters at the membrane wall and
collectively propel the vesicle, turning it into an active superstructure. To
further guide the motion of the superstructure, we add discontinuous features
to the rigid membrane boundary in the form of a kinked tip, which acts as a
steering component to direct the motion of the vesicle. We report that the
system’s geometrical and material properties, such as the aspect ratio and
Péclet number of the active rods as well as the kink angle and flexibility of
the membrane, determine the stacking of active particles close to the kinked
confinement and induce a diverse set of dynamical behaviors of the
superstructure, including linear and circular motion both in the direction of,
and opposite to, the kink. From a systematic study of these various behaviors,
we design vesicles with switchable and reversible locomotions by tuning the
confinement parameters. The observed phenomena suggest a promising mechanism
for particle transportation and could be used as a basic element to navigate
active matter through complex and tortuous environments.

##  I Introduction

The ongoing and intense interest in active matter inarguably involves its
promise as a precursor to microscopic colloidal robots. If the emergent,
collective motion of self-propelled particles central to active matter could
be directed and harnessed for work, colloidal machines that locomote, capture
and retrieve tiny objects, and other robotic functions could be achieved for
applications ranging from in vivo health diagnostics and drug delivery to
stealth. However, despite the rapid, ongoing development of new synthesis
techniques [1] and computational models [2, 3] of particles that convert
external energy into an internal driving force, their controlled and efficient
collective transport to mimic robotic entities still faces major technical
challenges. Such challenges include the robustness of active swarm navigation
against thermal noise and other perturbations, enabling switchable and
reversible mechanisms to direct transport, implementing programmable sensing,
and triggering responses to external stimuli and complex environments [4, 5].
A variety of external control strategies have been applied to realize sensing
and directed navigation [6], such as the use of external fields [7, 8, 9, 10]
or an external feedback loop [11, 12, 13, 14, 15]. In contrast to external
control strategies, autonomous control strategies could be advantageous for
many applications, but these strategies have received considerably less
attention in the literature.

One such autonomous approach might be to design active particle
superstructures that permit autonomous navigation as a result of their
intrinsic and self-emergent behavior. Here we can draw inspiration from
macroscopic robotic swarms. Particle robots with pre-programmed interactions
can perform robotic tasks such as targeted morphological changes [16] or
synchronized motion [17], but such high level of internal logic is not
necessary. For example, stochastic robotic collectives are capable of global
locomotion and are able to achieve complex goals by transforming task-
incapable single components into task-capable robotic swarms [18]. A promising
example is that of macroscopic, forward-propelling, rod-like robotic units
that self-organize when trapped inside a deformable confining boundary by
accumulating at the boundary [19, 20], turning the system into a moving
superstructure [21, 22]. In an effort to engineer cellular migration, similar
concepts have been reported in studies of shape modulation and transport
dynamics of vesicular membranes driven by encapsulated active components [23,
24, 25, 26, 27, 28, 29, 30, 22, 31, 32]. For instance, motile bacteria have
been used to drive droplets inside an emulsion [24] or to generate cell
extrusion in lipid vesicles before they propel the vesicle forward [23].

Inspired by these autonomous superstructure strategies, in this paper we
expand on the idea of encapsulation of active particles as a device to control
their collective motion by adding an asymmetric bias to the superstructure
model. Bias is a known requirement for achieving asymmetric behavior in active
matter systems [33, 34, 22, 31], and here we achieve asymmetry via the
geometry of the confining boundary. Using molecular dynamics simulation (see
methods section), we confine self-propelling rod-shaped particles inside rigid
vesicles whose discontinuous curvature creates kinks in the vesicle boundary.
We first show that the active rods tend to form clusters at the kink, similar
to clusters of active particles observed in channels [35, 36, 37], under
polygonal confinement [38] or near a chevron-shaped trap [39, 40]. The
resulting bias in the location of the emergent cluster leads to a tunable and
directed propulsion of the vesicular superstructure. Based on the properties
of the vesicle (rigidity and kink angle) and the active rods (Péclet number,
aspect ratio, particle density inside the vesicle) we observe linear motion
and circular motion, both in the direction of the kink (kink-forward) and
opposite to the direction of the kink (kink-backwards). We relate these four
different types of motion to the alignment and stacking of the active rods at
the kink and showcase how this mechanism can be used to create vesicles with
switchable and reversible dynamics.

##  II Method

Figure 1: Active rods confined in a vesicle. Sketch of a vesicle containing
rigid active rods (dark blue) with (a) aspect ratio $\alpha=3$ aligning
parallel to the long axis of the vesicle (beige) at a wide
kink($\theta_{kink}=5\pi/6$), and (b) $\alpha=2$ aligning at an angle with one
side of the vesicle wall (beige) at a narrow kink ($\theta_{kink}=\pi/2$). The
upper inset shows the (a) trimer ($\alpha=3$) rod, and (b) dimer ($\alpha=2$)
rod. The lower inset is a close-up view of the kinked part of the vesicle and
the stacking of active rods where (a) has parallel and (b) has angled
stacking.

We performed Brownian molecular dynamics simulations of self-propelled rod-
shaped particles confined inside a vesicle in a square box with length
$L_{b}=125\sigma$. We modeled the vesicles as chains of $N_{v}$=100 bonded
disks with diameter $\sigma=1$. The active rod-shaped particles are modeled as
a rigid body of $m={2,3}$ disks with diameter $\sigma$ that are linked
linearly end-to-end (see Fig 1) [41, 42]. The number of active rods per
vesicle perimeter varies between
$\rho_{N_{a}}=\frac{N_{a}}{N_{v}\sigma}\in\\{0.05,0.1,\dots,0.3\\}$. We chose
this definition of density to compare 2- and 3-rod systems that apply the same
amount of force to the vesicle wall. The equations of motion of each vesicle
disk and self-propelling rod is given by Brownian dynamics:

$\displaystyle\gamma\dot{\mathbf{r}_{i}}$
$\displaystyle=\sum_{j}\mathbf{F}_{ij}^{WCA}+F_{i}^{A}\hat{\mathbf{e}}_{i}+\mathbf{F}_{i}^{Ves}+\sqrt{2\gamma
k_{B}T}\eta_{i}(t)$ (1) $\displaystyle\gamma_{r}\dot{\omega_{i}}$
$\displaystyle=\sum_{j}\mathbf{T}_{ij}^{WCA}+\sqrt{2\gamma_{r}k_{B}T}\zeta_{i}(t)$
(2)

with translational drag coefficient $\gamma$ = 1, and rotational drag
coefficient $\gamma_{r}=\frac{\sigma^{3}\gamma}{3}$ per the Stokes-Einstein
relationship. The interparticle forces $\mathbf{F}_{ij}^{WCA}$ on each disk
and torques $\mathbf{T}_{ij}^{WCA}$ on each rod are derived from the purely
repulsive Weeks-Chandler-Anderson potential

$\displaystyle
U^{WCA}(r_{ij})=\begin{cases}4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]+\epsilon&r_{ij}<r_{c}=\sqrt[6]{2}\sigma\\\
0&r_{ij}>r_{c},\end{cases}$ (3)

where $r_{ij}$ is the distance between disk $i$ and $j$ and $r_{c}$ is the
cut-off distance. The functions $\eta_{i}(t)$ and $\zeta_{i}(t)$ are
normalized Gaussian white noise processes with zero mean
$\langle\eta_{i}(t)\rangle=0$, $\langle\zeta_{i}(t)\rangle=0$ and unit
variance
$\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$,
$\langle\zeta_{i}(t)\zeta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$.
The active force term $F_{i}^{A}\hat{\mathbf{e}}_{i}$ is only applied to the
rod center and points along the symmetry axis of the rod
$\hat{\mathbf{e}}_{i}=(\cos\omega_{i},\sin\omega_{i})$. The active force
magnitude $F_{i}^{A}=\text{Pe}\cdot k_{B}T$ is controlled by the Péclet number
$\text{Pe}\in\\{25,50,75\\}$ for the active rods and $F_{i}^{A}=0$ for the
vesicle particles. The final term $\mathbf{F}_{i}^{Ves}$ combines the vesicle
bond forces and, hence, is derived from the bond potentials

$\displaystyle U_{i}^{Ves}$
$\displaystyle=\sum_{j}U_{FENE}(r_{ij})+\sum_{k,l}U_{angle}(\theta_{ikl}),$
(4) $\displaystyle U_{FENE}(r_{ij})$
$\displaystyle=-\frac{1}{2}\kappa_{b}r_{0}^{2}\ln\left[1-\left(\frac{r_{ij}}{r_{0}}\right)^{2}\right],$
(5) $\displaystyle U_{angle}(\theta_{ikl})$
$\displaystyle=-\frac{1}{2}\kappa_{a}(\theta_{ikl}-\theta_{0})^{2}.$ (6)

Each set of neighboring disks $(i,j)$ comprising the vesicle boundary
(membrane) is bonded by the finitely extensible nonlinear elastic (FENE)
potential $U_{FENE}(r_{ij})$ [43] with bonding strength
$\kappa_{b}=100k_{B}T\sigma^{-2}$ and equilibrium distance $r_{0}=2\sigma$. To
apply a uniform bending rigidity to the vesicle membrane, we added a harmonic
angle potential with bending coefficient $\kappa_{a}=1000k_{B}T$ and rest
angle $\theta_{0}=\pi$ to every triplet $(ikl)$ of neighboring disks. We also
ran several simulations with smaller bending coefficients
$\kappa_{a}\in\\{50k_{B}T,125k_{B}T,250k_{B}T,500k_{B}T\\}$ (Fig. S1),
however, in the following all results are presented for rigid vesicle
membranes with $\kappa_{a}=1000k_{B}T$ unless explicitly stated otherwise.
Finally, we introduce a discontinuous element, which we refer to as a ”kink”,
by setting the rest angle $\theta_{0}$ at one of the triplets to
$\theta_{0}=\theta_{kink}\in\\{\frac{\pi}{3},\frac{5\pi}{12},\frac{\pi}{2},\frac{7\pi}{12},\frac{2\pi}{3},\frac{3\pi}{4},\frac{5\pi}{6}\\}$
(see Fig. 1).

We ran each simulation [44] at a temperature $k_{B}T=0.01$ for a time period
$t=2\times 10^{4}\tau$ with time step $\Delta t=1\times 10^{-4}\tau$ and unit
time $\tau=\sqrt{m\sigma^{2}/k_{B}T}$. To obtain proper statistics each point
in the parameter space is sampled over 200 independent replica simulations. We
used the open source molecular dynamics software HOOMD-blue [45] [v3.0.0] to
perform our simulations,the freud data analysis package[46] for cluster
analysis and the signac software package [47, 48] for data management.

##  III Results

To illustrate the range of dynamical behaviors observed for the vesicle
superstructure, we present their observed steady-state behavior as a function
of kink angle $\theta_{\text{kink}}$ and rod particle density $\rho_{N_{a}}$
for both dimer and trimer rod systems in Fig. 2 and SI Fig. S4. In all of our
simulation we observe that the microscopic driving forces of the active agents
are transferred to the vesicle in a coordinated manner, propelling the entire
vesicle superstructure. Flexible or uniformly curved vesicles move in random
and unpredictable directions (see Fig. 2) as seen in earlier studies [20].
However, at high bending rigidity $\kappa_{a}=1000k_{B}T$ and for non-uniform
membrane curvature, the activity of the rods is transferred to the vesicle in
a coordinated directed motion based on the position of the kink in the vesicle
membrane. Overall, we obtain four different patterns of motion for kinked
vesicles with different collective propelling mechanisms of the rod-shaped
particles: a motion in the direction of the kink (linear forward), a motion in
the opposite direction of the kink (linear backward), a kink forward circular
motion, in which the trajectories of the vesicles make loops, and a kink
backward circular motion (see Fig. 2). We describe each of the behaviors
below.

###  III.1 Linear motion

Vesicles driven by rigid dimers (Fig. 2a) or trimers (Fig. 2b) predominantly
exhibit linear motion at high $\theta_{kink}$. We identify that the linear
trajectories of vesicles are correlated to a parallel stacking (with respect
to the long axis of the vesicle) of the active particles inside the vesicle
near the kinked tip (see sketch in Fig. 1a). In this regime, the rod-shaped
particles first aggregate at the vesicle boundary individually and move along
the interface as reported in other studies of active matter in hard [19, 49]
and flexible confinement [50, 21]. Due to the discontinuity in curvature,
which is known to slow down and capture active particles [39], the rods
accumulate into dense stackings at the kink, align with the kink, and apply
collective local pressure to the vesicle boundary in the direction of the rod
alignment director. At low particle densities $\rho_{N_{a}}<0.15$ the density
distribution inside the vesicles indicates that the clusters are stable at the
kink in steady state (see Fig. 3a). Furthermore, the director of the aligned
clusters of both dimers and trimers is parallel to the symmetry axis of the
vesicle (see alignment angle distributions in Fig. 4), with particles
distributed roughly evenly on both sides of the kink (see Fig. 1a). As a
consequence, the vesicle superstructures are pushed forward with a small
variance in velocity and hardly any rotational component to its locomotion
(see velocity and angular velocity distribution in Fig. 4).

We find that the propulsion mechanism breaks down if the self-propelling force
is too weak ($Pe<50$) (see Fig. 3a), the vesicle membrane is too flexible
($\kappa_{a}<250k_{B}T$), or the perimeter of the vesicle exceeds 200$\sigma$.
In the former case, the rods are less likely to cluster, resulting in a more
uniform particle distribution inside the vesicle (see Fig. 3a) and a lower
yield of forward-moving vesicles. Additionally, clusters that do emerge at the
kink are prone to dissolving again due to the decrease in active particle
pressure. For flexible vesicles, we must consider shape fluctuations of the
membrane, which are otherwise negligible for the highest bending coefficient
$\kappa_{a}=1000k_{B}T$ (see Fig. S1). With decreasing $\kappa_{a}$, the
pressure applied to the vesicle walls by the active rods produces local
deformations and high curvature regions at the vesicle boundary. This creates
additional nucleation pockets for clustering that can capture more rods,
diminishing or eliminating the benefit of the kink. We also investigated
different vesicle perimeter lengths between $110\sigma$ and $500\sigma$ with
monomer diameter $\sigma$. We observed that the resulting circular and linear
motion of the vesicular superstructure guided by $\theta_{kink}$ gradually
start to fail when the perimeter exceeds 200$\sigma$. In these larger
vesicles, the effect that active particles spontaneously form crystallites
away from the kink and produce additional pockets of local high curvature to
the boundary cannot be neglected. Hence, the kink-induced-directionality is
lost.

The crossover from linear to circular motion (discussed below) occurs between
obtuse $\theta_{kink}=\frac{5\pi}{6}$ and acute angles
$\theta_{kink}=\frac{\pi}{2}$. While for most trimer systems the transition
shifts directly from linear forward to linear backward motion, we observe an
additional regime of predominantly linear backward-moving vesicles for dimer
systems and trimer systems with low activity. Here we distinguish between two
mechanisms. For $0.1<\rho_{N_{a}}<0.15$, the geometry of the kink, while still
inducing clustering, prevents the stability of large parallel stackings.
Hence, in some of our simulations this geometric frustration destabilizes rod
alignment either completely (see the symbol $\times$ in Fig.3) yet does not
break up the cluster completely, or partially shifts it to the opposite side
of the vesicle. Here the boundary features a second local curvature maximum
and the parallel rod packing is more stable. This causes the forward motion of
vesicles to become less reliable with decreasing kink angle and eventually to
flip to backward linear motion for kink angles close to
$\theta_{\text{kink}}=\frac{7\pi}{12}$. We observed a similar phenomenon in
the trimer-filled vesicles with increasing particle density, but the effect is
not as strong as in the dimer systems and only occurs at lower $Pe=25$ and
over a narrower region of
$\frac{\pi}{2}\leq\theta_{\text{kink}}\leq\frac{2\pi}{3}$.

The second mechanism that reverses the direction of motion of the dimer-filled
vesicle occurs for $\rho_{N_{a}}>0.15$. Here the dimers form two clusters, one
located at the kink and the other located at the opposite side of the vesicle.
For large obtuse $\theta_{\text{kink}}$ the cluster at the kink is larger,
effectively pushing the vesicle forward with high reliability. For smaller
$\theta_{\text{kink}}$ and less stable alignment of rods at the kink, the
cluster of rods at the opposite side of the vesicle becomes larger, eventually
dominating and reversing the direction of motion of the vesicle
superstructure. The two counteracting clusters also explain the decrease in
the velocity of the dimer-filled vesicles (see Fig. 3a) best observed in the
velocity distribution in Fig. 4a, where both fast vesicles with one dominant
cluster (peak at high velocity) and slow vesicles with two counteracting
clusters (peak close to $v=0$) are present. In contrast, trimer-filled
vesicles, where only one large cluster forms, move faster with more
contributing active agents.

###  III.2 Circular motion

In addition to linear motion, vesicle superstructures also exhibit circular
motion where they rotate with a mean angular velocity
$\langle\dot{\omega}\rangle$ around the superstructure’s center of mass while
moving linearly either forward or backward. While circular motion can also, in
principle, occur from an unbalanced distribution of parallel stacked particles
around the kink, that mechanism causes only a slight bend in the vesicle’s
trajectory. The radii of most circular trajectories with parallel stacked
particles are larger than 500$\sigma$, which are perceived as linear motion
over our simulation time. In most cases in which vesicles move with higher
angular velocity, we observe a different stacking behavior of the active
particles. Rather than parallel, the rods exhibit angled stacking near the
kinked tip. In circularly moving vesicles with acutely angled kinks
($\theta_{\text{kink}}>\frac{\pi}{2}$) the perpendicular velocity component
stems from an asymmetrically applied pressure at the kink (see sketch in Fig.
1b). As the likelihood of the rods unevenly occupying the space around the
kink increases with $\rho_{N_{a}}$ and decreases with $\theta_{\text{kink}}$,
the angular velocity increases in this region of phase space.

For vesicles with $\theta_{\text{kink}}<\frac{\pi}{2}$ the direction of the
rod alignment changes. Here, the angle is sufficiently narrow that the first
trapped rod is pushed to the boundary by all other incoming particles. As a
consequence, the first rod forms a new, effective ”kink” with the vesicle
boundary where another rod can settle. This mechanism cascades to produce an
array of rods that are all aligned with one side of the vesicle wall (see Fig.
1). Due to the resulting offset in the alignment angle
$\alpha\approx\frac{\theta_{\text{kink}}}{2}$ from the symmetry axis (see
third column in Fig. 4) we can decouple the applied forces at the kink into a
parallel and a perpendicular force component, resulting in a net circular
motion of the vesicle. The long stacking of particles at one side of the kink
is also apparent in the particle distribution in Fig. 3, where the peaks are
more stretched along the vesicle boundary for
$\theta_{\text{kink}}<\frac{\pi}{2}$ than for
$\theta_{\text{kink}}>\frac{\pi}{2}$. Compared to linearly moving vesicle
superstructures, however, both the velocity and angular velocity distribution
are wider and increase with $\rho_{N_{a}}$. The kink-forward and kink-backward
circular regime depends on the existence of a second cluster, which
counteracts the forward pushing component but not the rotational component of
the stacking. Hence, only dimer-filled vesicles can achieve a kink-backward
circular motion dominant regime as shown in Fig. 2.

To relate the stacking of the particles to the angular momentum of the vesicle
superstructure in more detail, we predict its angular velocity by analytically
calculating the rotation of different ideal stacking sequences in the vesicle.
For convenience we focus only on one specific vesicle system in the following
($\theta_{\text{kink}}=\pi/3$, $\rho_{N_{a}}=0.05$), although the same
approach also applies to all other vesicle shapes and rod densities. In this
system with five active rods we can build three different stackings that cause
a rotation (see Fig.5). Each stacking sequence represents the three different
ways of aligning the active rods with both sides of the kink, considering
symmetric equivalence. For each particle $i$ included in the vesicle transfers
its linear momentum to the vesicle such that,
$\sum_{i=1}^{n}|\vec{v}_{i,\bot}|\cdot|l_{i}|=I_{vesicle}\cdot\omega$, where
$l_{i}$ is the position vector of particle $i$ in relation of the center of
mass of the vesicle to the particle $i$, $\vec{v}_{i,\bot}$ is the active
velocity of particle $i$ perpendicular to $l_{i}$ and $I_{vesicle}$ is the
moment of inertia of the vesicle.

For comparison we obtain from the simulation the rotation angle of the vesicle
from its initial orientation at each time step (black bars in Fig. 5) and the
corresponding angular velocities (color bars in Fig. 5). We observe that the
angular velocities of each simulated vesicle fall into one of three groups.
Each group corresponds to one of the ideal stacking sequences, indicated by
the agreement in angular velocity for both dimer and trimer systems. As we
introduce more rod particles into the vesicle, additional types of stacking
sequences can be constructed, resulting in the splitting of the angular
velocities into even more groups and a wider distribution of the angular
velocity as shown in Fig.4.

###  III.3 switchable motion

All of the observed dynamical behaviors of the active vesicle superstructures
emanate from the spontaneous alignment of rods into clusters at the kink in
the vesicle boundary. If this propelling mechanism is to be used to navigate
and control the superstructure, we must be able to switch among the various
behaviors dynamically and reversibly. To test this, we perform three sets of
simulations that represent a horizontal (particle density), vertical (kink
angle) and loop (both density and kink angle in sequence) dynamical change in
regard to the steady-state diagram in Fig. 2.

To alter $\rho_{N_{a}}$ during our simulation we introduce two types of active
rods with similar shape that can be activated independently. By only turning
on one type of trimer confined within a vesicle with
$\theta_{\text{kink}}=\frac{5\pi}{6}$ we first observe the active particles
forming the parallel stacking that induces kink-forward linear motion while
the still passive trimers are dragged as cargo at the opposite end of the
vesicle. Once we activate the second rod type, the newly self-driven trimers
join the existing rod packing at the kink, adding to the propulsion and
accelerating the superstructure. The vesicle returns to its regular speed as
soon as one active rod type is turned off, breaks away from the cluster and
accumulates at the back of the vesicle as cargo again (see video 1). By
repeating this activation sequence with dimer particles we observe that the
vesicles slow down, stop their linear forward motion, or even exhibit backward
motion before they obtain their original speed again (see video 2). This
indicates that the speed of the vesicle is controllable and reversible through
the creation and destruction of counteractive clusters of rods.

Vesicles that continuously change their kink angle from
$\theta_{\text{kink}}=\frac{\pi}{2}$ to $\theta_{\text{kink}}=\frac{5\pi}{6}$
and back, dynamically switch between circular and linear motion (see videos 3
and 4). As predicted by the steady-sate diagram in Fig. 2, the parallel
stacking of active rods in the linear forward motion regime becomes
geometrically unstable or pushed to one side during the kick angle narrowing
process. Hence the dimers and trimer rearrange from aligning parallel to the
vesicle symmetry axis to aligning at an angle, causing subsequent circular
motion of the vesicle. By widening the angle again the particles can slide
from one side of the wall into the center and recover the original parallel
stacking and associated linear motion.

Although the steerability of the vesicle is most reliable in the low density
regime $\rho_{N_{a}}<0.25$, we can combine both the adjustable angle mechanism
and the multiple active particle type mechanism to trigger a variety of
preprogrammed superstructure locomotion. Video 5 shows a vesicle filled with
trimers that first moves linearly (one active/one passive type + obtuse kink
angle), increases its speed (two active types + obtuse kink angle), takes a
turn (two active types + acute kink angle), reduces its velocity to the
original speed (one active/one passive type + acute kink angle) and finally
exhibits a slow linear trajectory again. A similar loop in the phase space of
vesicles filled with dimers results in a sequence of forward circling (one
active/one passive type + acute kink angle), backward circling (two active
types + acute kink angle), backward linear motion (two active types + obtuse
kink angle), forward linear motion (one active/one passive type + obtuse kink
angle) and finally returning to forward circling. The trend in directionality
and spatial distribution of active particles that favors the kinked area of
the boundary weakens as the perimeter of the vesicle increases.

Figure 2:  Chart of directional motion of the vesicle with regards to
$\theta_{kink}$ and $\rho_{N_{a}}$ The right column is the legend for the left
four charts, and describes the statistically dominant motions of the vesicle
superstructures with and without a kink. The markers indicate unbiased motion
($\bigcirc$) and the four newly observed types of motion: linear forward (FW)
motion ($\triangle$), linear backward (BW) motion ($\bigtriangledown$),
circular forward motion (* ()), and circular backward motion (* )). The left
four steady-state charts are diagrams that indicate regions of the different
behaviors. The symbol $\times$ indicates that the majority of vesicles with
active rods have not clustered within the simulation time. Linear and circular
motion are sorted based on the stacking arrangement of the enclosed active
particles across the range of $\theta_{kink}$ and $\rho_{N_{a}}$ studied.
(Note that $\theta_{kink}=\pi$ corresponds to vesicles without a kink). The
color of the left and right halves of each marker indicates the velocity of
the vesicle normalized by the active velocity $v_{a}$ (left color bar), and
the angular velocity of the vesicle (right color bar), respectively. The
number on each marker denotes the percentage of vesicles that move according
to the marker. The remaining percentage of vesicles moves opposite to that
indicated. Vesicles with $\theta_{kink}=\pi$ have numberless markers because
there is no kink axis. Each data point (marker) is averaged over 200
independent replicas. (a) $\alpha=2$, (b) $\alpha=3$ Figure 3:  Spatial
distribution of the self-propelled rods confined in an anisotropic vesicle.
(a) Two-dimensional probability distributions of the rods in a vesicle. The
colors indicate the probability density per bin. (b) The average number of
particle clusters in the vesicle as a function of increasing number density
$\rho_{N_{a}}$ at different $\theta_{kink}$ for each aspect ratio $\alpha=2$
and $\alpha=3$. Each data point is averaged over 200 independent replicas.
Every vesicle in the per-vesicle map (a) and in (b) has bending rigidity
coefficient $\kappa_{a}=1000k_{B}T$. For other $\kappa_{a}$, see SI. Figure 4:
Histogram of the motion descriptors of a vesicle. The histograms describe the
distribution of two motion metrics (centroid velocity and angular velocity)
and alignment angles of clustered particles as a function of $\theta_{kink}$
and $\rho_{N_{a}}$. The set of histograms are plotted for (a) $\alpha=2$,
Pé=75 and (b) $\alpha=3$, Pé=50, respectively with
$\theta_{kink}=\pi/2,5\pi/6$. The first column indicates the velocity of the
vesicle centroid. The second column indicates the angular velocity of the
vesicle. The third column is the Kernel Density Estimation (KDE)[51]
approximation representing the distribution of alignment angles of clustered
particles between the vesicle symmetry axis and the long axis of each active
rod. Each row indicates the number density $\rho_{N_{a}}$ of active rods in
the vesicle. Solid lines at the top of the histogram are KDE approximations of
the corresponding $\theta_{kink}$ in matching colors. Each instance in the
histogram corresponds to a single simulation and the metrics are computed and
collected from the latter half of the simulation trajectory.

(a) $\alpha=2$, $Pe=75$

(b) $\alpha=3$, $Pe=50$

Figure 5:  Analytic calculation of angular velocities from different stacking
arrangements. From the simulation of vesicles with $\theta_{kink}=\pi/3$ and
$\rho_{N_{a}}=0.05$, the rotated angle of the vesicle from its initial
orientation is plotted against time. Each black line is from a single
simulation trajectory. On the right, histogram of angular velocities from each
instance is plotted in black. Sketches on the left are stacking arrangements
of active particles observed at the kinked tip of the vesicle. The analytical
angular velocities are calculated from each simulation replica, grouped by
each stacking sequence, and plotted in the histogram with matching color to
the border of the insets. Finally, the rotation angle plotted against time of
three different stackings of rods is interpolated from the analytical angular
velocity and drawn with matching color and line style to the insets. The
values squared in the histogram are the ratio of the analytical angular
velocity to the simulated angular velocity, averaged over replicas.

##  IV Discussion

In this paper, we used Brownian molecular dynamics simulations to study the
rich dynamical behavior of rigid kinked vesicles that contain self-propelling
rod-shaped particles. We identified that kinks in the vesicle membrane bias
the emergent clustering and alignment of the active agents. Based on the angle
of the kink, density of the active particles, and the active particle length,
the rods form different stacking sequences at or opposite to the kink and
collectively induce multiple types of directed motion by the vesicle
superstructure. Specifically, we identified a kink-forward linear motion, a
kink-backward linear motion, a kink-forward circular motion and a kink-
backward circular motion. By analytically calculating the rotation of
different rod stacking sequences in the vesicle, we showed how rod packing at
and around the kink correlates with the resulting superstructure motion.
Finally, we demonstrated that our model allows for dynamic interchangeability
between the different types of motion and adjustable superstructure speed by
changing the vesicle geometry or active particle properties.

Our findings suggest that both microscopic and macroscopic vesicular
superstructures could be controlled by just two steering parameters, the
geometry of the kinked vesicle and the number of active rods, bringing us
closer to realizing autonomous robotic entities. Motile bacteria or active
colloids encapsulated inside of a giant lipid vesicle are known to impart
motion to the entire vesicle [23]. Including anisotropic elements with
controllable position inside the vesicle membrane [52], or engulfing mixtures
of active particles that respond to different external stimuli [53, 54, 55],
could open the door for artificial cell systems that display directional
movement for targeted delivery or retrieval of micron-scale or smaller
objects. Similarly, in macroscopic machines, such as the superstructure
reported in [21], implementing adjustable bending stiffness to the confining
wall, either mechanically or via actuators [56], could be the key to guide the
location and collective alignment of the individual active machines and to
realize smarter, more adaptable swarm robots.

##  V Acknowledgements

This work was supported as part of the Center for Bio-Inspired Energy Science,
an Energy Frontier Research Center funded by the U.S. Department of Energy,
Office of Science, Basic Energy Sciences under Award # DE-SC0000989. This work
used the Extreme Science and Engineering Discovery Environment (XSEDE) [44],
which is supported by National Science Foundation grant number ACI-1548562
(XSEDE award DMR 140129); computational resources and services were also
provided by Advanced Research Computing (ARC), a division of Information and
Technology Services (ITS) at the University of Michigan, Ann Arbor.

## References

  * Wang  _et al._ [2020] W. Wang, X. Lv, J. L. Moran, S. Duan, and C. Zhou, Soft Matter 16, 3846 (2020).
  * Shaebani  _et al._ [2020] M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper, and H. Rieger, Nature Reviews Physics 2, 181 (2020).
  * Zöttl and Stark [2023] A. Zöttl and H. Stark, Annu. Rev. Condens. Matter Phys. 14, annurev (2023).
  * Gompper [2022] G. Gompper, Nat Mach Intell , 1 (2022), publisher: Nature Publishing Group.
  * Gompper  _et al._ [2020] G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, T. Kiørboe, E. Lauga, W. C. K. Poon, A. DeSimone, S. Muiños-Landin, A. Fischer, N. A. Söker, F. Cichos, R. Kapral, P. Gaspard, M. Ripoll, F. Sagues, A. Doostmohammadi, J. M. Yeomans, I. S. Aranson, C. Bechinger, H. Stark, C. K. Hemelrijk, F. J. Nedelec, T. Sarkar, T. Aryaksama, M. Lacroix, G. Duclos, V. Yashunsky, P. Silberzan, M. Arroyo, and S. Kale, J. Phys. , 68 (2020).
  * Nsamela  _et al._ [2022] A. Nsamela, A. I. Garcia Zintzun, T. D. Montenegro‐Johnson, and J. Simmchen, Small , 2202685 (2022).
  * Arlt  _et al._ [2018] J. Arlt, V. A. Martinez, A. Dawson, T. Pilizota, and W. C. K. Poon, Nat Commun 9, 768 (2018).
  * Merlitz  _et al._ [2020] H. Merlitz, H. D. Vuijk, R. Wittmann, A. Sharma, and J.-U. Sommer, PLoS ONE 15, e0230873 (2020).
  * Fernandez-Rodriguez  _et al._ [2020] M. A. Fernandez-Rodriguez, F. Grillo, L. Alvarez, M. Rathlef, I. Buttinoni, G. Volpe, and L. Isa, Nat Commun 11, 4223 (2020).
  * Singh  _et al._ [2020] D. P. Singh, A. Domínguez, U. Choudhury, S. N. Kottapalli, M. N. Popescu, S. Dietrich, and P. Fischer, Nat Commun 11, 2210 (2020).
  * Lavergne  _et al._ [2019] F. A. Lavergne, H. Wendehenne, T. Bäuerle, and C. Bechinger, Science 364, 70 (2019).
  * Bregulla  _et al._ [2014] A. P. Bregulla, H. Yang, and F. Cichos, ACS Nano 8, 6542 (2014).
  * Yang and Bevan [2020] Y. Yang and M. A. Bevan, Sci. Adv. 6, eaay7679 (2020).
  * Muiños-Landin  _et al._ [2021] S. Muiños-Landin, A. Fischer, V. Holubec, and F. Cichos, Sci. Robot. 6, eabd9285 (2021).
  * Yang  _et al._ [2022] L. Yang, J. Jiang, X. Gao, Q. Wang, Q. Dou, and L. Zhang, Nature Machine Intelligence 4, 480 (2022).
  * Rubenstein  _et al._ [2014] M. Rubenstein, A. Cornejo, and R. Nagpal, Science 345, 795 (2014).
  * Li  _et al._ [2019] S. Li, R. Batra, D. Brown, H.-D. Chang, N. Ranganathan, C. Hoberman, D. Rus, and H. Lipson, Nature 567, 361 (2019).
  * Savoie  _et al._ [2019] W. Savoie, T. A. Berrueta, Z. Jackson, A. Pervan, R. Warkentin, S. Li, T. D. Murphey, K. Wiesenfeld, and D. I. Goldman, Sci. Robot. 4, eaax4316 (2019).
  * Deblais  _et al._ [2018] A. Deblais, T. Barois, T. Guerin, P. Delville, R. Vaudaine, J. Lintuvuori, J. Boudet, J. Baret, and H. Kellay, Phys. Rev. Lett. 120, 188002 (2018).
  * Uplap  _et al._ [2023] S. Uplap, M. F. Hagan, and A. Baskaran, Design principles for transporting vesicles with enclosed active particles (2023), arXiv:2303.16095 [cond-mat].
  * Boudet  _et al._ [2021] J. F. Boudet, J. Lintuvuori, C. Lacouture, T. Barois, A. Deblais, K. Xie, S. Cassagnere, B. Tregon, D. B. Brückner, J. C. Baret, and H. Kellay, Sci. Robot. 6, eabd0272 (2021).
  * Gao and Li [2017] T. Gao and Z. Li, Phys. Rev. Lett. 119, 108002 (2017).
  * Le Nagard  _et al._ [2022] L. Le Nagard, A. T. Brown, A. Dawson, V. A. Martinez, W. C. K. Poon, and M. Staykova, Proc. Natl. Acad. Sci. U.S.A. 119, e2206096119 (2022).
  * Ramos  _et al._ [2020] G. Ramos, M. L. Cordero, and R. Soto, Soft Matter 16, 1359 (2020).
  * Paoluzzi  _et al._ [2016] M. Paoluzzi, R. Di Leonardo, M. C. Marchetti, and L. Angelani, Sci Rep 6, 34146 (2016).
  * Abaurrea-Velasco  _et al._ [2019] C. Abaurrea-Velasco, T. Auth, and G. Gompper, New J. Phys. 21, 123024 (2019).
  * Li and ten Wolde [2019] Y. Li and P. R. ten Wolde, Phys. Rev. Lett. 123, 148003 (2019).
  * Vutukuri  _et al._ [2020] H. R. Vutukuri, M. Hoore, C. Abaurrea-Velasco, L. van Buren, A. Dutto, T. Auth, D. A. Fedosov, G. Gompper, and J. Vermant, Nature 586, 52 (2020).
  * Iyer  _et al._ [2022] P. Iyer, G. Gompper, and D. A. Fedosov, Soft Matter 18, 6868 (2022).
  * Kokot  _et al._ [2022] G. Kokot, H. A. Faizi, G. E. Pradillo, A. Snezhko, and P. M. Vlahovska, Commun Phys 5, 91 (2022).
  * Sanchez  _et al._ [2012] T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Nature 491, 431 (2012).
  * Iyer  _et al._ [2023] P. Iyer, G. Gompper, and D. A. Fedosov, Soft Matter , 10.1039.D3SM00004D (2023).
  * Di Leonardo  _et al._ [2010] R. Di Leonardo, L. Angelani, D. Dell’Arciprete, G. Ruocco, V. Iebba, S. Schippa, M. P. Conte, F. Mecarini, F. De Angelis, and E. Di Fabrizio, Proc. Natl. Acad. Sci. U.S.A. 107, 9541 (2010).
  * Gentile  _et al._ [2020] K. Gentile, A. Somasundar, A. Bhide, and A. Sen, Chem 6, 2174 (2020).
  * Wensink and Löwen [2008] H. H. Wensink and H. Löwen, Phys. Rev. E 78, 031409 (2008).
  * Caprini  _et al._ [2020] L. Caprini, F. Cecconi, C. Maggi, and U. Marini Bettolo Marconi, Phys. Rev. Research 2, 043359 (2020).
  * Al Alam  _et al._ [2022] E. Al Alam, M. Brun-Cosme-Bruny, V. Borne, S. Faure, B. Maury, P. Peyla, and S. Rafaï, Phys. Rev. Fluids 7, L092301 (2022).
  * Frangipane  _et al._ [2019] G. Frangipane, G. Vizsnyiczai, C. Maggi, R. Savo, A. Sciortino, S. Gigan, and R. Di Leonardo, Nat Commun 10, 2442 (2019).
  * Kaiser  _et al._ [2012] A. Kaiser, H. H. Wensink, and H. Löwen, Phys. Rev. Lett. 108, 268307 (2012).
  * Reichhardt and Reichhardt [2017] C. O. Reichhardt and C. Reichhardt, Annu. Rev. Condens. Matter Phys. 8, 51 (2017).
  * Glaser  _et al._ [2020] J. Glaser, X. Zha, J. A. Anderson, S. C. Glotzer, and A. Travesset, Computational Materials Science 173, 109430 (2020).
  * Nguyen  _et al._ [2011] T. D. Nguyen, C. L. Phillips, J. A. Anderson, and S. C. Glotzer, Computer Physics Communications 182, 2307 (2011).
  * Warner [1972] H. R. Warner, Ind. Eng. Chem. Fund. 11, 379 (1972).
  * Towns  _et al._ [2014] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr, Comput. Sci. Eng. 16, 62 (2014).
  * Anderson  _et al._ [2020] J. A. Anderson, J. Glaser, and S. C. Glotzer, Computational Materials Science 173, 109363 (2020).
  * Ramasubramani  _et al._ [2020] V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson, and S. C. Glotzer, Computer Physics Communications 254, 107275 (2020).
  * Adorf  _et al._ [2018] C. S. Adorf, P. M. Dodd, V. Ramasubramani, and S. C. Glotzer, Computational Materials Science 146, 220 (2018).
  * Dice  _et al._ [2021] B. Dice, B. Butler, V. Ramasubramani, A. Travitz, M. Henry, H. Ojha, K. Wang, C. Adorf, E. Jankowski, and S. Glotzer (Austin, Texas, 2021) pp. 23–32.
  * Giomi  _et al._ [2013] L. Giomi, N. Hawley-Weld, and L. Mahadevan, Proc. R. Soc. A. 469, 20120637 (2013).
  * Nikola  _et al._ [2016] N. Nikola, A. P. Solon, Y. Kafri, M. Kardar, J. Tailleur, and R. Voituriez, Phys. Rev. Lett. 117, 098001 (2016).
  * noa [1992] in  _Multivariate Density Estimation_ (John Wiley & Sons, Ltd, 1992) pp. 125–193, section: 6 _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470316849.ch6.
  * Xin  _et al._ [2021] W. Xin, H. Wu, G. M. Grason, and M. M. Santore, Sci. Adv. 7, eabf1943 (2021).
  * Niese  _et al._ [2020] L. Niese, L. Wang, S. Das, and J. Simmchen, Soft Matter 16, 10585 (2020).
  * Salinas  _et al._ [2021] G. Salinas, K. Tieriekhov, P. Garrigue, N. Sojic, L. Bouffier, and A. Kuhn, J. Am. Chem. Soc. 143, 12708 (2021).
  * Ziepke  _et al._ [2022] A. Ziepke, I. Maryshev, I. S. Aranson, and E. Frey, Nat Commun 13, 6727 (2022).
  * Guo  _et al._ [2021] X. Guo, W. Li, and W. Zhang, Sci. China Technol. Sci. 64, 1663 (2021).

##  VI Supplementary

### Two dimensional probability density map of finding a rod in the vesicle

We use kernel density estimation (KDE) with a Gaussian kernel to estimate the
probability distribution of rods’ relative location in the vesicle in the
steady state. The probability to find a rod is measured per rectangular bin
located at $(w_{i},h_{j}),$ where $i,j\in[1,100]$. The bins are equally spaced
in the rectangular region that spans between the lower and upper limit of the
width of the vesicle along the $w$ direction and of the height of the vesicle
along the $h$ direction.

### Shape of the vesicle with respect to bending stiffness

The vesicle shape in steady state varies with vesicle stiffness. The vesicle
stiffness is set by the harmonic angle ”spring constant”, $\kappa_{a}$,
applied to every neighboring triplet of disks that composes the vesicle. The
shape elongation is determined by eigenvalues and eigenvectors of the vesicle
boundary gyration tensor. In two dimensions, the x-y component of the gyration
tensor is

$s_{xy}=\frac{1}{2N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}(x_{i}-x_{j})(y_{i}-y_{j})$

The square roots of the eigenvalues ($\sqrt{L1},\sqrt{L2}$) are the
characteristic principal-axis lengths (radii) of the ellipsoid that
describesthe shape of the vesicle. The metric for the shape descriptor is
defined as a ratio of the shortest gyration moment to the longest gyration
moment: $Ratio(L1,L2)=\frac{\sqrt{L1}}{\sqrt{L2}}$. The ratio is plotted in
Fig. S1 for each rod aspect ratio ($\alpha=2,3$) with various
$\kappa_{a}=50,125,250,500,1000k_{b}T$, $\theta_{kink}=\pi/2,5\pi/6$, and
number densities $\rho_{N_{a}}=0.1,0.15,0.2$. The value of the shape
descriptor decreases with decreasing vesicle stiffness, which means that the
vesicle is more elongated with more flexible boundaries. The reference points
marked with black squares are from the vesicles containing no propelling rods.
The shape descriptor deviates more from the reference point as vesicle
stiffness decreases. The higher number density $\rho_{N_{a}}$ leads to more
elongation.

Figure S1: Shape elongation of the vesicle. Shape descriptor measure of the
extent of elongation for different values of vesicle bending stiffness
$\kappa_{a}$ and different aspect ratios $\alpha$ of the rod. Color refers to
the number densities $\rho_{N_{a}}=0.1,0.15,0.2$. The marker is the averaged
descriptor value with error bar of matching color. The reference system for
each value of bending stiffness is from a vesicle with no active particles.
Figure S2: Ratio of parallel vs. angled stacking plotted over vesicles with
$\alpha=2$ for each data point with corresponding Pe, $\theta_{kink}$ and
$\rho_{N_{a}}$. For each moving vesicle the stacking was measured either
parallel or angled. Angled stacking corresponds to cases when the closest
active rod particle to the kink aligns with the vesicle boundary; otherwise
the rods are considered to be stacked in parallel. Figure S3: Ratio of
parallel vs. angled stacking plotted for vesicles with $\alpha=3$ for each
data point with corresponding Pe, $\theta_{kink}$ and $\rho_{N_{a}}$. The
measurement method is identical to that in Fig. S2.  Figure S4: Chart of
directional motion of the vesicle with regards to Pe=75, $\theta_{kink}$ and
$\rho_{N_{a}}$. The color scheme is the same as in Fig. 2 Figure S5:
Probability density plot of enclosed particle $\alpha=2$, $Pe=25$,
$\kappa_{a}=1000k_{B}T$. The color scheme is the same as in Fig. 3 Figure S6:
Probability density plot of enclosed particle $\alpha=2$, $Pe=50$,
$\kappa_{a}=1000k_{B}T$. The color scheme is the same as in Fig. 3 Figure S7:
Probability density plot of enclosed particle $\alpha=2$, $Pe=75$,
$\kappa_{a}=1000k_{B}T$. The color scheme is the same as in Fig. 3 Figure S8:
Probability density plot of enclosed particle $\alpha=3$, $Pe=25$,
$\kappa_{a}=1000k_{B}T$. The color scheme is the same as in Fig. 3 Figure S9:
Probability density plot of enclosed particle $\alpha=3$, $Pe=50$,
$\kappa_{a}=1000k_{B}T$. The color scheme is the same as in Fig. 3 Figure S10:
Probability density plot of enclosed particle $\alpha=3$, $Pe=75$,
$\kappa_{a}=1000k_{B}T$. The color scheme is the same as in Fig. 3 Figure S11:
Flipped trend in probability density plot of enclosed particle with various
$\theta_{kink}$ $\alpha=3$, $Pe=25,50$, $\kappa_{a}=1000k_{B}T$. The color
scheme is the same as in Fig. 3