††thanks: V.P.P. and H.T. contributed equally to this work.††thanks: V.P.P.
and H.T. contributed equally to this work.††thanks: Correspondence to
<EMAIL_ADDRESS><EMAIL_ADDRESS>Correspondence to
<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Supplementary Materials:
Ultrafast reversible self-assembly of living tangled matter

Vishal P. Patil  School of Humanities and Sciences, Stanford University, 450
Serra Mall, Stanford, CA 94305  Harry Tuazon  School of Chemical and
Biomolecular Engineering, Georgia Institute of Technology, Atlanta, GA 30318
Emily Kaufman  School of Chemical and Biomolecular Engineering, Georgia
Institute of Technology, Atlanta, GA 30318  Tuhin Chakrabortty  School of
Chemical and Biomolecular Engineering, Georgia Institute of Technology,
Atlanta, GA 30318  David Qin  Wallace H. Coulter Department of Biomedical
Engineering, Georgia Institute of Technology, Atlanta, GA 30332  Jörn Dunkel
Department of Mathematics, Massachusetts Institute of Technology, 77
Massachusetts Avenue, Cambridge, MA 02139  M. Saad Bhamla  School of Chemical
and Biomolecular Engineering, Georgia Institute of Technology, Atlanta, GA
30318

## Experimental methods

###  .1 Animals

We procured and reared blackworms as described in Tuazon, et al. [1]. We
placed another culture of worms in a 10-gallon aquarium filled with pool
filtered sand. Water is continuously filtered using an aquarium filter and
cooled to 12∘C using a water chiller. Studies with blackworms do not require
approval by an institutional animal care committee.

###  .2 Worm tangling and untangling experiments

We used a Leica MZ APO microscope (Heerbrugg, Switzerland) with an ImageSource
DFK 33UX264 camera (Charlotte, NC) at 30FPS to record all worm tangling
experiments placed in a 15mm confocal glass bottom petridish. To encourage
full entanglement, we placed worms that were kept in the aquarium into room
temperature water, which varied from 21.0±0.3∘C. A Chronos color (Burnaby,
British Columbia, Canada) high speed camera captured all disentangling
experiments at 1067FPS in a 35mm petridish. To stimulate rapid
disentanglement, we applied a quick 1-2s shock using a 9V battery. This method
does not cause harm to the worms. The worms were chosen from a population with
mean length $30\,$mm (2.s.f) and mean radius $0.3\,$mm (1.s.f). Untangling
data were collected from 5 trials each containing 5 worms. Every worm was
tracked yielding 25 trajectories. Tangling data were collected from 4 trials
each containing 5 worms. In total, 18 of the trajectories were tracked.

###  .3 Worm tracking

Recordings were converted into AVI format using Adobe© Premiere Pro and were
then imported into ImageJ [2] for image analysis. After enhancing the
contrast, we applied background subtraction and converted stacks to 8-bit. We
manually tracked the prostomium of the worms using TrackMate [3] for the
entirety of each stack or when we can no longer see the anterior segments.
This occurs when worms are entangled or escape outside the camera’s field of
view.

### Worm Ultrasound methods

Worms were fixed in a 16$\%$ gelatin solution using filtered water and Type A
gelatin (MP Biomedicals, Solon, OH). The solution was allowed to cool to 50∘C
before being placed in a desiccator vacuum chamber for 5 minutes to remove any
bubbles. The solution was then placed directly onto worms in a petridish and
was allowed to solidify for at least 1 hour. After immobilization, the sample
was transferred onto a flat gelatin disk; we then imaged them using a Vevo
2100 ultrasound system (VisualSonics, Toronto, Canada) and MS250 transducer
(21 MHz center frequency), acquiring 2D grayscale ultrasound images at a 30Hz
frame rate. The transducer was translated in the elevational direction to
yield a 3D image of the worm tangle (Fig. 1B). Contrast was adjusted within
the ultrasound system; the acquired data were then processed in MATLAB 2021b
(Mathworks, Natick, MA). The 3D images underwent denoising and manual
segmentation on MATLAB’s in-built GUI ‘Volume Segmenter’ to yield the final
worm tangle representation (Fig. 1C). The worm centerline curves were
extracted using the diffusion maps algorithm (Fig. S1A,B).

Figure S1: Segmented and smoothed ultrasound data (A) Segmented ultrasound
data for 3 different experiments with 12 worms each show tightly packed worms
exhibiting thigmotactic behavior. Markers correspond to markers in Fig. 1H,I.
Data set 3 (yellow diamond) is used as an initial condition for an untangling
simulation (Fig. 3E and movie S2). Scale bar 5mm. (B) Using the diffusion maps
algorithm, curves are fit to the point clouds representing each worm in (A).

## Mathematical modelling

### Geometry and topology of curves

In this section, we introduce the key geometrical and topological properties
of curves which are discussed in the main text, including link and contact
link.

#### Chirality

To define chirality, we introduce the Frenet frame. Let $\gamma(s)$ be a curve
parametrized by arc length with $\mathbf{t}=\gamma^{\prime}$, so
$|\mathbf{t}|=1$. The Frenet frame is an orthonormal frame,
$\\{\mathbf{t},\mathbf{n},\mathbf{b}\\}$, defined by

$\displaystyle\boldsymbol{t}^{\prime}(s)=\kappa\mathbf{n}\qquad\boldsymbol{n}^{\prime}(s)=-\kappa\mathbf{t}+\tau\mathbf{b}\qquad\boldsymbol{b}^{\prime}(s)=-\tau\mathbf{n}$

where primes denote $s$ derivatives. The quantities $\kappa(s)$ and $\tau(s)$
are the geometrical curvature and the torsion of the curve, respectively.
$\kappa$ is related to the bending strain of an elastic fiber. In particular,
if a cylindrical elastic fiber with radius $h$ has a centerline with curvature
$\kappa(s)$, then the bending strain is $\kappa h$ (Fig. 1F).

We define the chirality, $\chi$, of the 3D curve by

$\displaystyle\chi=\boldsymbol{t}\cdot(\boldsymbol{t}^{\prime}\times\boldsymbol{t}^{\prime\prime})$
$\displaystyle=\boldsymbol{t}\cdot\left(\kappa\boldsymbol{n}\times(\kappa^{\prime}\boldsymbol{n}+\kappa\boldsymbol{n}^{\prime})\right)$
$\displaystyle=\boldsymbol{t}\cdot\left(\kappa\boldsymbol{n}\times(\kappa^{\prime}\boldsymbol{n}-\kappa^{2}\boldsymbol{t}+\kappa\tau\boldsymbol{b})\right)$
$\displaystyle=\kappa^{2}\tau$

This can be written in terms of the tangent vectors close to a given point

$\displaystyle\boldsymbol{t}(s)\cdot\left(\boldsymbol{t}(s+h)\times\boldsymbol{t}(s-k)\right)$
$\displaystyle=\boldsymbol{t}\cdot\left((\boldsymbol{t}+h\boldsymbol{t}^{\prime}+\frac{1}{2}h^{2}\boldsymbol{t}^{\prime\prime})\times(\boldsymbol{t}-k\boldsymbol{t}^{\prime}+\frac{1}{2}k^{2}\boldsymbol{t}^{\prime\prime})\right)$
$\displaystyle=\boldsymbol{t}\cdot(\boldsymbol{t}^{\prime}\times\boldsymbol{t}^{\prime\prime})hk(h+k)/2$

In 2D, the signed curvature plays a similar role to chirality

$\displaystyle\chi_{2}=\mathbf{t}\wedge\mathbf{t}^{\prime}=\kappa$

#### Link

Figure S2: Link and contact link (A) The link of 2 closed curves is given by
the signed crossing number of a planar projection. (B,C) The link of open
curves measures topological obstruction. In (B), the green curve winds around
the purple curve, and the resulting configuration has nonzero link. In (C),
the green curve lies on top of the purple curve and the configuration has zero
link.

Here we define the linking number of two closed curves both diagrammatically
and in terms of the Gauss linking integral. We then show how these definitions
can be extended to give linking numbers for open curves, as in Ref. [4].

Link is a property of two closed curves. Consider two such closed curves
$\gamma_{1}$, $\gamma_{2}$

$\displaystyle\gamma_{i}:[0,L_{i}]\rightarrow\mathbb{R}^{3}$

The linking number can be defined in terms of the number of crossings, counted
with sign, between $\gamma_{1}$ and $\gamma_{2}$ in a projection of the curves
onto a plane (Fig. S2A). Consider a projection, $P_{\boldsymbol{n}}$ in the
direction $\boldsymbol{n}$

$\displaystyle P_{\boldsymbol{n}}=I-\boldsymbol{n}\boldsymbol{n}^{\top}$

where $\boldsymbol{n}$ is a unit vector. The resulting planar diagram will
have some number of crossings between the curves $\gamma_{1}$ and
$\gamma_{2}$. Suppose there are $k$ crossings, at points
$\gamma_{1}(s_{i}),\gamma_{2}(\sigma_{i})$, where $i=1,2,...,k$. These
crossing points satsify

$\displaystyle\Big{(}\gamma_{1}(s_{i})-\gamma_{2}(\sigma_{i})\Big{)}\parallel\boldsymbol{n}$

The sign of each crossing is

$\displaystyle\text{sgn}(i)=\text{sgn}\Big{(}\left[\gamma_{1}(s_{i})-\gamma_{2}(\sigma_{i})\right]\cdot\left[\gamma_{1}^{\prime}(s_{i})\times\gamma_{2}^{\prime}(\sigma_{i})\right]\Big{)}$

The total signed crossing number of the planar diagram obtained from the
projection $P_{\boldsymbol{n}}$ is then

$\displaystyle c_{\boldsymbol{n}}=\sum_{i}\text{sgn}(i)$

For closed curves, this signed crossing number is independent of the
projection vector $\boldsymbol{n}$, and so is equal to twice the linking
number

$\displaystyle Lk(\gamma_{1},\gamma_{2})=\frac{1}{2}c_{\boldsymbol{n}}$

The linking number can also be written as an integral by considering the
following smooth map

$\displaystyle\Gamma:[0,L_{1}]\times[0,L_{2}]$ $\displaystyle\rightarrow
S^{2}$ $\displaystyle(s,\sigma)$
$\displaystyle\mapsto\frac{\gamma_{1}(s)-\gamma_{2}(\sigma)}{|\gamma_{1}(s)-\gamma_{2}(\sigma)|}$

Link is the degree of the map $\Gamma$

$\displaystyle Lk(\gamma_{1},\gamma_{2})=\deg\Gamma$
$\displaystyle=\frac{1}{4\pi}\int
ds\,d\sigma\,\Gamma\cdot\left(\Gamma_{s}\times\Gamma_{\sigma}\right)$
$\displaystyle=\frac{1}{4\pi}\int
ds\,d\sigma\,\frac{(\gamma_{1}(s)-\gamma_{2}(\sigma))\cdot(\gamma_{1}^{\prime}(s)\times\gamma_{2}^{\prime}(\sigma))}{|\gamma_{1}(s)-\gamma_{2}(\sigma)|^{3}}$

Therefore, for closed curves, we have

$\displaystyle Lk(\gamma_{1},\gamma_{2})=\frac{1}{4\pi}\int
ds\,d\sigma\,\frac{(\gamma_{1}(s)-\gamma_{2}(\sigma))\cdot(\gamma_{1}^{\prime}(s)\times\gamma_{2}^{\prime}(\sigma))}{|\gamma_{1}(s)-\gamma_{2}(\sigma)|^{3}}=\frac{1}{2}c_{\boldsymbol{n}}$

where $c_{\boldsymbol{n}}$ is the signed crossing number defined above.

For open curves, we can still define link by the above integral, however its
interpretation in terms of signed crossing numbers is different. In
particular, the signed crossing number for open curves depends on the choice
of projection. As in Ref [4], the corresponding result is now

$\displaystyle Lk(\gamma_{1},\gamma_{2})=\frac{1}{4\pi}\int
ds\,d\sigma\,\frac{(\gamma_{1}(s)-\gamma_{2}(\sigma))\cdot(\gamma_{1}^{\prime}(s)\times\gamma_{2}^{\prime}(\sigma))}{|\gamma_{1}(s)-\gamma_{2}(\sigma)|^{3}}=\frac{1}{8\pi}\int
d^{2}\boldsymbol{n}\,c_{\boldsymbol{n}}$

The linking integral therefore gives the average signed crossing number over
all projections. This is a measure of the topological obstruction one curve
imposes on the motion of the other (Fig. S2B,C). For example, using
cylindrical polar coordinates, the configuration of Fig. S2B, in which one
curve winds around another, has

$\displaystyle(\gamma_{1}(s)-\gamma_{2}(\sigma))\cdot(\gamma_{1}^{\prime}(s)\times\gamma_{2}^{\prime}(\sigma))=-\mathbf{e}_{r}(\sigma)\cdot\left(\mathbf{e}_{z}\times\mathbf{e}_{\theta}(\sigma)\right)=1$

so $Lk\neq 0$. In Fig. S2C, the two curves can be easily separated from each
other

$\displaystyle(\gamma_{1}(s)-\gamma_{2}(\sigma))\cdot(\gamma_{1}^{\prime}(s)\times\gamma_{2}^{\prime}(\sigma))=\cos\sigma\,\mathbf{e}_{x}\cdot\left(\mathbf{e}_{z}\times\sin\sigma\,\mathbf{e}_{x}\right)=0$

and thus $Lk=0$.

#### Contact link

The contact link defined in the main text is a property of two curves
$\gamma_{1}$ and $\gamma_{2}$, and a thickness parameter, $r$

$\displaystyle
cLk(\gamma_{i},\gamma_{j}\,;r)=\begin{cases}|Lk(\gamma_{i},\gamma_{j})|\quad&\text{if
}\min_{s,t}|\gamma_{1}(s)-\gamma_{2}(t)|<2r\\\
0\quad&\text{otherwise}\end{cases}$ (1)

In other words, $cLk=|Lk|$ if tubes with radius $r$, centered around
$\gamma_{1}$ and $\gamma_{2}$ are in contact, and $cLk=0$ otherwise. For
worms, the value of $r$ is naturally given by the worm radius, so for
notational convenience, we often drop the explicit $r$-dependence of $cLk$.
However the change in contact link with $r$ is reminiscent of persistence
homology and so can also be a useful tangle quantity (Fig. 2F).

Given $N$ tubes with fixed $r$, and centerline curves,
$\gamma_{1},\gamma_{2},...\gamma_{N}$, we can construct the contact link
matrix, $C_{ij}$

$\displaystyle C_{ij}=cLk(\gamma_{i},\gamma_{j})$

The adjacency matrix $A_{ij}$ of the tangle graphs constructed in the main
text follows from thresholding the contact link matrix

$\displaystyle A_{ij}=\begin{cases}1\quad\text{if }C_{ij}\geq 1/2\\\
0\quad\text{if }C_{ij}<1/2\end{cases}$

We can obtain measures of the complexity of the tangle from the matrices
$A_{ij}$ and $C_{ij}$. In particular, we define the total contact link per
worm $T_{c}$, and the mean degree $\langle d\rangle$ as follows

$\displaystyle T_{c}(\gamma_{1},...,\gamma_{N})$
$\displaystyle=\frac{1}{N}\sum_{i,j}C_{ij}=\left\langle\sum_{j}C_{ij}\right\rangle$
(2) $\displaystyle\langle d\rangle$
$\displaystyle=\frac{1}{N}\sum_{i,j}A_{ij}$ (3)

$T_{c}$ is therefore a measure of the average amount of tangling between one
worm and the rest of the tangle. The mean degree can be thought of as the
thresholded analogue of the total contact link per worm. The radius dependence
of $T_{c}$ is shown in Fig. 2F. Since $T_{c}$ depends on the sharp contact
threshold in (1), we smooth our plots of $T_{c}$ against time (Fig. 3F and
Fig. 4E).

For larger values of $T_{c}$ and $\langle d\rangle$, the tangle is more
complex, and the tangle graph has more edges. When $T_{c}$ and $\langle
d\rangle$ are small, the $N$ curves do not form a coherent tangle. and the
tangle graph will be sparse and disconnected. We can therefore use $T_{c}$ and
$\langle d\rangle$ to map tangling to a percolation problem. Since a connected
graph on $n$ vertices has mean degree at least $2-1/n$, we expect $T_{c}>2$ to
correspond to tangled states, and $T_{c}<2$ to correspond to untangled states,
in the large $n$ limit. From our ultrasound reconstructions, we find (Fig.
2E,F)

$\displaystyle T_{c}$ $\displaystyle=1.8,\;2.1,\;1.7$ (4a)
$\displaystyle\langle d\rangle$ $\displaystyle=1.2,\;1.7,\;1.2$ (4b)

where the values correspond to the data sets in Fig. S1, from left to right.
We note that the low values of $\langle d\rangle$ are due to the fact that our
ultrasound worm tangles do not have connected tangle graphs (Fig. 2D,E).

### Elastic model

We model each worm using the Kirchhoff equations [5] together with internal
damping and a hydrodynamic friction [6, 7]. The worm centerline is given by
$\boldsymbol{x}(s,t)$ where $s\in[0,L]$ is an arc length parameter for the
unstretched fiber and $L$ is the worm length. The worm has material frame
$[\mathbf{d}_{1},\mathbf{d}_{2},\mathbf{d}_{3}]$ where
$\boldsymbol{x}^{\prime}\parallel\mathbf{d}_{3}$ and primes denote
$s$-derivtatives. The twist density is given by
$\theta^{\prime}=\mathbf{d}_{1}^{\prime}\cdot\mathbf{d}_{2}$. Similarly, we
define the angular velocity component
$\omega_{3}=\dot{\mathbf{d}}_{1}\cdot\mathbf{d}_{2}$, which is related to
$\theta^{\prime}$ by the identity

$\displaystyle\omega_{3}^{\prime}+\mathbf{d}_{3}\cdot\left(\mathbf{d}_{3}^{\prime}\times\dot{\mathbf{d}}_{3}\right)=\dot{\theta}^{\prime}$

The Kirchhoff equations, augmented by terms for damping, friction, contact and
activity, are

$\displaystyle\boldsymbol{f}^{\text{elast}}+\boldsymbol{f}^{\text{con}}+\boldsymbol{f}^{\text{fric}}+\boldsymbol{f}^{\text{act}}$
$\displaystyle=\rho A\dot{\boldsymbol{v}}+\gamma A\boldsymbol{v}-\eta
A\boldsymbol{v}^{\prime\prime}$ (5a)
$\displaystyle\tau^{\text{elast}}+\tau^{\text{fric}}$ $\displaystyle=\rho
I_{3}\dot{\omega}_{3}+\gamma I_{3}\omega_{3}-\eta
I_{3}(\omega_{3}\mathbf{t})^{\prime\prime}\cdot\mathbf{t}$ (5b)

where $\boldsymbol{v}=\dot{\boldsymbol{x}}$, $\rho$ is the density, $A$ is the
cross sectional area, $I_{3}$ is the perpendicular cross sectional moment of
inertia and $\gamma,\eta$ are damping coefficients, with dimensions

$\displaystyle[\gamma]=ML^{-3}T^{-1},\qquad\quad[\eta]=ML^{-1}T^{-1}$

We assume the fiber has circular cross section with radius $h$, so $A=\pi
h^{2}$ and $I_{3}=\pi h^{4}/2$. Here we focus on the case where the activity
is localized at the head of the worm, so $\boldsymbol{f}^{\text{act}}=0$, but
there is an active force at the boundary, $\mathbf{F}^{\text{act}}(L,t)$. To
write the boundary conditions, let $\mathbf{F}^{\text{tot}},T^{\text{tot}}$ be
the total force and torque on the rod, including the viscous dissipation terms

$\displaystyle\left(\mathbf{F}^{\text{tot}}\right)^{\prime}$
$\displaystyle=\boldsymbol{f}^{\text{elast}}+\boldsymbol{f}^{\text{con}}+\boldsymbol{f}^{\text{fric}}+\eta
A\boldsymbol{v}^{\prime\prime}$
$\displaystyle\left(T^{\text{tot}}\right)^{\prime}$
$\displaystyle=\tau^{\text{elast}}+\tau^{\text{fric}}+\eta
I_{3}(\omega_{3}\mathbf{t})^{\prime\prime}\cdot\mathbf{t}$

Here $\mathbf{F}^{\text{tot}},T^{\text{tot}}$ have units of force and torque,
whereas the terms $\boldsymbol{f}^{\text{elast}},\tau^{\text{elast}}$, etc.
have units of force density and torque density. Then the boundary conditions
describing a worm driven by an active force at the head are

$\displaystyle\mathbf{F}^{\text{tot}}(0,t)$
$\displaystyle=0,\qquad\quad\mathbf{F}^{\text{tot}}(L,t)=\mathbf{F}^{\text{act}}(t)$
$\displaystyle T^{\text{tot}}(0,t)$ $\displaystyle=0,\qquad\quad
T^{\text{tot}}(L,t)=0$

The forces in torques in equation (5) consist of an elastic force and torque
arising from the elastic energy of the fiber, a contact force,
$\boldsymbol{f}^{\text{con}}$, which prevents fibers from intersecting, a
frictional force $\boldsymbol{f}^{\text{fric}}$ and torque
$\tau^{\text{fric}}$. Following Ref. [5], our simulation framework stores a
discretized $\boldsymbol{x}$ and $\theta$, and the elastic forces and torques
are calculated by taking derivatives of a discretized elastic energy. The
discretized contact and friction forces are described further in Ref. [6] and
follow continuum expressions described below.

We can describe the elastic forces and torques on the left hand side of
equation (5) in more detail. The elastic forces and torques are derivatives of
the the elastic energy of the fiber, which is made up of the bending energy
$\mathcal{E}_{b}$, twisting energy, $\mathcal{E}_{\text{tw}}$ and stretching
energy $\mathcal{E}_{s}$

$\displaystyle\mathcal{E}_{b}=\frac{1}{2}E_{b}I\int_{0}^{L}ds\,\kappa^{2},\qquad\quad\mathcal{E}_{tw}=\frac{1}{2}\mu_{b}J\int_{0}^{L}ds\,\theta^{\prime
2},\qquad\quad\mathcal{E}_{s}=\frac{1}{2}EA\int_{0}^{L}ds\,\left(|\boldsymbol{x}^{\prime}|-1\right)^{2}$

where $I$ is the cross sectional moment of inertia, $J$ is the moment of
twist, $E$ is the Young’s modulus of the fiber, $E_{b}$ is the bending
modulus, and $\mu_{b}$ is the shear modulus. Under the assumption of circular
cross sections, $I=\pi h^{4}/4$ and $J=\pi h^{4}/2$. Owing to the complex
internal structure of worms, we take $E_{b}\neq E$. The Poisson’s ratio,
$\nu$, relates the shear modulus to the bending modulus,
$\mu_{b}=E_{b}/(2+2\nu)$. The elastic forces and torques are then

$\displaystyle\boldsymbol{f}^{\text{elast}}(s)$
$\displaystyle=-\frac{\delta\mathcal{E}_{b}}{\delta\boldsymbol{x}(s)}-\frac{\delta\mathcal{E}_{tw}}{\delta\boldsymbol{x}(s)}-\frac{\delta\mathcal{E}_{s}}{\delta\boldsymbol{x}(s)}$
$\displaystyle\tau^{\text{elast}}(s)$
$\displaystyle=-\frac{\delta\mathcal{E}_{tw}}{\delta\theta(s)}$

The remaining forces describe self-interactions within the fiber. There are
many different approaches to discretizing and implementing these interactions
[6, 8, 5, 9]. We begin with the contact forces, which prevent the fiber from
intersecting itself. Let
$\boldsymbol{r}(s,\sigma)=\boldsymbol{x}(s)-\boldsymbol{x}(\sigma)$ and
$\hat{\boldsymbol{r}}=\boldsymbol{r}/|\boldsymbol{r}|$. To quantify contact,
set

$\displaystyle p(s,\sigma)=1-\frac{\boldsymbol{r}(s,\sigma)}{2h_{\text{eff}}}$

where $h_{\text{eff}}-h$ is small. If $h_{\text{eff}}=h$, then $p=0$ when
$\boldsymbol{x}(s),\boldsymbol{x}(\sigma)$ are touching and $p=1$ when
$\boldsymbol{x}(s),\boldsymbol{x}(\sigma)$ overlap completely. Choosing
$h_{\text{eff}}\neq h$ has the desirable effect of smoothing out the contact
region and the associated contact forces [8]. When $h_{\text{eff}}>h$, we use
the threshold $h_{\text{eff}}$ to calculate the the contact link (1) from
simulation data. Using $p$ we can construct a contact potential energy, $V(p)$

$\displaystyle V(p)$
$\displaystyle=\frac{K}{2}p^{2}+Kp_{0}\left(\frac{1}{4}p^{4}+\frac{1}{6}p^{6}+\frac{1}{8}p^{8}\right)\quad\text{for
}p>0$ $\displaystyle V(p)$ $\displaystyle=0\quad\text{for }p<0$

where $K$ is an effective bulk modulus, and $p_{0}$ is a constant that sets
the rate at which the potential stiffens. Here, we take $p_{0}=100$. The
contact force is then

$\displaystyle\boldsymbol{f}^{\text{con}}(s)=-\int_{|\sigma-s|>2h_{\text{eff}}}d\sigma\,\boldsymbol{x}^{\prime}(s)\cdot\boldsymbol{x}^{\prime}(\sigma)\,\frac{dV}{dp}\bigg{|}_{p(s,\sigma)}\,\hat{\boldsymbol{r}}(s,\sigma)$

The integral excludes points $\boldsymbol{x}(\sigma)$ which are close to
$\boldsymbol{x}(s)$ along the fiber, i.e. points with $|\sigma-s|\leq
2h_{\text{eff}}$. This ensures that only points
$\boldsymbol{x}(s),\boldsymbol{x}(\sigma)$ which are legitimately in contact
are included in the contact force integral.

We implement the contact law by discretizing the rod [5, 6, 8] into vertices,
$\boldsymbol{x}_{i}=\boldsymbol{x}(s_{i})$, connected by links with midpoints
given by
$\bar{\boldsymbol{x}}_{i}=(\boldsymbol{x}_{i}+\boldsymbol{x}_{i+1})/2$. An
appropriate choice of $h_{\text{eff}}$ depends on the discretization length of
the rod, $|\bar{\boldsymbol{x}}_{i}|$. Contact forces as described above are
then implemented for link-link interactions and vertex-vertex interactions.
The length scale $h_{\text{eff}}$ is chosen for each interaction type so that
the average value of $h_{\text{eff}}$ is $h$. In our $30$-worm, $19$-worm
(Fig. 4D, movie S3) and $8$-worm (Fig. 3E movie S2) simulations, we have
$|\bar{\boldsymbol{x}}_{i}|=2h$ and we take $h_{\text{eff}}=1.25h$ for link-
link interactions, and $h_{\text{eff}}=0.75$ for vertex-vertex interactions.

Finally, we include hydrodynamic-type friction terms [6]

$\displaystyle\boldsymbol{f}^{\text{fric}}(s)$ $\displaystyle=-\zeta
A\int_{|\sigma-s|>2h}d\sigma\,\Theta\left[p(s,\sigma)\right]\frac{\dot{\boldsymbol{x}}(s)-\dot{\boldsymbol{x}}(\sigma)}{|\boldsymbol{r}(s,\sigma)|^{c}}$
$\displaystyle\tau^{\text{fric}}(s)$
$\displaystyle=\zeta_{tw}J\int_{|\sigma-s|>2h}d\sigma\,\Theta\left[p(s,\sigma)\right]\frac{\boldsymbol{r}(s,\sigma)\times\left(\dot{\boldsymbol{x}}(s)-\dot{\boldsymbol{x}}(\sigma)\right)\cdot\mathbf{d}_{3}(\sigma)}{|\boldsymbol{r}(s,\sigma)|^{c_{tw}}}$

Here, following Ref. [6] we take $c=2$ and $c_{tw}=4$, so the coefficients
$\zeta,\zeta_{tw}$ have the same dimensions

$\displaystyle[\zeta]=[\zeta_{tw}]=ML^{-2}T^{-1}$

We set $\zeta=\zeta_{tw}$.

### Worm head trajectories

In this section, we construct a stochastic model for worm head trajectories in
2D, which we then use to formulate a stochastic, 3D active force,
$\mathbf{F}^{\text{act}}(t)$, that drives worms. Coupling this active force to
the elastic model described above (5) can then be shown to produce tangling
and untangling dynamics (movies S2 and S3).

#### 2D Worm head trajectories

Our 2D stochastic model for the worm head trajectories is based on
experimental data (Fig. 3). We assume the worm head has average speed $v$, and
turns at average angular speed $\alpha$. The turning direction switches from
clockwise to anticlockwise at rate $\lambda$. These speeds are augmented by
stochastic terms governed by a translation diffusivity $D_{T}$ and rotational
diffusivity $D_{R}$

$\displaystyle d\mathbf{X}$
$\displaystyle=v\mathbf{N}_{\theta}\,dt+\sqrt{2D_{T}}*d\mathbf{B}(t)$ (6a)
$\displaystyle d\mathbf{\theta}$
$\displaystyle=(-1)^{S(t)}\alpha\,dt+\sqrt{2D_{R}}*dW(t)$ (6b)

where $\mathbf{B}(t)$ and $W(t)$ are independent Brownian processes, and
$S(t)=\hat{S}(t)+S_{0}$, where $\hat{S}(t)$ is a Poisson process of rate
$\lambda$, i.e. $\hat{S}(t)\sim Po(\lambda t)$, and $S_{0}\sim Ber(1/2)$ is
independent of $\hat{S}(t)$. The process $S(t)$ is then effectively a rate
$\lambda$ Poisson process with random initialization, so $S(t)-S(0)\sim
Po(\lambda t)$. In particular,

$\displaystyle\mathbb{E}\left[(-1)^{S(t)}\right]=\mathbb{E}\left[(-1)^{\hat{S}(t)+S_{0}}\right]=\mathbb{E}\left[(-1)^{\hat{S}(t)}\right]\mathbb{E}\Big{[}(-1)^{S_{0}}\Big{]}=0$

In Langevin notation, writing $\boldsymbol{x}=\mathbf{X}$ and
$\boldsymbol{n}_{\theta}=\mathbf{N}_{\theta}$, equation (6) becomes

$\displaystyle\dot{\boldsymbol{x}}$
$\displaystyle=v\boldsymbol{n}_{\theta}+\boldsymbol{\xi}_{T}(t)$ (7a)
$\displaystyle\dot{\theta}$ $\displaystyle=(-1)^{S(t)}\alpha+\xi_{R}(t)$ (7b)

where $\xi_{T},\xi_{R}$ are independent white noise processes, satisfying

$\displaystyle\langle\left(\xi_{T}\right)_{i}(t)\left(\xi_{T}\right)_{j}(t^{\prime})\rangle=2D_{T}\delta_{ij}\delta(t-t^{\prime}),\qquad\langle\xi_{R}(t)\xi_{R}(t^{\prime})\rangle=2D_{R}\delta(t-t^{\prime})$

where angle brackets denote expectations. Assuming the noise terms are small,
the trajectory shapes are determined by $\alpha,\lambda$ and $v$. In
particular, the chirality number, $\gamma=\alpha/2\pi\lambda$ sets the number
of full loops the head makes before switching direction, and $v/\alpha$ sets
the size of these loops. Although the noise terms in equation (6) are taken to
be Brownian, the nature of the noise in the real, biological system is
expected to be more complex. However, the choice of Brownian noise makes this
model more analytically tractable.

#### Rotation index of 2D trajectories

When $\theta(0)=0$, the value of $\theta$ is equal to the total winding of the
velocity vector around the origin, and is known as the rotation index,
$I_{R}$, of the curve

$\displaystyle I_{R}=\frac{1}{2\pi}\theta(t)$

Although this is distinct from the winding number of a trajectory about a
point, intuitively, we expect that trajectories with larger values of $I_{R}$
should tangle more readily. We can find the $I_{R}$ by calculating moments of
$\theta$. Under the assumption of Brownian noise.

$\displaystyle\mathbb{E}[\theta(t)]=\alpha\int_{0}^{t}dt^{\prime}\,\mathbb{E}\left[(-1)^{S(t^{\prime})}\right]=0$

Using Ito’s lemma

$\displaystyle
d(\theta^{2})=2\left((-1)^{S(t)}\alpha\theta+D_{R}\right)\,dt+2\theta\sqrt{2D_{R}}*dW(t)$

Thus

$\displaystyle\frac{d}{dt}\mathbb{E}[\theta^{2}]$
$\displaystyle=2D_{R}+2\alpha\mathbb{E}[(-1)^{S(t)}\theta]$
$\displaystyle=2D_{R}+2\alpha^{2}\mathbb{E}\left[(-1)^{S(t)}\int_{0}^{t}dt^{\prime}\,(-1)^{S(t^{\prime})}\right]$
$\displaystyle=2D_{R}+2\alpha^{2}\int_{0}^{t}dt^{\prime}\,\mathbb{E}\left[(-1)^{S(t)-S(t^{\prime})+2S(t^{\prime})}\right]$
$\displaystyle=2D_{R}+2\alpha^{2}\int_{0}^{t}dt^{\prime}\,e^{-2\lambda(t-t^{\prime})}$
$\displaystyle=2D_{R}+2\alpha^{2}e^{-2\lambda
t}\frac{1}{2\lambda}\left(e^{2\lambda t}-1\right)$
$\displaystyle=2D_{R}+\frac{\alpha^{2}}{\lambda}\left(1-e^{-2\lambda
t}\right)$

Integrating this gives

$\displaystyle\mathbb{E}\left[\theta^{2}\right]=\left(2D_{R}+\frac{\alpha^{2}}{\lambda}\right)t-\frac{\alpha^{2}}{2\lambda^{2}}\left(1-e^{-2\lambda
t}\right)$

We can simplify this expression by using the chirality number,
$\gamma=\alpha/2\pi\lambda$, which counts the typical number of loops of a
given handedness formed by the trajectory before it switches turning
direction. In particular, the rotation index satisfies

$\displaystyle\mathbb{E}\left[I_{R}^{2}\right]=\frac{1}{4\pi^{2}}\mathbb{E}\left[\theta^{2}\right]=\frac{1}{2\pi}\left(\frac{D_{R}}{\pi\alpha}+\gamma\right)\alpha
t-\frac{1}{2}\gamma^{2}\left(1-e^{-\alpha t/\pi\gamma}\right)$

Defining the dimensionless time $\tilde{t}=\alpha t$, and the we can simplify
this expression further

$\displaystyle\mathbb{E}[I_{R}^{2}]=\frac{1}{2\pi}\left(\frac{D_{R}}{\pi\alpha}+\gamma\right)\tilde{t}-\frac{1}{2}\gamma^{2}\left(1-e^{-\tilde{t}/\pi\gamma}\right)$
(8)

This is an increasing function of $\gamma$, suggesting that trajectories with
larger $\gamma$ have a greater propensity to tangle.

#### 3D Worm head dynamics

3D head dynamics can be built from this 2D model by adding a perturbation in
the $z$-direction. To mimic the trajectories described above, we introduce an
active force $\mathbf{F}^{\text{act}}$ at the head of each worm which rotates
at rate $\alpha$ and switches rotation direction at rate $\lambda$. More
concretely, setting $\mathbf{F}^{\text{act}}=q\boldsymbol{n}_{\theta}$ where
$\boldsymbol{n}_{\theta}=(\cos\theta,\sin\theta,0)$, we define the following
dynamics for the active head force

$\displaystyle\dot{\theta}(t)$ $\displaystyle=(-1)^{S(t)}\alpha+\xi_{R}(t)$
(9a) $\displaystyle\mathbf{F}^{\text{act}}(t)$
$\displaystyle=q\boldsymbol{n}_{\theta}(t)+\boldsymbol{\xi}_{F}(t)$ (9b)

The $\theta$ equation and the Poisson process $S(t)$ are as in (7), and
$\boldsymbol{\xi}_{F}$ in the force equation is a white noise process in 3D
satisfying

$\displaystyle\langle\left(\xi_{T}\right)_{i}(t)\left(\xi_{T}\right)_{j}(t^{\prime})\rangle=2D_{F}\delta_{ij}\delta(t-t^{\prime})$

where $D_{F}$ has units $[D_{F}]=M^{2}L^{2}T^{-3}$. The noise term therefore
acts as a perturbation, moving the filaments out of plane. In practice, a
combination of the noise term and contact effects produce the 3D tangled
structures we see in simulations (movies S2 and S3). The parameters $\alpha$
and $\lambda$ are explicit in these active force equations. The speed $v$
follows from considering the equations of motion for the elastic fiber, and
can be solved for numerically using our simulation framework.

### Mean-field tangling model

Using trajectory equations for the worm head, we can compute how much a single
worm tangles around a fixed set of obstacles, which represent the structure of
a fixed background tangle. This idea leads to a mean-field theory of tangling.
Although our numerical tests of this theory will necessarily use the 3D
dynamics described above (equation 9), for simplicity, we will build the mean-
field model in 2D, using the SDEs described in equation (6) and (7)

$\displaystyle\dot{\boldsymbol{x}}$
$\displaystyle=v\boldsymbol{n}_{\theta}+\boldsymbol{\xi}_{T}(t)$
$\displaystyle\dot{\theta}$ $\displaystyle=(-1)^{S(t)}\alpha+\xi_{R}(t)$

We further assume that $D_{T}$ is negligible. As a proxy for tangling, we can
construct a tangling parameter based on the winding of the curve
$\boldsymbol{x}(t)$ around a specified set of points
$\Lambda\subset\mathbb{R}^{2}$. The set $\Lambda$ introduces another length
scale into our model. We let $\ell$ be the characteristic spacing between
points in $\Lambda$. The winding number of $\boldsymbol{x}(t)$ about
$p\in\Lambda$ in the time taken to travel one worm length, $t=L/v$ is

$\displaystyle W_{p}(\boldsymbol{x})$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{L/v}dt\;\frac{d}{dt}\arg\left(\boldsymbol{x}(t)-p\right)$

Since tangling depends on geometrical contact as well as topological winding,
we keep track of both the winding and contact around obstacles $p\in\Lambda$.
We define the contact winding, $cW_{p}(\boldsymbol{x})$, in the same way as
contact link (1)

$\displaystyle
cW_{p}(\boldsymbol{x})=\begin{cases}|W_{p}(\boldsymbol{x})|\quad&\text{if
}\min_{t\in[0,L/v]}|\boldsymbol{x}(t)-p|<r\\\ 0&\text{otherwise}\end{cases}$
(10)

where $r$ is the contact threshhold, so the curve is said to be in contact
with a point $p\in\mathbb{R}^{2}$ if it comes within distance $r$ of $p$. In
addition, we threshhold the winding number, so that winding around a point
$p\in\Lambda$ is only counted if $cW_{p}\geq 1$. This is necessary to exclude
certain configurations where $W_{p}$ is non-negligible but the interaction
between $\boldsymbol{x}$ and $p$ is marginal. For example, for
$\boldsymbol{x}(t)=(\epsilon,vt-L/2)$, and $p=0$, we have

$\displaystyle
W_{p}(\boldsymbol{x})=\frac{1}{\pi}\arctan\left(\frac{L}{2\epsilon}\right)$

Thus $W_{p}$ can approach $1/2$ even though the curve $\boldsymbol{x}$ lies
adjacent to the point $p$. This motivates the following tangling parameter,
which combines the topological information of winding with the geometric
information of contact

$\displaystyle\mathcal{T}=\left\langle\sum_{p\in\Lambda}\Theta\left(cW_{p}-1\right)\right\rangle$
(11)

where $\Theta$ is the Heaviside step function. The obstacles in $\Lambda$ can
be thought of as other worms in the tangle. The tangling index therefore
counts the number of worms a given worm interacts with topologically. In this
picture, the tangling index can be viewed as the mean-field analogue of the
contact link based tangle quantifiers constructed earlier. In particular,
$\mathcal{T}$ is constructed similarly to the total contact link per worm
(equation 2) and the mean degree of the tangle graph (equation 3). This
motivates a connection with percolation theory, as discussed above. We assume
that a tangled state has a connected tangle graph and ask which values of
$\mathcal{T}$ correspond to dynamics (equation 9) capable of producing tangled
states. Since a connected graph on $n$ vertices has mean degree at least
$2-2/n$, we conjecture that $\mathcal{T}^{\ast}\approx 2$ is the critical
value separating tangled states from untangled states. This agrees with our
experimental data (Fig. 4C). In addition, our tomographic reconstructions have
total contact link per worm approximately equal to 2 (Fig. 2F, and equation
4), suggesting that living worm tangles are tuned for ease of untangling.

####  .3.1 Analysis of tangling index

Figure S3: Effect of noise and contact threshold on mean-field phase diagram
The phase diagrams depend weakly on noise and contact threshold. In Fig. 4C,
we choose $r/\ell=1/3$ and $D_{R}/\alpha=0.4$ (middle column).

We can write $\mathcal{T}$ as a function of dimensionless parameters by non-
dimensionalizing the the trajectory equations (7). Let $\tilde{t}=\alpha t$
and $\tilde{\boldsymbol{x}}=\boldsymbol{x}/\ell$, where $\ell$ is the
characteristic spacing between points in $\Lambda$. Then equation
$\eqref{langevin}$ becomes

$\displaystyle\frac{d\tilde{\boldsymbol{x}}}{d\tilde{t}}$
$\displaystyle=R\boldsymbol{n}_{\theta}$ (12a)
$\displaystyle\frac{d\theta}{d\tilde{t}}$
$\displaystyle=(-1)^{\tilde{S}(\tilde{t})}+\tilde{\xi}_{R}(\tilde{t})$ (12b)

where $R=v/\alpha\ell$ is the loop number,
$\tilde{S}(\tilde{t})=S(\tilde{t}/\alpha)$, the rotational noise term
$\tilde{\xi}_{R}(\tilde{t})=\alpha^{-1}\xi_{R}(\tilde{t}/\alpha)$ satisfies

$\displaystyle\left\langle\tilde{\xi}_{R}(\tilde{t})\tilde{\xi}_{R}(\tilde{t}^{\prime})\right\rangle=\frac{2D_{R}}{\alpha}\delta(\tilde{t}-\tilde{t}^{\prime})$

and we have assumed that the translational noise coefficient, $D_{T}$, is
negligible. The process $\tilde{S}(\tilde{t})$ is a Poisson process with
dimensionless rate $\lambda/\alpha=1/2\pi\gamma$ where $\gamma$ is the
chirality number

$\displaystyle\tilde{S}(\tilde{t})-\tilde{S}(0)\sim
Po\left(\frac{\tilde{t}}{2\pi\gamma}\right)$

We therefore find three dimensionless parameters from $\eqref{SDEnondim}$, the
chirality number $\gamma$, the loop number $R$ and a non-dimensional noise,
$D_{R}/\alpha$. To calculate $\mathcal{T}$, we run the process for time
$t=L/v$ and quantify contact using the contact threshold $r$ (equation 10).
The gives two additional dimensionless parameters, $L/\ell$ and $r/\ell$. For
a given obstacle set, $\Lambda$, the tangling index can therefore be written
as

$\displaystyle\mathcal{T}=\mathcal{T}\left(\gamma,R\,;\frac{D_{R}}{\alpha},\frac{r}{\ell},\frac{L}{\ell}\right)$
(13)

In this study, we have focused on the role of the dynamical dimensionless
parameters, $\gamma$ and $R$, which arise from equation (6). The effect of the
noise and contact parameters can be explored numerically (Fig. S3) and do not
change the underlying structure of the phase diagram. In Fig. 4C, we choose
$D_{R}/\alpha=0.4$ and $r/\ell=1/3$. The remaining degrees of freedom,
$L/\ell$ and $\Lambda$ govern the placement and spacing of winding obstacles.

The worm length, $L$, affects $\mathcal{T}$ by setting the time for which the
process $\mathbf{X}$ is run. We set $L/\ell=25$ in the phase diagrams shown
here (Figs. 4C, S3-S5). Larger $L$ will increase $\mathcal{T}$, since the
curve $\mathbf{X}$ will have more time to wind around points. This matches
with our intuition that longer filaments tangle more easily. From ultrasound
data, we find $\ell=3h$, and $L/h\approx 80$, where $h$ is worm radius, which
gives $L/\ell\approx 27$. Note that this implies that the chosen contact
threshold, $r$, is less than $2h$. This can be understood as a consequence of
the fact that our mean-field model is treating a dynamic 3D tangle as a static
2D point cloud. For example, a worm of radius $h$ cannot pass between two
parallel worms separated by $\ell=3h$, without causing some deformation. The
fact that $r<2h$ is due to the fact that our model does not explicitly account
for such deformation effects. In particular, our mean-field model does not
explicitly depend upon the worm radius $h$, although worm radius effects are
captured by the contact threshold $r$.

The set of obstacles, $\Lambda$ acts as the background tangle through which a
worm moves. Although this background could be amorphous, for simplicity, here
we assume that $\Lambda$ has a lattice structure. Numerical investigations
further reveal that the tangling index $\mathcal{T}$, does not strongly depend
on lattice type (Fig. S4A,B). This validates our assumption that $\Lambda$
affects $\mathcal{T}$ mostly through the lattice spacing $\ell$.

#### Mean-field trajectories through $\Lambda$

Figure S4: Effect of obstacle placement and spacing. (A,B) Choosing a square
(A) or triangular (B) lattice for the obstacle set $\Lambda$ does not affect
the structure of the tangling phase diagram. (C) Changing the $R$-axis reveals
a range of $R$-values for which tangling is mostly controlled by $\gamma$,
with small $\gamma$ producing untangled states and large $\gamma$ producing
tangled states. All phase diagrams (A-C) have
$D_{R}/\alpha=0.4,\,r/\ell=1/3,\,L/\ell=25$. (D-F) Sample trajectories run for
$0<t<0.4L$ corresponding to points on the phase diagrams (top row) marked by
black circles. Disks indicate obstacles in $\Lambda$, colored by the value of
their contact winding, $cW_{p}$, due to the displayed trajectory. (D,E) Near
the critical tangling index, $\mathcal{T}^{\ast}\approx 2$, trajectories
typically wind around points in both the clockwise and counterclockwise
directions. This can cause winding cancellation, where a point has low winding
despite being enclosed by a curve (D). These trajectories could be part of a
topological quick release mechanism, and are reminiscent of the “picture-
hanging puzzle” [10]. (F) Trajectories with large $R$ have very small
curvature, thus leading to untangled states. Friction effects often prevent
these linear motions from being a feasible untangling method. Trajectories in
(D-F) have low noise ($D_{R}/\alpha=0.2$) in order to illustrate their
topological properties.

The mean-field trajectories described by our phase diagrams (Fig. S4A-C)
capture a range of tangling and untangling gaits, including near-linear motion
at large $R$ (Fig. S4C). In particular these gaits often display interesting
topological motifs. For example, trajectories with tangling index close to
critical ($\mathcal{T}^{\ast}\approx 2$) typically wind around points in both
the clockwise and counterclockwise directions (Fig. S4D,E). This chirality
switching behavior can cause the trajectory to cancel out an earlier winding
by passing around the same point twice in opposite directions (Fig. S4D).
Physically if a worm winds clockwise around a point $p$, then clockwise around
another point $p^{\prime}$, and finally anticlockwise around $p$ once more,
the worm will still be entangled with $p$ due to winding around $p^{\prime}$.
These multi-worm interactions are not captured by our mean-field model,
however, the presence of such winding cancellation effects in critical worm
trajectories suggests that critical tangles are equipped with additional
topological quick release mechanisms for rapid untangling. We note that the
trajectories in Fig. S4D,E resemble the solution to the famous “picture-
hanging puzzle” [10], which is an example of one such quick release mechanism.

Expanding the $R$-axis of our phase diagrams demonstrates that extreme values
of $R$ produce untangling gaits (Fig. S4C, F). However, the small $R$
untangling gait requires producing very high curvatures, and is therefore
energetically unfavorable. Similarly, the large $R$ trajectory appears nearly
linear (Fig. S4F). However attempting linear motion from an initially tangled
state will incur a high friction cost, owing to capstan-like friction effects
[11]. Elastic relaxation to such a linear state will also require relatively
long timescales.

### Estimation of model parameters from data

Figure S5: Effect of smoothing on estimated worm head trajectory parameters
Worm head trajectory data are smoothed using a lowpass filter. The estimated
parameters depend only weakly on the choice of cutoff frequency. In Fig. 4C we
choose frequency cutoffs of $15$ms and $1.5$s (middle diagram). Background
phase diagrams have $D_{R}/\alpha=0.4,\,r/\ell=1/3,\,L/\ell=25$ (Fig. S3) as
in Fig. 4C.

Here we describe the estimation of the quantities $\alpha,\lambda,v,\ell$. The
parameters $\alpha,\lambda,v$ are obtained from trajectory data, and $\ell$,
which is a measure of tangle tightness, is obtained from the ultrasound data
sets.

We begin with the procedure for estimating $\alpha,\lambda,v$, by fitting
discrete trajectory data $\boldsymbol{x}_{i}=\boldsymbol{x}(t_{i})$, with
$t_{i+1}-t_{i}=\delta t$, to our SDE model (equation 6)

$\displaystyle d\mathbf{X}$
$\displaystyle=v\mathbf{N}_{\theta}\,dt+\sqrt{2D_{T}}*d\mathbf{B}(t)$
$\displaystyle d\mathbf{\theta}$
$\displaystyle=(-1)^{S(t)}\alpha\,dt+\sqrt{2D_{R}}*dW(t)$

As stated previously, we do not expect the noise in the real, biological
system to be Gaussian. Instead we proceed by smoothing the trajectory, and we
assume that this provides sufficient denoising to estimate $\alpha,\lambda,v$,
directly.

We first smooth $\boldsymbol{x}(t_{i})$ to obtain
$\tilde{\boldsymbol{x}}_{i}$. The smoothing operation used was the inbuilt
lowpass filter in MATLAB R2021a. The cutoff frequency chosen was
$1/\omega=15\,$ms for the untangling trajectories, and $1/\omega=1.5\,$s for
the tangling trajectories. However the estimated values of
$\gamma=\alpha/2\pi\lambda$ and $R=v/\alpha\ell$ do not depend strongly on the
choice of cutoff frequency (Fig. S5). From $\tilde{\boldsymbol{x}}$ we
reconstruct the velocity and angle

$\displaystyle\tilde{\boldsymbol{v}}_{i}$
$\displaystyle=\frac{\tilde{\boldsymbol{x}}_{i+1}-\tilde{\boldsymbol{x}}_{i-1}}{2\delta
t}$ $\displaystyle\theta_{i}$ $\displaystyle=\arg\tilde{\boldsymbol{v}}_{i}$

Finally, let $\tilde{\theta}_{i}$ be the result of passing $\theta_{i}$
through the same lowpass filter as $\boldsymbol{x}(t)$. We assume that
$\tilde{\boldsymbol{x}}_{i},\tilde{\boldsymbol{v}},\tilde{\theta}_{i}$ are
then smooth enough that $\alpha,\lambda$ and $v$ can be estimated directly.

$\displaystyle\alpha=\left\langle\left(\frac{\tilde{\theta}_{i+1}-\tilde{\theta}_{i}}{\delta
t}\right)^{2}\right\rangle^{1/2},\qquad
v=\left\langle\tilde{\boldsymbol{v}}_{i}\cdot\tilde{\boldsymbol{v}_{i}}\right\rangle^{1/2}$

This assumption is valid when the noise in the trajectory data is concentrated
in the high frequency sector.

The value of $\ell$ will typically depend on how loose or tight the tangle is.
We choose $\ell$ from a tangle follows. Let $\boldsymbol{x}_{i}(s)$, denote
the centerline curve of the $i$’th worm, where $s\in[0,L_{i}]$ and $L_{i}$ is
the length of the $i$’th worm. We take $\ell$ to be the average minimum
distance between a worm head and the tangle

$\displaystyle\ell=\frac{1}{N}\sum_{i}\min_{s,j\neq
i}|\boldsymbol{x}_{i}(L_{i})-\boldsymbol{x}_{j}(s)|$ (14)

where the curves are oriented so $\boldsymbol{x}_{i}(L_{i})$ is the worm head,
and $N$ is total number of worms. Since the worm heads and tails are
indistinguishable in our ultrasound data sets, in this case we measure $\ell$
by

$\displaystyle\ell=\frac{1}{2N}\sum_{i}\min_{s,j\neq
i}|\boldsymbol{x}_{i}(L_{i})-\boldsymbol{x}_{j}(s)|+\frac{1}{2N}\sum_{i}\min_{s,j\neq
i}|\boldsymbol{x}_{i}(0)-\boldsymbol{x}_{j}(s)|$

For our ultrasound reconstructions, this calculation gives $\ell\approx
1.5$mm, or $\ell\approx 3h$ in terms of the worm radius.

### Tangling and untangling simulations

Figure S6: Tangling and untangling simulations and their corresponding mean-
field parameters (A) Mean-field parameters corresponding to tangling (purple
markers) and untangling (green markers) gaits from experiments (circles) and
simulations (triangles, squares and diamonds). Circles represent average
values over data from 25 worms (green) and 18 worms (purple) shown in Fig. 4C.
Triangles denote simulations from Fig. 3E and movie S2; squares and diamonds
denote simulations from Fig. 4D and movie S3. The squares represent
simulations with low head speed (Fig. 4D, middle column) and diamonds
represent simulations with large head speed (Fig. 4D, right hand column). The
mean-field phase diagram is predictive of whether the simulations result in
tangled or untangled states (movies S2 and S3). (B) The tangle spacing $\ell$
fluctuates during a tangling simulation (Fig. 4D and movie S3). We measure
$\ell$ by taking the average value between $t=50/\alpha$ and $t=100/\alpha$.
Mean-field simulation parameters given by purple square in (A). (C) $19$-worm
artificial tangle obtained by taking the largest cluster of touching worms
generated by the simulation analyzed in (B) and shown in Fig. 4D (middle
column) and movie S3. This tangle is used as the initial condition for the
large speed simulations show in Fig. 4D (right hand column) and movie S3.
Scale bar 10$h$. (D) Tangle graph for the artificial tangle in (C). Edges are
present between worms which have contact link greater than $1/2$. The tangle
graph is not connected since the artificial tangle in (C) is chosen as a
cluster of touching, but not necessarily tangled, worms. (E) The relationship
between total contact link per worm and tube radius for this artificial tangle
is similar to those observed for real living tangles (Fig. 2F). The tube
radius dependence of total contact link also identifies an effective radius
worm radius, $h_{\text{eff}}=1.25h$ that arises from contact handling in
simulations.

Our tangling and untangling simulations (Figs. 3-4, movies S2-S3) agree with
mean-field predictions (Fig. S6A). These simulations use the elasticity
framework and the stochastic active head force (9) discussed above. The
parameters $\alpha,\lambda,v$ can therefore be directly input. On the other
hand, $\ell$ must be measured. For simulations which begin from a tangled
initial condition, a value of $\ell$ can be found by applying (14) to this
initial condition. On the other hand, during a tangling simulation, the
initial condition may not reflect the equilibrium value of $\ell$ (Fig. S6B).
To measure the tangle tightness in this case, we take the average value,
$\langle\ell\rangle$, between $t=50/\alpha$ and $t=100/\alpha$. Empirically,
this appears to give the value of $\ell$ time to equilibrate from the initial
condition, without including the decrease in $\ell$ that occurs once tangled
clusters begin to form (Fig. S6B).

In our simulations, we find values of $\ell$ that range from $7h$ to $10h$,
where $h$ is the worm radius. In contrast, our ultrasound reconstructions have
$\ell\approx 3h$. This discrepancy could arise from a variety of factors. For
example the real worms are soft and compressible, which could cause tighter
packings. The worms in the ultrasound data sets also have a higher variance in
length (Fig. S1) and smaller average length ($24\,$mm) than the simulated
worms (constant length $40\,$mm). The range of different lengths could enable
a collective to pack more densely. Finally, we note that the ratio $\ell/h$
could be sensitive to the precise contact handling method used in simulations,
and so part of the observed difference in the value of $\ell/h$ could be due
to the discretization error.

Overall, generated tangles have similarities to ultrasound reconstructions
(Fig. S6C-E). In particular, the radius dependence of total contact link per
worm exhibits a similar shape for the generated tangle (Fig. S6E) as for the
ultrasound data sets (Fig. 2F). The value of total contact link per worm in
the tangle shown in Fig. S6C is 3.6, which is larger than the observed values
for ultrasound tangles (1.7 to 2.1, Fig. 2F). We attribute this to the fact
that the dynamics used to generate this artificial tangle correspond to a
larger tangling index (Fig. S6A) than is observed for tangling worms (Fig.
4C).

####  .3.2 Simulation parameters

Our simulation parameters correspond to the following values for worm radius
$h$, worm length $L$, density $\rho$ and bending modulus $E_{b}$

$\displaystyle h=0.5\,\text{mm},\quad
L=40\,\text{mm},\quad\rho=10^{-3}\,\text{g mm}^{-3},\quad E_{b}=1\,\text{kPa}$

In addition we set an effective bulk modulus $K$, damping coefficients
$\gamma,\eta$, friction $\zeta$ and effective Young’s modulus $E$. The first
three parameters are fixed across simulations

$\displaystyle K=170\,\text{kPa},\quad\gamma=0.1\,\text{g
mm}^{-3}\text{s}^{-1},\quad\eta=3.7\,\text{g mm}^{-1}\text{s}^{-1}$

Our simulations reproduce tangling and untangling behavior (Fig. S6A) for a
range of values of the effective Young’s modulus $E$. The tangling and
untangling simulations in Fig. 3E (movie S2) and the low head speed
simulations in Fig. 4D (movie S3) have $E=10\,$kPa and the large head speed
simulations in Fig. 4D (movie S3) have $E=100\,$kPa. Note that the order of
magnitude of the values chosen for $E_{b},E$ and $K$ all fall within the
observed range of elastic moduli for C. elegans [12]. Finally, our simulations
are also robust to changes in the friction parameter $\zeta$. For the tangling
simulation in Fig. 3E and the low head speed simulations in Fig. 4D, we take
$\zeta=2.7\,\text{g mm}^{-2}\text{s}^{-1}$. The untangling simulation in Fig.
3E has $\zeta=2.2\,\text{g mm}^{-2}\text{s}^{-1}$ and the large head speed
simulations in Fig. 4D have $\zeta=2.0\,\text{g mm}^{-2}\text{s}^{-1}$.

The parameters of the active head forcing (equation 9) are indicated in Fig.
S6A, and given explicitly below. The force noise term is constant across all
simulations, $D_{F}=0.3\,\text{g}^{2}\text{mm}^{2}\text{s}^{-3}$. The
rotational noise strength varies from $D_{R}=5\,\text{s}^{-1}$ to
$D_{R}=10\,\text{s}^{-1}$ across simulations, which corresponds to
$D_{R}/\alpha=0.6$ for tangling simulations and $D_{R}/\alpha=0.04$ for
untangling simulations. The other parameters (to 2.s.f) are

$\displaystyle v=110\,\text{mm
s}^{-1},\quad\ell=5.0\,\text{mm},\quad\alpha=8.4\,\text{s}^{-1},\quad\lambda=0.7\,\text{s}^{-1}\qquad$
(Fig. S6A, purple triangle) $\displaystyle v=410\,\text{mm
s}^{-1},\quad\ell=1.7\,\text{mm},\quad\alpha=170\,\text{s}^{-1},\quad\lambda=130\,\text{s}^{-1}\qquad$
(Fig. S6A, green triangle) $\displaystyle v=160\,\text{mm
s}^{-1},\quad\ell=5.5\,\text{mm},\quad\alpha=8.4\,\text{s}^{-1},\quad\lambda=0.7\,\text{s}^{-1}\qquad$
(Fig. S6A, purple square) $\displaystyle v=160\,\text{mm
s}^{-1},\quad\ell=5.2\,\text{mm},\quad\alpha=8.4\,\text{s}^{-1},\quad\lambda=11\,\text{s}^{-1}\qquad$
(Fig. S6A, green square) $\displaystyle v=900\,\text{mm
s}^{-1},\quad\ell=3.3\,\text{mm},\quad\alpha=250\,\text{s}^{-1},\quad\lambda=20\,\text{s}^{-1}\qquad$
(Fig. S6A, purple diamond) $\displaystyle v=900\,\text{mm
s}^{-1},\quad\ell=3.3\,\text{mm},\quad\alpha=250\,\text{s}^{-1},\quad\lambda=200\,\text{s}^{-1}\qquad$
(Fig. S6A, green diamond)

Our phase diagram therefore predicts the outcome of multi-filament simulations
over a wide range of parameters, including for different values of noise.

## References

  * Tuazon  _et al._ [2022] H. Tuazon, E. Kaufman, D. I. Goldman, and M. S. Bhamla, Oxygenation-controlled collective dynamics in aquatic worm blobs, Integr. Comp. Biol. 10.1093/icb/icac089 (2022).
  * Schneider  _et al._ [2012] C. A. Schneider, W. S. Rasband, and K. W. Eliceiri, Nih image to imagej: 25 years of image analysis, Nat. Methods 9, 671 (2012).
  * Tinevez  _et al._ [2017] J.-Y. Tinevez, N. Perry, J. Schindelin, G. M. Hoopes, G. D. Reynolds, E. Laplantine, S. Y. Bednarek, S. L. Shorte, and K. W. Eliceiri, Trackmate: An open and extensible platform for single-particle tracking, Methods 115, 80 (2017).
  * Qinami and Grinspun [2017] K. Qinami and E. Grinspun, Quantifying geometric entanglement: the linking number of two open curves, [Unpublished manuscript] (2017).
  * Bergou  _et al._ [2008] M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, and E. Grinspun, Discrete elastic rods, ACM Trans. Graph. 27, 63 (2008).
  * Patil  _et al._ [2020a] V. P. Patil, J. D. Sandt, M. Kolle, and J. Dunkel, Topological mechanics of knots and tangles, Science 367, 71 (2020a).
  * Patil  _et al._ [2020b] V. P. Patil, Ž. Kos, M. Ravnik, and J. Dunkel, Discharging dynamics of topological batteries, Phys. Rev. Research 2, 043196 (2020b).
  * Tong  _et al._ [2022] D. Tong, A. Choi, J. Joo, and M. K. Jawed, A fully implicit method for robust frictional contact handling in elastic rods, arXiv preprint arXiv:2205.10309 (2022).
  * Bergou  _et al._ [2010] M. Bergou, B. Audoly, E. Vouga, M. Wardetzky, and E. Grinspun, Discrete viscous threads, ACM Trans. Graph. 29, 1 (2010).
  * Demaine  _et al._ [2014] E. D. Demaine, M. L. Demaine, Y. N. Minsky, J. S. Mitchell, R. L. Rivest, and M. Pǎtraşcu, Picture-hanging puzzles, Theory Comput. Syst. 54, 531 (2014).
  * Sano  _et al._ [2022] T. G. Sano, P. Johanns, P. Grandgeorge, C. Baek, and P. M. Reis, Exploring the inner workings of the clove hitch knot, Extreme Mech. Lett. 55, 101788 (2022).
  * Gilpin  _et al._ [2015] W. Gilpin, S. Uppaluri, and C. P. Brangwynne, Worms under pressure: bulk mechanical properties of c. elegans are independent of the cuticle, Biophys. J. 108, 1887 (2015).