1. Introduction

Experimental and theoretical models have been recently developed to examine fundamental aspects of collective motion exhibited by various biological systems. Following the seminal papers [18]-[17] and works cited therein we suggest that hydrodynamic interactions between the swimmers lead to collective motion when every bacterium interacts with other ones through the viscous environment. The implicit evidences of collective motion were confirmed by reduction in the effective viscosity [16]. The theoretical investigations of collective motion were based on the considering the motion in the framework of mechanical dynamical systems [13]-[20]. While the above presented theoretical and experimental results are related via the effective viscosity of suspension of swimmers, these relation lack simple and direct comparisons of viscosity for experimentally observed sets of bacteria and for simulated ensembles.

In the present paper, we propose a new quantitative criterion of collective behavior. The locations of bacteria are modeled by short segments having a small width randomly embedded in medium without overlapping. First, we theoretically simulate locations of particles called below by disordered sets of bacteria (DB sets for shortness) subjected to local viscous stresses and randomly reacted on hydrodynamic interactions. Second, we calculate the basic statistic moments of the constructed DB sets in terms of the generalized Eisenstein-Rayleigh sums ($e$-sums for shortness) introduced in [6]. Further, we construct the $e$-sums (2) for the observed experimental locations (31 film frames) of Bacillus subtilis in a very thin liquid film [18]. As the final step, we compare the $e$-sums of DB sets and of the experimental locations. The obtained results show a difference between these two sets of $e$-sums that demonstrates the collective motion of bacteria. We do not explain reasons for the collective motion and refer to [18]-[17], [13]-[20].

Our criterion of the collective behavior is based on a new RVE (representative volume element) theory proposed in [6], on the invariance of the effective transport properties on the conformal mappings [4] and on the algebraic dependence of the viscous lattice sums on the $e$-sums [6]. The main object of this theory is a set of $e$-sums systematically investigated in [10], [8], [2], [3], [4], [15]. According to the new RVE theory, the set of the $e$-sums (see (2) below and general formulas in [6]) determines macroscopically equivalent cells, i.e. cells having the same effective properties (conductivity, viscosity etc). The RVE is chosen as the minimal size cell from all the equivalent cells. The necessary justification of this theory can be found in the above cited works. It is worth noting that the usage of the $e$-sums implicitly takes into account high order correlation functions without their hard direct computations. In particular, we make simple implementation of the RVE theory [6] which demonstrates its numerical advantages in comparison with expensive computations based on the traditional statistical notion of the RVE.

It is worth noting that the criterion can be applied to various two-phase dispersed media (biological systems, composites etc).

---

2. Random location of segments

In the present section, we discuss theoretical simulations of elements embedded in a medium. Bacteria are modeled by short segments having non-overlapping thin $\delta$-security coatings. Random locations of segments are generated in the following way. The plane geometry is considered as the complex plane $\mathbb{C}$ of the complex variable $z=x+i y$ with standard designations accepted in complex analysis where $i$ denotes the imaginary unit, Re and Im the real and imaginary parts, the bar stands for the complex conjugation.

The complex numbers 1 and $i$ can be considered as the fundamental translation vectors on the complex plane $\mathbb{C}$. We introduce the $(0,0)$-cell as the square

$Q_{(0,0)}= \left\{z = t_{1}+i t_{2}: -\frac{1}{2}<t_{k}<\frac{1}{2}\;(k=1,2)\right\}.$

The square lattice $\mathcal{Q}$ consists of the cells $Q_{(m_{1},m_{2})}:=\{z\in \mathbb{C}:z-m_{1}-i m_{2}\in Q_{(0,0)}\}$, where $m_{1}, m_{2}$ run over integers.

Consider $N$ non-overlapping segments $\Gamma_{k}$ of length $l$ with the centers $b_k\in Q_{(0,0)}$ with the angle of inclination $\alpha_{k}\in [-\pi,\pi]$. (see Fig.1).

The centers $b_{k}$ are considered as random variables distributed in such a way that the segments

$\Gamma_{k} =\left\{z\in\mathbb{C}: z = b_k\pm\frac{l}{2}e^{i\alpha_k}t,\; 0\leq t\leq 1\right\}$

generate a set of uniformly distributed non-overlapping segments. Theoretically, this distribution can be introduced as the distribution of the variable $\mathbf{b}=(b_1, b_2, \cdots, b_N) \in Q_{(0,0)}^N$ with the restrictions $|b_m- b_n|\geq\delta$ for $m \neq n$ ($m, n =1,2,\ldots,N$). The separation parameter $\delta$ is taken equal to $\frac{l}{4}$. This value refers to the minimum distance between the centers of the bacteria. The choice of the constant $4$ is based on the experimentally observed widths of bacteria. It should be noted that the segments $\Gamma_k$ belong to the cell $Q_{(0,0)}$ in the double periodic topology when the opposite sides of $Q_{(0,0)}$ are glued by pairs.

The random variable $\mathbf{b}$ can be statistically realised for large $N$ by the Monte Carlo method to get numerical results. The following constructive procedure to generate random locations of segments is used. Let a random point $b_1$ is taken in accordance with the uniform distribution in $Q_{(0,0)}$. Next, a random angle $\alpha_1\in[-\pi,\pi]$ is chosen. Hence, the pair of points $\mathbf{z_1} = b_1 \pm \frac{l}{2}e^{i\alpha_1}$ is constructed. This is equivalent to the construction of the first segment $\Gamma_1$ determined by the point $b_1$ and angle $\alpha_1$. At the next step, we take a random point $b_2$ uniformly distributed in $Q\backslash H_{1}$, where $H_1 = \{z\in\mathbb{C}\colon|z-b_{1}|\leq\delta\}$. Further, a random angle $\alpha_2$ is selected and we check whether the segments with the ends $\mathbf{z_2} = b_2 \pm \frac{l}{2}e^{i\alpha_2}$ and $\mathbf{z_1}$ do not intersect. If it is true, we have the second random segment determined by the point $b_2$ and angle $\alpha_2$. If the segments intersect, we take a random point $b_2$ again in $Q\backslash H_{1}$ and randomly select a new random angle $\alpha_2$. In the same way, we take the next point $b_3$ uniformly distributed in $Q\backslash(H_{1}\cup H_{2})$, where $H_2 = \{z\in\mathbb{C}\colon|z-b_{2}|\leq\delta\}$ and so forth. The last random point $b_{N}$ is uniformly distributed in $Q\backslash \big(\bigcup_{k = 1}^{N-1}H_{k}\big)$ and determines a pair of points $\mathbf{z_N} = b_N \pm \frac{l}{2}e^{i\alpha_N}$ with the random angle $\alpha_N$.

Introduce the density segments associated to the conformally invariant conductivity (capacity) [4]

$\varrho(l,N) = N\left(\frac{l}{2}\right)^{2}, \label{dens}$

(1)

where $l$ is the length of the segment and $N$ the number of segments per representative cell. The algorithm described above generates a probability distribution $\mathcal{U}_{\varrho}$ depending on the density. This distribution models the DB sets. In the limit case $\varrho=0$, the distribution $\mathcal{U}_{0}$ becomes the well known Poisson distribution of points in the square.

In order to study the distributions $\mathcal{U}_{\varrho}$ for various densities $\varrho$ we will use the following $e$-sums [9]

$\begin{aligned} e_2 = \frac{1}{N^{2}}\sum\limits_{k=1}^{N}\sum\limits_{m=1}^{N}E_2(b_k -b_m),\\ e_{pp} = \frac{(-1)^p}{N^{p+1}}\sum\limits_{m=1}^{N}\left|\sum\limits_{k=1}^{N}E_p(b_m-b_k)\right|^{2}, \; p=2,3,4. \end{aligned} \label{eq:ee}$

(2)

Here, $E_p$ denotes the Eisenstein functions of order $p$ (see Appendix). The values of the basic sums $e_{2}$, $e_{22}$, $e_{33}$ and $e_{44}$ will be estimated using the Monte Carlo method.

---

3. Computation of $e$-sums for DB sets

The theoretical probabilistic distributions corresponding to disordered locations of bacteria are modeled in the previous section. Now, we propose an effective computational tool to properly describe $\mathcal{U}_{\varrho}$.

Let $M$ denote the number of simulated realizations of the unit cell (Monte Carlo experiments) with the random value $\mathbf{b}$. The parameters $M$ and $N$ must be chosen sufficiently large in order to obtain the stable averaged value of the $e$-sums. Theoretically, it is possible to consider only one experiment ($M=1$) and to take a huge number $N$. But computations are less expensive if $M$ is large and $N$ can be not so huge as in the case $M = 1$.

First, we determine the minimal $N$ for which the distribution is characterized in that the angles of the segments will be uniformly distributed over the interval $[-\pi,\pi]$. It is assumed that the segments directions are uniformly distributed if the following inequalities are fulfilled

$\left|\frac{1}{N}\sum\limits_{n=1}^{N}\rm{Re}\big[\exp(i\alpha_{n})\big]\right|\leq 0.15\;\;\rm{and}\;\; \left|\frac{1}{N}\sum\limits_{n=1}^{N}\rm{Im}\big[\exp(i\alpha_{n})\big]\right|\leq 0.15, \label{criterion}$

(3)

These conditions are satisfied for $N = 500$ (see Fig.2).

In order to estimate $M_0=M$ when the computations become stable we consider the dependence of the mean $\langle e_{44} \rangle$ on $M$ for $M\in[1,1500]$ experiments. All the sums $e_2$ and $e_{pp}$ ($p=2,3,4$) are estimated, but $e_{44}$ is characterized by major volatility. Here, $\langle e_{44} \rangle$ is equal to the mean value of $(e_{44})_m$ calculated for $\mathbf{b}$ generated in the $m$-th numerical experiment, more precisely,

$\langle e_{44} \rangle = \frac{1}{M} \sum\limits_{m=1}^M (e_{44})_m. \label{mean_e44}$

(4)

The results are shown in Fig.3. One can observe that errors do not exceed $2\%$ for $M > 700$. Therefore, the computations demonstrate that we have to simulate at least $M = 800$ cells each of them contains at least $N = 500$ segments.

The values $e_2, e_{22}, e_{33}, e_{44}$ are computed for the distributions $\mathcal{U}_{\varrho}$ for various densities $\varrho$. The sum $e_2$ must be equal to $\pi$ [9] for ideal macroscopically isotropic random locations of segments. This yields the first criterion equation for the macroscopic isotropy of structures. The sums $e_{pp}$ ($p=2,3,4$) describe high order basic terms of the distributions $\mathcal{U}_{\varrho}$. The next terms could describe more precisely $\mathcal{U}_{\varrho}$ for higher densities [8].

The average $e$-sums are selected in Table 1 for the densities $\varrho$ changing from $0.05$ to $0.6$ with the step $0.05$. The algorithm works too slowly for higher densities. The results from Table 1 can be extended by expensive computations to higher densities.

Table 1 contains the fundamental parameters of the uniform non-overlapping distribution $\mathcal{U}_{\varrho}$ of segments on the plane. The simulated theoretical distribution describes DB sets when bacteria may affect each other but these interactions are local and do not yield the collective behavior.

---

4. Collective motion of bacteria

In the present section, we pay attention to experimental results partially presented in [18]. The images of Bacillus subtilis in $31$ frames are used in computations. One of the typical frame is displayed in Fig.4.

We use algorithms of image processing and analysis to determine number, centers, angles of inclinations and length of bacteria. The density of bacteria is calculated by formula (1) and oscillates around value $0.15$. The results of the image processing and analysis are applied to computation of the values of the experimental $e$-sums. The results are selected in Table 2.

In order to compare the distributions of DB sets with the distribution of bacteria we have made theoretical calculation of $31$ samples. Because of the average number of bacteria in the frame is about $2050$, the minimum distance between the centers of the bacteria is $\frac{1}{4}$ of their length and density of bacteria equals about $0.15$, the theoretical calculations have been carried out for the following parameters $N = 2050$, $\varrho = 0.15$ and $\delta = \frac{l}{4}$. The length $l = 0.017108$ is normalized to the normalized unit area of the cell. The results of the calculations are shown in Table 3.

---

5. Conclusion

Comparing the results shown in the Tables 2 and 3 for the observed and theoretical distributions of bacteria, we can see that values of the corresponding $e$-sums differ. Therefore, the observed locations of bacteria are not fitted with the simulated disordered locations. We can conclude that the behavior of the bacteria is not disordered, hence, we suggest that it is collective. An exception is the values of $e_2$ which is close to $\pi \approx 3.14$ for the DB sets and the observed locations (see Table 4). This demonstrates averaged isotropy of the bacteria motion for the considered data.

In order to see that the theoretical and observed distributions are essentially different, we compare averaged basic sums (see Table 4) from Table 2 and Table 3. It is worth noting that changes of the parameters ($\delta$, $M$ etc) do not essentially impact onto deviations of the $e$-sums of DB sets. Moreover, the DB sets have significantly higher values $e_{pp}$ ($p=2,3,4$) than the observed ones (c.f., the bold line in Table 1 and the averaged data from Table 4).

The above analysis of collective behavior can be considered as the first application of the RVE theory [6] which will be extended in the future. In particular, we plan to study dynamical parameters of the bacteria distributions in time using the data displayed in Fig.5.

---

Appendix

Following [6], [7] we present constructive formulae for the Eisenstain-Rayleigh sums $S_m$ the Eisenstain functions $E_m$ corresponding to square lattice $\mathcal{Q}$ (see section 2).

The Eisenstein--Rayleigh lattice sums $S_{m}$ can be easily calculated through the rapidly convergent series (43) and recurrent formulae (44) from paper [7]. For the square array, we have $S_2 = \pi =3.14159$, $S_4 = 3.15121$, $S_8 =4.25577$, $S_{12} =3.93885$. Only the sums $S_m$ ($m=2$; $m=4,8,12,16, \ldots$) do not vanish for the square array.

The Eisenstein functions [21] are related to the Weierstrass function $\wp (z)$ [1] by the identities

$\begin{aligned} E_{2}(z) =\wp(z)+S_2,\\ E_{m}(z) =\frac{(-1)^{m}}{(m-1)!}\frac{d^{m-2} \wp(z)}{dz^{m-2}}, m=3,4,\ldots. \end{aligned} \label{Em}$

(5)

Every function (5) is doubly periodic and has a pole of order $m$ at $z=0$. For shortness, it is convenient to redefine the Eisenstein functions in (2) at zero as $E_p(b_k-b_m):=S_p$ for $k=m$.

---

Acknowledgments

We thank Alexei Sokolov who kindly presented the experimental results partially reported in [18] and Leonid Berlyand for stimulated discussion during his visit to Kraków and VM's visits to Penn State University.

The authors thank the anonymous reviewer for valuable comments and suggestions to improve the paper.

---