# Stability analysis for large-scale multi-agent molecular communication
systems

Taishi Kotsuka1, Yutaka Hori1 0000-0002-7897-3131 0000-0002-3253-4985

###### Abstract

Molecular communication (MC) is recently featured as a novel communication
tool to connect individual biological nanorobots in vivo. It is expected that
a large number of nanorobots can form large multi-agent MC systems through MC
to accomplish complex and large-scale tasks that cannot be achieved by a
single nanorobot. However, most previous models for MC systems assume a
unidirectional diffusion communication channel and cannot capture the feedback
between each nanorobot, which is important for multi-agent MC systems. In this
paper, we introduce the system theoretic model for large-scale multi-agent MC
systems using transfer functions, and then propose a method to analyze the
stability for multi-agent MC systems. The proposed method decomposes the
multi-agent MC system into multiple single-input and single-output (SISO)
systems, which facilitates to analyze the stability of the large-scale multi-
agent MC system. Finally, we demonstrate the proposed method by analyzing the
stability of a specific large-scale multi-agent MC system.

†† This work was supported in part by JSPS KAKENHI Grant Numbers JP21H05889,
JP22J10554, and JP23H00506, and in part by JST SPRING Grant Number JPMJSP2123.
1 Department of Applied Physics and Physico-Informatics, Keio University.
3-14-1 Hiyoshi, Kohoku-ku, Yokohama, Kanagawa 223-8522, Japan.
Email addresses<EMAIL_ADDRESS>(Taishi Kotsuka<EMAIL_ADDRESS>(Yutaka Hori).
TK, 0000-0002-7897-3131, YH, 0000-0002-3253-4985.

Keywords:  Molecular communication, Feedback control, Transfer function,
Biological control systems, Diffusion equation.

##  1 Introduction

In nature, bacteria are known to communicate with each other using diffusion-
based molecular communication (MC) to achieve cooperative functions such as
biofilm formation. Recently, many researchers attempted to use this MC as a
novel communication tool to connect individual biological nanorobots
(artificial cells) in vivo [1, 2, 3, 4, 5]. In particular, a large number of
dispersed nanorobots can form a large multi-agent system through MC, which is
expected to enable nanorobots to accomplish complicated and large-scale tasks,
leading to engineering applications such as targeted drug delivery [6, 7]. In
order to design such a multi-agent MC system that can achieve the desired
performance, it is important to construct a dynamical model for a large-scale
multi-agent MC system and analyze its fundamental properties.

An important feature that needs to be incorporated into the model of multi-
agent MC systems is the information sharing and the feedback of signals
between nanorobots. Dynamic models of one-to-one and one-to-$n$ MC systems
were constructed based on diffusion equations, and many communication
properties were analyzed [8, 9, 10, 11, 12, 13]. In [8], a system theoretic
model of a MC system which consists of a sender nanorobot and a receiver
nanorobot was constructed and the frequency response characteristics of the MC
channel were analyzed based on Green’s function of diffusion equation.
However, since these models assume unidirectional MC channels, they cannot be
easily extended to models for $n$-to-$n$ multi-agent systems, which require
bidirectional MC channels for feedback. A model capturing the feedback between
each nanorobot was proposed by the authors’ group [14], and its stability was
analyzed. However, the previously proposed stability analysis method [14] was
applicable only to the MC systems with two nanorobots. On the other hand,
there are several studies with large-scale multi-agent MC systems considering
the reaction in the nanorobots based on the reaction-diffusion (RD) equation
[15, 16, 17, 18]. However, RD-based MC models assume that the length between
the nanorobots is close enough so that the disruption of the signal by the MC
channel can be ignored, and thus, the model is not suitable for the analysis
and design of the multi-agent MC systems composed of a population of
distributed nanorobots for applications in vivo such as drug delivery. This
motivates us to develop a more versatile model and analysis method for large-
scale multi-agent MC systems for broader engineering applications of MC
systems.

In control engineering, many design/analysis methods for large-scale multi-
agent systems were developed [19, 20, 21, 22]. In [23], several models and
methods to analyze the stability of multi-agent systems based on graph theory
and control theory using transfer functions were introduced. However, since
these models consider high-speed communication such as electrical or optical
communication, the dynamics of communication channels are ignored while the
dynamics of MC channels need to be considered because of slow-speed
communication due to diffusion of signal molecules, and thus the
design/analysis methods previously developed in control engineering are not
applicable for multi-agent MC systems. From a control engineering perspective,
it is therefore a new challenge to construct system theoretic models and
methods to control large-scale multi-agent systems with communication channels
that have diffusion-based dynamics.

In this paper, we introduce a system theoretic model of large-scale multi-
agent MC systems that incorporates the transfer functions of bidirectional
communication and local reactions inside nanorobots, and then, propose a
method to analyze the stability of large-scale multi-agent MC systems as one
of the fundamental control properties. Specifically, we first approximately
model large-scale multi-agent MC systems as circulant MC systems that consist
of $n$ homogeneous nanorobots with periodic boundaries at the left/right ends
of the system. We then obtain a system theoretic model of the circulant MC
system, which contains a large-scale transfer function matrix based on the
transfer function of bidirectional MC channels previously derived by the
authors’ group [24]. The proposed method decomposes the transfer function
matrix of the circulant MC system into multiple single-input and single-output
(SISO) systems in the same spirit as the stability analysis for linear systems
with generalized frequency variables [25]. This decomposition allows for the
stability analysis of the large-scale multi-agent MC system by combining
simple analysis techniques for SISO systems. Moreover, since the decomposed
system captures the information for spatial frequencies of multi-agent MC
systems, the proposed method can analyze the response of multi-agent MC
systems based on the spatial frequency, which potentially leads to the
analysis for the condition for biological pattern formation such as Turing
pattern formation [26]. Finally, we demonstrate the proposed method by
analyzing the stability of a specific large-scale multi-agent MC system.

This paper is organized as follows. In the next section, we model large-scale
multi-agent MC systems as a circulant MC system with periodic boundaries based
on a diffusion equation. We then introduce system theoretic models of the
circulant MC system expressed as multi-inputs and multi-outpus (MIMO)
dynamical systems using a transfer function in Section 3. In Section 4, we
provide a method to analyze the stability of the circulant MC system by
decomposing the system into multiple SISO systems and show the relation
between the decomposed systems and spatial frequency. The proposed method is
demonstrated for the analysis of a specific large-scale multi-agent MC system
in Section 5. Finally, the paper is concluded in Section 6.

Notations: The following notations are used throughout this paper: The set of
real values is defined by $\mathbb{R}$, and the set of complex values is
defined by $\mathbb{C}$. The superscript is used to represent the dimension of
the vector space, e.g., $\mathbb{R}^{n}$. The identity matrix is defined by
$I_{n}\in\mathbb{R}^{n\times n}$. The Laplace transform of a function $z(t)$
is defined by
$Z(s)\coloneqq\mathcal{L}\\{z(t)\\}=\int^{\infty}_{0}z(t)e^{-st}dt$, where
$s=\sigma+j\omega$ is the complex variable with the real part $\sigma$ and the
imaginary part $\omega$. The determinant of a matrix $A$ is defined as
$\det(A)$.

##  2 Model of multi-agent MC systems and problem formulation

###  2.1 Mathematical models

We consider multi-agent MC systems with countless homogeneous nanorobots in
one-dimensional spatial domain, where the nanorobots could be bacterial cells
or vesicles that have reaction systems consisting of biomolecules such as
genes and proteins as shown in Fig. 1. Each nanorobot produces and emits
signal molecules, which can diffuse in the fluidic environment and are
absorbed by the other nanorobots to interact with the other molecules in the
nanorobot.

To model the dynamics of the multi-agent MC system, we consider a large number
of nanorobots, say $n$ nanorobots, in one-dimensional space. The left and the
right end of the space is connected as shown in Fig. 1 (blue dotted box),
i.e., the boundary is periodic, to approximately model the dynamics of
countless nanorobots using a small sub-unit in the space (see Fig. 1). We call
this system a circulant MC system. Let $\Sigma_{i}\,(i=1,2,\ldots,n)$ denote
the reaction system that captures the reaction of molecules in the $i$-th
nanorobot and the emission/absorption of the signal molecules and $\Gamma_{i}$
denote the $i$-th MC channel between the reaction systems $\Sigma_{i-1}$ and
$\Sigma_{i}$. In what follows, the subscript $i$ is taken by modulo $n$, i.e.,
$\Sigma_{0}:=\Sigma_{n}$ since periodic boundaries are applied.

Figure 1: The multi-agent MC system, where multiple nanorobots communicate
with each other via MC channels $\Gamma_{i}$. The reaction system $\Sigma_{i}$
captures the reaction of molecules in the nanorobot and the
emission/absorption of the signal molecule.

We denote the concentration of the signal molecules at position $r\in[0,L]$ in
the MC channel $\Gamma_{i}$ by $c_{i}(t,r)$, where $L$ is the communication
length of each channel $\Gamma_{i}$. In practice, the communication length
between nanorobots may not be uniform, but we here consider an average length
$L$ to approximately model the collective behavior of many nanorobots. The
dynamics of the concentration $c_{i}(t,r)$ can be then modeled by the
diffusion equation

$\frac{\partial c_{i}(t,r)}{\partial
t}=\mu\frac{\partial^{2}c_{i}(t,r)}{\partial r^{2}},$ (1)

where $\mu$ is the diffusion coefficient. The boundary conditions at $r=0$ and
$r=L$ are dynamic Dirichlet boundary conditions defined by

$\displaystyle c_{i}(t,0)$ $\displaystyle=y_{i-1}(t),$ (2) $\displaystyle
c_{i}(t,L)$ $\displaystyle=y_{i}(t),$ (3)

where $y_{i}(t)$ is the emitted concentration of the signal molecule from the
$i$-th reaction system $\Sigma_{i}$, which is thus determined by the dynamics
of the reaction system $\Sigma_{i}$.

For many practical examples [16, 27], the dynamics of the $i$-th reaction
system $\Sigma_{i}$ associated with $m-1$ species of molecules can be modeled
by nonlinear state-space models as

$\displaystyle\displaystyle\frac{d\bm{x}_{i}(t)}{dt}$
$\displaystyle=\bm{f}(\bm{x}_{i})+\bm{B}w_{i}(t),$ (4) $\displaystyle
y_{i}(t)$ $\displaystyle=\bm{C}\bm{x}_{i}(t),$ (5)

where the state $\bm{x}_{i}(t)=[x_{i,1},x_{i,2},\ldots,x_{i,m}]^{T}$ is the
concentrations of the molecules associated with reactions occurring inside and
outside of the $i$-th nanorobot, and $x_{i,l}$ is the concentration of the
$l$-the molecular species. The variable $x_{i,m}$ represents the concentration
of the signal molecule outside of the $i$-th nanorobot. The vector function
$\bm{f}(\cdot)$ represents the dynamics of the reactions in the nanorobot and
the membrane transport, and $\bm{B}=[0,\ldots,0,1]^{T}\in\mathbb{R}^{m}$ and
$\bm{C}=[0,\ldots,0,1]\in\mathbb{R}^{1\times m}$ are the input and the output
vectors, respectively. The variable $w_{i}(t)$ is the input fluxes from the MC
channels $\Gamma_{i}$ and $\Gamma_{i+1}$ to the $i$-th reaction system
$\Sigma_{i}$ following the Fick’s first law and the conservation law for the
total number of the signal molecule. Specifically, $w_{i}(t)$ is

$w_{i}(t)=\frac{\mu}{\Delta r}\left(\frac{\partial c_{i+1}(t,0)}{\partial
r}-\frac{\partial c_{i}(t,L)}{\partial r}\right),$ (6)

where $\Delta r$ is the size of the reaction system $\Sigma_{i}$. The output
$y_{i}(t)$ is the concentration of the signal molecule outside of the
nanorobot.

###  2.2 Problem formulation

A potential application of the multi-agent MC system is to synchronize the
states of nanorobots, which is important to make cooperative behaviors at a
population level. In fact, the circulant MC system can converge to a spatially
homogeneous equilibrium point at steady states by designing the reaction
systems and the MC channels appropriately. In other words, the population of
nanorobots can synchronize their state against small perturbation around an
equilibrium point at steady state. To analyze the convergence to the
synchronized states, an important concept is the stability of the system.
Specifically, the states of nanorobots converge to the homogeneous equilibrium
point if and only if the MC system is asymptotically stable around the
homogeneous equilibrium point. In what follows, we introduce the definitions
of the stability to develop mathematically rigorous stability analysis methods
in the following sections.

######  Definition 1.

Consider the system (1) – (6). A constant vector $\bm{x}^{*}\in\mathbb{R}^{m}$
is a spatially homogeneous equilibrium point if

  1. 1.

$\bm{f}(\bm{x}^{*})=0$, and

  2. 2.

$c_{i}(t,r)=x^{*}_{m}$ for all $i=1,2,\ldots,n$ and for all $r\in[0,L]$, where
$x^{*}_{m}$ is the $m$-th entry of $\bm{x}^{*}$.

######  Assumption 1.

The circulant MC system (1) – (6) has a spatially homogeneous equilibrium
point $\bm{x}^{*}$.

Next, we show the definition of the stability and the asymptotic stability of
the circulant MC system (1) – (6) around a spatially homogeneous equilibrium
point $\bm{x}^{*}$.

######  Definition 2.

Consider the system (1) – (6), and let $\bm{x}^{*}$ be a spatially homogeneous
equilibrium point. The spatially homogeneous equilibrium point is locally
stable if the following 1 and 2 hold:

  1. 1.

For all $\epsilon_{1}>0$, if there exists $\delta_{1}>0$ such that
$||\bm{x}_{i}(0)-\bm{x}^{*}||<\delta_{1}$, then
$||\bm{x}_{i}(t)-\bm{x}^{*}||<\epsilon_{1}$ for all $t\geq 0$ and
$i=1,2,\ldots,n$.

  2. 2.

For all $\epsilon_{2}>0$, if there exists $\delta_{2}>0$ such that
$\int^{L}_{0}{||c_{i}(0,r)-x_{m}^{*}||dr}<\delta_{2}$, then
$\int^{L}_{0}{||c_{i}(t,r)-x_{m}^{*}||dr}<\epsilon_{2}$ for $t\geq 0$ and all
$i=1,2,\ldots,n$.

######  Definition 3.

Consider the system (1) – (6), and let $\bm{x}^{*}$ be a spatially homogeneous
equilibrium point. The spatially homogeneous equilibrium point is locally
asymptotically stable if it is stable, and the following 1 and 2 hold:

  1. 1.

If there exists $\delta_{1}>0$ such that
$||\bm{x}_{i}(0)-\bm{x}^{*}||<\delta_{1}$, then
$\lim_{t\rightarrow\infty}{||\bm{x}_{i}(t)-\bm{x}^{*}||}=0$ for
$i=1,2,\ldots,n$.

  2. 2.

If there exists $\delta_{2}>0$ such that
$\int^{L}_{0}{||c_{i}(0,r)-x_{m}^{*}||dr}<\delta_{2}$, then

$\lim_{t\rightarrow\infty}{\int^{L}_{0}{||c_{i}(t,r)-x_{m}^{*}||dr}}=0$ (7)

for $i=1,2,\ldots,n$.

In the next section, we first express the circulant MC system as a multi-input
multi-output (MIMO) dynamical system using the transfer functions derived from
the diffusion equation (1) and the state-space model (5). We then show that
the MIMO circulant MC system can be decomposed into multiple single-input
single-output (SISO) systems to facilitate the stability analysis.

##  3 System theoretic model for circulant MC systems

Figure 2: The block diagram of the multi-agent MC system $\mathcal{M}$.

In this section, we first derive the transfer function of the circulant MC
system since the stability of systems can be analyzed by examining the roots
of the denominator of the transfer functions. We then show the MIMO
representation of the MC system based on the transfer functions.

Transfer functions represent the ratio of the Laplace transform of the
input/output signals and allow us to analyze and design the system in the
complex domain, where many tools in control theory can be applied. Let
$Y_{i}(s)\coloneqq\mathcal{L}[y_{i}(t)]$ and
$W_{i}(s)\coloneqq\mathcal{L}[w_{i}(t)]$ denote the Laplace transform of the
signals $y_{i}(t)$ and $w_{i}(t)$, respectievly. Specifically, $W_{i}(s)$ is
written as

$W_{i}(s)=\frac{\mu}{\Delta r}\left(\frac{\partial C_{i+1}(s,0)}{\partial
r}-\frac{\partial C_{i}(s,L)}{\partial r}\right),$ (8)

where $C_{i}(s,\cdot)\coloneqq\mathcal{L}[c_{i}(t,\cdot)]$. Then, the transfer
function of the reaction system $\Sigma_{i}$ is defined by

$h(s)\coloneqq\frac{Y_{i}(s)}{W_{i}(s)}=\bm{C}(sI_{m}-A)^{-1}\bm{B},$ (9)

where $A$ is the Jacobian matrix of $\bm{f(\bm{x})}$ at the homogeneous
equilibrium point $\bm{x}^{*}$. The transfer function $h(s)$ represents the
relation between the input $w_{i}(t)$ and the output $y_{i}(t)$ of the
reaction system $\Sigma_{i}$ in the complex domain.

We next consider the transfer function of the MC channel $\Gamma_{i}$. The
input-output relation between the concentrations $y_{i-1}(t)=c_{i}(t,0)$ and
$y_{i}(t)=c_{i}(t,L)$, and the concentration gradients $\partial
c_{i}(t,L)/\partial r$ and $\partial c_{i}(t,L)/\partial r$ can be expressed
using the Laplace transform of these variables as

$\displaystyle\frac{\partial C_{i}(s,0)}{\partial
r}=\frac{1}{2}g_{1}(s)C_{i}(s,0)+g_{2}(s)C_{i}(s,L),$ (10)
$\displaystyle\frac{\partial C_{i}(s,L)}{\partial
r}=-g_{2}(s)C_{i}(s,0)-\frac{1}{2}g_{1}(s)C_{i}(s,L),$ (11)

where

$\displaystyle g_{1}(s)$
$\displaystyle\coloneqq-\sqrt{\frac{s}{\mu}}\frac{2}{\tanh{\left(\frac{L}{\sqrt{\mu}}\sqrt{s}\right)}},$
(12) $\displaystyle g_{2}(s)$
$\displaystyle\coloneqq\sqrt{\frac{s}{\mu}}\frac{1}{\sinh{\left(\frac{L}{\sqrt{\mu}}\sqrt{s}\right)}}.$
(13)

Equations (10) – (11) are derived by the Laplace transform of the diffusion
equation (1) with the Dirichlet boundary conditions at $r=0$ and $r=L$ [24].

Based on the input-output relations (9) – (11), the circulant MC system
consisting of the MC channels $\Gamma_{i}$ and the reaction systems
$\Sigma_{i}$ can be represented as the MIMO system $\mathcal{M}$:

$\begin{split}\bm{Y}(s)&=H(s)\bm{W}(s),\\\
\bm{W}(s)&=G(s)\bm{Y}(s),\end{split}$ (14)

where $\bm{Y}(s):=[Y_{1}(s),Y_{2}(s),\ldots,Y_{n}(s)]^{T}$ and
$\bm{W}(s):=[W_{1}(s),W_{2}(s),\ldots,W_{n}(s)]^{T}$. The transfer function
matrix $H(s)$ is the diagonal matrix as

$H(s)=h(s)I_{n}$ (15)

representing the reaction systems. The transfer function matrix $G(s)$ is the
circulant tridiagonal matrix representing the MC channels defined by

$G(s)\coloneqq\frac{\mu}{\Delta
r}\left[\begin{matrix}g_{1}(s)&g_{2}(s)&0&\cdots&0&g_{2}(s)\\\
g_{2}(s)&g_{1}(s)&g_{2}(s)&\cdots&0&0\\\ 0&g_{2}(s)&\ddots&\ddots&0&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&0\\\
0&\ddots&\ddots&\ddots&\ddots&g_{2}(s)\\\
g_{2}(s)&0&\cdots&0&g_{2}(s)&g_{1}(s)\end{matrix}\right].$ (16)

The block diagram of the MIMO circulant MC system is shown in Fig. 2.

Based on the system theoretic model for the MIMO circulant MC system (14) with
the transfer functions, one can analyze the stability, robustness, and
performance of the system in the complex domain. However, these analyses
become difficult when the MIMO circulant MC system $\mathcal{M}$ consists of a
large number of nanorobots, i.e., $n$ is large since it requires computing the
determinant of a large non-rational matrix. In the next section, we propose
the method to analyze the stability of the MIMO circulant MC system
$\mathcal{M}$ by reducing the problem of computing the determinant of a large
matrix to that of computing multiple non-rational scalar equations.

##  4 Stability analysis for circulant MC systems

In this section, we show a method to analyze the stability of the circulant MC
system. Specifically, we first show that the stability analysis of the MIMO
circulant MC system with $n$ nanorobots can be reduced to that of $n$ SISO
systems by some linear transformation. We then discuss the physical
interpretation of the decomposed SISO system based on spatial frequency.

###  4.1 Stability analysis by decomposition of the MC system

We first show a necessary and sufficient stability condition for the closed-
loop system (14) based on the characteristic equation. The characteristic
equation $p(s)=0$ of the closed-loop system (14) is defined by

$p(s):=\det(I_{n}-G(s)H(s))=0.$ (17)

It is known in control engineering that the closed-loop (1) – (6) is locally
asymptotically stable if and only if all roots of the characteristic equation
$p(s)=0$ lie in the open left half plane (OLHP) of the complex plane (Theorem
5.7 of reference [28]).

The roots of the characteristic equation $p(s)=0$ are not easy to analyze
since $p(s)$ is defined by the determinant of a matrix with non-rational
functions. In what follows, we show a theorem that the stability analysis for
the MIMO system (14) can be reduced to that for $n$ SISO systems using the
property of the circulant matrix $G(s)$.

######  Theorem 1.

Consider the circulant MC system (1) – (6) and suppose Assumption 1 holds. Let
$\bm{x}^{*}$ denote a spatially homogeneous equilibrium point of the circulant
MC system. The following three statements are equivalent.

  1. 1.

The spatially homogeneous equilibrium point $\bm{x}^{*}$ is locally
asymptotically stable.

  2. 2.

All roots of $p(s)=0$ lie in the open left half plane of the complex plane,
i.e., $\mathrm{Re}[s]<0$.

  3. 3.

All roots of the characteristic equations

$\hat{p}_{i}(s):=1-\lambda_{i}(s)h(s)=0$ (18)

lie in the OLHP of the complex plane for all $i$, where

$\displaystyle i=\begin{cases}1,2,\ldots,\frac{n}{2}+1&\mathrm{when\
}n\mathrm{\ is\ even}\\\ 1,2,\ldots,\frac{n+1}{2}&\mathrm{when\ }n\mathrm{\
is\ odd}\end{cases}$ (19)

and

$\lambda_{i}(s)=g_{1}(s)+2g_{2}(s)\cos{\left(2\pi\frac{i-1}{n}\right)}.$ (20)

###### Proof.

The equivalence of the statement 1 and 2 is well-known in control engineering
(Theorem 5.7 of reference [28]), and thus we prove the equivalence of the
statements 2 and 3. Since the MC channel $G(s)$ is the circulant matrix,
$G(s)$ can be diagonalized by the discrete Fourier transform (DFT) matrix
$F\in\mathbb{C}^{n\times n}$ [29] as

$\Lambda(s):=FG(s)F^{*},$ (21)

where $\Lambda(s)$ is the diagonal matrix
$$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}$(\lambda_{1}(s),\lambda_{2}(s),\ldots,\lambda_{n}(s))$
and

$F=\frac{1}{\sqrt{n}}\left[\begin{matrix}1&1&1&\cdots&1\\\
1&\alpha&\alpha^{2}&\cdots&\alpha^{n-1}\\\
1&\alpha^{2}&\alpha^{4}&\cdots&\alpha^{2(n-1)}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
1&\alpha^{(n-1)}&\alpha^{2(n-1)}&\cdots&\alpha^{(n-1)(n-1)}\end{matrix}\right]$
(22)

with $\alpha=e^{-j2\pi/n}$. Using the DFT matrix $F$, the characteristic
equation (17) can be transformed as

$\displaystyle p(s)$
$\displaystyle=\det\left(F^{*}(I_{n}-FG(s)F^{*}FH(s)F^{*})F\right)$
$\displaystyle=\det(I_{n}-h(s)\Lambda(s))$
$\displaystyle=\prod_{i}^{n}[1-\lambda_{i}(s)h(s)]=\prod_{i}^{n}\hat{p}_{i}(s)=0,$
(23)

where we use $FH(s)F^{*}=h(s)I_{n}$ and $\det(F)\det(F^{*})=1$ because $F$ is
the unitary matrix. The theorem holds since Eq. (23) shows that the roots of
the characteristic equation $p(s)=0$ coincide with those of $\hat{p}_{i}(s)=0$
for all $i$. ∎

Theorem 1 shows that the problem of solving the determinant of the non-
rational matrix (17) can be reduced to that of solving non-rational scalar
equation (23). Theoretically, this means that MIMO system $\mathcal{M}$ with
the $n\times n$ circulant matrix $G(s)$ can be decomposed into $n$ SISO
systems $\hat{\mathcal{M}}_{i}$ as shown in Fig. 3, which facilitates the
analysis for the stability of the MC system using methods developed in control
engineering.

Figure 3: Decomposed multi-agent MC system $\hat{\mathcal{M}}_{i}$.

###  4.2 Physical interpretation of decomposed systems

Next, we discuss the physical interpretation of the decomposed system
$\hat{\mathcal{M}}_{i}$ shown in Fig. 3 based on spatial frequency. We
consider the decomposed subsystem $\hat{\mathcal{M}}_{i}$ shown in Fig. 3,
whose input-output relation is given by

$\begin{array}[]{cc}&\hat{Y}_{i}(s)=h(s)\hat{W}_{i}(s),\\\
&\hat{W}_{i}(s)=\lambda_{i}(s)\hat{Y}_{i}(s),\end{array}$ (24)

where

$\begin{array}[]{cc}&\displaystyle\hat{Y}_{i}(s)=\bm{F}_{i}Y(s)=\frac{1}{\sqrt{n}}\sum_{l=1}^{n}Y_{l}(s)\alpha^{(l-1)(i-1)},\\\
&\displaystyle\hat{W}_{i}(s)=\bm{F}_{i}^{T}W(s)=\frac{1}{\sqrt{n}}\sum_{l=1}^{n}W_{l}(s)\alpha^{(l-1)(i-1)},\end{array}$
(25)

and $\bm{F}_{i}\in\mathbb{C}^{1\times n}$ is the $i$-th row of the DFT matrix
$F$. Equation (25) corresponds to the discrete Fourier transform of the signal
$\bm{Y}(s)$ and $\bm{W}(s)$, and thus $\hat{Y}_{i}(s)$ and $\hat{W}_{i}(s)$
represent the Fourier component with a spatial frequency
$\xi_{i}=2\pi(i-1)/nL$ of the signals $\bm{Y}(s)$ and $\bm{W}(s)$,
respectively. This means that the decomposed subsystem $\hat{\mathcal{M}}_{i}$
consisting of $h(s)$ and $\lambda_{i}(s)$ expresses the dynamics of the
Fourier component of the signal with a specific spatial frequency
$\xi_{i}=2\pi(i-1)/nL$. Therefore, the stability of the subsystem
$\hat{\mathcal{M}}_{i}$ implies that the MC system converges to the spatially
homogeneous equilibrium state when a perturbation with a component of spatial
frequency $\xi_{i}$ is input to the system. Thus, the proposed method is able
to analyze for which Fourier component of perturbations the circulant MC
system converges to the spatially homogeneous equilibrium point by analyzing
the stability of the decomposed subsystems.

Next, we consider the circulant MC system with a large number of $n$ by
considering the limit of $n\rightarrow\infty$. For the infinite number of
nanorobots $n\rightarrow\infty$, Eq. (25) implies that the dynamics of the
subsystems $\hat{\mathcal{M}}_{i}$ is asymptotic to the sum of the Fourier
components for all continuous spatial frequencies $\xi$. We here show that the
stability analysis for decomposed subsystems $\hat{\mathcal{M}}_{i}$ when
$n\rightarrow\infty$ boils down to the problem of computing a non-rational
scalar equation for a certain range of the parameter.

######  Lemma 1.

Consider the characteristic equation (18) and define

$\mathbb{X}_{n}\coloneqq\left\\{s\in\mathbb{C}\mid
p_{i}(s)=0,\,i=1,2,\ldots,n\right\\}$ (26)

and

$\bar{\mathbb{X}}\coloneqq\left\\{s\in\mathbb{C}\mid\hat{p}(s)=0,\,\lambda(s)\coloneqq
g_{1}(s)+\beta g_{2}(s),\,\beta\in[-2,2]\right\\},$ (27)

where $\hat{p}(s):=1-\lambda(s)h(s)$. The following two statements hold:

  1. 1.

$\mathbb{X}_{n}\subset\bar{\mathbb{X}}$.

  2. 2.

In the limit of $n\rightarrow\infty$, $\mathbb{X}_{n}$ is asymptotic to
$\bar{\mathbb{X}}$.

Lemma 1 is derived since the cosine function in the second term of Eq. (20)
varies continuously for $i$ in the range from $-1$ to $1$ by taking the limit
$n\rightarrow\infty$. Lemma 1 shows that the characteristic equations (18) for
$i=1,2,\ldots,n$ is asymptotic to the characteristic equation $\hat{p}(s)=0$
for all $\beta$ by $n\rightarrow\infty$, and thus it implies that the
stability of the circulant MC system for $n\rightarrow\infty$ can be analyzed
by computing the roots of the characteristic equation $\hat{p}(s)=0$ for all
$\beta$.

######  Remark 1.

A potential application of Theorem 1 is the analysis of biological pattern
formation such as Turing pattern formation. The mechanism of Turing pattern
formation can be elucidated by the spatial frequency-based analysis for a
reaction-diffusion (RD) equation [26]. An RD system can be decomposed into
multiple subsystems, which express the dynamics of Fourier components for
spatial frequencies $\xi$. It is known that if the subsystem with the spatial
frequency $\xi=0$ is stable and one or more of the subsystems with nonzero
spatial frequencies $\xi_{u}$ is unstable, a spatial pattern with a period
corresponding to the spatial frequency $\xi_{u}$ is formed when spatial white
noise is input. This is because only the dynamics of the unstable Fourier
component with $\xi_{u}$ do not converge to an equilibrium point, and thus the
molecular concentration appears periodically in space corresponding to the
spatial frequency $\xi_{u}$. In the same way, considering the circulant MC
system, if the subsystem $\hat{M}_{1}(s)$ is stable and one or more of the
subsystems $\hat{M}_{i}(s)\,(i=2,3,\ldots,n)$ is unstable, the spatially
periodical pattern would emerge. By developing the proposed method, the
parameter conditions for the Turing pattern formation with the circulant MC
system might be derived.

###  4.3 Graphical stability test

In the previous section, we have shown the method to reduce the problem of
computing the determinant of a large matrix to that of computing multiple non-
rational scalar equations for the stability analysis of the MIMO circulant MC
system $\mathcal{M}$. However, finding all the roots of a non-rational scalar
equation containing exponential functions is not easy since there are an
infinite number of roots in some cases. We here show one of the useful methods
to graphically analyze if all the roots of the characteristic equation (18)
lie in the OLHP of the complex plane. To this end, we first show the
definition of the Nyquist contour, which is used for the graphical stability
test.

######  Definition 4.

A contour is the Nyquist contour if it consists of the imaginary axis from
$0-j\infty$ to $0+j\infty$ and a semicircular arc with a radius
$R\rightarrow\infty$ starting at $0+j\infty$ and traveling clockwise to
$0-j\infty$.

######  Lemma 2.

Consider the characteristic equation (18) and let $P$ denote the number of the
roots of the equation $\mathrm{den}(h(s))=0$ whose real part is positive,
where $\mathrm{den}(\cdot)$ is the denominator of the function. All roots of
the characteristic equation (18) lie in the OLHP of the complex plane if and
only if the Nyquist plot of the open-loop transfer function
$-\lambda_{i}(j\omega)h(j\omega)$ makes $P$ anti-clockwise encirclements of
the point -1 and does not pass through the point -1 in the complex plane,
where the Nyquist plot is the trajectory of a function along the Nyquist
contour.

Lemma 2 is based on the Nyquist stability criterion [30]. The proof is in A.
In summary, the procedure of the proposed method to analyze the stability of
the circulant MC system is as follows. First, we apply Theorem 1 to decompose
the stability analysis of the circulant MC system $\mathcal{M}$ into that of
$n$ subsystems $\hat{\mathcal{M}}_{i}$. We then use Lemma 2 to analyze the
stability of the decomposed subsystem $\hat{\mathcal{M}}_{i}$.

##  5 Numerical example

In this section, we demonstrate the stability analysis for a specific example
of the circulant MC system. In particular, we use Theorem 1 and Lemma 2 to
analyze the stability of the circulant MC system $\mathcal{M}$ for showing
that the stability of system $\mathcal{M}$ depends on the communication length
$L$ and the diffusion coefficient $\mu$. We then show that the state of the
circulant MC system either converges to the spatially homogeneous equilibrium
state or transitions to another equilibrium state depending on the spatial
frequency of small perturbations added to the system.

###  5.1 Stability analysis for the circulant MC system with an activator-
repressor-diffuser genetic circuit

Figure 4: The circulant MC system with ARD genetic circuit.

We consider the circulant MC system illustrated in Fig. 4, where each
nanorobot has an activator-repressor-diffuser (ARD) genetic circuit [16]. The
function $\bm{f}(\bm{x}_{i})$ of the reaction system is

$\bm{f}(\bm{x}_{i})=\left[\begin{matrix}\displaystyle-\delta_{a}x_{i,1}(t)+\gamma_{a}\frac{x_{i,1}(t)^{2}}{K_{a}^{2}+x_{i,1}(t)^{2}}\frac{K_{r}^{2}}{K_{r}^{2}+x_{i,2}(t)^{2}}\\\
\displaystyle-\delta_{r}x_{i,2}(t)+\gamma_{r}\frac{x_{i,1}(t)^{2}}{K_{a}^{2}+x_{i,1}(t)^{2}}\frac{K_{d}^{2}}{K_{d}^{2}+x_{i,3}(t)^{2}}\\\
\displaystyle-\delta_{d}x_{i,3}(t)+\gamma_{d}\frac{K_{r}^{2}}{K_{r}^{2}+x_{i,2}(t)^{2}}+k(x_{i,4}(t)-x_{i,3}(t))\\\
k_{p}(x_{i,3}(t)-x_{i,4}(t))\end{matrix}\right],$ (28)

where $x_{i,1}(t)$, $x_{i,2}(t)$, $x_{i,3}(t)$, and $x_{i,4}(t)$ are the
concentrations of activator, repressor, signal molecule in the $i$-th
nanorobot, and signal molecule outside of the nanorobot, respectively. The
parameters $\delta_{l}$, $\gamma_{l}$, and $K_{l}$ are the degradation rate,
the production rate, and the Michaelis Menten constant of the $l$-th molecular
species, respectively. The parameter $k_{p}$ is the membrane transport rate.
The average length of the communication channel is $L$, and the diffusion
coefficient is $\mu$. The parameter values are shown in Table 1, whose order
of magnitude is consistent with widely used values for numerical simulations
in synthetic biology [16, 31, 32]. In what follows, we analyze the stability
of the circulant MC system around the spatially homogeneous equilibrium point
$\bm{x^{*}}=[5.09,4.39,8.08,8.08]\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}$$
for different values of the communication length $L$ and the diffusion
coefficient $\mu$.

Table 1: Parameter sets for the analysis of the model in Fig. 1, where each nanorobot has an ARD genetic circuit. Parameter |  | Value
---|---|---
Production rate of $x_{i,1}(t)$ | $\gamma_{1}$ | $2.5\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
Production rate of $x_{i,2}(t)$ | $\gamma_{2}$ | $2.68\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
Production rate of $x_{i,3}(t)$ | $\gamma_{3}$ | $0.9\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
Degradation rate of $x_{i,1}(t)$ | $\delta_{1}$ | $0.07\,\mathrm{min^{-1}}$
Degradation rate of $x_{i,2}(t)$ | $\delta_{2}$ | $0.07\,\mathrm{min^{-1}}$
Degradation rate of $x_{i,3}(t)$ | $\delta_{3}$ | $0.009\,\mathrm{min^{-1}}$
Dissociation constant of $x_{i,1}(t)$ | $K_{1}$ | $11\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}$$
Dissociation constant of $x_{i,2}(t)$ | $K_{2}$ | $9\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}$$
Dissociation constant of $x_{i,3}(t)$ | $K_{3}$ | $11\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}$$
Membrane transport rate | $k_{p}$ | $0.5\,$\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$

We consider $n=100$ nanorobots. The transfer function matrix of the MC channel
$G(s)$ is obtained as the $100\times 100$ circulant matrix, and the reaction
system is $H(s)=h(s)I$, where $h(s)=\bm{C}(sI-A)^{-1}\bm{B}$. The Jacobian
matrix $A$ around the spatially homogeneous equilibrium point
$\bm{x^{*}}=[5.09,4.39,8.08,8.08]\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}$$ is
obtained as

$A=\left[\begin{matrix}0.0453&-0.0312&0&0\\\ 0.0993&-0.0700&-0.0266&0\\\
0&-0.0636&-0.5900&0.5000\\\ 0&0&0.5000&-0.5000\end{matrix}\right],$ (29)

which leads to the transfer function

$h(s)=\frac{10s^{3}+6.1469s^{2}+0.128s+0.0003}{s^{4}+1.1147s^{3}+0.0701s^{2}+0.0003s}.$
(30)

The roots of the equation $\mathrm{den}(h(s))=0$ are computed as $-1.0480$,
$0.0018+0.0217j$, $0.0018-0.0217j$, and $-0.0704$. Since the equation
$\mathrm{den}(h(s))=0$ has two roots whose real part is positive, the reaction
system $h(s)$ is unstable. Using Theorem 1, the stability of the circulant MC
system can be analyzed by computing the roots of the characteristic equation
(18) for $i=1,2,\ldots,51$.

We then use Lemma 2 to analyze whether all the roots of the characteristic
equation (18) lie in the OLHP of the complex plane. Since the reaction system
$h(j\omega)$ has two roots with positive real part, namely $P=2$, if the
Nyquist plot of $-\lambda_{i}(j\omega)h(j\omega)$ makes two anti-clockwise
encirclements of the point -1 and does not pass through the point $-1$ in the
complex plane as the frequency $\omega$ increases for $i=1,2,\ldots,51$, the
circulant MC system $\mathcal{M}$ is locally asymptotically stable, otherwise
unstable.

Figure 5: (A) The parameter map for the stability of the circulant MC system
$\mathcal{M}$ for different communication lengths $L$ and diffusion
coefficients $\mu$. (B) The behavior of the concentration of the repressor
molecule in the reaction system $\Sigma_{i}$ around the spatially homogeneous
equilibrium point for $i=1,2,\ldots,100$ when
$\mu=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
and $L=70\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$, and (C) when
$L=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$.

Fig. 5 (A) shows the parameter map for the stability of the circulant MC
system for $n=100$ for different communication lengths $L$ and diffusion
coefficients $\mu$. Fig. 5 (B) shows the concentration behavior of the
repressor $x_{i,1}(t)$ in 100 nanorobots for
$L=70\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ and
$\mu=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
when a small perturbation $\bm{\eta}_{r}\in\mathbb{R}^{100}$ generated by
random values is added to
$\bm{\bar{x}}_{4}(0)\coloneqq[x_{1,4}(0),x_{2,4}(0),\ldots,x_{100,4}(0)]$,
whose $i$-th entry is
$x_{i,4}(0)=8.08\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}$$. The other
molecular concentrations $x_{i,1}(0)$, $x_{i,2}(0)$, and $x_{i,3}(0)$ are not
perturbed around $\bm{x}^{*}$ for all $i$. We can see that all molecular
concentrations converge to the spatially homogeneous equilibrium point
$\bm{x}^{*}$ corresponding to the analytical result in Fig. 5. On the other
hand, Fig. 5 (C) shows the concentration behavior of the repressor
$x_{i,1}(t)$ in 100 nanorobots for
$L=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ and
$\mu=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
when a small perturbation $\bm{\eta}_{r}$ inputs to $\bm{\bar{x}}_{4}(0)$. We
can see that the molecular concentrations transition to a different
equilibrium point than $\bm{x}^{*}$. Thus, Fig. 5 verifies the proposed
stability analysis method, which is helpful for the design of stable multi-
agent MC systems, where the nanorobot population can achieve a cooperative
behavior by synchronizing their state.

Figure 6: The Nyquist plot of the open-loop transfer function
$-\lambda_{i}(j\omega)h(j\omega)$ of the subsystems $\hat{\mathcal{M}}_{i}$
for $i=1,2$ with the communication length
$L=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ and the diffusion coefficient
$\mu=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$.
The cross mark represents the point $-1$.

As seen in Fig. 5, the circulant MC system becomes stable as the communication
length $L$ is greater. This result implies that the variation of the molecular
concentration from the spatially homogeneous equilibrium point induced by a
perturbation flows out into the large MC channels, thereby facilitating the
convergence of the variation of the molecular concentration to zero. Also, as
the diffusion coefficient $\mu$ increases, the MC system becomes unstable.
This is because, when the diffusion coefficient $\mu$ is large, the variation
of the molecular concentration can reach the neighboring nanorobots, which
amplifies the variation of the molecular concentration in the nanorobot. On
the other hand, for extremely small $\mu$, the variation of the molecular
concentration in the nanorobot cannot be released into the MC channel, and
thus the variation increases gradually, which leads to the unstable MC system.
These properties are consistent with those of the quorum sensing mechanism
[33], such as the property that a greater $L$ stabilizes the MC system,
meaning that the system with a lower density of nanorobots can easily
transition to another state.

###  5.2 Stability analysis based on spatial frequency

We next analyze the stability of the circulant MC system based on the spatial
frequency to show that the convergence of the state of the circulant MC system
to the spatially homogeneous equilibrium point depends on the Frequency
component of the small perturbation added to the system.

For demonstration purpose, we consider the circulant MC system with $n=3$
nanorobots illustrated in Fig. 4 with the parameter values in Table 1. Figure
6 depicts the Nyquist plot of the open-loop transfer function
$-\lambda_{i}(j\omega)h(j\omega)$ of the subsystem $\hat{\mathcal{M}}_{i}$ for
$i=1,2$ when the diffusion coefficient
$\mu=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
and the communication length $L=10\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$.
The figure shows that the Nyquist plot of $-h(j\omega)\lambda_{1}(j\omega)$
(red line) does not encircle the point $-1$ while the Nyquist plot of
$-h(j\omega)\lambda_{2}(j\omega)$ (blue line) encircles the point $-1$ twice.
Thus the subsystem $\hat{\mathcal{M}}_{1}$ is unstable while the subsystem
$\hat{\mathcal{M}}_{2}$ is stable, which leads to the circulant MC system
$\mathcal{M}$ unstable. This result implies that the state of the circulant MC
system converges to the spatially homogeneous equilibrium state although the
system $\mathcal{M}$ is unstable when the perturbation that does not include
the Fourier component with the spatial frequency $\xi_{1}$.

Figure 7: (A) (Left) The spatial distribution of the perturbation
$\bm{\eta}_{r}$, which inputs to $\bm{\tilde{y}}(0)$, and (Right) the time
series data of the distribution of the concentration variation of the signal
molecule $\tilde{c}_{i}(t,r)$ for all $i$. (B) (Left) The spatial distribution
of the perturbation $\bm{\eta}_{c}$, which inputs to $\bm{\tilde{y}}(0)$, and
(Right) the time series data of the distribution of the concentration
variation of the signal molecule $\tilde{c}_{i}(t,r)$ for all $i$.

Fig. 7 shows the time series data for the concentration variation of the
signal molecule around the spatially homogeneous equilibrium point
$x^{*}_{4}$, which is $\tilde{c}_{i}(t,r)=c_{i}(t,r)-x_{4}^{*}$, for all
$i=1,2,3$ when the perturbation $\bm{\eta}_{r}\in\mathbb{R}^{3}$ and the
perturbation $\bm{\eta}_{c}\in\mathbb{R}^{3}$, whose $i$-th entry is
$\eta_{ri}=\cos{(\xi_{2}iL)}$, respectively input to the concentration
variation
$\bm{\tilde{x}}_{4}(0)\coloneqq[\tilde{x}_{1,4}(0),\tilde{x}_{2,4}(0),\ldots,\tilde{x}_{100,4}(0)]$,
where the $i$-th entry of $\bm{\tilde{x}}_{4}(0)$ is
$\tilde{x}_{i,4}(0)\coloneqq x_{i,4}-x_{4}^{*}$. The other molecular
concentrations $x_{i,1}(0)$, $x_{i,2}(0)$, and $x_{i,3}(0)$ are not perturbed
around $\bm{x}^{*}$ for all $i$. When the perturbation $\bm{\eta}_{r}$ inputs,
the concentration variation $\tilde{c}_{i}(t,r)$ diverges as seen in Fig. 7
(A), implying that the molecular concentration transitions to another
equilibrium point. On the other hand, when the perturbation $\bm{\eta}_{c}$ is
added to the system, the concentration variation $\tilde{c}_{i}(t,r)$
converges to 0 for all $i$ as seen in Fig. 7 (B). These results arise because
the perturbation $\bm{\eta}_{r}$ includes the Fourier component with the
spatial frequency $\xi_{1}$ while the perturbation $\bm{\eta}_{c}$ does not
include the Fourier component with $\xi_{1}$, which corresponds to the
analytical result. Figures 6 and 7 confirm that the molecular concentration
either converges to the spatially homogeneous equilibrium point or transitions
to another state depending on the spatial frequency of the perturbation. This
phenomenon is very similar to the mechanism of Turing pattern formation,
suggesting that the proposed method can be used to analyze spatial frequency-
based biological phenomena such as Turing pattern formation.

##  6 Conclusion

In this paper, we have modeled large-scale multi-agent MC systems as circulant
MC systems using transfer functions based on a diffusion equation. We have
then proposed the method to analyze the stability for the circulant MC systems
by decomposing the transfer function matrix of the circulant MC system into
multiple single-input and single-output (SISO) systems. The proposed method
can analyze the response of multi-agent MC systems for the spatial pattern of
perturbation, which potentially leads to the analysis of spatial frequency-
based biological phenomena such as Turing pattern formation. Finally, we have
demonstrated the proposed method by analyzing the stability of a specific
large-scale multi-agent MC system.

##  Appendix A Proof of Lemma 2

For the proof of Lemma 2, we first introduce the Nyquist stability criterion
and then show the equivalence between the Nyquist stability criterion and
Lemma2. Consider the characteristic equation (18), and let $\hat{P}$ denote
the number of the roots with positive real part of the equation
$\mathrm{den}(h(s))\mathrm{den}(\lambda_{i}(s))=0$, where

$\mathrm{den}(\lambda_{i}(s))=\sqrt{\frac{\mu}{s}}\sinh{\left(\frac{L}{\sqrt{\mu}}\sqrt{s}\right)}.$
(31)

The Nyquist stability criterion [30] states that all roots of the
characteristic equation (18) lie in the OLHP of the complex plane if and only
if the Nyquist plot of the open-loop transfer function
$-\lambda_{i}(j\omega)h(j\omega)$ makes $\hat{P}$ anti-clockwise encirclements
of the point -1 and does not pass through the point -1 in the complex plane.

Next, we show $\hat{P}=P$. Let $\mathbb{Y}$ denote the solution set of the
equation $\mathrm{den}(h(s))\mathrm{den}(\lambda_{i}(s))=0$. The solution set
$\mathbb{Y}$ is equivalent to the set of the solution satisfying the equation
$\mathrm{den}(h(s))=0$ or the equation $\mathrm{den}(\lambda_{i}(s))=0$. Since
$\mathrm{den}(\lambda_{i}(s))=0$ has no roots, the solution set $\mathbb{Y}$
is equivalent to the solution set of the equation $\mathrm{den}(h(s))=0$.
Therefore, the number of the roots with positive real part of the equation
$\mathrm{den}(h(s))\mathrm{den}(\lambda_{i}(s))=0$ equals the number of the
roots with positive real part of the equation $\mathrm{den}(h(s))=0$, i.e.,
$\hat{P}=P$.

## References

  * [1] T. Suda and T. Nakano, “Molecular communication : a personal perspective,” IEEE Transactions on NanoBioscience, vol. 17, no. 4, pp. 424–432, 2018.
  * [2] D. Bi, A. Almpanis, A. Noel, Y. Deng, and R. Schober, “A Survey of Molecular Communication in Cell Biology: Establishing a New Hierarchy for Interdisciplinary Applications,” IEEE Communications Surveys and Tutorials, vol. 23, no. 3, pp. 1494–1545, 2021.
  * [3] C. A. Söldner, E. Socher, V. Jamali, W. Wicke, G. S. Member, A. Ahmadzadeh, H.-g. Breitinger, A. Burkovski, K. Castiglione, R. Schober, and H. Sticht, “A Survey of Biological Building Blocks for Synthetic Molecular Communication Systems,” IEEE Communications Surveys and Tutorials, vol. 22, no. 4, pp. 2765–2800, 2020.
  * [4] S. Lotter, L. Brand, V. Jamali, M. Schäfer, H. M. Loos, H. Unterweger, S. Greiner, J. Kirchner, C. Alexiou, D. Drummer, G. Fischer, A. Buettner, and R. Schober, “Experimental research in synthetic molecular communications – part ii,” IEEE Nanotechnology Magazine, vol. 17, no. 3, pp. 54–65, 2023.
  * [5] N. Farsad, H. B. Yilmaz, A. Eckford, C. B. Chae, and W. Guo, “A comprehensive survey of recent advancements in molecular communication,” IEEE Communications Surveys and Tutorials, vol. 18, no. 3, pp. 1887–1919, 2016.
  * [6] M. Femminella, G. Reali, and A. V. Vasilakos, “Molecular communications model for drug delivery,” IEEE Transactions on NanoBioscience, vol. 14, no. 7, pp. 935–945, 2015.
  * [7] W. Gao and J. Wang, “Synthetic micro/nanomotors in drug delivery,” Nanoscale, vol. 6, pp. 10486–10494, 2014.
  * [8] M. Pierobon and I. Akyildiz, “A physical end-to-end model for molecular communication in nanonetworks,” IEEE Journal on Selected Areas in Communications, vol. 28, no. 4, pp. 602–611, 2010.
  * [9] U. A. Chude-Okonkwo, R. Malekian, and B. T. Maharaj, “Diffusion-controlled interface kinetics-inclusive system-theoretic propagation models for molecular communication systems,” EURASIP Journal on Advances in Signal Processing, vol. 2015, no. 1, pp. 1–23, 2015.
  * [10] Y. Huang, F. Ji, Z. Wei, M. Wen, X. Chen, Y. Tang, and W. Guo, “Frequency Domain Analysis and Equalization for Molecular Communication,” IEEE Transactions on Signal Processing, vol. 69, pp. 1952–1967, 2021.
  * [11] S. Lotter, A. Ahmadzadeh, and R. Schober, “Channel modeling for synaptic molecular communication with re-uptake and reversible receptor binding,” in Proceedings of the International Conference on Communications, pp. 1–7, 2020.
  * [12] M. Schäfer, W. Wicke, W. Haselmayr, R. Rabenstein, and R. Schober, “Spherical diffusion model with semi-permeable boundary: A transfer function approach,” in Proceedings of the International Conference on Communications, pp. 1–7, 2020.
  * [13] B. C. Akdeniz, A. E. Pusane, and T. Tugcu, “2-d channel transfer function for molecular communication with an absorbing receiver,” in Proceedings of the IEEE International Conference on Intelligent Computer Communication and Processing, pp. 938–942, IEEE, 2017.
  * [14] S. Hara, T. Kotsuka, and Y. Hori, “Modeling and stability analysis for multi-agent molecular communication systems : a case study for two agents,” in Proceedings of the SICE Annual Conference, pp. 659–662, 2021.
  * [15] T. Kotsuka and Y. Hori, “Spatial Frequency-Based Characterization of Disturbance Rejection in Molecular Communication Systems,” IEEE Transactions on Molecular, Biological, and Multi-Scale Communications, vol. 8, no. 1, pp. 36–43, 2022.
  * [16] Y. Hori, H. Miyazako, S. Kumagai, and S. Hara, “Coordinated spatial pattern formation in biomolecular communication networks,” IEEE Transactions on Molecular, Biological, and Multi-Scale Communications, vol. 6, no. 2, pp. 111–121, 2015.
  * [17] J. Hsia, W. J. Holtz, D. C. Huang, M. Arcak, and M. M. Maharbiz, “A Feedback Quenched Oscillator Produces Turing Patterning with One Diffuser,” PLOS Computational Biology, vol. 8, no. 1, p. e1002331, 2012.
  * [18] K. Kashima, T. Ogawa, and T. Sakurai, “Selective pattern formation control: Spatial spectrum consensus and Turing instability approach,” Automatica, vol. 56, pp. 25–35, 2015.
  * [19] J. Qin, Q. Ma, Y. Shi, and L. Wang, “Recent Advances in Consensus of Multi-Agent Systems: A Brief Survey,” IEEE Transactions on Industrial Electronics, vol. 64, no. 6, pp. 4972–4983, 2017.
  * [20] Y. Cao, W. Yu, W. Ren, and G. Chen, “An overview of recent progress in the study of distributed multi-agent coordination,” IEEE Transactions on Industrial Informatics, vol. 9, no. 1, pp. 427–438, 2013.
  * [21] J. Fax and R. Murray, “Information Flow and Cooperative Control of Vehicle Formations,” IEEE Transactions on Automatic Control, vol. 49, no. 9, pp. 1465–1476, 2004.
  * [22] T. Kim, S. Hara, and Y. Hori, “Cooperative control of multi-agent dynamical systems in target-enclosing operations using cyclic pursuit strategy,” International Journal of Control, vol. 83, no. 10, pp. 2040–2052, 2010.
  * [23] R. Olfati-Saber, J. A. Fax, and R. M. Murray, “Consensus and Cooperation in Networked Multi-Agent Systems,” Proceedings of the IEEE, vol. 95, no. 1, pp. 215–233, 2007.
  * [24] T. Kotsuka and Y. Hori, “A control-theoretic model for bidirectional molecular communication systems,” IEEE Transactions on Molecular, Biological and Multi-Scale Communications, vol. 9, no. 3, pp. 274–285, 2023.
  * [25] S. Hara, T. Hayakawa, and H. Sugata, “LTI systems with generalized frequency variables: A unified framework for homogeneous multi-agent dynamical systems,” SICE Journal of Control, Measurement, and System Integration, vol. 2, no. 5, pp. 299–306, 2009.
  * [26] A. M. Turing, “The chemical basis of morphogenesis,” Philosophical Transactions of the Royal Society B: Biological Sciences, vol. 237, no. 641, pp. 37–72, 1952.
  * [27] T. S. Gardner, C. R. Cantor, and J. J. Collins, “Construction of a genetic toggle switch in Escherichia coli,” Nature, vol. 403, no. 6767, pp. 339–342, 2000.
  * [28] K. Zhou, J. Doyle, and K. Glover, Robust and optimal control.  Englewood Cliffs, N.J: Prentice Hall, 1996.
  * [29] K. R. Rao and P. Yip, The transform and data compression handbook.  The electrical engineering and signal processing series, Boca Raton, Fla: CRC Press, 2001.
  * [30] S. Skogestad and I. Postlethwaite, Multivariable Feedback Control: Analysis and Design.  Hoboken, NJ, USA: John Wiley & Sons, Inc., 2005.
  * [31] S. Basu, Y. Gerchman, C. H. Collins, F. H. Arnold, and R. Weiss, “A synthetic multicellular system for programmed pattern formation,” Nature, vol. 434, no. 7037, pp. 1130–1134, 2005.
  * [32] X. Li, J. Jin, X. Zhang, F. Xu, J. Zhong, Z. Yin, H. Qi, Z. Wang, and J. Shuai, “Quantifying the optimal strategy of population control of quorum sensing network in Escherichia coli,” npj Systems Biology and Applications, vol. 7, no. 1, p. 35, 2021.
  * [33] C. M. Waters and B. L. Bassler, “Quorum sensing: cell-to-cell communication in bacteria,” Annual Review of Cell and Developmental Biology, vol. 21, pp. 319–346, 2005.


#  The Shannon-McMillan Theorem Proves Convergence to Equiprobability of
Boltzmann’s Microstates

Arnaldo Spalvieri∗
∗Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di
Milano

###### Abstract

The paper shows that, for large number of particles and for distinguishable
and non-interacting identical particles, convergence to equiprobability of the
$W$ microstates of the famous Boltzmann-Planck entropy formula $S=k\log(W)$ is
proved by the Shannon-McMillan theorem, a cornerstone of information theory.
This result further strengthens the link between information theory and
statistical mechanics.

##  I Introduction

Entropy plays a key role in communication and information theory and in
thermodynamics. Starting from the concept of entropy, many authors built in
the past bridges between information and communication theory and
thermodynamics, so that all the textbooks of statistical mechanics and thermal
physics have chapters or sections where the basics of information theory are
introduced. A recent survey about the links between entropy in communications
and in thermodynamics is [1]. Making a comprehensive review of the
bibliography on this topic is out of the scope of the present paper. We limit
ourselves to mention the pioneering work of Jaynes [2] and the large attention
that has been attracted in the past by the connections between information
theoretic inequalities and the irreversibility of certain thermodynamical
processes, see [3], Chapter 4 of [4], [5].

The aim of this paper is to put the information-theoretic concept of typical
set inside the framework of statistical mechanics, thus strengthening the link
between the two disciplines. Typicality is not a new word in statistical
mechanics. It has been introduced in [6] for system that are not at the
equilibrium to mean that these systems, in their spontaneous evolution,
typically tend to the maximum entropy. The difference with [6] is that (non)
equilibrium is not our concern and that the typicality we refer to is the one
that is technically established in information theory by the Shannon-McMillan
theorem. Specifically, this paper shows that the Shannon-McMillan theorem
proves that the microstates of the famous Boltzmann-Plank entropy formula are
equiprobable. This result is surprising because, since the times of Boltzmann,
the entire scientific community has postulated equiprobability of microstates.
Today, all the most influential textbooks of statistical mechanics and thermal
physics unanimously say that microstates’ equiprobability is the central
postulate of statistical mechanics, see for instance 4.2 of [7], chapter 6 of
[8], 16.1 of [9], 1.1 of [10]. In 3.4 of [11], the author not only postulates
equiprobability, but also expresses, in the passage that we hereafter quote,
strong scepticism about the possibility of proving it.

It is, unfortunately, impossible to prove with mathematical rigor that the
principle of equal a priori probabilities applies to many-particle systems.
Over the years, many people have attempted this proof, and all have failed.
Not surprisingly, therefore, statistical mechanics was greeted with a great
deal of scepticism when it was first proposed in the late 1800s. One of the
its main proponents, Ludwig Boltzmann, became so discouraged by all of the
criticism that he eventually committed suicide. Nowadays, statistical
mechanics is completely accepted into the cannon of physics - quite simply
because it works.

The outline of the paper is as follows. In Section II we introduce our model
of the thermodynamical system, where the only assumption that we make is that
the system is made by a large number of non-interacting and identical
particles. We say our model because it is slightly different from the commonly
accepted one even if, in the thermodynamical limit of infinite number of
particles, our model becomes compatible with the standard one. Section III
shows that, with the further assumption of distinguishable particles, the set
of microstates accessible to the thermodynamical system is the information-
theoretic typical set of the Shannon-McMillan theorem. A crucial property of
the typical set is that the probability distribution of its element converges
to uniformity, therefore the identification between the typical set and the
set of accessible microstates proves convergence to equiprobability of the
accessible microstates. Finally, in Section IV we recap the main points of the
paper and draw the conclusions.

##  II Notation and System Model

Let the uppercase calligraphic character, e.g. ${\cal X}$, denote discrete
random variables, let the set $\\{{\cal X}\\}$ be the support set of the
random variable, e.g.

$\\{{\cal X}\\}=\\{x_{0},x_{1},\cdots\\}$

and let $|{\cal X}|$ be the number of elements of the support set. The
probability $Pr({\cal X}=x)$ of the event ${\cal X}=x$ is denoted with the
shorthand $p(x)$,

$Pr({\cal X}=x)=p(x),\ \ \forall\ x\in\\{{\cal X}\\}.$

The probability distribution of ${\cal X}$ is the deterministic function
$p(\cdot)$ of the random ${\cal X}$ with

$\\{p({\cal X})\\}=\\{p(x_{0}),p(x_{1}),\cdots\\}$

and

$Pr(p({\cal X})=p(x))=p(x),\ \ \forall\ \ x\in\\{{\cal X}\\}.$

Consider a thermodynamical system made by a fixed number $N$ of non-
interacting identical particles with one degree of freedom each. The
thermodynamical system is modelled here by the random vector

${\cal E}_{1}^{N}=({\cal E}_{1},{\cal E}_{2},\cdots,{\cal E}_{N})$

of the energies of the $N$ particles. We assume that energy is quantized and
that the random entries of the random vector are independent (non interacting)
identically distributed (identical) (i.i.d.) random variables (particles),
that is

$p({\cal E}_{1}^{N})=\prod_{i=1}^{N}p({\cal E}_{i}),$ (1)

with

$p({\cal E}_{1})=p({\cal E}_{2})=\cdots=p({\cal E}_{N}).$

Since particles are i.i.d., one of them represents all. The random energy of
the generic particle, that we call ${\cal E}$, takes its values in the set
$\\{{\cal E}\\}$ of the allowed energy levels

$\\{{\cal E}\\}=\\{\epsilon_{0},\epsilon_{1},\cdots\\}.$ (2)

The expected energy of the particle is

$E=\sum_{\epsilon\in\\{{\cal E}\\}}\epsilon p(\epsilon).$

The total random energy of the system is

$\sum_{i=1}^{N}{\cal E}_{i}.$

By the i.i.d. assumption and by the law of large numbers, when
$N\rightarrow\infty$, the system’s random energy divided by the number of
particles, that is, the random mean energy, tends to the deterministic
expected energy:

$\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{i=1}^{N}{\cal E}_{i}=E.$ (3)

Convergence of the limit can be in various senses. In this paper, we will
focus our attention on the weak law of large numbers (Khinchin law), hence we
consider convergence of the limit (3) in probability, meaning that, for every
$\eta>0$,

$\lim_{N\rightarrow\infty}Pr\left(\left|\frac{1}{N}\sum_{i=1}^{N}{\cal
E}_{i}-E\right|>\eta\right)=0.$ (4)

The limit in (4) is interpreted by saying that, for every $\eta$, if $N$ is
sufficiently large, then the total energy of the system lies in a narrow
interval $\pm N\eta$ around $NE$ with high probability.

A comment is in order about the approach pursued in this paper, that is based
on the expected energy per particle, in comparison to the classical approach
of thermal physics and statistical mechanics where the total system energy is
considered. The assumption that system’s total energy lies in a narrow
interval is the basis of the classical microcanonical ensemble approach, see
for instance chapters 1 and 2 of [10]. In the microcanonical ensemble
approach, the system is assumed to have a large and fixed number of particles
and to be closed. Since $N$ is large and fixed, basically there is no
difference between the expected energy per particle and the total energy
divided by the number of particles. When the number of particles is fixed but
not so large, the classical approach in statistical mechanics and thermal
physics is that of the canonical ensemble, see for instance chapter 3 of [10].
In the canonical ensemble approach, the thermodynamical system is at the
thermal equilibrium at a temperature $T$ with a heat bath and interacts with
it to maintain the equilibrium condition. One of the consequences of the
interaction between system and heat bath is that the total energy of the
system fluctuates, hence it becomes substantially random. Due to the random
nature of total system’s energy, in the canonical ensemble approach it can be
no more used to characterize the system, which is characterized instead by the
temperature of the heat bath. While system’s total energy becomes
substantially random and must be abandoned in favor of the temperature of the
heat bath, the expected energy, which is deterministic, can be related in a
deterministic way to the temperature of the heat bath by the appropriate
energy-temperature relation. For instance, the high temperature approximation
for wholly kinetic energy is

$E=\frac{kT}{2},$

where $T$ is in Kelvin degrees and $k$ is Boltzmann’s constant. Hence,
although in the following we will limit ourselves to $N\rightarrow\infty$, we
see that the expected energy per particle can be used to characterize both the
microcanonical and the canonical ensembles.

##  III Equiprobability of Accessible Microstates

When we become aware of the result of a random experiment, in our case, the
energy level of a particle, the surprise that we experience is the following
deterministic function of the random result:

$\displaystyle-\log(p({\cal E})),$

where $\log(x)$ indicates the logarithm of $x$ and the base of the logarithm
depends on the context. In information theory the base is 2 while in physics
the base is Euler’s constant $e$. Being a function of random variable, the
surprise is itself a random variable. Its expectation is Shannon’s entropy,
that we call $H$:

$H=-\sum_{\epsilon\in\\{{\cal E}\\}}p(\epsilon)\log(p(\epsilon)).$ (5)

If distinguishable particles, or indexed particles, are considered, an outcome
of ${\cal E}_{1}^{N}$ is an energy microstate, or, simply, microstate. For
instance, a microstate is

${\cal E}_{1}^{N}:\ {\cal E}_{1}=\epsilon_{5},\ {\cal E}_{2}=\epsilon_{31},\
\cdots,\ {\cal E}_{N}=\epsilon_{18},$ (6)

meaning that the energy of the first particle is $\epsilon_{5}$, the energy of
the second particle is $\epsilon_{31}$, $\cdots$, the energy of the $N$-th
particle is $\epsilon_{18}$. This definition of microstate matches the
commonly accepted one, see e.g. 1.1 of [10], where the microstates of the
system are the independent solutions of the Schröedinger equation of the
system whose eigenvalue is the total energy of the system. Due to the
assumption of independency between particles, the system decouples into $N$
independent equations whose individual solutions are the wave functions whose
energy eigenvalues are the energy levels (2), the sum of the individual energy
eigenvalues being equal to system’s total energy. By now we settle for the
definition of microstate, we will return soon on the concept of accessible
microstate and on the energy constraint that the system must comply with.

When we become aware of the microstate visited by the system, the surprise
that we experience is

$\displaystyle-\log(p({\cal E}_{1}^{N})).$

The expectation of the above surprise is the Shannon entropy of the system,
hence the Shannon entropy of the random vector ${\cal E}_{1}^{N}$:

$-\sum_{\epsilon_{1}^{N}\in\\{{\cal
E}_{1}^{N}\\}}p(\epsilon_{1}^{N})\log(p(\epsilon_{1}^{N})).$ (7)

In statistical mechanics, (7) is called Gibbs entropy. By the i.i.d.
assumption (1), the entropy (7) is equal to $n$ times the Shannon entropy of
one particle:

$-\sum_{\epsilon_{1}^{N}\in\\{{\cal
E}_{1}^{N}\\}}p(\epsilon_{1}^{N})\log(p(\epsilon_{1}^{N}))=NH,$

see 2.6.6 of [4]. One immediately recognizes that

$\displaystyle\lim_{N\rightarrow\infty}-\frac{1}{N}\log(p({\cal E}_{1}^{N}))$
$\displaystyle=\lim_{N\rightarrow\infty}-\frac{1}{N}\sum_{i=1}^{N}\log(p({\cal
E}_{i}))$ (8) $\displaystyle=H,$ (9)

where (8) is the assumption of independency between particles and (9) follows
from the assumption of identically distributed random variables and from the
law of large numbers. The equality between the leftmost term and the rightmost
one in the above two equalities shows that, as $N\rightarrow\infty$, the
distribution of the microstates that the system can visit converges to the
uniform distribution, because $H$ is a deterministic and fixed quantity,
independent of the specific microstate ${\cal E}_{1}^{N}$. The convergence of
the limit (8) to (9) is in various senses. Convergence in probability, that is

$\lim_{N\rightarrow\infty}Pr\left(\left|\frac{1}{N}\log(p({\cal
E}_{1}^{N})))+H\right|>\eta\right)=0,$ (10)

is referred to as the Asymptotic Equipartition Property (AEP) in chapter 3 of
[4],111Here it is the probability that is equally partitioned, not the energy,
hence this AEP has nothing to do with the classical energy equipartition
property. while many authors refer to (10) as to the Shannon-McMillan theorem.
Almost everywhere convergence, that strengthens convergence in probability, is
known as the Shannon-McMillan-Breiman theorem [12]. Exactly as it happens with
total system’s energy, here, for sufficiently large but finite $N$, the
surprise about system’s microstate lies in a narrow range $\pm N\eta$ around
$NH$ with high probability and, exactly as it happens with the expected energy
per particle $E$, $H$, which is per particle, can be used to characterize the
entropy both in the microcanonical and the in canonical ensemble approaches
even if, in what follows, we will limit ourselves to $N\rightarrow\infty$.

###  III-A Typical Set and Accessible Microstates

Starting from (10) we arrive to the definition of typical set, which is the
subset $\\{{\cal T}_{1}^{N}(\eta)\\}$ of $\\{{\cal E}_{1}^{N}\\}$ made by the
vectors ${\cal E}_{1}^{N}$ whose probability, for any $\eta$ and for
sufficiently large $N$, is in the narrow range

$e^{-N(H+\eta)}\leq p({\cal E}_{1}^{N})\leq e^{-N(H-\eta)}.$ (11)

The properties of the typical set are that the number of its elements is in
the narrow range

$e^{N(H+\eta)}\geq|{\cal T}_{1}^{N}(\eta)|\geq(1-\eta)e^{N(H-\eta)}$ (12)

and that the probability that the system visits one element of the typical set
is

$Pr({\cal E}_{1}^{N}\in\\{{\cal T}_{1}^{N}(\eta)\\})>1-\eta,$ (13)

see again [4]. Formulas (11)-(13) show that with high probability the system
visits, or, in the language of statistical mechanics, makes access to, one
among $|{\cal T}_{1}^{N}(\eta)|$ microstates. Also, the distribution of
microstates converges to the uniform one as $N\rightarrow\infty$, because
$|{\cal T}_{1}^{N}(\eta)|$ is basically the inverse of $p({\cal E}_{1}^{N})$,
which is constrained to lie in the narrow range (11). Note that convergence to
equiprobability of accessible microstates is a consequence of the i.i.d.
assumption and of $N\rightarrow\infty$, hence, in the end, of the law of large
numbers, not a postulate.

The AEP divides the support set of energy microstates into two subsets: the
typical set, that is the set of accessible microstates, and its complement,
that is the set of microstates that the system cannot access. The division can
be appreciated by observing that

$|{\cal E}|-e^{H}\geq 0,$

with equality when the energy levels that the single particle can visit are
equiprobable, hence

$|{\cal E}_{1}^{N}|=|{\cal E}|^{N}\geq e^{NH}.$ (14)

Writing, with some abuse of mathematics, (12) as

$|{\cal T}_{1}^{N}(\eta)|\approx e^{NH},$

we see that, when inequality (14) fits with equality, the size of the typical
set is the entire support set made of $|{\cal E}|^{N}$ elements, while, when
$|{\cal E}|>e^{H}$, the number of elements of the typical set is lower than
$|{\cal E}|^{N}$, leading to the division of the support set into the two
subsets mentioned above. In this case, the system will visit with high
probability the typical set and will visit with vanishingly small probability
the complement of the typical set to the support set.

We now return on accessibility and on the energy constraint to comment again
on the key difference between our approach and the standard one. In the
standard approach, for a microstate to be accessible, the energy eigenvalue of
the wave function that solves system’s Schröedinger equation of the system
must be equal to the constraint of total energy that is imposed on the system,
while here it is the probability distribution $\\{p({\cal E})\\}$ that must
have a prescribed expected energy per particle. Compatibility between our
approach and the standard one is guaranteed because, when
$N\rightarrow\infty$, the energy eigenvalue of the system divided by the
number of particles is forced by the law of large numbers to be equal to the
expected energy of one particle, see (3).

###  III-B Boltzmann-Plank entropy

The equality between the Boltzmann-Plank entropy and $N$ times the Shannon
entropy of one particle is widely accepted in textbooks of statistical
mechanics and thermal physics, see e.g. eqn. 18.14 of [9]. However, while in
the standard approach both the Boltzmann-Plank entropy and the Shannon entropy
are purely deterministic quantities, in our approach the Boltzmann-Planck
entropy is the random quantity that we are going to introduce in (18), that,
divided by the number of particles, becomes equal to the deterministic Shannon
entropy when the number of particles tends to infinity. We hereafter discuss
the key passage, the passage where our model becomes compatible with the
standard combinatorial approach to the evaluation of the number of accessible
microstates.

Let ${\cal N}(\epsilon)$ be the random number of particles that visit the
generic energy level $\epsilon$, with

$\sum_{\epsilon\in\\{{\cal E}\\}}{\cal N}(\epsilon)=N.$ (15)

Compatibility of the accessible microstates with constraints imposed on the
system will be guaranteed by imposing compatibility of the set

$\\{{\cal N}(\epsilon_{0}),{\cal N}(\epsilon_{1}),\cdots\\}$ (16)

with the distribution $\\{p({\cal E})\\}$. The random number of microstates is
the multinomial coefficient

${\cal W}=\frac{N!}{\prod_{\epsilon\in\\{{\cal E}\\}}{\cal N}(\epsilon)!}$
(17)

and the random Boltzmann-Plank entropy is

${\cal S}=k\log({\cal W})$ (18)

Randomness of the set (16) and, as a consequence, of the number of microstates
(17) and of the Boltzmann-Plank entropy (18), marks the difference between our
approach and the standard, deterministic, one. Exactly as in the case of the
energy, once again the law of large numbers bridges the gap between randomness
and determinism, thus making our random approach compatible with the
deterministic standard approach, hence the key passage is222Specifically, step
by step one has  $\displaystyle\lim_{N\rightarrow\infty}\frac{1}{kN}{\cal
S}=\lim_{N\rightarrow\infty}\frac{1}{N}\log({\cal W})$
$\displaystyle=\lim_{N\rightarrow\infty}\frac{1}{N}\left(\log(N!)-\sum_{{\epsilon\in\\{{\cal
E}\\}}}\log({\cal N}(\epsilon)!)\right)$
$\displaystyle=\lim_{N\rightarrow\infty}\frac{1}{N}\left(N\log\left(\frac{N}{e}\right)-\sum_{{{\epsilon\in\\{{\cal
E}\\}}}}{\cal N}(\epsilon)\log\left(\frac{{\cal
N}(\epsilon)}{e}\right)\right)$ (19)
$\displaystyle=\lim_{N\rightarrow\infty}\frac{1}{N}\left(N\log(N)-\sum_{{{\epsilon\in\\{{\cal
E}\\}}}}{\cal N}(\epsilon)\log({\cal N}(\epsilon))\right)$ (20)
$\displaystyle=\lim_{N\rightarrow\infty}\sum_{{{\epsilon\in\\{{\cal
E}\\}}}}\frac{{\cal N}(\epsilon)}{N}\left(\log(N)-\log({\cal
N}(\epsilon))\right)$ (21)
$\displaystyle=\lim_{N\rightarrow\infty}-\sum_{{{\epsilon\in\\{{\cal
E}\\}}}}\frac{{\cal N}(\epsilon)}{N}\log\left(\frac{{\cal
N}(\epsilon)}{N}\right)$ $\displaystyle=-\sum_{\epsilon\in\\{{\cal
E}\\}}p(\epsilon)\log(p(\epsilon))$ (22) $\displaystyle=H,$ where in (20) and
(21) we use (15), (22) is the law of large numbers and (19) is Stirling’s
formula that, in the big ${\cal O}$ notation, is
$\log(N!)=N\log\left(\frac{N}{e}\right)+{\cal O}(\log(N)),$ and we put in (19)
$\lim_{N\rightarrow\infty}\frac{{\cal
O}(\log(N))}{N}=\lim_{N\rightarrow\infty}\frac{{\cal O}(\log({\cal
N}(\epsilon)))}{N}=0,\ \ \forall\ \epsilon\in\\{{\cal E}\\}.$ The law of large
number is also invoked in [13] in a passage that is basically the same as our
(22) to show the equivalence of the phase space formulation of Gibbs and
Boltzmann-Plank entropies.

$\lim_{N\rightarrow\infty}\frac{1}{N}\log({\cal W})=H.$ (23)

From (12) and (23) we conclude that, for every small $\eta>0$,

$\lim_{N\rightarrow\infty}\frac{1}{N}\log({\cal
W})=\lim_{N\rightarrow\infty}\frac{1}{N}\log(|{\cal T}_{1}^{N}(\eta)|).$ (24)

Therefore, not only our definition of energy microstates exemplified in (6) is
the commonly accepted one, but also the logarithm of the number of accessible
microstates that we obtain by our approach is, for $N\rightarrow\infty$, equal
to the logarithm of the standard combinatorial definition (17) of accessible
microstates.

##  IV Conclusion

The paper has demonstrated that there is no need of assuming or postulating
equiprobability of accessible microstates, because convergence of microstates’
distribution to uniformity is a consequence of the large number of particles,
of the i.i.d. assumption and of the assumption of distinguishable particles.
Specifically, equations (11)-(13) show that convergence to equiprobability is
a consequence of the AEP and equation (24) shows that the microstates we are
dealing with are the microstates that are counted by the standard
combinatorial approach.

It is worth observing that equiprobability has nothing to do with the thermal
equilibrium, hence with entropy maximization. Actually, the only conditions
for convergence to equiprobability are the assumptions that we made, which do
not include thermal equilibrium and/or entropy maximization, therefore the
distribution converges to the uniform one also if the system is not at the
thermal equilibrium. What can change between a system at the equilibrium and a
system that is not at the equilibrium when both are subject to the same
constraints is the probability distribution of energy levels and, with it,
Shannon’s entropy and, with it again, the number of accessible microstates,
which, however, tend to equiprobability. Since the accessible microstates
remain virtually equiprobable also in systems that are not at the equilibrium,
the entropy of these systems is still given by the Boltzmann-Plank entropy
formula, provided that the number of accessible microstates is conveniently
expressed through the Shannon entropy of the specific (non-equilibrium)
probability distribution of the energy of the individual particle by (5),
(12).

When $N$ is not large, e.g. $N=1$, microstates are in general not
equiprobable, independently of the entropy of the system.

To summarize, the main result of this paper is:

  * •

for large $N$ and i.i.d. particles, the accessible microstates are always
virtually equiprobable, also for systems that are not at the thermal
equilibrium;

  * •

for small $N$, the accessible microstates are never equiprobable, also for
systems that are at the thermal equilibrium with a heat bath.

The important case of indistinguishable particles is not treated in the paper
and is left to future research.

## References

  * [1] Martinelli, M. ”Photons, Bits and Entropy: From Planck to Shannon at the Roots of the Information Age.” Entropy, 2017, 19, 341.
  * [2] Jaynes, E. T. ”Information theory and statistical mechanics.” Physical review, 1957, 4, 620–630.
  * [3] Merhav, N. ”Physics of the Shannon Limits,” IEEE Trans. on Inform. Theory, 2010, 9, 4274–4285.
  * [4] Cover, T. M.; and Thomas, J. A. ”Elements of information theory, 2nd edition.” Wiley series in telecommunications and signal processing, 2006.
  * [5] Landauer, R. ”Irreversibility and heat generation in the computing process.” IBM journal of research and development 5.3 (1961): 183-191.
  * [6] Goldstein, Sheldon, and Joel L. Lebowitz. ”On the (Boltzmann) entropy of non-equilibrium systems,” 2004, Physica D: Nonlinear Phenomena 193.1-4, 53–66.
  * [7] Kardar, M. ”Statistical Physics of Particles.” Cambridge University Press, 2007.
  * [8] Huang, K. ”Statistical Mechanics.” Wiley, 1987.
  * [9] Sekerka, R. F. ”Thermal Physics: Thermodynamics and Statistical Mechanics for Scientists and Engineers.” Elsevier, 2015.
  * [10] Pathria, R. K.; Beale, P. D. ”Statistical Mechanics, 3rd edition.” Elsevier, 2011.
  * [11] Fitzpatrick, R. ”Thermodynamics and Statistical Mechanics.” World Scientific, 2020.
  * [12] Breiman, L. ”The individual ergodic theorem of information theory.” The Annals of Mathematical Statistics 28.3 (1957): 809-811.
  * [13] Frigg, R.; Werndl, C. ”A guide for the perplexed.” Probabilities in physics 2 (2011): 115.