Spatial Confinement Affects the Heterogeneity and Interactions Between Shoaling Fish

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 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 ± 0.2 (mean ± stdev) (Figure 1).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 ρ g = N/(π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.
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) by comparing positions ⃗ r(t) as a function of elapsed time (τ ) 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 τ , the MSD turnover and plateau are set by the finite size of the arena.For small τ , we observe power-law scaling (i.e.MSD ∼ τ α ) where α characterizes the type of motion.We find that α depends on confinement, demonstrating that fish motion is super-diffusive (1 < α < 2) and approaches ballistic motion (α − → 2) as containers get larger (Figure 2e).The MSD of the shoal's center of mass has similar arena-size dependent values of α (Figure 2e, Supplemental Figure 1).
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 ⃗ 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 ⟨A H ⟩ /π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 ⃗ 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.
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 A/N int (Figure 4c).
We note that the time courses of the A H (Figure 3b), A (Figure 4b), and 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 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.
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 ∼ 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 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).

Simulations
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: where k and θ are the shape and scale parameters, respectively, while Γ( * ) 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 The scale parameter θ as a function of arena radii R.
(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.25cmor 44.25cm, 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.To elucidate the influence of the arena radius on fish interactions and the interplay between heterogeneity and dynamical activities, we implement a Monte Carlo that simulates the fish positions with different arena radii and β = 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 + σδX i,n , where X i,n represents the position of fish i at step n, σ is the step size, and δ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 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 ∆E is accepted with probability min(1, e −β∆E ).We conducted the simulation with a step size σ = 0.4 over 10 7 iterations, recording the internal area metrics every 100 step.As β 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 ∼ 3.75 ± 0.1 and a scale parameter θ(R) = 0.029 R 2 , which leads to a quadratic relationship between ⟨A⟩ or A 0 and R. To find out the parameters that best describe the experiments, we systematically sweep the parameter space of r ranging from 0.1 to 20 and β 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 β 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 β display a non-monotonic variation as a function of R, with a transition point located between R = 14cm and 34.25cm.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 β 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 β = 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 β(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 ρ 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 v 2 ∼ 1/β.
Our simulations indeed show that an increase in β 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 ≈ 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 = 14cm and R = 34.25cmsuch 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.

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.25cmarena to 18.8s in the smallest R = 8cm arena.Scaling exponent α is calculated via a power-law fit for τ ≤ 1s.

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 .

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

Figure 2 :
Figure 2: Confinement affects motion.(a) Images of arenas with radii R = 8cm, 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 (τ ) for fish in five different arenas.(e) The short-time power-law slope of MSD (α) as a function of radii (R).Experimental duplicates of α are plotted in grey (N = 3 for each), with the mean plotted in black.α 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).

Figure 3 :
Figure 3: Shoal area fluctuates in time.(a) Convex hull of 25 fish in R = 44.25cmarena.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 (πR 2 ) as a function of time for 25 fish in R=34.25cm arena.(c) Time average convex hull area ⟨A H ⟩ normalized by arena area for arenas of different radii.Error bars are one standard deviation.

Figure 4 :Figure 5 :
Figure 4: Confinement affects local fish packing.(a) Voronoi tessellation (polygons) of 25 fish (grey markers) in R = 44.25cmarena.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 (πR 2 ) in R = 34.25cmarena.(c) Sum of internal areas (Σ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.

Figure 6 :
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.25and 44.25 cm radii arenas.f Γ (k, θ) denotes the PDF of Gamma distribution defined in Equation (4).(d,e) At R = 34.25cmand 44.25cm the Gamma distribution fit is done for A < 1.5A 0 and a normalization factor is applied to f Γ (k, θ) to align the modes.Insets in (d,e) show exponential tails (black dashed lines) with decay parameters of 0.0098 and 0.0038cm −2 , respectively.(f) The shape parameter k as a function of arena radii R. (f-inset)

Figure 7 :
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 = 8cm, 11cm, and 34.25cm.(d) The preferred distance r 0 and (d-inset) parameter β as a function of arena radii R from Monte Carlo simulations.