# Stability analysis for circulant structured multi-agent molecular
communication systems

Taishi Kotsuka, Yutaka Hori This work was supported in part by JSPS KAKENHI
Grant Numbers JP21H05889 and JP22J10554, and in part by JST SPRING Grant
Number JPMJSP2123.T. Kotsuka and Y. Hori are with the Department of Applied
Physics and Physico-Informatics, Keio University, Kanagawa 223-8522 Japan.
Correspondence should be addressed to Y. Hori (email<EMAIL_ADDRESS>yhori@appi.keio.ac.jp).

###### Abstract

In this paper, we introduce the system theoretic model for the multi-agent MC
systems represented by multi-input and multi-output (MIMO) systems using the
transfer functions, and then propose a method to analyze the stability for the
special case of the circulant structured multi-agent MC systems. The proposed
method decomposes the MIMO MC system into multiple single-input and single-
output (SISO) systems, which facilitates to analyze of the stability of the
large-scale multi-agent MC system. Finally, we demonstrate the proposed method
to analyze the stability of a specific MC system.

###### Index Terms:

Molecular communication, Feedback control, Modeling, Biological control
systems, Diffusion equation.

##  I Introduction

Recently, many studies for molecular communication (MC) systems were conducted
to control the behavior of dispersed nanorobot populations toward engineering
applications such as targeted drug delivery [1, 2, 3, 4, 5]. In general MC
systems, the dispersed nanorobots transfer signals to each other via MC
channels to achieve cooperative performance, and the fundamental properties of
such MC channels were analyzed based on diffusion equation [6, 7, 8, 9]. One
of the key roles of MC in the control of the multi-agent system is to
stabilize and synchronize the reactions inside the individual nanorobot. To
control the behavior of nanorobot populations, the dynamical reaction system
in each nanorobot needs to be stabilized via the MC channels. To this end, it
is crucial to model the multi-agent MC systems that incorporates the dynamics
of reactions in the nanorobots and develop tools for stability analysis.

There are several studies of multi-agent MC systems considering the reaction
in the nanorobots based on the reaction-diffusion (RD) equation [10, 11, 12].
However, RD-based MC models assume that the distance between the nanorobots
are close enough to ignore disruption of the signal by the MC channel, and
thus, the model is not suitable for the analysis and design of the MC systems
composed of a population of distributed nanorobots for applications in vivo
such as drug delivery. This limitation was somewhat relaxed by introducing a
system theoretic model of multi-agent dynamical nanorobots that incorporates a
detailed model of reactions and diffusion by the authors’ group [13]. However,
the previously proposed stability analysis method [13] was applicable only to
the MC systems with two nanorobots. This motivates us to develop a more
versatile stability analysis method for large-scale multi-agent MC systems for
broader engineering applications of MC systems.

In this paper, we introduce a system theoretic model of multi-agent MC systems
that incorporates the transfer functions of bidirectional MC channels and
local reactions inside nanorobots, and then, propose a stability analysis
method for a special case of the circulant structured homogeneous multi-agent
MC systems. The proposed method decomposes the transfer function matrix of the
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 [14]. This decomposition allows for the stability analysis
of the large-scale multi-agent MC system by combining simple analysis
techniques for SISO systems. Finally, we demonstrate the proposed method to
analyze the stability of a specific MC system.

##  II Multi-agent MC systems

We consider one-dimensional multi-agent MC systems consisting of $n$
nanorobots and $(n+1)$ MC channels shown in Fig. 1. Let $\Sigma_{i}$ denote
the $i$-th reaction system that captures the reaction in the $i$-th nanorobot
and the diffusion of the signal molecules, and $\Gamma_{i}=[0,L_{i}]$ denote
the $i$-th MC channel between the reaction systems or between the reaction
system and the left/right boundaries of the MC system with $L_{i}$ being the
communication distance.

Figure 1: The multi-agent MC system

We denote the location and the concentration of the signal molecules in the MC
channel $\Gamma_{i}$ by $r_{i}$ and $c_{i}(t,r_{i})$, respectively. The
dynamics of the concentration $c_{i}(t,r_{i})$ can be modeled by the diffusion
equation

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

where $\mu$ is the diffusion coefficient. The boundary conditions at $r_{i}=0$
and $r_{i}=L_{i}$ are Dirichlet boundary conditions as $c_{i}(t,0)=y_{i-1}(t)$
and $c(t,L_{i})=y_{i}(t)$ for $i=2,3,\cdots,n$, respectively. Dirichlet or
Neumann boundary conditions can be applied to the left end boundary of
$\Gamma_{1}$ and the right end boundary of $\Gamma_{n+1}$ depending on the
situations. The initial condition is

$c_{i}(0,r_{i})=0;\ \forall r_{i}\in[0,L_{i}].$ (2)

For many practical examples, the $i$-th reaction system $\Sigma_{i}$ can be
modeled by non-linear state-space models as

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

where the state $\bm{x}(t)\in\mathbb{R}^{m}$ is the concentrations of the
molecules associated with reactions occurring inside and outside of the
nanorobot, and the $m$-th entry of $\bm{x}(t)$ represents the concentration of
the signal molecule outside of the nanorobot. The vector function
$\bm{f}(\cdot)$ represents the dynamics of the reactions in the nanorobot and
the membrane transport, $\bm{B}=[0,\cdots,0,1]^{T}\in\mathbb{R}^{m}$ and
$\bm{C}=[0,\cdots,0,1]\in\mathbb{R}^{m}$ are the input and 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}$
represented by

$w_{i}(t)=\frac{\mu}{\Delta r_{i}}\left(\frac{\partial c_{i+1}(t,0)}{\partial
r_{i+1}}-\frac{\partial c_{i}(t,L_{i})}{\partial r_{i}}\right),$ (4)

where $\Delta r_{i}$ 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 .

The dynamical multi-agent MC system can synchronously converge to a spatially
homogeneous equilibrium point
$\bm{x}_{1}(t)=\bm{x}_{2}(t)=\cdots=\bm{x}_{n}(t)$ at steady state if the
reactions are appropriately designed. In what follows, we are interested in
the local stability analysis of the multi-agent MC system. Specifically, we
first express the multi-agent 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 (3). We then show that the circulant
structured MIMO MC system can be decomposed into multiple single-input single-
output (SISO) systems under certain conditions to facilitate the stability
analysis.

##  III Stability analysis for multi-agent MC systems

In this section, we first introduce the MIMO representation of the MC system.
We then show that the stability analysis of the $n$-dimensional MIMO MC system
can be reduced to that of $n$ SISO systems for a specific class of the MC
channel.

###  III-A System theoretic model

The MIMO representation of the MC system consisting of the MC channels
$\Gamma_{i}$ and the reaction systems $\Sigma_{i}$ shown in Fig. 2 (A) can be
expressed as

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

where $\bm{Y}(s):=[Y_{1}(s),Y_{2}(s),\cdots,Y_{n}(s)]^{T}$ and
$\bm{W}(s):=[W_{1}(s),W_{2}(s),\cdots,W_{n}(s)]^{T}$ with $Y_{i}(s)$ and
$W_{i}(s)$ being the Laplace transform of $y_{i}(t)$ and $w_{i}(t)$ as
$Y_{i}(s):=\mathcal{L}[y_{i}(t)]$ and $W_{i}(s):=\mathcal{L}[w_{i}(t)]$. The
transfer function matrix $H(s)$ is the diagonal matrix representing the
reaction systems, where the $(i,i)$-th entry is $h(s)=\bm{C}(sI-A)^{-1}\bm{B}$
and $A$ is the Jacobian matrix of $\bm{f(\bm{x})}$ at the homogeneous
equilibrium. The transfer function matrix $G(s)$ is the symmetric matrix
representing the MC channels defined by

$G(s)=\left[\begin{matrix}g_{1}^{b}(s)&g_{1}^{s}(s)&0&\cdots&0&g_{p}^{s}(s)\\\
g_{1}^{s}(s)&g_{2}^{t}(s)&g_{2}^{s}(s)&\cdots&0&0\\\
0&g_{2}^{s}(s)&\ddots&\ddots&0&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&0\\\
0&\ddots&\ddots&\ddots&\ddots&g_{n-1}^{s}(s)\\\
g^{s}_{p}(s)&0&\cdots&0&g^{s}_{n-1}(s)&g^{b}_{n}(s)\end{matrix}\right],$ (6)

where

$\displaystyle g^{t}_{i}(s)$ $\displaystyle=$
$\displaystyle-\hat{s}\left(\frac{1}{\tanh{\left(L_{i}\hat{s}\right)}}+\frac{1}{\tanh{\left(L_{i+1}\hat{s}\right)}}\right),$
(7) $\displaystyle g^{s}_{i}(s)$ $\displaystyle=$
$\displaystyle\frac{\hat{s}}{\sinh{\left(L_{i+1}\hat{s}\right)}}$ (8)

with $\hat{s}=\sqrt{s/\mu}$. The functions $g^{b}_{1}(s)$ and $g^{b}_{n}(s)$
are determined by the boundary conditions for the left end boundary of
$\Gamma_{1}$ and the right end boundary of $\Gamma_{n}$, respectively. The
function $g^{s}_{p}(s)$ represents the MC channel between the reaction systems
$\Sigma_{1}$ and $\Sigma_{n}$ with distance $L_{0}$, which appears when the
periodic boundary $c_{1}(t,0)=y_{n}(t)$ is applied as

$\displaystyle
g^{s}_{p}(s)=\begin{cases}\displaystyle\frac{\hat{s}}{\sinh{\left(L_{0}\hat{s}\right)}}&(\mathrm{Periodic\,boundary})\\\
\displaystyle 0&(\mathrm{Otherwise}).\end{cases}$ (9)

The functions $g^{t}_{i}(s)$ and $g^{s}_{i}(s)$ are obtained by the Laplace
transform of the diffusion equation (1) with Dirichlet boundary conditions.

Figure 2: (A) The MIMO representation of the multi-agent MC system. (B) The
decomposed multi-agent MC system.

Example 1. Consider an MC system, where countless homogeneous nanorobots are
spaced with same distance $L$ as shown in Fig. 3. The MC system can be
approximately modeled by the unit MC system composed of $n$ nanorobots with
the periodic boundary $c_{1}(t,0)=y_{n}(t)$, which leads to
$g^{b}_{1}(s)=g^{t}_{1}(s)$ and $g^{b}_{n}(s)=g^{t}_{n}(s)$. The reaction
system is $H(s)=h(s)I$, where $h(s)=\bm{C}(sI-A)^{-1}\bm{B}$. The MC channel
$G(s)$ of the unit MC system is modeled by a circulant transfer matrix

$G(s)=\mathrm{circ}(g^{t}(s),g^{s}(s),0,\cdots,0,g_{p}^{s}(s)),$ (10)

where $g^{t}(s)$ and $g^{s}(s)$ are defined by (7) and (8) with
$L_{i}=L_{i+1}=L$, respectively.

Figure 3: The multi-agent MC system in Example 1 and the unit MC system.

###  III-B Stability condition for periodic MC systems

We here show a necessary and sufficient stability condition for the closed-
loop system (5) based on the characteristic equation. The characteristic
polynomial $p(s)$ of the closed-loop system (5) is

$p(s):=\det(I-H(s)G(s)).$ (11)

The closed-loop is asymptotically stable if and only if the real part of all
the roots of the characteristic equation $p(s)=0$ are negative.

The roots of the characteristic equation $p(s)$ are not easy to analyze since
it involves the computation of the determinant. In what follows, we show a
theorem that the stability analysis for the MIMO system (5) with the circulant
matrix $G(s)$ can be reduced to that for $n$ SISO systems.

Theorem 1. Consider the MC system with the homogeneous reaction systems $h(s)$
and the MC channel (10) with the periodic boundary $c_{1}(t,0)=y_{n}(t)$. The
closed-loop MC system is asymptotically stable if and only if the real part of
all the roots of the characteristic equations

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

are negative for all $i\,(i=1,2,\cdots,n)$, where

$\lambda_{i}(s)=g^{t}_{i}(s)+g^{s}_{i}(s)\left(\alpha^{i-1}+\alpha^{(i-1)(n-1)}\right)$
(13)

with $\alpha=e^{-j2\pi/n}$.

Proof. 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}$ as

$\Lambda(s):=F^{*}G(s)F,$ (14)

where $\Lambda(s)$ is the diagonal matrix 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],$
(15)

with $\alpha=e^{-j2\pi/n}$. Note that $FF^{*}=I$ and $F^{-1}=F^{*}$, where
$F^{*}$ is the conjugate transpose of $F$, since $F$ is the unitary matrix.
Using the DFT matrix $F$ the characteristic polynomial (11) can be transformed
as

$\displaystyle p(s)$ $\displaystyle=$
$\displaystyle\det\left(F(I-h(s)F^{*}G(s)F)F^{*}\right)$ (16) $\displaystyle=$
$\displaystyle\det(I-h(s)\Lambda(s))$ $\displaystyle=$
$\displaystyle\prod_{i}^{n}[1-h(s)\lambda_{i}(s)]=\prod_{i}^{n}\hat{p}_{i}(s).$

The theorem holds since Eq. (16) 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 $n$-dimensional MIMO system (5) with the circulant matrix
$G(s)$ can be decomposed into $n$ SISO systems as shown in Fig. 2 (B), which
facilitates the analysis for the stability of the MC system using methods such
as Nyquist stability criterion developed in control engineering.

##  IV Numerical example

In this section, we demonstrate the proposed method to analyze the stability
of a specific multi-agent MC system. In particular, we use Theorem 1 to
decompose the MIMO MC system into the multiple SISO systems, and then analyze
the stability of the decomposed SISO systems using Nyquist stability criterion
to show the stable and unstable cases.

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

$\bm{f}(\bm{x})=\left[\begin{matrix}\displaystyle-\delta_{a}x_{a}+\gamma_{a}\frac{x_{a}^{2}}{K_{a}^{2}+x_{a}^{2}}\frac{K_{r}^{2}}{K_{r}^{2}+x_{r}^{2}}\\\
\displaystyle-\delta_{r}x_{r}+\gamma_{r}\frac{x_{a}^{2}}{K_{a}^{2}+x_{a}^{2}}\frac{K_{d}^{2}}{K_{d}^{2}+x_{d}^{2}}\\\
\displaystyle-\delta_{d}x_{d}+\gamma_{d}\frac{K_{r}^{2}}{K_{r}^{2}+x_{r}^{2}}+k(x_{e}(t)-x_{d}(t))\\\
k(x_{d}(t)-x_{e}(t))\end{matrix}\right],$ (17)

where $x_{a}$, $x_{r}$, $x_{d}$, and $x_{e}$ are the concentrations of
activator, repressor, signal molecule in the nanorobot, and signal molecule
outside of the nanorobot, respectively. $\delta_{i}$ and $\gamma_{i}$ are the
degradation rate and the production rate of the corresponding molecular
species, and $K_{i}$ is the Michaelis Menten constant. The parameters values
are referred from Case A of Sec. $\mathrm{I}\mathrm{I}\mathrm{I}$ in [11] with
$k=0.05\,$\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$. Each MC channel $\Gamma_{i}$
has the same distance $L=50\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ and the
diffusion coefficient
$\mu=83\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\cdot\mathrm{s}^{-1}$$. In
what follows, we approximate the MC system by the unit MC system with 4
nanorobots, and analyze the stability of the unit MC system around the
equilibrium point based on the transfer functions using the proposed method.

The transfer function matrix of the MC channel $G(s)$ is obtained as the
$4\times 4$ circulant matrix shown in Eq. (10), and the reaction system is
$H(s)=h(s)I$, where $h(s)=\bm{C}(sI-A)^{-1}\bm{B}$ with $A$ being the Jacobian
matrix of Eq. (17). Using Theorem 1 the characteristic polynomial $p(s)$ of
the MC system can be decomposed into $\hat{p}_{i}(s)$ with $n=4$, where
$\hat{p}_{1}(s)=\hat{p}_{4}(s)$ and $\hat{p}_{2}(s)=\hat{p}_{3}(s)$.

We draw Nyquist plots to find the number of the roots $Z$ of $\hat{p}_{i}=0$
whose real parts are non-negative. The number $Z$ can be calculated by
$Z=N+P$, where $N$ is the number that the trajectory of
$-h(j\omega)\lambda_{i}(j\omega)$ encircles clockwise around the point
$(-1,j0)$, and $P$ is the number of the non-negative roots of
$-h(j\omega)\lambda_{i}(j\omega)$. Fig. 4 (A) depicts the Nyquist plot of the
MC channel for the production rate
$\gamma_{a}=2.5\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$,
which shows each trajectory of $-h(j\omega)\lambda_{i}(j\omega)$ does not
encircle around $(-1,j0)$. Since $-h(j\omega)\lambda_{i}(j\omega)$ has no
unstable poles, i.e. $P=0$, the number of the non-negative roots of the
closed-loop is $Z=0$, and thus the system is stable. On the other hand, Fig. 4
(B) depicts the Nyquist plot of the MC channel for
$\gamma_{a}=3.0\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$,
which shows each trajectory of $-h(j\omega)\lambda_{i}(j\omega)$ does not
encircle around $(-1,j0)$. Since $P=2$, the number of the non-negative roots
of the closed-loop is $Z=2$, which leads that the closed-loop of the MC system
is unstable.

Fig. 5 shows the concentration behavior of each molecular concentration in the
ARD genetic circuit for different production rate $\gamma_{a}$ when
perturbation inputs to each molecular concentration around the equilibrium
point $\bm{x}^{*}$. The molecular concentration converges to the equilibrium
point when
$\gamma_{a}=2.5\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$,
while the molecular concentration oscillates when
$\gamma_{a}=3.0\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$.
Thus, Fig. 5 verifies the proposed stability analysis method in Theorem 1,
which is helpful for the analysis and design of multi-agent dynamical MC
system.

Figure 4: Nyquist plots of the MC system (A) with the production rate
$\gamma_{a}=2.5\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
and (B) with the production rate
$\gamma_{a}=3.0\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$.
Figure 5: The behavior of each molecular concentration in the ARD genetic
circuit around the equilibrium point. (A) The production rate
$\gamma_{a}=2.5\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
and the equilibrium point $\bm{x}^{*}=[7.6,15.6,14.5,14.5]$. (B) The
production rate
$\gamma_{a}=3.0\,$\mathrm{\SIUnitSymbolMicro}\mathrm{M}\cdot\mathrm{m}\mathrm{i}\mathrm{n}^{-1}$$
and the equilibrium point $\bm{x}^{*}=[7.1,17.4,12.4,12.4]$.

##  V Conclusion

In this paper, we have formulated the system theoretic model for multi-agent
MC systems expressed by MIMO systems using the transfer functions based on a
diffusion equation. We have then proposed the method to analyze the stability
for the special case of the circulant structured multi-agent MC systems by
decomposing the $n$-dimensional MIMO MC system into $n$ SISO systems. Finally,
we have demonstrated the proposed method to analyze the stability of a
specific MC system.

## References

  * [1] T. Suda and T. Nakano, “Molecular communication : a personal perspective,” _IEEE Trans. Nanobiosci._ , 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 Commun. Surveys Tuts._ , 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 Commun. Surveys Tuts._ , vol. 22, no. 4, pp. 2765–2800, 2020.
  * [4] M. Femminella, G. Reali, and A. V. Vasilakos, “Molecular communications model for drug delivery,” _IEEE Trans. Nanobiosci._ , vol. 14, no. 7, pp. 935–945, 2015.
  * [5] W. Gao and J. Wang, “Synthetic micro/nanomotors in drug delivery,” _Nanoscale_ , vol. 6, pp. 10 486–10 494, 2014. [Online]. Available: http://dx.doi.org/10.1039/C4NR03124E
  * [6] M. Pierobon and I. Akyildiz, “A physical end-to-end model for molecular communication in nanonetworks,” _IEEE J. Select. Areas Commun._ , vol. 28, no. 4, pp. 602–611, 2010. [Online]. Available: http://ieeexplore.ieee.org/document/5452953/
  * [7] U. A. Chude-Okonkwo, R. Malekian, and B. T. Maharaj, “Diffusion-controlled interface kinetics-inclusive system-theoretic propagation models for molecular communication systems,” _Eurasip J. Adv. Signal Process._ , vol. 2015, no. 1, pp. 1–23, 2015. [Online]. Available: http://dx.doi.org/10.1186/s13634-015-0275-1
  * [8] Y. Huang, F. Ji, Z. Wei, M. Wen, X. Chen, Y. Tang, and W. Guo, “Frequency Domain Analysis and Equalization for Molecular Communication,” _IEEE Trans. Signal Processing_ , vol. 69, pp. 1952–1967, 2021. [Online]. Available: https://ieeexplore.ieee.org/document/9380963/
  * [9] T. Kotsuka and Y. Hori, “Frequency response of diffusion-based molecular communication channels in bounded environment,” in _Proc. Eur. Control Conf. (ECC)_ , 2022, pp. 327–332.
  * [10] ——, “Spatial Frequency-Based Characterization of Disturbance Rejection in Molecular Communication Systems,” _IEEE Trans. Mol. Biol. Multi-Scale Commun._ , vol. 8, no. 1, pp. 36–43, 2022. [Online]. Available: https://ieeexplore.ieee.org/document/9441010/
  * [11] Y. Hori, H. Miyazako, S. Kumagai, and S. Hara, “Coordinated spatial pattern formation in biomolecular communication networks,” _IEEE Trans. Mol. Biol. Multi-Scale Commun._ , vol. 6, no. 2, pp. 111–121, 2015.
  * [12] 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. [Online]. Available: https://dx.plos.org/10.1371/journal.pcbi.1002331
  * [13] S. Hara, T. Kotsuka, and Y. Hori, “Modeling and stability analysis for multi-agent molecular communication systems : a case study for two agents,” in _Proc. SICE Annual Conf._ , 2021, pp. 659–662.
  * [14] 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. [Online]. Available: https://doi.org/10.9746/jcmsi.2.299


11institutetext: Institute for Risk and Uncertainty, University of Liverpool
11email<EMAIL_ADDRESS>United
Kingdom Atomic Energy Authority 33institutetext: Centre for Intelligent
Infrastructure, University of Strathclyde
33email<EMAIL_ADDRESS>
# Correlated Boolean Operators for Uncertainty Logic

Enrique Miralles-Dolz 1122 Ander Gray 1122 Edoardo Patelli 33 Scott Ferson 11

###### Abstract

We present a correlated and gate which may be used to propagate uncertainty
and dependence through Boolean functions, since any Boolean function may be
expressed as a combination of and and not operations. We argue that the and
gate is a bivariate copula family, which has the interpretation of
constructing bivariate Bernoulli random variables following a given Pearson
correlation coefficient and marginal probabilities. We show how this copula
family may be used to propagate uncertainty in the form of probabilities of
events, probability intervals, and probability boxes, with only partial or no
knowledge of the dependency between events, expressed as an interval for the
correlation coefficient. These results generalise previous results by Fréchet
on the conjunction of two events with unknown dependencies. We show an
application propagating uncertainty through a fault tree for a pressure tank.
This paper comes with an open-source Julia library for performing uncertainty
logic.

###### Keywords:

imprecise probability uncertainty logic Boolean functions uncertainty
propagation copula

##  1 Introduction

00footnotetext: This version of the article has been accepted for publication,
after peer review and is subject to Springer Nature’s AM terms of use, but is
not the Version of Record and does not reflect post-acceptance improvements,
or any corrections. The Version of Record is available online at:
https://doi.org/10.1007/978-3-031-08971-8_64

The logical conjunction ($\land$) is a function
$f:\\{0,1\\}^{2}\rightarrow{\\{0,1\\}}$ that returns a value of 1 if and only
if both of the inputs are 1. The logical values assigned to the inputs
represent the truth value of certain propositions or events, and in Boolean
algebra it is required absolute certainty about these truth values (i.e. they
are either true (1) or false (0)). That is, variables take the form
$p\in\\{0,1\\}$. However, this requirement is often too restrictive to be used
in practical applications where truth values have some degree of uncertainty,
and therefore an extension of classical Booleans to uncertain Booleans would
be desirable. We define uncertain Booleans as probabilities $p\in[0,1]$
defining a precise Bernoulli distribution (an event), interval probabilities
$p\subseteq[0,1]$ defining a set of Bernoulli distributions (a credal set),
and probability boxes with $\text{range}(p)\subseteq[0,1]$. This paper
presents a method to perform uncertainty logic with these structures.

Note that it is possible to build any binary Boolean operation from two
primitive operations, for example from and and not. Table 1 shows the 16
($2^{4}$) possible binary and unary Boolean operations, written in terms of
two primitives operations $\land$ and $\neg$. Note that some of these
operations (e.g. Identity and Zero) are trivial, and with the other operations
being written in terms of previously derived operations for brevity. Since all
of these operations may be written in terms of an $\land$ and a $\neg$
operation, and since the not operation is $not(A)=\neg A=1-A$, it follows that
we therefore only need to describe a correlated and operator, and all other
binary operators, and further more complicated Boolean functions, may be
calculated in terms of these two operators.

Table 1: Summary of the 16 possible binary and unary Boolean operations written in terms of two primitive operations $\land$ and $\neg$. $A$ | $1100$ | Adopted Name | Adopted Symbol | Expansion
---|---|---|---|---
$B$ | $1010$ |  |  |
$w_{1}$ | $1100$ | Identity | $A$ | trivial
$w_{2}$ | $1010$ | Identity | $B$ | trivial
$w_{3}$ | $0000$ | Zero | 0 | trivial
$w_{4}$ | $1111$ | One | 1 | trivial
$w_{5}$ | $1000$ | And | $A\land B$ | primitive
$w_{6}$ | $0011$ | Not | $\neg A$ | $1-A$
$w_{7}$ | $0101$ | Not | $\neg B$ | $1-B$
$w_{8}$ | $1110$ | Or |  $A$ $\lor$ $B$ | $\neg((\neg A)\land(\neg B))$
$w_{9}$ | $0111$ | Nand |  $A$ nand $B$ | $\neg(A\land B)$
$w_{10}$ | $0001$ | Nor |  $A$ nor $B$ |  $\neg(A$ $\lor$ $B)$
$w_{11}$ | $0110$ | Exclusive Or |  $A$ xor $B$ |  $(A$ $\lor$ $B)\land(A$ nand $B)$
$w_{12}$ | $1001$ | Equivalence | $A\equiv B$ |  $\neg(A$ xor $B)$
$w_{13}$ | $1011$ | Implication | $A\implies B$ |  $\neg A$ $\lor$ $B$
$w_{14}$ | $1101$ | Implication | $B\implies A$ |  $A$ $\lor$ $\neg B$
$w_{15}$ | $0100$ | Inhibition | $A\Longmapsto B$ | $\neg(A\implies B)$
$w_{16}$ | $0010$ | Inhibition | $B\Longmapsto A$ | $\neg(B\implies A)$

When events $A$ and $B$ are independent, their logical conjunction is
calculated as $\mathbb{P}(A\land B)=\mathbb{P}(A)\mathbb{P}(B)$. However, the
assumption of independence has a significant consequence quantitatively, as
shown in [3]. Therefore, after extending the mathematical structures to
describe events $A$ and $B$ from classical Booleans to uncertain Booleans, the
next desirable extension would be on this assumption of independence between
their probabilities. For example, consider the two following random bit-
vectors each with the same marginal probabilities $P(A)=P(B)=0.5$ (a sequence
of two fair coin tosses)

$\displaystyle A$ $\displaystyle=\\{0,1,1,1,1,0,0,1,1,1,0,0,1\\},$
$\displaystyle B$ $\displaystyle=\\{0,0,0,1,1,1,1,1,1,1,0,0,1\\}.$ (1)

Although the individual coin tosses is fair (the sample mean for these 12
tosses is $\sim 0.6$ for each vector), the vectors are corrected, that is, the
outcome of one throw can influence the other, with $\rho_{AB}=0.3$ (sample
correlation is $\sim 0.35$).

In [9], it is derived a model for the conjunction employing the Pearson
correlation coefficient to capture dependence (referred to as the Lucas
model). Unlike for continuous distributions, two marginals and a correlation
coefficient is sufficient to completely define a bivariate Bernoulli random
variable [7]. The Lucas model is defined as

$\mathbb{P}(A\land
B)=\mathbb{P}(A)\mathbb{P}(B)+\rho_{AB}\sqrt{\mathbb{P}(A)\mathbb{P}(\neg
A)\mathbb{P}(B)\mathbb{P}(\neg B)}\;,$ (2)

where $\rho_{AB}$ is the Pearson correlation coefficient of $A$ and $B$.

However, this model can return misleading results when certain combinations of
probabilities of events and correlations are employed. For example,
considering $\mathbb{P}(A)=0.3$ and $\mathbb{P}(B)=0.2$ with $\rho_{AB}=-1$
(opposite dependence), the Lucas model returns $\mathbb{P}(A\land B)=-0.123$,
which is obviously erroneous. The fact that Equation 2 is returning a negative
probability is simply because the probability assigned for events A and B
cannot have a correlation of -1. This means that, for some probabilities of
events, the Pearson correlation coefficient cannot take any value in [-1, 1],
but in some subset $S\subseteq[-1,1]$.

This subset $S$ can be found through the Fréchet inequalities, which define
the lower and upper bounds of the logical conjunction given the probability of
its events [5], and are written as

$\mathbb{P}(A\land
B)\in[\max(\mathbb{P}(A)+\mathbb{P}(B)-1,0),\min(\mathbb{P}(A),\mathbb{P}(B))]\;.$
(3)

Substituting the bounds in Equation 2, and rearranging for $\rho_{AB}$, the
subset $S=[\underline{\rho}_{AB},\overline{\rho}_{AB}]\subseteq[-1,1]$ can be
found as

$\begin{cases}\underline{\rho}_{AB}=\frac{\max(\mathbb{P}(A)+\mathbb{P}(B)-1,0)-\mathbb{P}(A)\mathbb{P}(B)}{\sqrt{\mathbb{P}(A)\mathbb{P}(\neg
A)\mathbb{P}(B)\mathbb{P}(\neg B)}}\\\ \\\
\overline{\rho}_{AB}=\frac{\min(\mathbb{P}(A),\mathbb{P}(B))-\mathbb{P}(A)\mathbb{P}(B)}{\sqrt{\mathbb{P}(A)\mathbb{P}(\neg
A)\mathbb{P}(B)\mathbb{P}(\neg B)}}\;.\\\ \end{cases}$ (4)

With these definitions for the lower and upper bound of the Pearson
correlation coefficient, the subset $S$ for the previous example with
$\mathbb{P}(A)=0.3$ and $\mathbb{P}(B)=0.2$ is [-0.327, 0.763], and not [-1,1]
as previously guessed.

The proposed correlated and operation combines the Lucas model and the Fréchet
inequalities to restrict the former to return probabilities in [0,1] for any
specified Pearson correlation. If the introduced correlation is greater than
$\overline{\rho}$, then the upper Fréchet bound is returned. On the other
hand, if it is lower than $\underline{\rho}$, then the model gives the lower
Fréchet bound. Lastly, when the introduced correlation is in $S$, then the
probability is calculated following the Lucas model. Therefore, the correlated
and operator is written as

$\mathbb{P}(A\land B)=\\\
\begin{cases}\max(\mathbb{P}(A)+\mathbb{P}(B)-1,0)\qquad\qquad\text{ if
}\rho_{AB}\leq\underline{\rho}_{AB}\\\
\min(\mathbb{P}(A),\mathbb{P}(B))\qquad\qquad\qquad\qquad\text{ if
}\rho_{AB}\geq\overline{\rho}_{AB}\\\
\mathbb{P}(A)\mathbb{P}(B)+\rho(A,B)\sqrt{\mathbb{P}(A)\mathbb{P}(\neg
A)\mathbb{P}(B)\mathbb{P}(\neg B)}\text{ otherwise .}\end{cases}$ (5)

We present an extension for the Boolean operators allowing for uncertainty not
only in the inputs, but also in the dependence, which now can be specified for
any interval correlation $\rho\subseteq[-1,1]$, having the Fréchet bounds as a
special case when $\rho=[-1,1]$. Functions of Bernoulli random variables with
uncertainty in dependence generally yield interval probabilities, we thus show
how intervals may also be propagated through the derived operations.

##  2 Correlated and as a Copula

In this section we argue that the above derived correlated and is a bivariate
copula (2-copula) family, parameterised by a correlation coefficient $\rho$.
Rewritten in a more standard copula notation:

$C_{\rho}(u,v)=\begin{cases}W(u,v)&\text{if }\rho\leq\underline{\rho}_{uv}\\\
M(u,v)&\text{if }\rho\geq\overline{\rho}_{uv}\\\
uv+\rho\sqrt{u(1-u)v(1-v)}&\text{otherwise}\;,\end{cases}$

where $W(u,v)=\max(u+v-1,0)$ and $M(u,v)=\min(u,v)$ are the Fréchet-Hoeffding
copula bounds. A 2-copula $C$ is any function $C:[0,1]^{2}\rightarrow[0,1]$
with the following properties:

  1. 1.

Grounded: $C(0,v)=C(u,0)=0$,

  2. 2.

Uniform margins: $C(u,1)=u;\;C(1,v)=v$,

  3. 3.

2-increasing:
$C(u_{2},v_{2})-C(u_{2},v_{1})-C(u_{1},v_{2})+C(u_{1},v_{1})\geq 0$
for all $0\leq u_{1}\leq u_{2}\leq 1$ and $0\leq v_{1}\leq v_{2}\leq 1$.

It is easy to see that the first two properties hold. The third property is
harder to demonstrate, and in this work we provide no proof, despite the Lucas
model $uv+\rho\sqrt{u(1-u)v(1-v)}$ is non-decreasing in $u$ and $v$. However,
non-decreasing is a necessary but not sufficient criteria for 2-increasing
[12]. Yet Durante and Jaworski [2] prove that $C$ is a copula iff it satisfies
criteria 1. and 2. and if the partial derivatives are increasing (Corollary
2.4). That is, for every $u\in[0,1]$,

$v\mapsto\frac{\delta C(u,v)}{\delta u}$

is increasing on $[0,1]$. The partial derivatives of the above and gate is

$\frac{\delta C_{\rho}(u,v)}{\delta u}=\begin{cases}0&\text{if
}\rho_{uv}\leq\underline{\rho}_{uv}\\\ 1&\text{if
}\rho_{uv}\geq\overline{\rho}_{uv}\\\
v+\rho\frac{v(1-u)(1-v)-uv(1-v)}{2\sqrt{uv(1-u)(1-v)}}&\text{otherwise}\end{cases}$

which from observation $C_{\rho}$ follows. Figure 1 shows $C_{\rho}(u,v)$ for
various values of $\rho$. If $C_{\rho}(u,v)$ is a copula family, then it is a
complete copula family, as it includes the two Fréchet-Hoeffding bounds $W$
and $M$, corresponding to minimal (when $\rho=-1$) and maximal (when $\rho=1$)
correlation respectively, and the independence copula $\Pi(u,v)=uv$ when
$\rho=0$.

Figure 1: $C_{\rho}(u,v)$ for $\rho=\\{-1,-0.8,-0.4,0.4,0.8,1\\}$

For the main results of this paper it is not required that the presented
operator is a copula, only that it is non-decreasing. However it is
interesting since t-norms, which are similar functions to copulas, are widely
used in fuzzy logic to model and operations [8]. The interpretation of t-norms
in fuzzy logic is often unclear, whilst the presented $C_{\rho}$ has a clear
probabilistic interpretation. Given two marginal Bernoulli random variables
with $\mathbb{P}(A)$ and $\mathbb{P}(B)$, and Pearson correlation coefficient
$\rho$, $C_{\rho}(\mathbb{P}(A),\mathbb{P}(B))$ returns the joint probability
of events $A$ and $B$ occurring, i.e. $\mathbb{P}(A=1,B=1)$, which is one
element of the joint probability of the bivariate Bernoulli distribution, with
the other elements being $\mathbb{P}(A=1,B=0)$, $\mathbb{P}(A=0,B=1)$, and
$\mathbb{P}(A=0,B=0)$.

##  3 Interval probabilities

Since $C_{\rho}(u,v)$ is non-decreasing in $u$ and $v$, interval values for
$u$ and $v$ can be simply evaluated with endpoints. Moreover, in
$C_{\rho}(u,v)$ the correlation $\rho$ only has a single occurrence, and so it
can be evaluated exactly with interval arithmetic. However, a useful
observation is that $C_{\rho}(u,v)$ is also non-decreasing in $\rho$, and an
interval value for $\rho$ induces an imprecise copula. The concept of
imprecise copulas, which are a bounded set of copulas, has been discussed in
[10]. From Figure 1 it can be seen that increasing $\rho$ yields larger or
equal values of joint probabilities, i.e. $C_{\rho_{1}}\leq C_{\rho_{2}}$ for
$\rho_{1}\leq\rho_{2}$, and therefore interval uncertainty in $\rho$ may also
be evaluated with endpoints. Therefore, given two events with interval
probabilities $u=[\underline{u},\overline{u}]$ and
$v=[\underline{v},\overline{v}]$, and partially known correlation
$\rho=[\underline{\rho},\overline{\rho}]$, the interval probability of their
conjunction can be evaluated as

$\displaystyle\underline{\mathbb{P}}(u=1,v=1)$
$\displaystyle=C_{\underline{\rho}}(\underline{u},\underline{v}),$
$\displaystyle\overline{\mathbb{P}}(u=1,v=1)$
$\displaystyle=C_{\overline{\rho}}(\overline{u},\overline{v}).$

As an example, for $u=[0.2,0.3]$, $v=[0.45,0.5]$ and $\rho=[-0.2,0.4]$,
rigorous bounds on their conjunction can be calculated as
$\mathbb{P}(u=1,v=1)=[0.0502,0.2417]$.

The full joint probability distribution of the bivariate Bernoulli can also be
found by noticing that, for example,
$\mathbb{P}(u=1,v=0)=\mathbb{P}(u\land\neg v)$, with the not operator defined
as $\neg v=1-v$, which may be evaluated with $C_{\rho}$. However, some careful
consideration is required regarding the correlation coefficient. If
probabilities $u$ and $v$ have correlation $\rho_{uv}$, then $u$ and $1-v$
will have $\rho_{u\neg v}=-1*\rho_{uv}$ [11], i.e. when complementing an
event, the correlation must be reversed. The other elements of the joint
distribution of the bivariate Bernoulli can be calculated as

$\displaystyle\mathbb{P}(u=1,v=0)$ $\displaystyle=C_{-1*\rho}(u,1-v),$
$\displaystyle\mathbb{P}(u=0,v=1)$ $\displaystyle=C_{-1*\rho}(1-u,v),$
$\displaystyle\mathbb{P}(u=0,v=0)$ $\displaystyle=C_{\rho}(1-u,1-v).$

Notice that when the probabilities are complemented twice, the correlation
stays the same, since it has been negated twice: $\rho_{\neg u,\neg
v}=-(-\rho_{uv})=\rho_{uv}$. Table 2 shows the computed joint distribution
using the previous example of $u=[0.2,0.3]$, $v=[0.45,0.5]$ and
$\rho=[-0.2,0.4]$.

Table 2: Interval bounds on the joint probability of a bivariate Bernoulli
with marginals $\mathbb{P}(u=1)=[0.2,0.3]$ , $\mathbb{P}(v=1)=[0.45,0.5]$ and
interval correlation $\rho=[-0.2,0.4]$.

$\mathbb{P}(u,v)$ | 0 | 1 | $\mathbb{P}(u)$
---|---|---|---
0 | $[0.304,0.52]$ | $[0.223,0.44]$ | $[0.7,0.8]$
1 | $[0.02,0.2106]$ | $[0.0502,0.2417]$ | $[0.2,0.3]$
$\mathbb{P}(v)$ | $[0.5,0.55]$ | $[0.45,0.5]$ |

###  3.1 Other Boolean operations

Once a correlated $and$ operation has been constructed, other correlated
Boolean operations can be defined in terms of this operation and a not, shown
in Table 1. For example, or can be defined as

$A\text{ or }B=\mathbb{P}(A\lor B)=1-\mathbb{P}((1-A)\land(1-B)),$

which may be written in terms of $C_{\rho}$ is

$A\lor B=1-C_{\rho}(1-A,1-B).$ (6)

The other operations can be similarly expanded, taking care to negate $\rho$
appropriately.

##  4 Probability Boxes

The previous section describes a method to perform logical operations on
uncertain Booleans characterised by interval probabilities $p\subseteq[0,1]$,
which define a bounded set of Bernoulli distributions. A possible
generalisation of this is to have a distributional or imprecise distributional
(p-box) characterisation, e.g., any p-box $p$ whose range is a subset of the
unit interval $\text{range}(p)\subseteq[0,1]$. The binary events involved with
Boolean operations can be expressed in the form of Bernoullis (precise), set
of Bernoullis (interval), or distributional Bernoullis (p-box). Note that here
we are not describing the arithmetic of real functions ($f:R^{m}\to R^{n}$),
but are describing an extension of Boolean functions with p-box inputs, i.e.,
events or uncertain Booleans with uncertainty characterised by p-boxes. This
is particularly relevant for c-boxes [4], which are p-box shaped confidence
distributions for binomial inference with limited data. That is, they are a
confidence characterisation of an uncertain Boolean given some sample set,
e.g., the data for $A$ in Equation 1.

${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\rho\ =\
-1}$${\color[rgb]{0.15,0.43,0.76}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.43,0.76}\rho\
=\
-0.5}$${\color[rgb]{0.31,0.54,0.07}\definecolor[named]{pgfstrokecolor}{rgb}{0.31,0.54,0.07}\rho\
=\
0}$${\color[rgb]{0.77,0.49,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.77,0.49,0}\rho\
=\
0.5}$${\color[rgb]{0.56,0.07,1}\definecolor[named]{pgfstrokecolor}{rgb}{0.56,0.07,1}\rho\
=\ 1}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\rho\
=\ [0.5,\ 1]}$ Figure 2: (Left) Show logical and between the two K-out-of-N
c-boxes $X\sim\text{KN}(5,6)$, $Y\sim\text{KN}(16,20)$ for different values of
precise correlation. (Right) shows logical and with the interval correlation
$\rho=[0.5,1]$, and is the envelope of the orange and purple c-boxes.
Independence is used in the upper level throughout $C_{XY}=\Pi$.

Since $C_{\rho}$ is a non-decreasing binary operator, it can also be readily
evaluated with the convolutions used in p-box arithmetic. Given two random
variables with distribution functions $F_{X}$ and $F_{Y}$, correlated by
copula $C_{XY}$, a binary operation $Z=f(X,Y)$ can be evaluated with the
following Lebesgue-Stieltjes integral

$F_{Z}(z)=\int_{f\\{z\\}}dC_{XY}(F_{X}(x),F_{Y}(y)),$

where the integration domain is the set
$f\\{z\\}=\\{(x,y)|x,y\in\mathbb{R},f(x,y)<z\\}$ (the set of all $x$ and $y$
for which $f(x,y)<z$). Note that the copula $C_{XY}$ is not the same copula
used to define the correlated and operation $C_{\rho}$, where $C_{XY}$ defines
the dependence between the random variables $X$ and $Y$ (which defines the
uncertainty we have about the events), and $C_{\rho}$ defines the correlation
between the events themselves. We elaborate on this difference in Section 4.1.

Inserting $f=C_{\rho}$ into the above equation, and since the operation is
non-decreasing we have

$\displaystyle\underline{F_{Z}}(z)$
$\displaystyle=\int_{C_{\underline{\rho}}}dC_{XY}\left(\underline{F}_{X}(x),\underline{F}_{Y}(y)\right),$
(7) $\displaystyle\overline{F_{Z}}(z)$
$\displaystyle=\int_{C_{\overline{\rho}}}dC_{XY}\left(\overline{F}_{X}(x),\overline{F}_{Y}(y)\right).$
(8)

for two p-boxes $X=[\underline{F}_{X},\overline{F}_{X}]$,
$Y=[\underline{F}_{Y},\overline{F}_{Y}]$, and interval correlation
$\rho=[\underline{\rho},\overline{\rho}]$. Software for performing rigorous
correlated p-box arithmetic is readily available for bounding the above
integrals efficiently [6]. The above integrals are those usually used in p-box
arithmetic, except that the binary operation $C_{\rho}$ is parameterised by
$\rho$, which could be an interval. In the interval case, the envelope of the
two end points yields the bounds on the output p-box. Figure 2 shows Equations
7 and 8 evaluated for two K-out-of-N c-boxes [4] $X\sim\text{KN}(5,6)$,
$Y\sim\text{KN}(16,20)$ for different values of event correlations, shown in
different colours. The right figure shows logical and with an interval
correlation $\rho=[0.5,1]$, which is the envelope of the purple and orange
c-boxes. Note in Figure 2 independence was used for the upper level $C_{XY}$,
but independence is not necessarily the only choice for $C_{XY}$. In fact, any
copula $C_{XY}$ can be used in Equations (7) and (8), which can greatly
influence results, as much as the choice of correlation $\rho$ does. We
explore this difference in the next section.

###  4.1 Two levels of dependence

The copula $C_{XY}$ in the and operation (7) and (8) plays a different role to
the correlation $\rho$. The correlation $\rho$ is the dependence between the
events $\mathbb{P}(A)=a$ and $\mathbb{P}(B)=b$, i.e., the correlation between
the two random bit-vectors (1), and is the only dependence that plays a role
when two uncertain Booleans are characterised by precise probabilities as real
values. When considering two interval probabilities, two distributions, or two
p-boxes, $\rho$ still correlates the events as before, however one could worry
that the uncertainty characterising the marginal probabilities might share
some bivariate information. That is, the p-boxes have a dependence which the
copula $C_{XY}$ characterises, distinctly from the event correlation $\rho$.

${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}C_{XY}=W}$${\color[rgb]{0.15,0.43,0.76}\definecolor[named]{pgfstrokecolor}{rgb}{0.15,0.43,0.76}C_{X}{}_{Y}=\Phi_{-0.5}}$${\color[rgb]{0.31,0.54,0.07}\definecolor[named]{pgfstrokecolor}{rgb}{0.31,0.54,0.07}C_{XY}=\Pi}$${\color[rgb]{0.77,0.49,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.77,0.49,0}C_{XY}=\Phi_{0.5}}$${\color[rgb]{0.56,0.07,1}\definecolor[named]{pgfstrokecolor}{rgb}{0.56,0.07,1}C_{XY}=M}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}C_{XY}\
=\
}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}[}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}W,\
M}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}]}$ Figure
3: Shows the conjunction of two p-boxes with varying upper-level dependence.
(Left) Show logical and between the two K-out-of-N c-boxes
$X\sim\text{KN}(5,6)$, $Y\sim\text{KN}(16,20)$ for constant independent event
correlation $\rho$, but different values of upper-level dependence. $\Phi_{r}$
is the Gaussian copula with parameter value $r$. (Right) shows logical and
with independence for the lower-level $\rho$, but unknown dependence in the
upper-level.

We call the event correlation $\rho$ as “lower-level dependence” and the
copula $C_{XY}$ as the “upper-level dependence”. This upper-level dependence
can have a profound impact on the conjunction of two uncertain Booleans. The
left of Figure 3 shows the variation in the same two c-boxes as Figure 2, but
with a constant lower-level independence and a varying upper-level $C_{XY}$.
The copula $\Phi_{r}$ is a Gaussian copula with parameter $r$. Note that the
green p-box, showing independence on both levels, is the same in both Figures
2 and 3.

As for the lower-level, the dependence on the upper-level could be imprecisely
known, that is, it may not be possible to know $C_{XY}$. Using p-box
arithmetic, unknown dependence can be propagated through a non-decreasing
binary operation using the following convolutions [3] which has been adapted
using the imprecise probabilistic conjunction $C_{\rho}$

$\displaystyle\underline{F_{Z}}(z)$
$\displaystyle=\inf_{C_{\underline{\rho}}(x,y)=z}\left[W\left(\underline{F}_{X}(x),\underline{F}_{Y}(y)\right)\right],$
(9) $\displaystyle\overline{F_{Z}}(z)$
$\displaystyle=\sup_{C_{\overline{\rho}}(x,y)=z}\left[W^{d}\left(\overline{F}_{X}(x),\overline{F}_{Y}(y)\right)\right],$
(10)

where $C^{d}$ is the dual copula of $C$: $C^{d}(u,v)=u+v-C(u,v)$. The right of
Figure 3 shows the result of the p-box conjunction with unknown upper-level
dependence and with lower-level independence.

The p-box bounds from unknown upper-level dependence are not only rigorous but
are also best-possible in the sense that they cannot be made tighter without
introducing additional dependence assumptions. One may find the breath (grey
shaded area) of this p-box surprising, and how strong the independence
assumption is. One may also wonder why the p-box on the right is not a simple
envelope of the p-boxes on the left. This is because the left p-boxes only
consider Gaussian dependencies, whilst the right p-box considers all possible
copulas $C_{XY}$, of which there are an infinite number.

Risk analysts tell us that for a fault tree analysis to be probabilistic the
event probabilities should be characterised by distributions [1]. If so, then
the issue of the two levels of dependence immediately arises. As far as we are
aware, this issue has not been widely addressed. Generally, the upper-level
independence is assumed. This section has shown that this assumption has a
strong effect on the results of logical operations involving uncertainty.

##  5 Application

In fault tree analysis, Boolean operators are employed to calculate a system’s
probability of failure. These are useful to, for example, understand what are
the most likely routes of failure, show compliance with the reliability
requirements, or designing monitoring strategies. To carry out these analyses
it is required to define an event tree, where the connections between events
are made with Boolean operators, with their respective event probabilities and
dependencies. Then, the Boolean operators are calculated backwards from the
top event to find its probability of failure.

Figure 4 represents a fault tree for a pressure tank system, derived in [13,
3]. It displays the Boolean operators $\land$ (AND) and $\lor$ (OR) connecting
the failure events $E_{i}$ (with $E_{1}$ being the top event), and the system
components tank (T), relay (K2), pressure switch (S), on-switch (S1), timer
relay (R), and relay (K1). For the sake of brevity, the function of these
components will not be explained in this paper.

Figure 4: Fault tree diagram for the pressure tank system.

To measure the probability of failure of the system (event $E_{1}$), it is
necessary to estimate the failure rates for its components T, K2, S, S1, R and
K1, and the dependence of the events in the Boolean operations. The most
common method of calculation, as in [13], is to calculate the probability of
failure assuming independence and known rates of failure for the components.
In [3], the assumption of lower-level independence was relaxed including some
known/unknown dependence, and extending the point probabilities to intervals,
showing that the assumptions in [13] underestimate the probability of failure.
To demonstrate the capabilities of the method derived in this paper, two
scenarios on the components probability of failure will be calculated, with
three lower-level dependence cases for each (independence, mixed dependence,
and unknown dependence):

  1. 1.

An interval scenario, where the probability of failure of the components is
modeled as intervals.

  2. 2.

A p-box scenario, which generalises [3] combining p-boxes and intervals. The
probabilities of failure of the relays (K1, K2, and R) have been extended to
follow a K-out-of-N c-box [4].

The specific values for each scenario and dependence are indicated in Tables 3
and 4 respectively.

Table 3: Probability of failure for the tank pressure system components. In the interval scenario, all probabilities are modeled in the form of intervals. In the p-box scenario, tank and switches have an interval probability, whilst the relays follow a K-out-of-N c-box. Component | Interval scenario | P-box scenario
---|---|---
T | ${[}4.5\times 10^{-6},5.5\times 10^{-6}{]}$ | ${[}4.5\times 10^{-6},5.5\times 10^{-6}{]}$
K2 | ${[}2.5\times 10^{-5},3.5\times 10^{-5}{]}$ | $KN(3,10^{5})$
S | ${[}0.5\times 10^{-4},1.5\times 10^{-4}{]}$ | ${[}0.5\times 10^{-4},1.5\times 10^{-4}{]}$
K1 | ${[}2.5\times 10^{-5},3.5\times 10^{-5}{]}$ | $KN(3,10^{5})$
R | ${[}0.5\times 10^{-4},1.5\times 10^{-4}{]}$ | $KN(1,10^{4})$
S1 | ${[}2.5\times 10^{-5},3.5\times 10^{-5}{]}$ | ${[}2.5\times 10^{-5},3.5\times 10^{-5}{]}$
Table 4: Dependence in the system’s events for the three different cases: independence, mixed dependence, and unknown dependence. Event | Independence | Mixed Dependence | Unknown Dependence
---|---|---|---
$E_{1}$ | $\rho=0$ | $\rho=0$ | $\rho=[-1,1]$
$E_{2}$ | $\rho=0$ | $\rho=[-1,1]$ | $\rho=[-1,1]$
$E_{3}$ | $\rho=0$ | $\rho=0.15$ | $\rho=[-1,1]$
$E_{4}$ | $\rho=0$ | $\rho=[-0.2,0.2]$ | $\rho=[-1,1]$
$E_{5}$ | $\rho=0$ | $\rho=1$ | $\rho=[-1,1]$

###  5.1 Interval scenario

On the left of Figure 5 it is shown the probability of failure $E_{1}$ for the
interval scenario, where the red area belongs to the independence case, blue
to the mixed dependence case, and green to the unknown dependence. It is
possible to see how the uncertainty increases as the dependence assumptions
are removed from the Boolean operations. Intervals for the probability of
$E_{1}$ are indicated in Table 5.

Table 5: Probability of $E_{1}$ in the interval scenario. Case | $\mathbb{P}(E_{1})$
---|---
Independence | ${[}2.950\times 10^{-5},4.053\times 10^{-5}{]}$
Mixed Dependence | ${[}2.949\times 10^{-5},6.551\times 10^{-5}{]}$
Unknown Dependence | ${[}2.499\times 10^{-5},1.905\times 10^{-4}{]}$

###  5.2 P-box scenario

The results in the P-box scenario are more complicated to interpret since the
calculations are probabilities of the probability of the event $E_{1}$.
However, meaningful analyses can be drawn from them. For example, assuming the
required failure probability of the tank system is no more than $10^{-4}$, a
possible question could be what is the probability of it being lower or equal
than $10^{-4}$, or how likely is the system to comply with that requirement.
Table 6 includes the results of such inquiry. In the case of all the events
being independent, the probability of fulfilling the requirements goes from
0.969 to 1, so one can infer the tank pressure system is likely to fulfil the
requirements. When the independence assumption is removed, and the mixed
dependence scenario is adopted, it is possible to see how the uncertainty on
the assessment increases. Finally, if all the dependence assumptions are
relaxed, the probability goes from $0$ to $1$, meaning that no guarantees can
be given on the reliability of the tank pressure system. These results are
illustrated in Figure 5 (right), which depicts how dramatic the consequences
can be when assuming certain dependencies for the events.

Table 6: Probability of $\mathbb{P}(E_{1})<=10^{-4}$ in the P-box scenario. Case | $\mathbb{P}(E_{1})<=10^{-4}$
---|---
Independence | [0.969, 1]
Mixed Dependence | [0.88, 1]
Unknown Dependence | [0, 1]

Figure 5: (Left) Interval of $\mathbb{P}(E_{1})$ for different dependence
assumptions. (Right) Probability box of $\mathbb{P}(E_{1})$ for different
dependence assumptions.

##  6 Conclusion

In this paper, classical Boolean functions have been generalised to be able to
operate with precise probabilities, intervals, p-boxes, and with any input
correlation. We show an application of the generalisation of Boolean
operations calculating the probability of failure of a pressure tank system.
The uncertainty is propagated through the fault tree with different
combinations of probabilities of the events, intervals, and p-boxes, under
different dependence assumptions. The results suggest that assumptions on the
probability or dependence of the events have a strong impact on the outcome of
the analysis, and these should be carefully addressed in any serious
assessment.

Also, the issue of the two levels of dependence has been introduced, and shown
that it can have a dramatic effect on the results of a probabilistic risk
assessment.

The computational resources to perform correlated Boolean operations are
available in the following open-source Julia package:
https://github.com/Institute-for-Risk-and-Uncertainty/UncLogic.jl

## Acknowledgements

We thank William (Bill) Huber (from Analysis & Inference) for his advice in
the early stages of this research. This research was funded by the EPSRC and
ESRC CDT in Risk and Uncertainty (EP/L015927/1), established within the
Institute for Risk and Uncertainty at the University of Liverpool. This work
has been carried out within the framework of the EUROfusion Consortium, funded
by the European Union via the Euratom Research and Training Programme (Grant
Agreement No 101052200 — EUROfusion). Views and opinions expressed are however
those of the author(s) only and do not necessarily reflect those of the
European Union or the European Commission. Neither the European Union nor the
European Commission can be held responsible for them.

## References

  * [1] Atwood, C., LaChance, J., Martz, H., Anderson, D., Englehardt, M., Whitehead, D., Wheeler, T.: Handbook of parameter estimation for probabilistic risk assessment. US Nuclear Regulatory Commission (2003)
  * [2] Durante, F., Jaworski, P.: A new characterization of bivariate copulas. Communications in Statistics—Theory and Methods 39(16), 2901–2912 (2010)
  * [3] Ferson, S., Nelsen, R.B., Hajagos, J., Berleant, D.J., Zhang, J., Tucker, W.T., Ginzburg, L.R., Oberkampf, W.L.: Dependence in probabilistic modeling Dempster—Shafer theory and probability bounds analysis. Tech. rep., Sandia National Lab.(SNL-NM), Albuquerque, NM (United States) (2015)
  * [4] Ferson, S., O’Rawe, J., Balch, M.: Computing with confidence: imprecise posteriors and predictive distributions. In: Vulnerability, Uncertainty, and Risk: Quantification, Mitigation, and Management, pp. 895–904 (2014)
  * [5] Fréchet, M.: Généralisation du théoreme des probabilités totales. Fundamenta mathematicae 1(25), 379–387 (1935)
  * [6] Gray, A., Ferson, S., Patelli, E.: Probabilityboundsanalysis.jl: Arithmetic with sets of distributions. Submitted to the Proceedings of JuliaCon (2021)
  * [7] Joe, H.: Multivariate models and multivariate dependence concepts. CRC press (1997)
  * [8] Klir, G.J., Folger, T.A.: Fuzzy sets, uncertainty, and information. Prentice-Hall, Inc. (1987)
  * [9] Lucas, D.J.: Default correlation and credit analysis. The Journal of Fixed Income 4(4), 76–87 (1995)
  * [10] Montes, I., Miranda, E., Pelessoni, R., Vicig, P.: Sklar’s theorem in an imprecise setting. Fuzzy Sets and Systems 278, 48–66 (2015)
  * [11] Nelsen, R.B.: An introduction to copulas. Springer Science & Business Media (2007)
  * [12] Schweizer, B., Sklar, A.: Probabilistic metric spaces. Courier Corporation (2011)
  * [13] Vesely, W.E., Goldberg, F.F., Roberts, N.H., Haasl, D.F.: Fault tree handbook. Tech. rep., Nuclear Regulatory Commission Washington DC (1981)