# Flocking and Collision Avoidance for a Dynamic Squad of Fixed-Wing UAVs
Using Deep Reinforcement Learning

Chao Yan, Xiaojia Xiang, Chang Wang*, and Zhen Lan This work was supported in
part by the Science and Technology Innovation 2030-Key Project of New
Generation Artificial Intelligence under Grant 2020AAA0108200, in part by the
National Natural Science Foundation of China under Grant 61825305, Grant
61803377, and Grant 61906203.College of Intelligence Science and Technology,
National University of Defense Technology, Changsha 410073, China, {yanchao17,
xiangxiaojia, wangchang07<EMAIL_ADDRESS>*Corresponding author

###### Abstract

Developing the flocking behavior for a dynamic squad of fixed-wing UAVs is
still a challenge due to kinematic complexity and environmental uncertainty.
In this paper, we deal with the decentralized flocking and collision avoidance
problem through deep reinforcement learning (DRL). Specifically, we formulate
a decentralized DRL-based decision making framework from the perspective of
every follower, where a collision avoidance mechanism is integrated into the
flocking controller. Then, we propose a novel reinforcement learning algorithm
PS-CACER for training a shared control policy for all the followers. Besides,
we design a plug-n-play embedding module based on convolutional neural
networks and the attention mechanism. As a result, the variable-length system
state can be encoded into a fixed-length embedding vector, which makes the
learned DRL policy independent with the number and the order of followers.
Finally, numerical simulation results demonstrate the effectiveness of the
proposed method, and the learned policies can be directly transferred to semi-
physical simulation without any parameter finetuning.

##  I Introduction

Due to the capability limitation of a single unmanned aerial vehicle (UAV),
multi-UAV collaboration has attracted increasing attention [1, 2, 3]. One of
the fundamental and challenging problems is the flocking control of UAVs
without collision [4]. Traditional methods such as model predictive control
[5] and consensus theory [6] usually depend on precise physical models, which
are complex and difficult to obtain in practice.

As an alternative, reinforcement learning (RL) [7], [8] can be used for the
flocking control problem. For example, Hung et al. used the
Dyna-Q($\lambda{}$) algorithm [9] and the Q($\lambda{}$) algorithm [10] to
learn flocking control policies for fixed-wing UAVs. Speck et al. [11]
combined the SARSA algorithm with object-focused learning to implement the
formation control of UAV swarms. However, the above methods discretized the
state and action spaces, which was inappropriate for controlling the UAVs in
more realistic environments. Besides, the UAVs were assumed to fly at
different fixed-altitudes to simplify the collision avoidance problem.

Deep reinforcement learning (DRL) has been proved effective for high-
dimensional and continuous control problems in robotics [12, 13, 14]. In our
previous work, we proposed a DRL algorithm, i.e., continuous actor-critic with
experience replay (CACER), for flocking with fixed-wing UAVs in continuous
state and action spaces [15, 16]. Following the previous work [9, 10], we also
simplified the collision avoidance problem by assuming that the UAVs
maintained the same constant speed but different altitudes. However, some
real-world applications such as surveying and mapping require the UAVs to
flock at the same altitude. This requirement makes flocking control more
challenging, because the collision avoidance problem must be considered. In
contrast to the previous work [9], [10], [15], this paper solves a more
challenging flocking control problem. Specifically, we not only allow the UAVs
to change speed as needed, but also design a collision avoidance mechanism and
integrate it into the flocking controller.

Besides, it is difficult to construct a leader-follower flocking control model
for a dynamic squad of UAVs. The main reason is that the length of the system
state is relevant with the number of followers, while the DRL-based control
policies typically require a fixed-length input for deep neural networks. Sui
et al. [17] followed the method in [18] where long short-term memory network
(LSTM) [19] was applied to process the states of other followers sequentially
in the reverse order of their distances to the decision-making follower, under
the assumption that the nearest neighbor had the biggest effect on the
follower. However, this assumption is not always true due to the influence of
other factors such as speed and heading. Instead, we design a customized
network module to deal with the collision avoidance problem for a dynamic
squad of fixed-wing UAVs.

The main contributions of this paper are as follows:

  * •

A decentralized DRL-based framework is developed to address the flocking
control and collision avoidance problem for a dynamic squad of fixed-wing
UAVs.

  * •

A novel DRL algorithm is proposed for training the flocking controller, where
a plug-n-play embedding module based on convolutional neural networks and the
attention mechanism is designed to handle variable-length system state.

  * •

The proposed method can be directly transferred from numerical simulation to
semi-physical simulation without any parameter finetuning.

The rest of this paper is organized as follows. Section II formulates the
flocking and collision avoidance problem, followed by the proposed PS-CACER
algorithm and the designed SEMP embedding module in Section III. Section IV
discusses the simulation results. Finally, Section V concludes this paper.

##  II Problem Formulation

In this section, we describe the flocking problem and the kinematics of fixed-
wing UAV. Then, we formulate this problem as a Markov decision process (MDP)
in the RL framework.

###  II-A Problem Description

In our flocking scenario, a unique leader is followed by a variable number of
followers. These followers are homogeneous and fly at the same altitude [10,
20]. The leader is remotely controlled by a human operator via a ground
control station. We assume that each UAV is able to obtain the state of other
UAVs through the inner communication channel [21]. Each follower has to select
its steering commands independently to maintain a certain distance $\rho$ from
the leader ( $d_{1}<\rho<d_{2}$ ) while avoiding collisions with other
followers simultaneously, as shown in Figure 1.

Figure 1: The top view of the relationship between the leader and followers.
Note that the number of followers, i.e., $n$, is a variable.

###  II-B Kinematics of Fixed-Wing UAVs

We use the same numerical kinematic model of fixed-wing UAVs with stochastic
disturbances as in [16] to generate simulated samples for training the DRL-
based control policy. We rewrite the established kinematics as follows:

${\rm{\dot{\xi}}}=\frac{d}{{dt}}\left\\{{\begin{array}[]{*{20}{l}}x\\\ y\\\
\psi\\\ \phi\\\
v\end{array}}\right\\}=\left\\{{\begin{array}[]{*{20}{c}}{v\cos\psi+{\eta_{x}}}\\\
{v\sin\psi+{\eta_{y}}}\\\
{-\left({{\alpha_{g}}/v}\right)\tan\phi+{\eta_{\psi}}}\\\
{f_{\phi}\left({\phi,{\phi_{d}}}\right)}\\\
{f_{v}\left({v,{v_{d}}}\right)}\end{array}}\right\\},$ (1)

where $\left({x,\;y}\right)$ is the planar position; $\psi$ is the heading
angle; $\phi$ is the roll-angle; $v$ is the airspeed; ${\alpha_{g}}$ is the
gravitational acceleration;
$\left({{\eta_{x}},\;{\eta_{y}},\;{\eta_{\psi}}}\right)$ are disturbance terms
that follow normal distributions; the roll dynamics
$f_{\phi}\left({\phi,\;{\phi_{d}}}\right)$ and the airspeed dynamics
$f_{v}\left({v,\;{v_{d}}}\right)$ describe the response relationship between
(i) the desired airspeed ${v_{d}}$ and the actual airspeed $v$, and (ii) the
roll-angle setpoint ${\phi_{d}}$ and the actual roll-angle $\phi$,
respectively.

###  II-C MDP of Flocking with Collision Avoidance

####  II-C1 State Representation

As mentioned above, there are one leader and $n$ followers. From the
perspective of an arbitrary follower, e.g., follower-$i$ ($i=1,\ 2,\ \ldots,\
n$), the flock of UAVs can be divided into three groups: the follower-i itself
(termed as the ego-follower), the leader, and the other followers
{follower-$j$ $|$ $j=1,\ 2,\ldots,\ n,\ j\not=i\\}$. Accordingly, the system
state is composed of three parts: the ego-follower’s state
${{\rm{\xi}}_{e}}={\rm{\xi}}_{f}^{i}$, the leader’s state ${{\rm{\xi}}_{l}}$,
and the other followers’ states ${\rm{\xi}}_{o}=\\{{\rm{\xi}}_{f}^{j}$ $|$
$j=1,\ 2,\ \ldots,\ n,\ j\not=i\\}$, where
${\rm{\xi}}:=(x,\;y,\;\psi,\;\phi,\;v)$ is a tuple that describes the state of
a single fixed-wing UAV according to its kinematics.

To reduce the redundancy of the representation of the system state, the ego-
follower and leader joint state
${s^{e}}:=\left({s_{1}^{e},\;s_{2}^{e},\;s_{3}^{e},\;s_{4}^{e},\;s_{5}^{e},\;s_{6}^{e},\;s_{7}^{e},\;s_{8}^{e},\;s_{9}^{e}}\right)$
is defined as:

$\begin{aligned} \left[{\begin{array}[]{*{20}{l}}{s_{1}^{e}}\\\
{s_{2}^{e}}\end{array}}\right]&=\left[{\begin{array}[]{*{20}{c}}{\cos{\psi_{l}}}&{\sin{\psi_{l}}}\\\
{-\sin{\psi_{l}}}&{\cos{\psi_{l}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{l}}{{x_{e}}-{x_{l}}}\\\
{{y_{e}}-{y_{l}}}\end{array}}\right]\\\ s_{3}^{e}&={\psi_{e}}-{\psi_{l}}\\\
\vspace{0.5ex}s_{4}^{e}&={\phi_{e}}\\\ \vspace{0.5ex}s_{5}^{e}&={\phi_{l}}\\\
\vspace{0.5ex}s_{6}^{e}&=\phi_{d}^{l}\\\ \vspace{0.5ex}s_{7}^{e}&={v_{e}}\\\
\vspace{0.5ex}s_{8}^{e}&={v_{l}}\\\
\vspace{0.5ex}s_{9}^{e}&=v_{d}^{l}\end{aligned},$ (2)

where $(s_{1}^{e},s_{2}^{e})$ is the planar position of the ego-follower
relative to the leader; $s_{3}^{e}\in[-\pi,\pi)$ is the difference in the
heading between the leader and the follower. Besides, $\phi_{d}^{l}$ and
$v_{d}^{l}$ are the roll-angle and velocity setpoints of the leader,
respectively.

The ego-follower and other followers joint state
$s^{o}=\\{s_{j}^{o}:=\left({s_{1,j}^{o},\;s_{2,j}^{o},\;s_{3,j}^{o},\;s_{4,j}^{o},\;s_{5,j}^{o}}\right)$$|{}$
$j=1,\ 2,\ \ldots,\ n,\ j\not=i\\}$ is defined as:

$\begin{aligned} \left[{\begin{array}[]{*{20}{l}}{s_{1,j}^{o}}\\\
\vspace{0.5ex}{s_{2,j}^{o}}\end{array}}\right]&=\left[{\begin{array}[]{*{20}{c}}{\cos{\psi_{e}}}&{\sin{\psi_{e}}}\\\
\vspace{0.5ex}{-\sin{\psi_{e}}}&{\cos{\psi_{e}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{l}}{{x_{f}^{j}}-{x_{e}}}\\\
\vspace{0.5ex}{{y_{f}^{j}}-{y_{e}}}\end{array}}\right]\\\
s_{3,j}^{o}&={\psi_{f}^{j}}-{\psi_{e}}\\\
\vspace{0.5ex}s_{4,j}^{o}&={\phi_{f}^{j}}\\\
\vspace{0.5ex}s_{5,j}^{o}&=v_{f}^{j}\\\ \end{aligned},$ (3)

where $\left({s_{1,j}^{o},s_{2,j}^{o}}\right)$ represents the planar position
of another follower-$j$, relative to the ego-follower;
$s_{3,j}^{o}\in[-\pi,\pi)$ denotes the difference in the heading between the
ego-follower and follower-$j$.

Consequently, concatenating the two joint states above, we construct the state
representation of the system $s:=\left({s^{e},\;s^{o}}\right)$. We note that
its dimensionality depends on the number of followers.

####  II-C2 Action Space

The followers are maneuvered by executing the roll and velocity actions. In
this paper, we define the action $a:=\left({a_{r},\;a_{v}}\right)$ in
continuous spaces, where the roll action
${a_{r}}\in[-\frac{\pi}{18},\frac{\pi}{18}]$ and the velocity action
${a_{v}}\in\left[{-1,\;+1}\right]$. The current roll-angle and airspeed of a
follower are denoted by $\phi$ and $v$, respectively. The next roll-angle
setpoint ${\phi_{d}}$ can be calculated by:

$\phi_{d}=\left\\{\begin{array}[]{ll}\vspace{0.5ex}r_{\mathrm{bd}}&\mathrm{if}\
\phi+a_{r}>r_{\mathrm{bd}}\\\ \vspace{0.5ex}-r_{\mathrm{bd}}&\mathrm{if}\
\phi+a_{r}<-r_{\mathrm{bd}}\\\
\vspace{0.5ex}\phi+a_{r}&\mathrm{otherwise}\end{array}\right.,$ (4)

where $[-{r_{{\rm{bd}}}},\;{r_{{\rm{bd}}}}]$ is the allowed range of the roll-
angle setpoint.
Similarly, the next velocity setpoint ${v_{d}}$ is defined by:

$v_{d}=\left\\{\begin{array}[]{ll}\vspace{0.5ex}{v_{\max}}&\mathrm{if}\
v+a_{v}>{v_{\max}}\\\ \vspace{0.5ex}{v_{\min}}&\mathrm{if}\
v+a_{v}<{v_{\min}}\\\
\vspace{0.5ex}v+a_{v}&\mathrm{otherwise}\end{array}\right.,$ (5)

where the ${v_{\max}}$ and ${v_{\min}}$ denote the maximum velocity and the
minimum velocity of the follower, respectively.

####  II-C3 Reward Function

In this paper, the followers aim to avoid collision between each other while
flocking with the leader. Therefore, the reward function consists of two
pieces: the flocking reward and the collision penalty. Specifically, in order
to facilitate the agent (i.e., the ego-follower) to maintain a suitable
distance from the leader, the flocking reward ${r_{l}}$ is designed as follows
[10]:

$r_{l}=-\max\left\\{{{d_{e}},\;\frac{{{d_{1}}\left|{s_{3}^{e}}\right|}}{{\pi\left({1+\omega{d_{e}}}\right)}}}\right\\},$
(6)

where
${d_{e}}=\max\left\\{{m\left({{d_{1}}-\rho}\right),0,\rho-{d_{2}}}\right\\}$.
$d_{1}$ and $d_{2}$ are the inner and outer radius of the desired annulus
shown in Figure 1, respectively.
$\rho=\sqrt{{{\left|{s_{1}^{e}}\right|}^{2}}+{{\left|{s_{2}^{e}}\right|}^{2}}}$
denotes the distance between the leader and the agent. Both $\omega$ and $m$
are tuning parameters.

Besides, the collision penalty $r_{c}^{j}$ is defined to prevent the agent
from colliding with the follower-$j$, as follows:

$r_{c}^{j}=-\max\left\\{{m\left({{d_{1}}-\rho_{j}^{o}}\right),0}\right\\},$
(7)

where
$\rho_{j}^{o}=\sqrt{{{\left|{s_{1,j}^{o}}\right|}^{2}}+{{\left|{s_{2,j}^{o}}\right|}^{2}}}$
is the distance from the agent to follower-$j$.

Lastly, the final reward function $r$ is specified as:

$r={r_{l}}+\sum\nolimits_{j}{r_{c}^{j}}.$ (8)

##  III Approach

###  III-A PS-CACER

In our previous work [15], [16], we proposed an algorithm called continuous
actor-critic with experience replay (CACER). In contrast to other actor-critic
algorithms, e.g., DDPG [22], one distinctive feature of CACER is in its
positive-temporal difference (TD) update scheme [23] for training the policy
(actor). In other words, CACER updates its policy only when the TD-error is
positive [24]. In this paper, we propose a novel algorithm PS-CACER by
extending CACER to a multi-agent scenario (see Algorithm 1). Specifically, PS-
CACER has a plug-n-play embedding module (see Section III-B) based on
convolutional neural networks and the attention mechanism to handle variable-
length system state, while CACER only deals with fixed-length system state.

Algorithm 1 PS-CACER

0: $N_{s}\ -$ maximum time steps; $N_{b}\ -$ training batch size; $M\ -$
desired number of training episodes

1: Empty replay memory $D$ with capacity $N$

2: Initialize actor $Act_{\theta^{A}}(s)$ and critic $V_{\theta^{V}}(s)$
randomly

3: for $episode=1$ to $M$ do

4: Initialize the number of followers $n$ randomly

5: Represent initial state $s_{i}$ for each follower $i$

6: for $t=1$ to $N_{s}$ do

7: for follower $i=1$ to $n$ do

8: Select action ${a_{i}}\leftarrow
Act_{\theta^{A}}\left(s_{i}\right)+\mathcal{N}$ with the current policy and
Gaussian exploration noise

9: Execute the selected action ${a_{i}}$

10: end for

11: for follower $i=1$ to $n$ do

12: Represent new system state ${s^{\prime}_{i}}$

13: Calculate immediate reward ${r_{i}}$

14: Store tuple $({s_{i}},{a_{i}},{r_{i}},s^{\prime}_{i})$ in $D$

15: $s_{i}\leftarrow s_{i}^{\prime}$

16: end for

17: Sample a minibatch of $N_{b}$ tuples
$({s_{k}},{a_{k}},{r_{k}},s^{\prime}_{k})$ from $D$

18: Empty temporal buffer $D^{\prime}=\varnothing$

19: for tuple $k=1$ to $N_{b}$ do

20: Calculate TD-error: ${\delta_{k}}={r_{k}}+\gamma\cdot
V_{\theta^{V}}\left({s^{\prime}}_{k}\right)-V_{\theta^{V}}\left(s_{k}\right)$

21: Store tuple $k$ to $D^{\prime}$ if ${\delta_{k}}>0$

22: end for

23: Optimize actor by minimizing the loss:
$\frac{1}{{\left\|{D^{\prime}}\right\|}}{\sum\nolimits_{k^{\prime}}{\left\|{{a_{k^{\prime}}}-Act_{\theta^{A}}\left({s_{k^{\prime}}}\right)}\right\|}^{2}}$

24: Optimize critic by minimizing the loss of
$\frac{1}{{{N_{b}}}}{\sum\nolimits_{k}{\left\|{{\delta_{k}}}\right\|}^{2}}$

25: end for

26: end for

As mentioned in Section II-A, the followers are homogeneous. Thus, we adopt
the parameter sharing (PS) approach [25] and allow all the followers to share
the parameters of a single policy. This shared policy can be optimized more
efficiently with the data of experiences collected by all the followers. We
note that the PS-CACER algorithm follows the centralized-learning and
decentralized-execution fashion [26]. In other words, during the learning
phase, experiences obtained by all the followers simultaneously are stored
into a shared experience replay memory, and the shared policy in terms of deep
neural networks are trained with the stored experiences centrally
(centralized-learning). Nevertheless, during the execution phase, each
follower independently selects and executes its action by following the shared
policy based on its perceived state (decentralized-execution).

###  III-B Network Architecture

Inspired by [27] and [28], we design a customized network module based on
convolutional neural networks and attention mechanisms, which enables our
algorithm to adapt to the changes in the number of followers.

As illustrated in Figure 2, the Embedding module is composed of two
convolutional layers (Conv), two Squeeze-and-Excitation (SE) blocks [28], a
transpose layer (Transpose), a max-pooling layer (MaxPooling), and a flatten
layer (Flatten). First, the variable-length input (i.e., the ego-follower and
other followers joint state ${s^{o}}$) is passed by two convolutional layers
(i.e., Conv1 and Conv2) successively. The filter size of Conv1 is equal to the
length of $s_{j}^{o}$ , while the filter size of Conv2 is equal to the number
of filters of Conv1. With this specialized structure, each feature map
extracted by convolutional layers only depends on one follower (one of the
other followers). Besides, benefitting from the properties of CNNs, i.e.,
weight sharing and shift invariance [28], the sequencing operations according
to the distances between the leader and the followers as in [17, 18] are no
longer required. In other words, the output are invariant to the indexing of
the followers.

Figure 2: Embedding module. Each convolutional layer (Conv) is featured by the
number of filters, filter size, stride, and activation mode. Other layers are
represented by their type/name and output size. Figure 3: Network
architecture. Each fully-connected layer (Dense) is featured by its type,
number of neurons, and activation mode. Other layers are represented by their
type/name and output size.

After each convolutional layer, a SE block is added to improve the capacity
for feature extraction of networks. The SE block is a channel attention
mechanism. Each SE block uses a global average pooling layer
(GlobalAvgPooling) in the squeeze phase and two fully-connected layers (Dense)
with different activation functions in excitation phase, followed by a
channel-wise scaling operation (Multiply). By explicitly modelling the
relationship between channels, this architectural unit can enhance its
selectivety to the more informative features while suppressing the less useful
ones, boosting the representational power of the network [28]. After the
second Multiply layer, a max-pooling layer is added to create a fixed-size
output, independent of the input order and length. Lastly, the output of the
max-pooling layer is flattened by a flatten layer to create a fixed-length
embedding vector. In this way, the variable-length input can be embedded into
a fixed-length embedding vector.

We term our scheme of the Embedding module as SE-MaxPooling (SEMP). The
designed SEMP embedding module can handle inputs with arbitrary length and its
output is independent of the input order. We emphasize that this module is
plug-n-play, which means the SEMP embedding module can be easily integrated
with other network architectures, as well as any reinforcement learning
algorithm.

The entire architecture of the actor network and the critic network is
depicted in Figure 3. We process the ego-follower and leader joint state
${s^{e}}$ and the ego-follower and other followers joint state ${s^{o}}$
separately [17]. Specifically, a dense layer with ReLU activation function is
used to extract the features of ${s^{e}}$, and the Embedding module is used to
encode the ${s^{o}}$ into a fixed-length vector. After that, the two outputs
above are merged and then subsequently fed into 3 dense layers with different
activation functions. Note that both the actor and the critic use the same
network architecture up to the last output layer. The output layer of the
critic uses a linear activation function, while the output payer of the actor
uses a hyperbolic tangent (tanh) activation function.

##  IV Simulation Results

In this section, we evaluate the proposed PS-CACER algorithm with the SEMP
scheme in both numerical simulation and semi-physical simulation.

###  IV-A Simulation Setup

In this paper, the PS-CACER algorithm was trained with a total of 30000
episodes ($M=$ 30000), in which each episode had a maximal number of 60 time
steps ($N_{s}=$ 60). All network parameters were updated with 64 batch size
($N_{b}=$ 64), Adam optimizer with 0.001 (for actor) and 0.0001 (for critic)
learning rates. The exploration parameter was annealed exponential from 0.5 to
0.05 over a period of 2000 episodes, and then fixed at 0.05 thereafter. The
empirical values of essential parameters are listed in Table I [15, 16].

TABLE I: Parameter Settings Name |  Value |  Name |  Value
---|---|---|---
d1 |  40 |  d2 |  65
$\omega$ |  0.05 |  m |  2
${\alpha_{g}}$ |  9.8 |  ${r_{{\rm{bd}}}}$ |  $\frac{\pi}{6}$
${v_{\max}}$ |  18 |  ${v_{\min}}$ |  12
$N$ |  100000 |  $\gamma$ |  0.95

###  IV-B Numerical Simulation

In order to verify the effectiveness of the proposed SEMP scheme, three
existing state-of-the-art approaches, LSTM [17], SA (social attentive pooling)
[29], and CNNMP (convolutional neural networks with max pooling)[27], were
selected as the benchmarks. To make a fair comparison, we kept their network
architectures the same except for the embedding module. Besides, the above
solutions also used the same learning algorithm, i.e., the PS-CACER algorithm,
to update their network parameters. We note that, i) the number of LSTM cells
was set to 64, in accordance with the length of the fixed-length vector
encoded by SEMP. ii) the network and its parameters used by SA were the same
as [29]. iii) The network structure of CNNMP was the same as SEMP, except that
SEMP merged a channel attention mechanism by adding two SE blocks.

In the training phase, we evaluated the performance of the proposed algorithm
using the average reward $G_{{\rm{Avg}}}$ obtained by one agent within a
certain number of $n$e episodes:

${G_{{\rm{Avg}}}}=\frac{1}{{n{\kern
1.0pt}{N_{e}}{N_{s}}}}\sum\limits_{i=1}^{n}{\sum\limits_{p=1}^{{N_{e}}}{\sum\limits_{t=1}^{{N_{s}}}{{}^{i}r_{t}^{p}}}},$
(9)

where ${}^{i}r_{t}^{p}$ is the immediate reward obtained by agent $i$ at the
time step $t$ of the episode $p$ according to (8). Without loss of generality,
we set $N_{e}=$ 100\. At each episode, the number of followers $n$ was
randomly selected from 3 to 10. At each time step, the initial state of the
leader and followers were initialized randomly, and the steer commands of the
leader were generated randomly from the action spaces defined in Section
II-C2.

Figure 4: Learning curves of LSTM, SA, CNNMP, and SEMP. TABLE II: Comparison Results of LSTM, SA, CNNMP, and SEMP Method | Average Reward [Avg / Var] | Collision Rate (%) [Avg / Var]
---|---|---
$n$ = 4 |  $n$ = 6 |  $n$ = 8 |  $n$ = 10 |  $n$ = 4 |  $n$ = 6 |  $n$ = 8 |  $n$ = 10
LSTM [17] | -70.99/330.43 | -76.07/292.45 | -81.89/97.87 | -92.70/152.85 | 0.179/0.019 | 0.300/0.058 | 0.325/0.026 | 0.429/0.029
SA [29] | -49.25/97.86 | -59.01/104.62 | -67.20/101.59 | -79.72/155.76 | 0.282/0.120 | 0.272/0.016 | 0.375/0.025 | 0.554/0.038
CNNMP [27] | -48.88/127.72 | -57.51/133.02 | -63.39/112.21 | -72.98/122.63 | 0.215/0.048 | 0.282/0.026 | 0.294/0.025 | 0.404/0.036
SEMP (Ours) | -36.95/58.32 | -45.65/79.99 | -51.91/60.34 | -60.54/49.67 | 0.140/0.011 | 0.192/0.016 | 0.285/0.021 | 0.317/0.021

The learning curves of the four above approaches are shown in Figure 4. As can
be seen, the LSTM scheme has the worst performance in terms of both learning
efficiency and the obtained average reward. One reasonable explanation is that
LSTM has too many parameters to optimize. Besides, the learning curves of SA
and CNNMP are almost overlapped, indicating that the two schemes have similar
performance. Additionally, the average reward obtained by SEMP grows as fast
as CNNMP in the early stage. However, with the increase of training episodes,
the average reward of SEMP still gradually increases, and finally obtains a
higher reward. This comparison results validate that the attention mechanism
with SE blocks used by SEMP is effective. Overall, SEMP finally obtains the
highest reward with a similar learning efficiency. This result demonstrates
that our SEMP scheme has advantage over the three benchmarks.

After training, we tested the control policy learned by the PS-CACER algorithm
with the SEMP scheme (hereafter referred to as the learned SEMP policy for
short) in a flocking task with eight followers lasting for 200 seconds.
Initially, the leader was followed by four followers (Follower #1 $\sim$ #4 in
Figure 5). Since the time step 100, the other four followers (Follower #5
$\sim$ #8 in Figure 5) joined the flock simultaneously. The testing results
are visualized in Figure 5. The generated trajectories illustrate that the
eight followers can successfully keep up with the leader while avoiding
collisions with each other. This result indicates that the learned SEMP policy
can handle a variable number of followers.

Figure 5: The visualized results of the eight followers following the leader
using the learned SEMP policy, i.e., trajectories (left), the distances
($\rho$) to the leader (top-right), and the minimum distance ($mindis$) among
the followers (bottom-right).

To further quantitatively evaluate the performance of the learned SEMP policy,
we defined collision rate as a metric aside from the average reward. We
considered that a collision would happen if the distance between two followers
was less than a certain value, e.g., 2 meters. Thus, the collision rate meant
the percentage of being too close (i.e., less than 2 meters) among the
followers during the testing episode. We compared the learned SEMP policy with
three baselines by averaging the results of 200 episode. The parameter
settings are the same as the training phase, but the maximum time step
($N_{s}$) for each episode was set to 180. Table II and Figure 6 compares the
results of the four methods.

As can be seen from Table II, the learned policies with one set of network
parameters can adapt to a dynamic squad of UAVs, without the need of
retraining in a new environment with different number of followers. When the
number of followers increases, the average reward decreases and the collision
rate increases. The reason is that increasing the number of followers leads to
a higher probability of collision. However, compared with the other three
baselines, our SEMP method enables the follower to obtain the largest average
reward, the lowest collision rate, and the lowest variance, regardless of the
number of followers. This result shows that the SEMP method outperforms the
three existing state-of-the-art methods.

Figure 6: Box-and-whisker plot for comparison results of LSTM, SA, CNNMP, and
SEMP. Left: Average reward. Right: Collision rate.

###  IV-C Semi-Physical Simulation

In addition to the numerical simulation, we also tested the generalization
performance of the learned control policy by conducting a hardware-in-loop
(HIL) experiment in a high fidelity semi-physical simulation system [15],
[16]. Instead of the kinematic model of fixed-wing UAVs with stochastic
disturbances used in the training phase, we selected the professional tool
X-Plane 10111https://www.x-plane.com/manuals/desktop/ as the flight simulator
for the testing. X-Plane 10 can simulate complex environmental conditions,
such as weather changes and wind disturbances.

In this experiment, the five followers used the learned SEMP policy to flock
with the leader lasting for 200 time steps (seconds). At each time step, the
leader selected its roll-angle setpoint according to the pre-planned paths and
changed its velocity randomly, then broadcasted its state information to all
the followers. Each follower used the actor network with the trained
parameters to determine its steering commands according to its own state, the
leader’s state, and the other followers’ states.

Figure 7: The trajectory results in the semi-physical simulation. Note that
the three representative snapshots displayed on the top are captured from the
ground control station. The top-left and top-right snapshots show a third
person perspective of the UAV flock, while the top-middle snapshot is top
view. Figure 8: The distances between the leader and five followers ($\rho$)
and the minimum distance among five followers ($mindis$) in the semi-physical
simulation.

The trajectory results are depicted in Figure 7; the distances between the
leader and the five followers ($\rho$) as well as the minimum distance among
the five followers ($mindis$) are shown in Figure 8. As can be seen, the
distances between the leader and the followers were maintained around
75$\,\rm{m}$ most of the time (the average distance is 74.44$\,\rm{m}$). This
means that the five followers were able to keep up with the leader steadily,
even if the leader changed its heading sharply. More importantly, the minimum
value of $mindis$ is larger than 2$\,\rm{m}$, which means that there were no
collisions among the followers during this experiment. The above results
demonstrate that the proposed SEMP method enables the followers to avoid
collisions between each other while flocking with the leader. We note that the
control policy employed by the five followers in this semi-physical simulation
is the SEMP policy learned from the previous numerical simulation, without any
parameter finetuning. This result demonstrates that the learned SEMP policy
can be directly generalized to new situations.

##  V Conclusion

In this paper, we have proposed a deep reinforcement learning algorithm to
solve the leader-follower flocking control and collision avoidance problem for
a dynamic squad of fixed-wing UAVs. Specifically, we have designed a
decentralized DRL-based framework where the collision avoidance policy is
integrated into the flocking controller for a variable number of UAVs. Then,
we have proposed the PS-CACER algorithm for training multiple agents, in which
a customized embedding module, SEMP, is integrated to handle the variable-
length input for deep neural networks. With this module, the learned policies
are adaptive to the changes of the number of UAVs without the need of
retraining. Numerical simulation results have demonstrated that the proposed
method outperforms three existing state-of-the-art methods, i.e., LSTM, SA,
and CNNMP. Finally, the learned policy can be directly transferred to a semi-
physical simulation without any parameter finetuning. In the future, we will
further evaluate our method with fixed-wing UAVs in real-world environments.

## References

  * [1] Y. Tang, Y. Hu, J. Cui, F. Liao, M. Lao, F. Lin, and R. S. Teo, “Vision-aided multi-UAV autonomous flocking in GPS-denied environment,” IEEE Transactions on Industrial Electronics, vol. 66, no. 1, pp. 616-626, Jan. 2019.
  * [2] C. H. Liu, X. Ma, X. Gao and J. Tang, “Distributed Energy-Efficient Multi-UAV Navigation for Long-Term Communication Coverage by Deep Reinforcement Learning,” in IEEE Transactions on Mobile Computing, vol. 19, no. 6, pp. 1274-1285, 1 June 2020.
  * [3] P. C. Lusk, X. Cai, S. Wadhwania, A. Paris, K. Fathian and J. P. How, “A Distributed Pipeline for Scalable, Deconflicted Formation Flying,” in IEEE Robotics and Automation Letters, vol. 5, no. 4, pp. 5213-5220, Oct. 2020.
  * [4] S. Chung, A. A. Paranjape, P. Dames, S. Shen, and V. Kumar, “A Survey on Aerial Swarm Robotics,” IEEE Transactions on Robotics, vol. 34, no. 4, pp. 837-855, Aug. 2018.
  * [5] Y. Kuriki, and T. Namerikawa, “Formation control with collision avoidance for a multi-UAV system using decentralized MPC and consensus-based control,” in European Control Conference (ECC), Linz, 2015, pp. 3079-3084.
  * [6] O. Saif, I. Fantoni, and A. Zavala-Río, “Distributed integral control of multiple UAVs: precise flocking and navigation,” IET Control Theory & Applications, vol. 13, no. 13, pp. 2008-2017, 2019.
  * [7] R.S. Sutton, and A.G. Barto, Reinforcement learning: An introduction. Cambridge, USA: MIT Press, 1998.
  * [8] C. Yan, and X. Xiang, “A path planning algorithm for UAV based on improved Q-learning,” in International Conference on Robotics and Automation Sciences (ICRAS), Wuhan, 2018, pp. 1-5.
  * [9] S. Hung, S. N. Givigi, and A. Noureldin, “A Dyna-Q ($\lambda{}$) approach to flocking with fixed-wing UAVs in a stochastic environment,” in IEEE International Conference on Systems, Man, and Cybernetics, Kowloon, 2015, pp. 1918-1923.
  * [10] S. Hung, and S. N. Givigi, “A Q-learning approach to flocking with UAVs in a stochastic environment,” IEEE Transactions on Cybernetics, vol. 47, no. 1, pp. 186-197, Jan. 2017.
  * [11] S. C. Speck, and D. J. Bucci, “Distributed UAV swarm formation control via object-focused, multi-objective SARSA,” in Annual American Control Conference (ACC), Milwaukee, WI, 2018, pp. 6596-6601.
  * [12] C. Yan, X. Xiang, and C. Wang, “Towards Real-Time Path Planning through Deep Reinforcement Learning for a UAV in Dynamic Environments,” Journal of Intelligent & Robotic Systems, vol. 98, no. 2, pp. 297-309, 2020.
  * [13] L. Tai, G. Paolo, and M. Liu, “Virtual-to-real deep reinforcement learning: Continuous control of mobile robots for mapless navigation,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Vancouver, BC, 2017, pp. 31-36.
  * [14] A. Singla, S. Padakandla, and S. Bhatnagar, “Memory-based deep reinforcement learning for obstacle avoidance in UAV with limited environment knowledge,” in IEEE Transactions on Intelligent Transportation Systems, vol. 22, no. 1, pp. 107-118, Jan. 2021.
  * [15] C. Wang, C. Yan, X. Xiang, and H. Zhou, “A continuous actor-critic reinforcement learning approach to flocking with fixed-wing UAVs,” in Asian Conference on Machine Learning (ACML), Nagoya, Japan, 2019, pp. 64-79.
  * [16] C. Yan, X. Xiang, and C. Wang, “Fixed-Wing UAVs flocking in continuous spaces: A Deep reinforcement learning approach,” Robotics and Autonomous Systems, vol. 131, p. 103594, Sep. 2020.
  * [17] Z. Sui, Z. Pu, J. Yi, and T. Xiong, “Formation control with collision avoidance through deep reinforcement learning,” in International Joint Conference on Neural Networks (IJCNN), Budapest, Hungary, 2019, pp. 1-8.
  * [18] M. Everett, Y. F. Chen, and J. P. How, “Motion planning among dynamic, decision-making agents with deep reinforcement learning,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Madrid, 2018, pp. 3052-3059.
  * [19] S. Hochreiter, and J. Schmidhuber, “Long short-term memory,” Neural Computation, vol. 9, no. 8, pp. 1735-1780, 1997.
  * [20] C. Wang, J. Wang, and X. Zhang, “A deep reinforcement learning approach to flocking and navigation of UAVs in large-scale complex environments,” in IEEE Global Conference on Signal and Information Processing (GlobalSIP), Anaheim, CA, USA, 2018, pp. 1228-1232.
  * [21] C. Kownacki, D. Odziej, J. Bordeneuve-Guibé, A. Drouin, and C. Roos, “Flocking algorithm for fixed-wing unmanned aerial vehicles,” in Advances in Aerospace Guidance Navigation and Control, Springer, pp. 415-431, 2015.
  * [22] T. Lillicrap, J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” in International Conference on Learning Representations (ICLR), IEEE, 2016.
  * [23] X. B. Peng, G. Berseth, K. Yin, and M. Van De Panne, “DeepLoco: Dynamic locomotion skills using hierarchical deep reinforcement learning,” ACM Transactions on Graphics, vol. 36, no. 4, pp. 1-13, 2017.
  * [24] H. van Hasselt, “Reinforcement learning in continuous state and action spaces,” In Reinforcement Learning, Springer, 2012, pp. 207-251.
  * [25] K. G. Jayesh, E. Maxim, and K. Mykel, “Cooperative multi-agent control using deep reinforcement learning”, in International Conference on Autonomous Agents and Multiagent Systems (AAMAS), Cham, 2017, pp. 66-83.
  * [26] R. Lowe, Y. Wu, A. Tamar, J. Harb, P. Abbeel, and I. Mordatch, “Multi-agent actor-critic for mixed cooperative-competitive environments,” in International Conference on Neural Information Processing Systems (NIPS), 2017, pp. 6382-6393.
  * [27] C. Hoel, K. Wolff, and L. Laine, “Automated speed and lane change decision making using deep reinforcement learning,” in International Conference on Intelligent Transportation Systems (ITSC), Maui, HI, 2018, pp. 2148-2155.
  * [28] J. Hu, L. Shen, and G. Sun, “Squeeze-and-excitation networks,” in IEEE Conference on Computer Vision and Pattern Recognition (CVPR), Salt Lake City, USA, 2018, pp. 7132-7141.
  * [29] C. Chen, Y. Liu, S. Kreiss, and A. Alahi, “Crowd-robot interaction: Crowd-aware robot navigation with attention-based deep reinforcement learning,” in International Conference on Robotics and Automation (ICRA), Montreal, Canada, 2019, pp. 6015-6022.


# Spatial Confinement Affects the Heterogeneity and Interactions between
Shoaling Fish

Gabriel Kuntz  Seattle University, Department of Physics, Seattle WA, 98122,
USA  Junxiang Huang  Northeastern University, Department of Physics, Boston,
MA, 02115, USA  Mitchell Rask  Seattle University, Department of Physics,
Seattle WA, 98122, USA  Alex Lindgren-Ruby  Seattle University, Department of
Physics, Seattle WA, 98122, USA  Jacob Y. Shinsato  Seattle University,
Department of Physics, Seattle WA, 98122, USA  Dapeng Bi  Northeastern
University, Department of Physics, Boston, MA, 02115, USA  A. Pasha Tabatabai
Seattle University, Department of Physics, Seattle WA, 98122, USA  California
Polytechnic State University, Physics Department, San Luis Obispo CA, 93410,
USA<EMAIL_ADDRESS>
###### Abstract

Living objects are able to consume chemical energy and process information
independently from others. However, living objects can coordinate to form
ordered groups such as schools of fish. This work considers these complex
groups as living materials and presents imaging-based experiments of
laboratory schools of fish to understand how this non-equilibrium activity
affects the mechanical properties of a group. We use spatial confinement to
control the motion and structure of fish within quasi-2D shoals of fish. Using
image analysis techniques, we make quantitative observations of the
structures, their spatial heterogeneity, and their temporal fluctuations.
Furthermore, we utilize Monte Carlo simulations to replicate the
experimentally observed area distribution patterns which provide insight into
the effective interactions between fish and confirm the presence of a
confinement-based behavioral preference transition. In addition, unlike in
short-range interacting systems, here structural heterogeneity and dynamic
activities are positively correlated as a result of complex interplay between
spatial arrangement and behavioral dynamics in fish collectives.

## Introduction

Nature provides fantastic examples of complex collective behaviors on many
length scales in order to accomplish certain tasks. For example, cells within
tissues coordinate to successfully close wounds [1, 2], ants build structures
to overcome obstacles [3, 4], and fish form cohesive groups to improve
computations about their environment [5, 6]. In each of these examples, the
interactions between individuals lead to function on a larger scale.

Understanding the details of the interactions between individuals within these
complex groups is an active area of research [7, 8]. Previously, it has been
shown that metric interactions, where constituents within a certain distance
interact, qualitatively capture the collective behaviors seen in flocking [9].
However, closer inspections in a variety of species suggest that the true
interactions are more likely visual, topological [10, 11, 12], or more
complicated [13]. The pursuit of understanding these interactions is valuable
to understanding fundamental problems in complex systems.

However, this connection between constituent interactions and bulk behavior
parallels the language used to describe and design materials. Examples of this
include tuning the interaction strength between colloids to influence
colloidal gel rheology [14, 15, 16, 17], the interaction specificity within
DNA hybridized colloids and material structure [18, 19, 20], and the relative
physical parameters within models of epithelial tissue and tissue fluidity
[21, 22]. Therefore, this search for the relationship between interactions and
bulk behavior is also critical for defining collective systems as living
materials in-and-of themselves. In particular, these living materials are
fundamentally non-equilibrium due to the local consumption of energy by each
entity, and definitions of the mechanical [23, 24, 25] or thermodynamic [26,
27] properties of these systems will help determine the material possibilities
of these types of systems.

In this paper, we aim to understand the mechanical properties of a quasi-2D
living material - groups of aquarium fish within the lab; for simplicity, we
confine fish to thin volumes of water. We make quantitative observations of
groups of swimming fish using image analysis techniques to identify fish
positions and trajectories. We track fluctuations of structures at the local
and group level as self-generated deformations. We find that we can control
the motile behaviors of these fish by varying the level of spatial confinement
and that this change in individual motion is correlated with a changing
heterogeneity of the group. In addition, we employ Monte Carlo simulations to
recreate the complex dynamics of fish interactions and examine the effects of
varying spatial constraints on the group’s mechanical properties. By using
simulations to replicate experimentally measured distributions, we infer
effective interactions between individuals. Through a wide variety of metrics,
we observe a behavioral transition as a function of spatial confinement which
correlates with changes in structural heterogeneity. In contrast to colloidal
gels involving short-range interactions, the simulations uncover a positive
correlation between structural heterogeneity and dynamic activities.

## Results

Cardinal tetra fish are imaged moving freely within quasi-2D cylindrical
arenas with a depth of 1.5cm$\pm$0.1cm (Methods). Fish positions $\vec{r}$ are
identified using the open-source software Trex [28] which identifies fish body
cross-sections through image contrast. We also use this software to connect
fish positions in time to build trajectories $\vec{r}(t)$ and calculate
instantaneous velocities $\vec{v}(t)$ for each fish (Methods). Since we now
have the projected area of each fish, we use the mid-line length to
characterize fish size. We find fish have a length of $L=2.0\pm 0.2$ (mean
$\pm$ stdev) (Figure 1).

Figure 1: (a) Image of 50 fish within arena with radius $R=44.25$cm. Arena
boundary outlined for clarity. (b) Probability distribution of fish mid-line
length ($L$) for 50 fish in (a).

We record videos of 25 fish within arenas of different radii ($R$) to
investigate the effects of confinement (Figure 2a). Altering the radius $R$
adjusts the global area density $\rho_{g}=N/(\pi R^{2})$. However, it is
crucial to note that fish can exhibit significant local density fluctuations
since they do not uniformly occupy the entire space [29, 30]. We observe that
arena size influences the probability distribution of speeds ($v$), where
decreasing $R$ biases the distributions towards slower speeds (Figure 2b).
This broad distribution of speeds is consistent with the stop-start motion
associated with fish motility [31]. We fit each speed probability distribution
to a modified Rayleigh distribution.

$P(v)=\frac{v+b}{a}e^{-(v+b)^{2}/(2a)}$ (1)

where $a$ and $b$ are fitting parameters associated with the width and shift
of the distribution accordingly. We find that both of these parameters
increase with $R$ (Figure 2c). This functional form was chosen to resemble the
2D Maxwell-Boltzmann distribution in an attempt to make parallels between the
motion of molecules at thermodynamic equilibrium and the motion of fish out-
of-equilibrium; the connection is elaborated on further in the Discussion
Section.

While Figure 2b indicates that the magnitude of motion is affected by $R$, it
does not describe the persistence of motion. As such, we calculate the mean
squared displacement (MSD)

$\mathrm{MSD}(\tau)=\left<(\vec{r}(t+\tau)-\vec{r}(t))^{2}\right>$ (2)

by comparing positions $\vec{r}(t)$ as a function of elapsed time ($\tau$) for
each fish; we average the MSD from all fish within an experiment to generate a
single ensemble-averaged MSD (Figure 2d, Methods). For large $\tau$, the MSD
turnover and plateau are set by the finite size of the arena. For small
$\tau$, we observe power-law scaling (i.e. $\mathrm{MSD}\sim\tau^{\alpha}$)
where $\alpha$ characterizes the type of motion. We find that $\alpha$ depends
on confinement, demonstrating that fish motion is super-diffusive
($1<\alpha<2$) and approaches ballistic motion ($\alpha\xrightarrow{}2$) as
containers get larger (Figure 2e). The MSD of the shoal’s center of mass has
similar arena-size dependent values of $\alpha$ (Figure 2e, Supplemental
Figure 1).

Figure 2: Confinement affects motion. (a) Images of arenas with radii $R=8$cm,
11cm, 14cm, 34.25cm & 44.25cm containing 25 fish. Arena boundaries were added
for clarity. The red arrow denotes radius $R$. The scale bar is 10cm. (b)
Example probability distributions of fish speed ($v$) for 25 fish in the
differently sized arenas (markers) and best fit Gaussian functions (lines).
(c) Fitting parameters $a$ (left-blue) and $b$ (right-red) for the lines in
(b). (d) Mean squared displacement (MSD) as a function of elapsed time
($\tau$) for fish in five different arenas. (e) The short-time power-law slope
of MSD ($\alpha$) as a function of radii ($R$). Experimental duplicates of
$\alpha$ are plotted in grey ($N=3$ for each), with the mean plotted in black.
$\alpha$ for the shoal center of mass motion (blue). Error bars are one
standard deviation. The plotted color darkens as $R$ increases in (b) and (d).

We next asked how these differences in motion affected the organization of
fish in groups. To investigate the effects of confinement on structural and
material properties of the shoal we use fish positions $\vec{r}(t)$ to
calculate the time-varying convex hull (Figure 3a, Methods) which we use to
define the overall geometric size of the group. The area of the convex hull
$A_{H}$ fluctuates considerably over time, signifying that the group is
exploring different fractions of the space (Figure 3b). While $A_{H}$
fluctuates in time, the time-averaged fraction of space occupied by the shoal
$\left<A_{H}\right>/\pi R^{2}$ is consistent across different confinements
(Figure 3c); on average, the shoal will fill the space to an equal extent
regardless of the amount of space it has available to it (Supplemental Figure
2).

The size of a shoal characterized by $A_{H}$ is prone to bias by fish that do
not move with the group. Therefore, we define a local measurement of the space
occupied by individual fish by calculating the Voronoi tessellation using the
fish positions $\vec{r}$. The Voronoi cell for a particular fish is the space
that is most proximal to a fish; this acts as an amorphous unit cell. We
restrict our structural analysis at each frame to fish that have Voronoi cells
completely enclosed within the convex hull (Figure 4a); these ‘internal’ fish
have Voronoi areas which are both closed and do not drastically change with
small neighbor movements. We calculate the areas of each of these Voronoi
cells $A$ for individual fish and find that they fluctuate through time as
well. In Figure 4b, we show an example Voronoi area that fluctuates by an
order of magnitude in area over the plotted observation window. We also note
that the data point frequency is not constant over the observation window; no
Voronoi area is calculated for this fish if it fails to be an ‘internal’ fish
or if we cannot uniquely identify all fish in a particular frame.

Figure 3: Shoal area fluctuates in time. (a) Convex hull of 25 fish in
$R=44.25$cm arena. Fish positions are grey markers. Convex hull is defined by
the black dashed line. The image width is $2R$ with the arena boundary
outlined for clarity. (b) Example of convex hull area ($A_{H}$) normalized to
arena area ($\pi R^{2}$) as a function of time for 25 fish in R=34.25cm arena.
(c) Time average convex hull area $\left<A_{H}\right>$ normalized by arena
area for arenas of different radii. Error bars are one standard deviation.

The number of these internal fish $N_{int}$ varies widely over time (Figure 4c
inset) and each internal fish has an average of approximately six topological
neighbors, defined as neighbors that share a Voronoi edge (Supplemental Figure
3). Therefore, to further understand the relationship between local
fluctuations in $A$ and shoal level fluctuations, we estimate the time-varying
net size of all internal fish $\sum A/N_{int}$ (Figure 4c). We note that the
time courses of the $A_{H}$ (Figure 3b), $A$ (Figure 4b), and $\sum A/N_{int}$
(Figure 4c) all come from the same experiment and are plotted over the same
range of time. The relative size of fluctuations within $A_{H}$ and $\sum
A/N_{int}$ are qualitatively similar. However, these are not equivalent to the
fluctuations of $A$ over the same time period. Therefore, the group and the
individual area fluctuations are not mirrored, and fish areas do not all
uniformly expand or shrink collectively.

Figure 4: Confinement affects local fish packing. (a) Voronoi tessellation
(polygons) of 25 fish (grey markers) in $R=44.25$cm arena. Internal fish (red
polygons) are a subset of all fish and have all vertices within the convex
hull (black dashed line). The image width is $2R$ with the arena boundary
outlined for clarity. (b) Example time evolution of a single fish Voronoi area
($A$) normalized by arena area ($\pi R^{2}$) in $R=34.25$cm arena. (c) Sum of
internal areas ($\Sigma A$) normalized by number of internal fish ($N_{int}$)
as a function of time for data in (b). (c-inset) Number of internal fish
($N_{int}$) as a function of time for data in (b). (d) Probability
distributions ($P(A)$) of internal areas ($A$) for 25 fish for different arena
sizes. The vertical dashed line indicates the peak of a distribution $A_{0}$.
(e) $A_{0}$ as a function of arena radii ($R$). Bold circles are averages
across different fish groups, small data points are individual experiments and
error bars are one standard deviation. (f) Scaled probability distributions
from (d). Colors in (d) and (f) darken with increasing $R$.

Since shoals occupy larger spaces in larger arenas (Figure 4c), the Voronoi
region associated with each fish must also vary with $R$. Indeed, we see this
dependence in the probability distributions of $A$ for all internal fish
within an experiment, where increasing $R$ biases the distribution to larger
areas $A$ (Figure 4d). We define the modal area $A_{0}$ which increases with
arena size, however, it does not increase $\sim R^{2}$ as would be expected
for a 2D gas (Figure 4e).

By comparing the normalized probability distributions ($P(A)*A_{0}$ vs
$A/A_{0}$), we show that the fluctuations observed are statistically similar
between small arena sizes yet vary significantly for larger arenas (Figure
4f); this is consistent with the arena size dependence of MSD scaling in
Figure 2. We also show that these distributions are not equivalent to
distributions made from randomly generated points, consistent with the fact
that fish are not randomly occupying space (Supplemental Figure 4).

Upon inspection, the probability distribution of observed internal Voronoi
areas is not symmetric about the mode $A_{0}$ (Figure 5a). Here, we take an
approach that is similar to the computational modeling of cells in tissues via
the Vertex and Self-Propelled Voronoi models where deviations of a cell from a
modal area are associated with an energy cost for that cell [21, 22]. To
understand the underlying dynamics that result in the asymmetric distributions
in Figure 4d, we fit two separate parabolic functions to each distribution
such that

$P(A_{0})-P(A)\sim\begin{cases}k_{c}(A-A_{0})^{2},&\text{if }A<A_{0}\\\
k_{e}(A-A_{0})^{2},&\text{if }A>A_{0}\end{cases}$ (3)

with constants $k_{c}$ and $k_{e}$ associated with compression and expansion,
respectively, of Voronoi areas away from the modal area $A_{0}$ (Figure 5a,
Methods). The ratio $k_{e}/k_{c}$ of these constants is a signature of the
effective interactions between fish, which varies with $R$ (Figure 5b).

Figure 5: Confinement affects area distribution asymmetry.(a) Example
probability distribution of internal fish voronoi areas ($A$) for $R=14$cm
arena (143,963 instances). Parabolic fits to $A<A_{0}$ (red) and $A>A_{0}$
(cyan) associated with compression and expansion, respectively. (b) Ratio of
coefficients of expansion $k_{e}$ and compression $k_{c}$ found from fit in
(a) as a function of arena radius ($R$). Small markers are individual
experiments, and the large black markers are mean with one standard deviation
error bar.

## Simulations

Figure 6: The heterogeneity of internal Voronoi areas increases as the density
is decreased. (a-e) Experimental area $A$ distributions are fit to Gamma
distribution (red) for $R=8,11,14,34.25\mathrm{and\ }44.25$ cm radii arenas.
$f_{\Gamma}(k,\theta)$ denotes the PDF of Gamma distribution defined in
Equation (4). (d,e) At $R=34.25\mathrm{cm}$ and $44.25\mathrm{cm}$ the Gamma
distribution fit is done for $A<1.5A_{0}$ and a normalization factor is
applied to $f_{\Gamma}(k,\theta)$ to align the modes. Insets in (d,e) show
exponential tails (black dashed lines) with decay parameters of $0.0098$ and
$0.0038\mathrm{cm}^{-2}$, respectively. (f) The shape parameter $k$ as a
function of arena radii $R$. (f-inset) The scale parameter $\theta$ as a
function of arena radii $R$.

Previous literature has established that the area distribution of cells within
two-dimensional random Voronoi networks adheres to a Gamma distribution [32].
The probability density function (PDF) for such a distribution is
mathematically represented as:

$f_{\Gamma}(x)=\frac{x^{k-1}e^{-x/\theta}}{\theta^{k}\Gamma(k)},$ (4)

where $k$ and $\theta$ are the shape and scale parameters, respectively, while
$\Gamma(*)$ denotes the gamma function. For random Voronoi networks, a shape
parameter $k=3.63$ is reported to yield the best fit to the observed cell area
distribution [33]. In addition, for hard disks, the distribution of Voronoi
free area, which is the difference between the actual Voronoi cell area and
the minimum cell area at close packing, is well described by a Gamma
distribution with $k$ between $3$ to $4$ [34].

Intriguingly, despite the non-equilibrium and non-random nature of the fish
collective, the internal area distributions at smaller radii $R$ conform
closely to a Gamma distribution (Figure 6a-c). The shape parameter $k$, which
governs the distribution asymmetry and tail behavior, is greater than the non-
interacting limit ($k=3.6$). Large $k$ values represent a more symmetric bell-
shaped curve indicative of low heterogeneity and a tendency for cell areas to
aggregate around the mean. This pattern suggests that interactions between
fish within small arenas lead to a more homogeneous structural arrangement;
this reduces variability and allows each fish to navigate and occupy space
more effectively. As the arena radius $R$ increases, there is a clear increase
in the Voronoi area heterogeneity with two distinct signatures: a decrease in
the $k$ value and the emergence of an exponential tail in the distribution
(Figure 6d-f). At large arenas such as $R=34.25$cm or $44.25$cm, the internal
areas initially adhere to a Gamma distribution up to a threshold around
$1.5A_{0}$. Beyond this point, the distribution exhibits strong exponential
tails indicative of highly heterogeneous shoal structure and significant
probabilities of finding large cells and large local density fluctuations.
This phenomenon bears resemblance to the behavior observed in granular
aggregates with capillary interactions, while the fish aggregate is unique in
its ability to adapt the $k$ value, which is obtained by fitting the bulk of
the distribution to a Gamma distribution and is observed to change across a
broad range [35].

The structural heterogeneity and dynamics are frequently interlinked. In
colloidal gels formed through short-range attractive interactions, an increase
in interaction strength leads to an increase in structural heterogeneity and
dynamical arrest [36, 37]. However, our observations in fish collectives
present a contrasting scenario. Due to the long-range nature of their
interactions, both structural heterogeneity and dynamical activities escalate
with a greater radius $R$. In these expansive arenas, local densities
experience more pronounced fluctuations which result in varied behavioral
patterns: some fish form tightly-knit compact shoals while others
simultaneously and independently navigate the container.

Figure 7: Interaction strength and length scale are influenced by confinement.
(a-c) Sample comparisons between experimental internal area distributions $A$
and Monte Carlo simulation (red line) for radii $R=8$cm, 11cm, and 34.25cm.
(d) The preferred distance $r_{0}$ and (d-inset) parameter $\beta$ as a
function of arena radii $R$ from Monte Carlo simulations.

To elucidate the influence of the arena radius on fish interactions and the
interplay between heterogeneity and dynamical activities, we implement a Monte
Carlo model that simulates the fish positions with different arena radii and
$\beta=1/k_{B}T$. At each simulation step $n$, we randomly select a fish,
indexed as $i$, and propose its subsequent position as
$X_{i,n+1}=X_{i,n}+\sigma\delta X_{i,n}$, where $X_{i,n}$ represents the
position of fish $i$ at step $n$, $\sigma$ is the step size, and $\delta
X_{i,n}$ is Gaussian white noise with unit variance. The Morse potential is
adapted to emulate the complex dynamics of fish schooling by depicting the
attractive and repulsive inter-fish forces. The potential is expressed as

$U(r)=\left(1-e^{1-r/r_{0}}\right)^{2},$ (5)

where $r$ is the actual distance between two fish, and $r_{0}$ is the
preferred distance. A proposed move that results in an energy change $\Delta
E$ is accepted with probability $min(1,e^{-\beta\Delta E})$. We conducted the
simulation with a step size $\sigma=0.4$ over $10^{7}$ iterations, recording
the internal area metrics every $100$ step. As $\beta$ approaches zero, the
influence of interactions among fish vanishes, corresponding to the ideal gas
limit. Within this limit, the statistical properties of fish under different
confinements are essentially the same, except for a scaling factor determined
by the arena’s radius. For example, the internal Voronoi areas adhere to a
Gamma distribution, characterized by a constant shape parameter $k\sim 3.75\pm
0.1$ and a scale parameter $\theta(R)=0.029~{}R^{2}$, which leads to a
quadratic relationship between $\left<A\right>$ or $A_{0}$ and $R$. To find
out the parameters that best describe the experiments, we systematically sweep
the parameter space of $r_{0}$ ranging from $0.1$ to $20$ and $\beta$ ranging
from $0$ to $2$ for different simulations, and the resulting area
distributions are then statistically analyzed using the chi-squared method to
ascertain the optimal values for $r_{0}$ and $\beta$ corresponding to each
radius $R$. As shown in Figure 7, despite its simplified nature, the Monte
Carlo simulations successfully reproduce the experimentally observed area
distributions. This indicates that the core principles and rules embedded in
the Monte Carlo model are effective in capturing the essential dynamics of
fish interactions.

Interestingly, the preferred distance $r_{0}$ and $\beta$ display a non-
monotonic variation as a function of $R$, with a transition point located
between $R=14$cm and $34.25$cm. This non-monotonic variation suggests a
complex interplay between individual space requirements and the benefits of
social interactions. In small arenas, the high density compels fish to
maintain a small $r_{0}$. When more space is available, fish tend to increase
their preferred distance $r_{0}$ to avoid overcrowding and reduce stress, with
a decreasing $\beta$ indicative of increased activities. With sufficient
space, however, behavior changes and a decrease in $r_{0}$ imply a shift
toward preserving the advantages of schooling such as enhanced communication
and collective vigilance. Such a transition explains the dramatic decrease of
$A_{0}/R^{2}$ for the big radii in Figure 4e, in contrast to a constant ratio
in the non-interacting limit or $\beta=0$. The presence of a transition point
indicates a threshold at which the fish alter their spacing behavior, possibly
to balance the conflicting needs for individual space and group affiliation as
a strategic response to maximize the evolutionary benefits of schooling.

## Discussion

In this manuscript, we have treated a quasi-2D shoal of fish as a living
material, and we have demonstrated the ability to control the average motion
of individuals as well as the structures present within the group simply by
changing the extent of confinement while keeping the number of individuals
within the group constant.

In describing the distribution of speeds in Figure 2, we found that a Rayleigh
function (Equation 1) described our data well. This is reminiscent of a 2D
Maxwell-Boltzmann distribution and suggests that this experimental system can
teach us something about energy usage within non-equilibrium systems. To make
this analogy, we consider the traditional Maxwell-Boltzmann distribution which
describes the speed $v$ of molecules with mass $m$ of a gas at a given
temperature $T$. In this form, $k_{B}T$ is the thermal energy scale where
$k_{B}$ is the Boltzmann constant. If we consider $a=k_{B}T/m$ and that $b$ is
a fitting parameter for offset speed in Equation 1, then $a$ is analogous to
the amount of energy inputted. Traditionally this energy would be via thermal
fluctuations per particle, however, $a$ is not thermal in origin. Instead,
this energy input comes from the energy usage of the fish towards swimming and
is therefore a measure of non-equilibrium activity. This trend of changing
effective energy with confinement is also observed through the interaction
energy within Monte Carlo simulations, where $\beta(R)$ in Figure 7d is
similar to $1/a$ in Figure 2c.

Extending this non-equilibrium thermodynamic analogy further, the increase in
non-equilibrium activity $a$ with system size $R$ is correlated with a
decrease in global density $\rho_{g}$. This interpretation of Figure 2b,c
suggests that fish consume more energy while swimming in the larger arenas.
Therefore, our observed trend is similar to an ideal gas at constant pressure:
gas molecules must have more thermal energy to maintain a constant pressure if
there are fewer molecules.

Drawing on statistical mechanics, we make the analogy that the dynamics of a
material $\left<v^{2}\right>\sim 1/\beta$. Our simulations indeed show that an
increase in $\beta$ which represents a decrease in motion, is concurrent with
decreases in structural heterogeneity (larger shape parameter $k$). Therefore,
we again suggest that an increase in structural heterogeneity is correlated
with an increase in dynamics within fish shoals.

We observe a decoupling of the local fluctuations in fish areas with the group
level fluctuations in size, indicating that expansions and contractions of the
group are not homogeneously distributed amongst individuals (Figures 3b, 4b,
and 4c). We find that the distribution of internal Voronoi areas is asymmetric
about the mode (Figure 5) and is robust to the method of measurement
(Supplemental Figures 5, 6, and 7). However, the area distribution of non-
interacting cells within a two-dimensional random Voronoi networks follows a
Gamma distribution with $k=3.6$, which is equivalent to $k_{e}/k_{c}\approx
0.6$. As such, deviations from this non-interacting result in Figure 5b are a
direct result of changing interactions between fish as a result of
confinement.

We find that changes in structural and dynamic trends occur between $R=14$cm
and $R=34.25$cm such as fish motion (Figure 2e), modal Voronoi area (Figure
4e), distributions of Voronoi areas (Figure 4f), the shapes of these
distributions (Figures 5b, 6), and the inferred radii and energy of
interaction between fish (Figure 7). As such, we conclude that spatial
confinement is a robust method to control the dynamical and mechanical
properties of this non-equilibrium material.

## Methods

### Fish Care

Cardinal tetra (Paracheirodon axelrodi) fish are purchased from the Aquarium
Co-Op in Edmonds WA. Fish are housed in a $\sim 50$ gallon living aquarium at
a maximum density of 1 fish/gallon. Lights are set to a 12h-12h day-night
cycle; all experiments are done during the daytime setting. Water is kept at a
temperature between $77^{\circ}$ and $79^{\circ}$F. Water pH is monitored and
adjusted around 7.2-7.4 and ammonia, nitrate, and nitrite levels are monitored
using the API Freshwater Test Kit and kept at undetectable levels. Fish are
fed once daily. All fish housing and handling is approved by an Institutional
Animal Care and Use Committee.

### Experimental Details

For experimental observations, 25 fish are chosen at random from sets of 30-60
fish in all data except Figure 1 where 50 fish were used. The fish are
transferred to shallow cylindrical arenas made for experiments. The
temperature of the observation tank is kept between $77^{\circ}$ and
$79^{\circ}$F and the water used for this tank is directly taken from the
living aquarium. The water depth is kept at 1.5cm with a tolerance of $\pm
0.1$cm across the arena. This minimizes 3D fish crossings that affect our fish
identification algorithm. Arenas are made from either custom acrylic or white
PVC. Clear acrylic arenas are lined with white tape to match PVC. Arena sizes
are radii $R=8$cm, 11 cm, 14cm, 34.25cm, and 44.25cm.

Shallow arenas are submerged in a large $\sim 200$ gallon water bath with
active heating and water circulation which acts as a thermal reservoir but
avoids generating any flow in the observation arenas. Water within the shallow
container containing the fish is static except when perturbed by the fish
within the observation arenas. Fish are left undisturbed in the imaging arena
for a minimum of one hour before imaging for acclimatization.

Fish are back-lit by submerged broadband visible light. A diffusive acrylic
layer separates the light source and the imaging aquarium base which helps to
homogeneously illuminate the field of view. Room lights are turned off during
acclimatization and experiments. Videos are recorded from overhead with a
Pixelink PL-D7620 machine vision camera at 10 frames per second for up to 60
minutes.

### Image analysis

Videos are taken with lighting optimized to ensure shadows, bubbles, and any
other visual noise are minimized before using the open-source software Trex
[28] to threshold the videos and determine position, velocity, and orientation
for individual fish. The pixel-to-centimeter conversion is found by taking a
photograph of a ruler at the bottom of the arena after each experiment without
disturbing the camera setup.

### MSD Analysis

The mean-squared-displacement MSD is calculated for each fish and averaged for
all fish in an experiment. When we lose continuity in fish trajectories due to
tracking errors, we ensure that the MSD is only calculated for consecutively
tracked frames. The average consecutive track lengths for any given fish range
from 59s in the largest $R=44.25$cm arena to 18.8s in the smallest $R=8$cm
arena. Scaling exponent $\alpha$ is calculated via a power-law fit for
$\tau\leq 1$s.

### Fitting Probability Distributions

Two parabolas are fit to a probability distribution smoothed with a Gaussian
filter and forced through a common peak $A_{0}$ in data such as Figure 5a.
Each parabola is either fit to values less than $A_{0}$ for compression or
greater than $A_{0}$ for expansion. The range of values around $A_{0}$ that
the parabolas are fit to be the same for both and is defined separately for
each experiment. This range falls between 1/3 to 2/3 the value of $A_{0}$.

## Acknowledgements

APT was supported by the M.J. Murdock Charitable Trust Award Number
FSU-201913717. APT, GK, MR, and JS were supported by the National Science
Foundation Award Number 2137509. JH, DB, and APT were supported by National
Science Foundation Award Number DMR-2046683.

## Author contributions statement

GK and MR collected experimental data. GK, MR, and ALR contributed analytical
tools. GK, JH, MR, and JS analyzed data. JH and DB conducted simulations. APT,
GK, JH, DB, and MR wrote the manuscript. APT and GK edited the manuscript. APT
designed and conceived the work.

## Additional information

The authors have no competing interests.

## References

  * [1] Visar Ajeti, A. Pasha Tabatabai, Andrew J. Fleszar, Michael F. Staddon, Daniel S. Seara, Cristian Suarez, M. Sulaiman Yousafzai, Dapeng Bi, David R. Kovar, Shiladitya Banerjee, and Michael P. Murrell.  Wound healing coordinates actin architectures to regulate mechanical work.  Nature Physics, 15(7):696–705, July 2019.  Number: 7 Publisher: Nature Publishing Group.
  * [2] Agustí Brugués, Ester Anon, Vito Conte, Jim H. Veldhuis, Mukund Gupta, Julien Colombelli, José J. Muñoz, G. Wayne Brodland, Benoit Ladoux, and Xavier Trepat.  Forces driving epithelial wound healing.  Nature physics, 10(9):683–690, September 2014.
  * [3] Hungtang Ko, Mathias Hadgu, Keyana Komilian, and David L. Hu.  Small fire ant rafts are unstable.  Physical Review Fluids, 7(9):090501, September 2022.
  * [4] Jason M. Graham, Albert B. Kao, Dylana A. Wilhelm, and Simon Garnier.  Optimal construction of army ant living bridges.  Journal of Theoretical Biology, 435:184–198, December 2017.
  * [5] Andrew Berdahl, Colin J. Torney, Christos C. Ioannou, Jolyon J. Faria, and Iain D. Couzin.  Emergent Sensing of Complex Environments by Mobile Animal Groups.  Science, 339(6119):574–576, February 2013.
  * [6] Andrew M Hein, Sara Brin Rosenthal, George I Hagstrom, Andrew Berdahl, Colin J Torney, and Iain D Couzin.  The evolution of distributed sensing and collective computation in animal populations.  eLife, 4:e10955, December 2015.  Publisher: eLife Sciences Publications, Ltd.
  * [7] Y. Katz, K. Tunstrom, C. C. Ioannou, C. Huepe, and I. D. Couzin.  Inferring the structure and dynamics of interactions in schooling fish.  Proceedings of the National Academy of Sciences, 108(46):18720–18725, November 2011.
  * [8] Vivek H. Sridhar, Liang Li, Dan Gorbonos, Máté Nagy, Bianca R. Schell, Timothy Sorochkin, Nir S. Gov, and Iain D. Couzin.  The geometry of decision-making in individuals and collectives.  Proceedings of the National Academy of Sciences, 118(50):e2102157118, December 2021.
  * [9] Tamás Vicsek, András Czirók, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet.  Novel Type of Phase Transition in a System of Self-Driven Particles.  Physical Review Letters, 75(6):1226–1229, August 1995.
  * [10] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic.  Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study.  Proceedings of the National Academy of Sciences, 105(4):1232–1237, January 2008.
  * [11] Marcelo Camperi, Andrea Cavagna, Irene Giardina, Giorgio Parisi, and Edmondo Silvestri.  Spatially balanced topological interaction grants optimal cohesion in flocking models.  Interface Focus, 2(6):715–725, December 2012.
  * [12] Vijay Kumar and Rumi De.  Flocking: Influence of Metric Versus Topological Interactions.  arXiv:2109.10088 [cond-mat, physics:nlin, physics:physics], September 2021.  arXiv: 2109.10088.
  * [13] Colin J. Torney, Myles Lamont, Leon Debell, Ryan J. Angohiatok, Lisa-Marie Leclerc, and Andrew M. Berdahl.  Inferring the rules of social interaction in migrating caribou.  Philosophical Transactions of the Royal Society B: Biological Sciences, 373(1746):20170385, May 2018.
  * [14] S. A. Shah, Y.-L Chen, K. S. Schweizer, and C. F. Zukoski.  Viscoelasticity and rheology of depletion flocculated gels and fluids.  The Journal of Chemical Physics, 119(16):8747–8760, October 2003\.
  * [15] Peter J. Lu, Emanuela Zaccarelli, Fabio Ciulla, Andrew B. Schofield, Francesco Sciortino, and David A. Weitz.  Gelation of particles with short-range attraction.  Nature, 453(7194):499–503, May 2008.
  * [16] Lilian C. Hsiao, Richmond S. Newman, Sharon C. Glotzer, and Michael J. Solomon.  Role of isostaticity and load-bearing microstructure in the elasticity of yielded colloidal gels.  Proceedings of the National Academy of Sciences, 109(40):16029–16034, October 2012.  Publisher: Proceedings of the National Academy of Sciences.
  * [17] Mohammad Nabizadeh and Safa Jamali.  Life and death of colloidal bonds control the rate-dependent rheology of gels.  Nature Communications, 12(1):4274, July 2021.  Bandiera_abtest: a Cc_license_type: cc_by Cg_type: Nature Research Journals Number: 1 Primary_atype: Research Publisher: Nature Publishing Group Subject_term: Colloids;Gels and hydrogels;Rheology Subject_term_id: colloids;gels-and-hydrogels;rheology.
  * [18] W. Benjamin Rogers and John C. Crocker.  Direct measurements of DNA-mediated colloidal interactions and their quantitative modeling.  Proceedings of the National Academy of Sciences, 108(38):15687–15692, September 2011.  Publisher: Proceedings of the National Academy of Sciences.
  * [19] Yu Wang, Yufeng Wang, Xiaolong Zheng, Étienne Ducrot, Jeremy S. Yodh, Marcus Weck, and David J. Pine.  Crystallization of DNA-coated colloids.  Nature Communications, 6(1):7253, June 2015.  Number: 1 Publisher: Nature Publishing Group.
  * [20] Daichi Hayakawa, Thomas E. Videbaek, Douglas M. Hall, Huang Fang, Christian Sigl, Elija Feigl, Hendrik Dietz, Seth Fraden, Michael F. Hagan, Gregory M. Grason, and W. Benjamin Rogers.  Geometrically programmed self-limited assembly of tubules using DNA origami colloids.  Proceedings of the National Academy of Sciences, 119(43):e2207902119, October 2022.  Publisher: Proceedings of the National Academy of Sciences.
  * [21] Dapeng Bi, J. H. Lopez, J. M. Schwarz, and M. Lisa Manning.  A density-independent glass transition in biological tissues.  Nature Physics, 11(12):1074–1079, December 2015.  arXiv:1409.0593 [cond-mat, physics:physics].
  * [22] Xingbo Yang, Dapeng Bi, Michael Czajkowski, Matthias Merkel, M. Lisa Manning, and M. Cristina Marchetti.  Correlating cell shape and cellular stress in motile confluent tissues.  Proceedings of the National Academy of Sciences of the United States of America, 114(48):12663–12668, November 2017.
  * [23] Rachel M. Lee, Douglas H. Kelley, Kerstin N. Nordstrom, Nicholas T. Ouellette, and Wolfgang Losert.  Quantifying stretching and rearrangement in epithelial sheet migration.  New journal of physics, 15(2):025036, February 2013.
  * [24] Rui Ni, James G. Puckett, Eric R. Dufresne, and Nicholas T. Ouellette.  Intrinsic Fluctuations and Driven Response of Insect Swarms.  Physical Review Letters, 115(11):118104, September 2015.
  * [25] Kasper van der Vaart, Michael Sinhuber, Andrew M. Reynolds, and Nicholas T. Ouellette.  Mechanical spectroscopy of insect swarms.  Science Advances, 5(7):eaaw9305, July 2019.
  * [26] S. C. Takatori, W. Yan, and J. F. Brady.  Swim Pressure: Stress Generation in Active Matter.  Physical Review Letters, 113(2):028103, July 2014.
  * [27] Mark J. Bowick, Nikta Fakhri, M. Cristina Marchetti, and Sriram Ramaswamy.  Symmetry, Thermodynamics and Topology in Active Matter.  arXiv:2107.00724 [cond-mat], August 2021.  arXiv: 2107.00724.
  * [28] Tristan Walter and Iain D Couzin.  TRex, a fast multi-animal tracking system with markerless identification, and 2D estimation of posture and visual fields.  eLife, 10:e64000, February 2021.  Publisher: eLife Sciences Publications, Ltd.
  * [29] N. C. Makris.  Fish Population and Behavior Revealed by Instantaneous Continental Shelf-Scale Imaging.  Science, 311(5761):660–663, February 2006.
  * [30] Ch. Becco, N. Vandewalle, J. Delcourt, and P. Poncin.  Experimental evidences of a structural and dynamical transition in fish school.  Physica A: Statistical Mechanics and its Applications, 367:487–493, 2006.
  * [31] Gen Li, Dmitry Kolomenskiy, Hao Liu, Ramiro Godoy-Diana, and Benjamin Thiria.  Intermittent versus continuous swimming: An optimization tale.  Physical Review Fluids, 8(1):013101, January 2023.  Publisher: American Physical Society.
  * [32] T Kiang.  Random fragmentation in two and three dimensions.  Zeitschrift für Astrophysik, Vol. 64, p. 433, 64:433, 1966.
  * [33] D Weaire, JP Kermode, and J Wejchert.  On the distribution of cell areas in a voronoi network.  Philosophical Magazine B, 53(5):L101–L105, 1986.
  * [34] Nitin Kumar, Harsh Soni, Sriram Ramaswamy, and A. K. Sood.  Flocking at a distance in active granular matter.  Nature Communications, 5(1):4688, September 2014.  Number: 1 Publisher: Nature Publishing Group.
  * [35] Michael Berhanu and Arshad Kudrolli.  Heterogeneous structure of granular aggregates with capillary interactions.  Physical review letters, 105(9):098002, 2010.
  * [36] Clare J Dibble, Michael Kogan, and Michael J Solomon.  Structure and dynamics of colloidal depletion gels: Coincidence of transitions and heterogeneity.  Physical Review E, 74(4):041403, 2006.
  * [37] Lilian C Hsiao, Benjamin A Schultz, Jens Glaser, Michael Engel, Megan E Szakasits, Sharon C Glotzer, and Michael J Solomon.  Metastable orientational order of colloidal discoids.  Nature Communications, 6(1):8507, 2015.