Two-Way Molecular Communications

TL;DR

This paper investigates two-way molecular communication systems for nano devices, proposing models and interference cancellation techniques to improve performance when using a single molecule type.

Contribution

It introduces a novel two-way molecular communication system using a single molecule type and develops models and cancellation techniques to enhance performance.

Findings

Proposed a channel model for single-molecule two-way communication.

Analyzed bit error rate and throughput of the system.

Demonstrated performance improvements with interference cancellation.

Abstract

For nano-scale communications, there must be cooperation and simultaneous communication between nano devices. To this end, in this paper we investigate two-way (a.k.a. bi-directional) molecular communications between nano devices. If different types of molecules are used for the communication links, the two-way system eliminates the need to consider self-interference. However, in many systems, it is not feasible to use a different type of molecule for each communication link. Thus, we propose a two-way molecular communication system that uses a single type of molecule. We develop a channel model for this system and use it to analyze the proposed system's bit error rate, throughput, and self-interference. Moreover, we propose analog- and digital- self-interference cancellation techniques. The enhancement of link-level performance using these techniques is confirmed with both numerical and analytical results.

Figures

Tables

Table 1. Table I: Simulation Parameters for Channel Model Verification

System Parameter	Notation	Values
Distance (${\text{Tx}}_i$ to ${\text{Rx}}_i$)	$d_1$	$1.5\mu m$
Distance (${\text{Tx}}_i$ to ${\text{Rx}}_j$)	$d_{{\text{Tx}}_1}^{{\text{Rx}}_2}$	$5\mu m$-$d_1$
Number of Molecules for Bit 1	$N_1$	50000
Diffusion Coefficient	D	$100\mu m^2/s$
Radius of Receiver	$r_{r1}$=$r_{r2}$	$5\mu m$
Simulation Time Step	$\Delta t$	${10}^{-5}s$
Simulation Duration		$0.1s$
Simulation Replication		10
Number of the terms in (13) and (33)	m, n	100000

Table 2. Table II: Simulation Parameters for BER and Throughput Comparison

System Parameter	Notation	Values
Distance (${\text{Tx}}_i$ to ${\text{Rx}}_i$)	$d_1$=$d_2$	$1.5\mu m$
Distance (${\text{Tx}}_i$ to ${\text{Rx}}_j$)	$d_{{\text{Tx}}_2}^{{\text{Rx}}_1}$=$d_{{\text{Tx}}_1}^{{\text{Rx}}_2}$	$5\mu m$-$d_1$
Number of Molecules for Bit 1	$N_1$	300, 400, 500
Diffusion Coefficient	D	$100\mu m^2/s$
Radius of Receiver	$r_r$	$5\mu m$
Simulation Time Step	$\Delta t$	${10}^{-5}s$
Molecular Noise Variance	${\sigma }_{\text{noise}}^2$	$100$
Considered ISI Period		0.6s
Number of the terms in (13) and (33)	m, n	100000
Replication		10

Table 3. Table III: Throughput ratio (FD/HD), t s subscript 𝑡 𝑠 t_{s} and BER values of the half-duplex system using BCSK and the full-duplex system using BCSK.

		$N_1=300$				$N_1=400$				$N_1=500$
	$t_s^{\text{HD}}=2t_s^{\text{FD}}$	$0.100s$	$0.200s$	$0.300s$	$0.400s$	$0.100s$	$0.200s$	$0.300s$	$0.400s$	$0.100s$	$0.200s$	$0.300s$	$0.400s$
Case 1	BER (FD)	0.1155	0.0229	0.0083	0.0045	0.0887	0.0097	0.0022	0.0009	0.0721	0.0045	$6.7\times {10}^{-4}$	$1.9\times {10}^{-4}$
	BER (HD)	0.0931	0.0086	0.0022	0.0010	0.0660	0.0017	0.0002	$5\times {10}^{-5}$	0.0492	$3.1\times {10}^{-4}$	$1.5\times {10}^{-5}$	$1.5\times {10}^{-6}$
	THP Ratio	1.9506	1.9711	1.9876	1.9929	1.9515	1.9840	1.9959	1.9983	1.9517	1.9916	1.9987	1.9996
	$t_s^{\text{FD}}$=$t_s^{\text{HD}}$	$0.100s$	$0.200s$	$0.300s$	$0.400s$	$0.100s$	$0.200s$	$0.300s$	$0.400s$	$0.100s$	$0.200s$	$0.300s$	$0.400s$
Case 2	BER (FD)	0.0229	0.0045	0.0023	0.0014	0.0097	0.0009	0.0003	0.0001	0.0045	$1.9\times {10}^{-4}$	$4.7\times {10}^{-5}$	$1.1\times {10}^{-5}$
	BER (HD)	0.0931	0.0086	0.0022	0.0010	0.0660	0.0017	0.0002	$5\times {10}^{-5}$	0.0492	$3.1\times {10}^{-4}$	$1.5\times {10}^{-5}$	$1.5\times {10}^{-6}$
	THP Ratio	1.0774	1.0041	0.9999	0.9997	1.0603	1.0008	0.9999	0.9999	1.0469	1.0001	1.0000	1.0000
	$t_s^{\text{HD}}$	$0.100s$	$0.200s$	$0.300s$	$0.400s$	$0.100s$	$0.200s$	$0.300s$	$0.400s$	$0.100s$	$0.200s$	$0.300s$	$0.400s$
Case 3	$t_s^{\text{FD}}$	$0.056s$	$0.147s$	$0.320s$	$0.630s$	$0.056s$	$0.165s$	$0.390s$	$2.785s$	$0.056s$	$0.178s$	$0.398s$	N/A
	THP Ratio	1.7747	1.3603	0.9375	0.6349	1.7853	1.2122	0.7692	0.1436	1.7820	1.1236	0.7538	N/A

Table 4. Table IV: Throughput ratio (FD/HD) of the half-duplex system using QCSK and the full-duplex system using BCSK with optimized D-SIC and A-SIC

123 $N_1$ $1^{2^3}$	$t_s^{\text{HD}}\left(s\right)$	$t_s^{\text{FD}}\left(s\right)$	Throughput Ratio (FD/HD)
300	0.200	0.100	1.222
	0.300	0.150	1.295
	0.400	0.200	1.274
400	0.200	0.100	1.134
	0.300	0.150	1.188
	0.400	0.200	1.162
500	0.200	0.100	1.078
	0.300	0.150	1.125
	0.400	0.200	1.098

Equations

$$\frac{\partial p(r,t|r_0)}{\partial t}=D{\nabla }^{2}p(r,t|r_0),$$

$$\begin{array}{rl}{G}_1^{\text{Tx}}(t)& ={\int }_0^t{g}_1^{\text{Tx}}(t^{\prime })d{t}^{\prime }\\ & ={\int }_0^t\frac{r_r}{r_r+d}\frac{d}{\sqrt{4\pi D{t}^{\prime 3}}}{e}^{-d^2/4D{t}^{\prime }}d{t}^{\prime }\\ & =\frac{r_r}{r_r+d}\text{erfc}\left(\frac{d}{\sqrt{4Dt}}\right),\end{array}$$

$$P_{ij}[k]\triangleq P_{ij}(k{t}_s,(k+1)t_s),k\in {\mathbb{N}}_0,$$

$$N_{ij}[k]\sim \mathcal{B}(N_1,P_{ij}[k]),$$

$$y_{{\text{Rx}}_j}[n]\triangleq \sum _{k=0}^{n-1}(N_{ij}[k]\cdot x_i[n-k]+\underset{\text{self-interference}}{\underbrace{N_{jj}[k]\cdot x_j[n-k]}})+s[n].$$

$$N_{ij}[k]\approx \mathcal{N}(N_1{P}_{ij}[k],N_1{P}_{ij}[k](1-P_{ij}[k])),$$

$$\begin{array}{rl}{y}_{{\text{Rx}}_j}[n]& \sim \mathcal{N}({\mu }_{n,\text{total}},{\sigma }_{n,\text{total}}^2)\\ {\mu }_{n,\text{total}}& =\sum _{k=0}^{n-1}{N}_1(P_{ij}[k]x_i[n-k]+P_{jj}[k]x_j[n-k])\\ {\sigma }_{n,\text{total}}^2& ={\sigma }_{\text{noise}}^2+\sum _{k=0}^{n-1}{N}_1(P_{ij}[k](1-P_{ij}[k])x_i[n-k]\\ & +P_{jj}[k](1-P_{jj}[k])x_j[n-k]).\end{array}$$

$$P_{ij}^{\phi }[k]=P_{ij}(k{t}_s+T_c,(k+1)t_s).$$

$$N_{ij}^{\phi }[k]\sim \mathcal{B}(N_1,P_{ij}^{\phi }[k]).$$

$$\begin{array}{rl}{y}_{{\text{Rx}}_j}^{\phi }[n]& \triangleq \sum _{k=0}^{n-1}(N_{ij}^{\phi }[k]\cdot x_i[n-k]+\underset{\text{self-interference}}{\underbrace{N_{jj}^{\phi }[k]\cdot x_j[n-k])}}\\ & -\underset{\text{when D-SIC applied}}{\underbrace{\mathbb{E}\left[N_{jj}^{\phi }[0]\cdot x_j[n]\right]}},\\ \text{where}& \mathbb{E}\left[N_{jj}^{\phi }[0]\cdot x_j[n]\right]=N_1{P}_{jj}^{\phi }[0]x_j[n].\end{array}$$

$$\begin{array}{rl}{y}_{{\text{Rx}}_j}^{\phi }[n]& \sim \mathcal{N}({\mu }_{n,\text{total}}^{\phi },({\sigma }_{n,\text{total}}^{\phi }{)}^2)\\ {\mu }_{n,\text{total}}^{\phi }& =\sum _{k=0}^{n-1}{N}_1(P_{ij}^{\phi }[k]x_i[n-k]+P_{jj}^{\phi }[k]x_j[n-k])\\ & -\underset{\text{when D-SIC applied}}{\underbrace{N_1{P}_{jj}^{\phi }[0]x_j[n]}}\\ ({\sigma }_{n,\text{total}}^{\phi }{)}^2& ={\sigma }_{\text{noise}}^2+\sum _{k=0}^{n-1}{N}_1(P_{ij}^{\phi }[k](1-P_{ij}^{\phi }[k])x_i[n-k]\\ & +P_{jj}^{\phi }[k](1-P_{jj}^{\phi }[k])x_j[n-k]).\end{array}$$

$$\begin{array}{rl}\mu & ={\tanh }^{-1}\left(\frac{fz}{x^2+y^2+z^2+(f/2{)}^2}\right),\\ \eta & ={\tan }^{-1}\left(\frac{f(x^2+y^2{)}^{1/2}}{x^2+y^2+z^2-(f/2{)}^2}\right),\\ \varphi & ={\tan }^{-1}\left(\frac{y}{x}\right).\end{array}$$

$$\begin{array}{rl}\underset{t\to \infty }{\lim }& F_1^{{\text{Tx}}_1}(t)=\sqrt{2(\cosh {\mu }_0-\cos {\eta }_0)}\\ & \times \left(\sum _{m=0}^{\infty }\frac{e^{-(m+\frac{1}{2}){\mu }_1}\sinh [(m+\frac{1}{2})({\mu }_0+{\mu }_2)]}{\sinh [(m+\frac{1}{2})({\mu }_1+{\mu }_2)]}{P}_m(\cos {\eta }_0)\right),\\ \underset{t\to \infty }{\lim }& F_2^{{\text{Tx}}_1}(t)=\sqrt{2(\cosh {\mu }_0-\cos {\eta }_0)}\\ & \times \left(\sum _{m=0}^{\infty }\frac{e^{-(m+\frac{1}{2}){\mu }_2}\sinh [(m+\frac{1}{2})({\mu }_1-{\mu }_0)]}{\sinh [(m+\frac{1}{2})({\mu }_1+{\mu }_2)]}{P}_m(\cos {\eta }_0)\right),\end{array}$$

$$f=r_{r2}\sinh \left({\cosh }^{-1}\left(\frac{a}{r_{r2}}\right)\right).$$

$${\nabla }^2{k}_1(\vec{r})=0,$$

$${\nabla }^2{k}_2(\vec{r})=0.$$

$$k_1(\vec{r})=1,k_2(\vec{r})=0\text{when}{d}_{\vec{r}}^{{\text{Rx}}_1}=0,$$

$$k_1(\vec{r})=0,k_2(\vec{r})=1\text{when}{d}_{\vec{r}}^{{\text{Rx}}_2}=0,$$

$$k_1(\vec{r}),k_2(\vec{r})\to 0\text{as}{d}_{\vec{r}}^{{\text{Rx}}_1},d_{\vec{r}}^{{\text{Rx}}_2}\to \infty ,$$

$$a=\frac{{\ell }^2+r_{r2}^2-r_{r1}^2}{2\ell },$$

$$\begin{array}{rl}\underset{t\to \infty }{\lim }{F}_1^{{\text{Tx}}_1}(t)=& \sqrt{\frac{\cosh {\mu }_0-\cos {\eta }_0}{\cosh ({\mu }_0-2{\mu }_1)-\cos {\eta }_0}}\\ =& \frac{r_{r1}}{r_{r1}+d_{\vec{r}}^{{\text{Rx}}_1}},\text{when}{r}_{r2}\to 0,\end{array}$$

$$F_2^{{\text{Tx}}_1}(t)={\int }_0^t{f}_2^{{\text{Tx}}_1}(t^{\prime })d{t}^{\prime },$$

$$f_2^{{\text{Tx}}_1}(t^{\prime })=g_2^{{\text{Tx}}_1}(t^{\prime })-\alpha (t^{\prime }),$$

$$f_2^{{\text{Tx}}_1}(t^{\prime })=g_2^{{\text{Tx}}_1}(t^{\prime })-{\int }_0^{t^{\prime }}{\int }_{{\Omega }_1}{f}_1^{{\text{Tx}}_1}(\tau ,s_{1,{\text{Tx}}_1})g_2^{s_{1,{\text{Tx}}_1}}(t^{\prime }-\tau )d{s}_{1,{\text{Tx}}_1}d\tau ,$$

$$f_1^{{\text{Tx}}_1}(t^{\prime })=g_1^{{\text{Tx}}_1}(t^{\prime })-{\int }_0^{t^{\prime }}{\int }_{{\Omega }_2}{f}_2^{{\text{Tx}}_1}(\tau ,s_{2,{\text{Tx}}_1})g_1^{s_{2,{\text{Tx}}_1}}(t^{\prime }-\tau )d{s}_{2,{\text{Tx}}_1}d\tau ,$$

$$f_2^{{\text{Tx}}_1}(t^{\prime })=g_2^{{\text{Tx}}_1}(t^{\prime })-{\int }_0^{t^{\prime }}{g}_2^{s_1^{\prime }(t^{\prime },\tau )}(t^{\prime }-\tau ){\int }_{{\Omega }_1}{f}_1^{{\text{Tx}}_1}(\tau ,s_{1,{\text{Tx}}_1})d{s}_{1,{\text{Tx}}_1}d\tau ,$$

$$f_2^{{\text{Tx}}_1}(t^{\prime })=g_2^{{\text{Tx}}_1}(t^{\prime })-{\int }_0^{t^{\prime }}{g}_2^{s_1^{\prime }(t^{\prime },\tau )}(t^{\prime }-\tau )f_1^{{\text{Tx}}_1}(\tau )d\tau .$$

$$f_1^{{\text{Tx}}_1}(t^{\prime })=g_1^{{\text{Tx}}_1}(t^{\prime })-{\int }_0^{t^{\prime }}{g}_1^{s_2^{\prime }(t^{\prime },\tau )}(t^{\prime }-\tau )f_2^{{\text{Tx}}_1}(\tau )d\tau ,$$

$$\begin{array}{r}{F}_2^{{\text{Tx}}_1}(t)=G_2^{{\text{Tx}}_1}(t)-{\int }_0^t{g}_2^{s_1^{\ell }}(t^{\prime })\ast f_1^{{\text{Tx}}_1}(t^{\prime })d{t}^{\prime }\\ F_1^{{\text{Tx}}_1}(t)=G_1^{{\text{Tx}}_1}(t)-{\int }_0^t{g}_1^{s_2^{\ell }}(t^{\prime })\ast f_2^{{\text{Tx}}_1}(t^{\prime })d{t}^{\prime }.\end{array}$$

$$\begin{array}{r}{F}_2^{{\text{Tx}}_1}(t)=G_2^{{\text{Tx}}_1}(t)-g_2^{s_1^{\ell }}(t)\ast F_1^{{\text{Tx}}_1}(t)\\ F_1^{{\text{Tx}}_1}(t)=G_1^{{\text{Tx}}_1}(t)-g_1^{s_2^{\ell }}(t)\ast F_2^{{\text{Tx}}_1}(t).\end{array}$$

$$\begin{array}{rl}\underset{t\to \infty }{\lim }{F}_2^{{\text{Tx}}_1}(t)& =\underset{t\to \infty }{\lim }{G}_2^{{\text{Tx}}_1}(t)-\underset{t\to \infty }{\lim }{G}_2^{s_1^{\ell }}(t)\underset{t\to \infty }{\lim }{F}_1^{{\text{Tx}}_1}(t)\\ \underset{t\to \infty }{\lim }{F}_1^{{\text{Tx}}_1}(t)& =\underset{t\to \infty }{\lim }{G}_1^{{\text{Tx}}_1}(t)-\underset{t\to \infty }{\lim }{G}_1^{s_2^{\ell }}(t)\underset{t\to \infty }{\lim }{F}_2^{{\text{Tx}}_1}(t),\end{array}$$

Full text

I Introduction

Over the past decade, developments in the field of nano robotics have enabled the use of nano devices in various technologies and especially in those used by the bio-medical industry [2, 3, 4, 5, 6]. Well-organized clusters of nano devices can be used for drug delivery applications and artificial immune systems. Each cluster is responsible for a single task, e.g., discovering or destroying of pathogens. Since a nano device can only perform simple tasks, it is important to have a communication system among nano devices. Radio frequency (RF)-based communication is not suitable for nano devices because of physical limitations such as the size of the antenna, which is typically proportional to the wavelength of the electromagnetic (EM) wave in order to maximize efficiency [7, 8]. Furthermore, EM waves—especially at high frequencies—do not propagate well in the body [9].

Thus, researchers have focused on molecular communication as an alternative to RF-based communication, where information is transmitted via molecules. One such system is that of molecular communication via diffusion (MCvD). Here molecules are propagated in an environment by diffusion [10]. An MCvD system mainly consists of the following: transmitter nodes capable of emitting and modulating information through molecules, receiver nodes capable of receiving and demodulating molecular signals, information molecules to transfer information, and a fluid environment to host nodes and molecules. Physical realization of the transceiver is also an important issue. The authors in [11, 12] used Quorum Sensing [13] to synchronize the activities of a cluster of engineered bacteria implemented on a chip. Using this orchestrated bacteria cluster as a node, the authors in [14] proposed a diffusion-based molecular communication network.

One of the main challenges in MCvD is to establish channel models for representing the molecular received signal (i.e., a time-dependent solution of the fraction of received molecules by time $t$). Some known channel models assume that the arrival time of molecules are a first-passage time process (i.e., information molecules are absorbed whenever they hit a receiver) [15, 16, 17, 18, 19]. Compared to the case of transparent receivers, using the absorbing receiver makes it hard to derive the channel model due to the additional boundary condition. In this paper, we consider the absorbing receivers. The authors in [18] modeled the molecular received signal in a three-dimensional (3-D) environment– a point source represented a transmitter, and an absorbing sphere represented a receiver. For this basic topology, it is possible to acquire an analytical closed form of the channel model representing the molecular received signal due to spherical symmetry. If the system has more than one absorbing sphere (receiver), however, such symmetry disappears. It is then difficult to model the arrival times mathematically.

Other challenges in MCvD include low transmission rates due to severe inter-symbol interference (ISI). In MCvD, ISI occurs when the molecules of a previous symbol are absorbed by the target receiver in the current or future symbol slots. The heavy tail nature of impulse responses in MCvD causes severe ISI. Thus, several researchers [20, 21] have suggested ISI mitigation techniques, including enzymatic degradation of ISI using different molecule types. All of these works assumed a one-way MCvD system whereby molecules are transmitted in one direction from the transmitter to the receiver. While this simplifies the design, recent work in full-duplex radio communications [22, 23, 24] indicates that data rate gains along with other performance advantages may be obtained from two-way communication.

In this paper, we propose a two-way MCvD system that uses a single type of molecule for simultaneous communication between two nano devices. If each of the directional communication links uses a different type of molecule, then self-interference111Self-interference in MCvD refers to a phenomenon in which a molecule emitted from the transmitter is absorbed by its own rather than the intended receiver. (SI) will not occur [25, 26, 27]. However, this is not a feasible solution. First, the nano devices are too simple to perform complex tasks and, second, the number of molecule types will increase rapidly, on the order of $\binom{\ell}{2}$ where $\ell$ is the number of communication links. Therefore, we propose a two-way MCvD system that uses a single type of molecule for each link. Two communication models—a half-duplex system and a full-duplex system—are considered. In the half-duplex system, each paired transceiver (i.e., for $\text{Rx}_{1}$, the paired transmitter is $\text{Tx}_{2}$) operates alternately with respect to time. In the full-duplex system, each paired transceiver operates simultaneously with respect to time (i.e., it is not necessary to wait until the other transceiver ends its communication). However, the analytical modeling of such a system is not without difficulty—two absorbing spheres are to use a single type of molecule. Unfortunately, to implement simultaneous communication between two nano devices, we cannot use previous studies on MCvD channel models as such studies assumed a one-way MCvD system with a single absorbing sphere. Therefore, in this paper, we analytically model the molecular received signal in the case of two absorbing spheres, propose SI cancellation (SIC) techniques, and analyze the proposed system’s performance in terms of bit error rate (BER) and throughput. The main contributions of this paper are as follows:

• •

The channel model of a two-way MCvD system with two absorbing sphere receivers is derived. To the best of our knowledge, it is the first work that derives the channel model for multi absorbing receiver case. We improve the channel model derived in [1]. More precisely, we analytically derive the asymptotic capture probability (i.e., the probability that the molecules eventually being absorbed by the receivers). By using the derived asymptotic capture probability, we also derive an approximation of the time-dependent solution of the fraction of the received molecules by time $t$. The derived channel model is validated through particle-based simulations.

• •

The paper derives the BER expression for two two-way MCvD systems—one with a half-duplex system and one with a full-duplex system. The theoretical BERs are then validated through particle-based simulations.

• •

The paper proposes two SIC techniques—analog SI cancellation (A-SIC) and digital SI cancellation (D-SIC)—as analytical results show that it is impossible to achieve reliable data transmission without SIC in the full-duplex system. Using the derived channel model and the BER expression, we investigate the optimal values for the normalized detection threshold and the discarding time (i.e., initial part of the received signal is ignored by discarding time $T_{c}$) in order to minimize the BER. The throughputs and BERs of both systems are analyzed. We compare the SIC-adopted full-duplex system with binary concentration shift keying (BCSK) to the half-duplex system with quadrature concentration shift keying (QCSK), which is not provided in [1].

The rest of this paper is organized as follows. In Section II, we discuss the conventional one-way MCvD model. In Section III, we introduce the proposed two-way MCvD system and SIC techniques. In Section IV, we investigate the channel model of a two-way MCvD system with two absorbing sphere receivers. Then we present the channel model verifications and BER formulations for the two considered systems. In Section V, we present performance analysis results in terms of BER and throughput for the two considered systems. Finally, in Section VI, we present our conclusions.

II One-Way Molecular Communication

We start by providing details of the conventional one-way MCvD system. Also, we present details of the single absorbing receiver channel model so that the reader may understand the differences and challenges with respect to the corresponding model for two-way molecular communication system.

II-A Conventional one-way MCvD system

The conventional one-way MCvD system consists of a point source (point transmitter) and an absorbing sphere (receiver). In Fig. 1, the point source (point transmitter), Tx, is separated from the absorbing sphere (receiver), $\text{Rx}_{1}$, the radius of which is denoted by $r_{r}$, by a distance of $d$. Here we focus on three molecular processes: emission, propagation, and reception. The emission process is related to the modulation of the data bits onto the physical properties of the molecules or the emission time[28]. The propagation process is governed by diffusion and flow[29, 30, 31]. In this paper, we only consider diffusion. The reception process is related to the acquisition of the molecules at the receiver and the demodulation of the data bits.

Regarding the propagation process, the interactions between diffusing molecules are ignored since the messenger molecules are assumed to be chemically stable. We assume that the transmitter and the receiver are fully synchronized, which can be achieved by the method introduced in [32].

II-B *Channel Model for One-way Molecular Communication with a Single Receiver *

In diffusion-based systems, propagation of the molecules is governed by Fick’s second law in a 3-D environment; that is

[TABLE]

where $\nabla^{2}$, $p(r,t|r_{0})$, and $D$ are the Laplacian operator, the molecule distribution function at time $t$ and distance $r$ given the initial distance $r_{0}$, and the diffusion constant, respectively. The value of $D$ depends on the temperature, viscosity of the fluid, and the Stokes radius of the molecule [33].

Fig. 1 illustrates a simple topology of one-way MCvD. In [18], the expected channel response of one-way MCvD is presented and analyzed from a channel characteristics perspective. Also, a time-dependent solution for a fraction of molecules hitting a single absorbing sphere ($\text{Rx}_{1}$) by time $t$ is presented, as follows:

[TABLE]

where Tx, $g^{\text{Tx}}_{1}(t)$, $r_{r}$, $d$, and $\text{erfc}(\cdot)$ represent the location of the point transmitter, the instantaneous hitting probability density (i.e., the arrival time distribution for $\text{Rx}_{1}$ when the molecules are emitted from Tx), the radius of the receiver, the distance from Tx to $\text{Rx}_{1}$, and the complementary error function, respectively.

III Two-Way Molecular Communication

In this section, we introduce two modes of operation for the proposed two-way MCvD system, i.e., full-duplex and half-duplex based on time-division. In the former case, severe SI is observed when both transceivers simultaneously modulate signals using the same type of molecule. Therefore, in this case, it is necessary to apply an SIC technique.

III-A Topology

We consider a 3-D environment with two point sources (transmitters) and two absorbing spheres (receivers). Each transmitter emits molecules without directionality. The propagation of molecules is governed by (1). The molecules are immediately absorbed when they reach the surface of any one of the receivers. Since the absorbed molecules are removed from the system, each molecule is detected at most once.

Fig. 2 shows the model of the proposed system. The blue pipes in Fig. 2 represent the connection between the transmitting part and the receiving part. They are assumed to be transparent to the released molecules. $\text{Tx}_{i}$ releases molecules that are intended to be absorbed by $\text{Rx}_{j}$ ($i\neq j$). If the molecules that are released from $\text{Tx}_{i}$ are absorbed by $\text{Rx}_{i}$ (not the desired result), then we call this SI. The distance between $\text{Tx}_{i}$ and $\text{Rx}_{i}$ is denoted by $d_{i}$, and the distance from a point $p$ to $\text{Rx}_{j}$ is denoted by $d_{p}^{\text{Rx}_{j}}$.

The system model illustrated in Fig. 2 is similar to full duplex radios with RF technologies that have antenna separation as passive analog cancellation and digital self interference cancellation. We leave the fabrication issues for the transceivers for our future work and will focus on fundamental analyses of the proposed system.

III-B Communication Model & Modulation

Consider the following two modes of operation for the proposed two-way MCvD system:

• •

Half-duplex system: $\text{Tx}_{1}$ and $\text{Tx}_{2}$ release molecules alternately (i.e., when $\text{Tx}_{1}$ emits molecules, $\text{Rx}_{1}$ and $\text{Rx}_{2}$ receive the molecules but $\text{Rx}_{1}$ does not count the molecules).

• •

Full-duplex system: $\text{Tx}_{1}$ and $\text{Tx}_{2}$ release molecules simultaneously that are intended for $\text{Rx}_{2}$ and $\text{Rx}_{1}$, respectively. Receiver $\text{Rx}_{i}$ receives and counts the molecules that are emitted by $\text{Tx}_{i}$ and $\text{Tx}_{j}$.

In the half-duplex system, at least half of the elements of the bit sequences are not used (i.e., $\text{Tx}_{1}$ does not emits molecules when $\text{Tx}_{2}$ emits molecules), which is not the case for the full-duplex system. Therefore, the ISI and the SI are much more severe in the full-duplex system. For the $n\text{th}$ symbol period, the molecular received signal is composed of $2n$ bits including the current symbols and the previous $2n\!-\!2$ symbols sent from the two transmitters. The bit sequences for the transmitters are denoted by $x_{1}[1\!:\!n]$ and $x_{2}[1\!:\!n]$.

For the modulation, we use binary/quadrature concentration shift keying (BCSK, QCSK)  [28]. We let $N_{1}$ denote the number of molecules for encoding bit-1, and we define that there will be no emission in the case of bit-0 for BCSK. Each of the transmitters has its bit sequences $x_{i}$ to encode, where $x_{i}[k]$ denotes the symbol in the $k\text{th}$ symbol duration for $\text{Tx}_{i}$. We define $P_{ij}[k]$ as the probability that molecules emitted from $\text{Tx}_{i}$ hit $\text{Rx}_{j}$ in the $k\text{th}$ symbol duration after the emission, which is formulated as follows:

[TABLE]

where $t_{s}$ and ${P}_{ij}$(${t}_{1}$,${t}_{2}$) denote the symbol duration and the probability that molecules emitted from $\text{Tx}_{i}$ are absorbed by $\text{Rx}_{j}$ but not $\text{Rx}_{i}$ between time $t_{1}$ and $t_{2}$ after the emission. In the rest of this paper we call $P_{ij}[k]$ the channel coefficient.

In (3), $P_{ij}[0]$ corresponds to the probability of being absorbed in the current symbol slot. We let $y_{\text{Rx}_{j}}[n]$ denote the number of molecules that are absorbed by $\text{Rx}_{j}$ in the $n\text{th}$ symbol slot. Note that $y_{\text{Rx}_{j}}[n]$ can be affected by the number of molecules released from (i) a pair source at the current symbol slot, (ii) a pair source at the previous time slots, (iii) a non pair source at the current symbol slot, and (iv) non pair source at the previous time slots. To formulate $y_{\text{Rx}_{j}}[n]$, we define $N_{ij}[k]$ as follows:

[TABLE]

where $\mathscr{B}(m,p)$ is a binomial distribution with $m$ trials and success probability $p$. Then, $y_{\text{Rx}_{j}}[n]$ can be formulated as follows:

[TABLE]

To consider the misoperations of the receiver, we add the noise term which is assumed to be a Gaussian random variable, $s[n]\sim\mathscr{N}(0,\sigma_{\text{noise}}^{2})$. For the sake of tractability, we approximate the binomial distribution as follows:

[TABLE]

where $\mathscr{N}(\mu,\sigma^{2})$ represents a Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$. Hence, $y_{\text{Rx}_{j}}[n]$ can be expressed as follows (i.e., as a Gaussian random variable, where mean and variance values are dependent upon transmitted bit sequences):

[TABLE]

III-C Self-Interference Cancellation

Since the proposed two-way MCvD system comprises two transceivers that use the same type of molecule, the system’s receivers are unable to distinguish molecules in terms of the transmitting source. For example, if $\text{Tx}_{1}$ sends bit-1 and $\text{Tx}_{2}$ sends bit-0, the molecules are released only from $\text{Tx}_{1}$. However, those molecules can also be absorbed by $\text{Rx}_{1}$, which is not desired. Then $\text{Rx}_{1}$ may decode the received signal as bit-1, even though its paired transmitter $\text{Tx}_{2}$ sends bit-0. In fact, most of the molecules released from $\text{Tx}_{1}$ will be absorbed by $\text{Rx}_{1}$ because $\text{Tx}_{1}$ is much closer to $\text{Rx}_{1}$ than $\text{Rx}_{2}$. Hence, in this case, the number of received molecules is mostly dependent on the transmitted symbol from the unpaired transmitter, which makes for infeasible communication. Therefore, we propose the following two SIC techniques:

• •

Analog self-interference cancellation (A-SIC): the initial part (i.e., between time 0 and $T_{c}$) of the molecular received signal for each symbol slot is ignored (see Fig. 3 for the state diagram). In the rest of this paper, we call $T_{c}$ as discarding time.

• •

Digital self-interference cancellation (D-SIC): we predict the number of SI molecules (i.e., the number of absorbed molecules originating from the unpaired transmitter) from the current bit and subtract it from the molecular received signal

Fig. 4 shows the full-duplex system with A-SIC and D-SIC. Note that the two SIC techniques can be applied separately. The channel coefficient of the system with A-SIC, $P^{\varphi}_{ij}[k]$, is given as follows:

[TABLE]

Hence, $N_{ij}[k]$ of the proposed two-way MCvD system with A-SIC is denoted by $N^{\varphi}_{ij}[k]$ and becomes

[TABLE]

Furthermore, $y_{\text{Rx}_{j}}[n]$ with A-SIC and D-SIC becomes

[TABLE]

$\mathbb{E}\left[\cdot\right]$ is the expectation operation. By applying these two SIC techniques, (7) is transformed into

[TABLE]

IV Channel Model & BER Formulation of Two-Way Molecular Communication

We formulate the BER as a function of detection threshold $\tau$, the number of molecules for encoding bit-1 ($N_{1}$), $P^{\varphi}_{ij}[k]$ ($P^{\varphi}_{ij}[k]$ when SIC is applied), and the symbol duration $t_{s}$. In this section, we first derive the channel model to obtain the channel coefficients (i.e., $P_{ij}[k]$ and $P^{\varphi}_{ij}[k]$) and then utilize these channel coefficients in the BER calculations.

IV-A Channel Model

As we mentioned, the channel model is a time-dependent solution of the fraction of received molecules by time $t$. In the proposed two-way MCvD system, we cannot use (2) directly to derive the channel model since we have to consider the events of molecules being absorbed by each $\text{Rx}_{i}$, which are not independent of each other. The proposed BER formula and the SIC techniques are based on the channel model. However, there is no analytical closed-form solution in the literature for the case of two absorbing spherical receivers. In prior work on molecular MIMO [34, 35], researchers considered two pairs of point transmitters and fully absorbing receivers. To obtain the channel models, the authors in [34] and [35] utilized a one-way MCvD system channel model for a single receiver [18] and fitted the coefficients accordingly. In this paper, we propose a new approach to analytically derive the multi-receiver channel model.

Fig. 5 describes the notations and the cartesian coordinates in the derivation. The centers of two receivers are located on $z$-axis and the distance between them is denoted by $\ell$ (i.e., $(0,0,a)$ and $(0,0,a-\ell)$). Our goal is to obtain the fraction of received molecules by time $t$ for the receivers when $\text{Tx}_{1}$ is the emitter (i.e., $F_{1}^{\text{Tx}_{1}}(t)$ and $F_{2}^{\text{Tx}_{1}}(t)$ for the receivers $\text{Rx}_{1}$ and $\text{Rx}_{2}$, respectively).

In the rest of this subsection, we first derive a closed form of $\lim_{t\to\infty}F_{1}^{\text{Tx}_{1}}(t)$ and $\lim_{t\to\infty}F_{2}^{\text{Tx}_{1}}(t)$, which correspond to the probability that the molecule eventually being absorbed by $\text{Rx}_{1}$ and $\text{Rx}_{2}$, respectively. This so-called asymptotic capture probability was studied in [36]. The authors in [36] derive an exact solution of the asymptotic capture probability for the case of two absorbing spheres, which has an identical radius. In this paper we generalize this solution to release the equal-radius constraint since our system adopts a general topology.

In the derivation of the asymptotic capture probability, we use bispherical coordinate $(\mu,\eta,\phi)$ [37, p. 1298]). Note that bispherical coordinate is obtained by rotating bipolar axes about the line between the two foci, $(0,0,\pm f)$ in cartesian coordinates:

[TABLE]

Theorem 1**.**

The asymptotic capture probabilities in the cases of two receivers, i.e., $\lim_{t\to\infty}F_{1}^{\text{Tx}_{1}}(t)$ and $\lim_{t\to\infty}F_{2}^{\text{Tx}_{1}}(t)$, are as follws:

[TABLE]

where $\mu_{1}\!\!=\!\!\cosh^{-1}\left(\!\frac{\ell-a}{r_{r1}}\!\right)$ and $\mu_{2}\!\!=\!\!\cosh^{-1}\left(\!\frac{a}{r_{r2}}\!\right)$, $P_{m}$ is the $m$-th order Legendre polynomial, and $(\mu_{0},\eta_{0},\phi_{0})$ is the bispherical representation of the cartesian coordinate of $\text{Tx}_{1}$ (i.e., $(x_{0},y_{0},z_{0})$) with two foci $(0,0,\pm f)$,

[TABLE]

Proof.

Let $k_{i}(x_{0},y_{0},z_{0})$ denotes an asymptotic capture probability for $\text{Rx}_{i}$, where $\vec{r}\!\!=\!\!(x_{0},y_{0},z_{0})$ is the initial position of the molecule. The asymptotic capture probability satisfies the Laplace equation [38],

[TABLE]

The boundary conditions in this case are

[TABLE]

[TABLE]

[TABLE]

where $d_{\vec{r}}^{\text{Rx}_{1}}$ and $d_{\vec{r}}^{\text{Rx}_{1}}$ is the distance from the initial point $\vec{r}$ to $\text{Rx}_{1}$ and $\text{Rx}_{2}$, respectively. Without loss of generality, we set $a$ as follows:

[TABLE]

By setting $f$ and $a$ as (14) and (19), the surface of $\text{Rx}_{1}$ and $\text{Rx}_{2}$ are represented as $\mu\!=\!\mu_{1}$ and $\mu=-\mu_{2}$, respectively. This enables us to utilize the typical solutions of the Laplace equation [37, p. 1299] and derive the solution of (15) as (13). ∎

It can be easily shown that when $r_{r2}\to 0$, the asymptotic capture probability for $\text{Rx}_{2}$ (i.e., $\lim_{t\to\infty}F_{2}^{\text{Tx}_{1}}(t)$ in (13)) goes to 0 as $\mu_{2}\to\infty$. Moreover, by using [37, eqs. (10.3.70)], we get the following:

[TABLE]

which is the asymptotic capture probability for the single absorbing receiver case [18].

Based on the derived asymptotic capture probability, we provide a derivation of the approximation of the channel model. We investigate the multi-receiver channel model by considering the possible diffusion paths in the case of a two-way MCvD system. In Fig. 6, three possible traversal paths (diffusion paths) for a molecule emitted from $\text{Tx}_{1}$ are shown. $F_{2}^{\text{Tx}_{1}}(t)$ can be expressed as follows:

[TABLE]

where $f_{2}^{\text{Tx}_{1}}(t)$ denotes the instantaneous hitting probability density for the molecules following Path 2. Note that $\text{Tx}_{2}$ will not be considered in the derivation because we can obtain the whole case by superposition.

Remark 1**.**

$\text{Rx}_{1}$* will be regarded as a non-intended receiver on the path to $\text{Rx}_{2}$ and vice versa.*

Remark 2**.**

To obtain the probability corresponding to Path 2, we subtract the probability corresponding to Path 3 from the channel model of the one-way MCvD with a single point transmitter, which corresponds to the probability of $\{\text{Path\leavevmode\nobreak\ 2}\cup\text{Path\leavevmode\nobreak\ 3}\}$.

Hence, we consider the instantaneous hitting probability densities for $\text{Rx}_{2}$ as follows:

[TABLE]

where $\alpha(t^{\prime})$ corresponds to the instantaneous hitting probability densities for the molecules following Path 3. Additionally, $g^{\text{Tx}_{1}}_{2}(t^{\prime})$ corresponds to the single receiver case that is given in the literature as (2), i.e., the union of Path 2 and Path 3 as if $\text{Rx}_{1}$ is not in the environment. Therefore, if we calculate $\alpha(t^{\prime})$, then we can calculate $f_{2}^{\text{Tx}_{1}}(t^{\prime})$, which in turn will lead us to $F_{2}^{\text{Tx}_{1}}(t^{\prime})$.

By definition, each of the molecules that is moving through Path 3 visits $\text{Rx}_{1}$ at least once. Therefore, we can segment Path 3 into two parts. For those molecules originating from $\text{Tx}_{1}$, we denote the first hitting point on the surface of $\text{Rx}_{1}$ as $s_{{1},{\text{Tx}_{1}}}$ and the corresponding first hitting time as $\tau$. Note that $s_{{1},{\text{Tx}_{1}}}$ can be an arbitrary point on the surface of $\text{Rx}_{1}$ and $\tau$ can be any real value less than $t^{\prime}$.

The instantaneous hitting probability density for $s_{{1},{\text{Tx}_{1}}}$ on the surface of $\text{Rx}_{1}$ (i.e., the arrival time distribution for $s_{{1},{\text{Tx}_{1}}}$ when the molecules are emitted from Tx) is denoted by $f_{1}^{\text{Tx}_{1}}(t^{\prime},\!s_{{1},{\text{Tx}_{1}}}\!)$. Since $s_{{1},{\text{Tx}_{1}}}$ is the first hitting point on the surface of $\text{Rx}_{1}$, it can be regarded as the starting point of the successive path to $\text{Rx}_{2}$. Hence, (22) can be rewritten as

[TABLE]

where $\Omega_{1}$ indicates the points on the surface of $\text{Rx}_{1}$. To obtain $f_{2}^{\text{Tx}_{1}}(t^{\prime})$, we also need to consider $f_{1}^{\text{Tx}_{1}}(t^{\prime})$ in a similar way. Thus, we have the following:

[TABLE]

where $s_{{2},{\text{Tx}_{1}}}$ is the first hitting point on $\text{Rx}_{2}$ that is analogous to $s_{{1},{\text{Tx}_{1}}}$ for $\text{Rx}_{1}$.

When we apply the mean value theorem for integration to the surface integration in (23), we get

[TABLE]

where $s^{\prime}_{1}(t^{\prime},\tau)$ is a point $\in\Omega_{1}$ . After the surface integration in (25), we obtain

[TABLE]

In the same way, we can express $f_{1}^{\text{Tx}_{1}}(t^{\prime})$ as

[TABLE]

where $s^{\prime}_{2}(t^{\prime},\tau)$ is a point $\in\Omega_{2}$.

Note that to derive $s^{\prime}_{i}(t^{\prime},\tau)$ requires the closed forms of $f_{i}^{\text{Tx}_{1}}(t)$, therefore, it is unattainable to solve the simultaneous equations (26) and (27). Therefore, we set $s^{\prime}_{i}(t^{\prime},\tau)$ as a constant $s^{\ell}_{i}\!\in\!\Omega_{i}$, and then derive the approximation of $f_{i}^{\text{Tx}_{1}}(t)$ by solving (26) and (27). We determine the two constants $s^{\ell}_{1}$ and $s^{\ell}_{2}$ by matching the asymptotic values of approximation forms and (13). By substituting $s^{\prime}_{i}(t^{\prime},\tau)$ with $s^{\ell}_{i}$ and integrating (26) and (27), we obtain the following:

[TABLE]

Then, we get

[TABLE]

Next, we take $t\to\infty$ in (28), which gives us

[TABLE]

In (30), $G^{s^{\ell}_{1}}_{2}(t)$ and $G^{s^{\ell}_{1}}_{2}(t)$ are expressed as follows:

[TABLE]

where $d_{s^{\ell}_{1}}^{\text{Rx}_{2}}$ is the distance from $s^{\ell}_{1}$ to $\text{Rx}_{2}$, and $d_{s^{\ell}_{2}}^{\text{Rx}_{1}}$ is the distance from $s^{\ell}_{2}$ to $\text{Rx}_{1}$. By solving (30) simultaneously, we obtain

[TABLE]

We obtain $d_{s^{\ell}_{1}}^{\text{Rx}_{2}}$ and $d_{s^{\ell}_{2}}^{\text{Rx}_{1}}$ from the fact that (13) is equal to (32) (See Appendix A).

Fianlly, we derive the approximation forms of $F_{1}^{\text{Tx}_{1}}(t)$ and $F_{2}^{\text{Tx}_{1}}(t)$ based on (29), as follows (see Appendix B):

[TABLE]

where $A=\frac{(r_{r1}+d_{s^{\ell}_{2}}^{\text{Rx}_{1}})(r_{r2}+d_{s^{\ell}_{1}}^{\text{Rx}_{2}})}{r_{r1}r_{r2}}$, $B=-\left(d_{s^{\ell}_{1}}^{\text{Rx}_{2}}+d_{s^{\ell}_{2}}^{\text{Rx}_{1}}\right)/\sqrt{D}$, $b_{1}=\left(d_{s^{\ell}_{2}}^{\text{Rx}_{1}}+d_{\text{Tx}_{1}}^{\text{Rx}_{2}}\right)/\sqrt{D}$, $b_{2}=\left(d_{s^{\ell}_{1}}^{\text{Rx}_{2}}+d_{\text{Tx}_{1}}^{\text{Rx}_{1}}\right)/\sqrt{D}$, $c_{1,1}=\frac{Ar_{r1}}{(A-1)(r_{r1}+d_{\text{Tx}_{1}}^{\text{Rx}_{1}})}$, $c_{1,2}=\frac{Ar_{r1}r_{r2}}{(A-1)(r_{r2}+d_{\text{Tx}_{1}}^{\text{Rx}_{2}})(r_{r1}+d_{s^{\ell}_{2}}^{\text{Rx}_{1}})}$, $c_{2,1}=\frac{Ar_{r2}}{(A-1)(r_{r2}+d_{\text{Tx}_{1}}^{\text{Rx}_{2}})}$, $c_{2,2}=\frac{Ar_{r1}r_{r2}}{(A-1)(r_{r1}+d_{\text{Tx}_{1}}^{\text{Rx}_{1}})(r_{r2}+d_{s^{\ell}_{1}}^{\text{Rx}_{2}})}$, and \big{[}H(u)\big{]}^{a}_{b}=H(a)-H(b).

IV-B Channel Model Verification

In this Section, we verify our channel model by particle-based simulation. We repeat the simulation for a number of times, and take the mean value with respect to the number of trials. The number of trials is presented in Tables II and III. In each simulated trial, 50000 molecules are released. The received molecules are distinguished according to the transmitter that emits them. In the simulations, only $\text{Tx}_{1}$ releases molecules which is sufficient to verify the theoretical channel model.

To implement Brownian motion for the emitted molecules, our simulator records and updates the position of each molecule at each time step $\Delta t$. The position of the emitted molecules, $\mathbf{X}_{p}(t)$, changes by $\Delta\mathbf{X}_{p}$ after simulation time step $\Delta t$ as in (34)[39]. The simulation parameters used for verification of the channel model are given in Table I.

[TABLE]

Through extensive simulations, we obtain the fraction of received molecules for each receiver (i.e., $\text{Rx}_{1}$ and $\text{Rx}_{2}$) at each time step during the simulation time, i.e., $F_{2,\text{sim}}^{\text{Tx}_{1}}(t)$ and $F_{1,\text{sim}}^{\text{Tx}_{1}}(t)$, respectively. Then, we use (33) to calculate $F_{2}^{\text{Tx}_{1}}(t)$ and $F_{1}^{\text{Tx}_{1}}(t)$. We compare the analytical results with the simulation results, as shown in Fig. 7. It is worth noting that the asymptotic capture probabilities (i.e., $\lim_{t\to\infty}F_{1}^{\text{Tx}_{1}}(t)$ and $\lim_{t\to\infty}F_{2}^{\text{Tx}_{1}}(t)$ in (13)) for Fig. 7 are $0.6414$ and $0.2932$, respectively. By the long-term simulation (i.e., $t=1000s$), we observe that the simulation results and (33) converge to those values. We depict the number of received molecules in Fig. 8, which corresponds to the impulse response of the channel. Theoretical values in Fig. 8 at $t=k\Delta t$ are obtained by calculating

[TABLE]

Using the derived channel model, we substitute $P_{ij}[k]$ and $P^{\varphi}_{ij}[k]$ as follows:

[TABLE]

In Fig. 9, we depict the comparison of the channel coefficients $P_{ij}[k]$.

IV-C BER Formula for Two-Way Molecular Communication

In this section, we formulate the BER in terms of the Q-function (i.e., the tail probability of the standard normal distribution). For the receiver $\text{Rx}_{j}$, an error occurs when the result of decoding is different from the bit transmitted from the $\text{Tx}_{i}$. If $\text{Tx}_{i}$ encodes bit-1, an error occurs when $y_{\text{Rx}_{j}}[n]$ is less than the detection threshold $\tau_{d}$. If $\text{Tx}_{i}$ encodes bit-0, then an error occurs when $y_{\text{Rx}_{j}}[n]$ is greater than the detection threshold $\tau_{d}$.

Considering the transmitted bit sequences $x_{i}$ and $x_{j}$, we obtain the error probabilities at the $n\text{th}$ symbol slot as (37) where x_{i}[1\!:\!n\!\scalebox{0.75}[1.0]{,-,}\!1] denotes the bits transmitted previously from $\text{Tx}_{i}$:

[TABLE]

where $\mu_{n,\text{total}}$ and $\sigma^{2}_{n,\text{total}}$ are defined in (7), $Q(\cdot)$ is the Q-function, and

[TABLE]

Note that \sum_{x_{i}[1:n\!\scalebox{0.75}[1.0]{,-,}\!1]} and $\sum_{x_{j}}$ in (37) indicate the summation over all possible bit sequences x_{i}[1\!\!:\!\!n\!\scalebox{0.75}[1.0]{,-,}\!1] and $x_{j}[1\!\!:\!\!n]$, respectively. The BER expression for the full-duplex system with the proposed SIC techniques is obtained by substituting $y_{\text{Rx}_{j}}[n]$, $\mu_{n,\text{total}|x_{i},x_{j}}$ and $\sigma^{2}_{n,\text{total}|x_{i},x_{j}}$ in (38) with $y^{\varphi}_{\text{Rx}_{j}}[n]$, $\mu^{\varphi}_{n,\text{total}|x_{i},x_{j}}$, and $(\sigma^{\varphi}_{n,\text{total}})^{2}$, respectively.

V Particle-Based Simulation Results

In this section, we compare the proposed half-duplex and full-duplex systems in terms of BER and throughput. Throughput of the systems are evaluated as follows:

[TABLE]

where $M$ is the number of bits transmitted in one symbol, ${P}_{e}$ is the BER of the system, and $t_{s}$ denotes the symbol duration. The system parameters for the rest of our work are summarized in Table II. For convenience, we denote $t_{s}$ of the half-duplex and full-duplex systems as $t^{\text{HD}}_{s}$ and $t^{\text{FD}}_{s}$, respectively. In Table III and Table IV, throughput ratio is calculated as a ratio of the throughput of full-duplex system and the half-duplex system.

In the half-duplex system, each receiver operates only when the paired transmitter releases the molecular signal. Hence, the operating time of the receiver (i.e., detection period) is half of the symbol duration. In the full-duplex system, the detection period of each receiver is equal to the symbol duration. Roughly, we can expect faster but less accurate communications in the full-duplex system if we use the same modulation technique and detection period for both systems. By the theoretical and simulation BER analysis, we confirm that the proposed SIC techniques are necessary in the full-duplex system. Therefore we analyze the BER improvement in the full-duplex system with SIC to find, numerically, the optimal values for the normalized detection threshold and the discarding time of the SIC to minimize BER. Through the numerical parameter optimization, the throughput of the full-duplex system with SIC is compared to the half-duplex system. For a fair comparison, we evaluate the throughput in the following three cases, considering that the throughput is a function of M, $t_{s}$, and $P_{e}$:

Half-duplex system (BCSK) vs. full-duplex system (BCSK) with SIC, where $t^{\text{FD}}_{s}\!=\!t^{\text{HD}}_{s}/2$. 2. 2.

Half-duplex system (BCSK) vs. full-duplex system (BCSK) with SIC. We set the same $t_{s}$ for both systems (i.e., $t^{\text{FD}}_{s}\!=\!t^{\text{HD}}_{s}$). 3. 3.

Half-duplex system (BCSK) vs. full-duplex system (BCSK) with SIC. We empirically adjust $t^{\text{FD}}_{s}$ to make the BER of both systems close. 4. 4.

Half-duplex system using quadrature concentration shift keying (QCSK) vs full-duplex system (BCSK) with SIC, where $t^{\text{FD}}_{s}\!=\!t^{\text{HD}}_{s}/2$.

V-A BER Analysis

Fig. 10(a) depicts the simulation and theoretical BERs of the half-duplex system using BCSK. The x-axis is the normalized threshold ($\tau_{m}$), which is $\tau_{d}/N_{1}$. First of all, the simulation and theoretical values match each other well. Since the half-duplex system is not susceptible to SI, we do not need to apply the proposed SIC techniques to this system. On the other hand, we observe from Fig. 11 that the BER of the full-duplex system is nearly 0.3 if we do not apply the SIC techniques. In Fig. 10(a), we can see the optimal normalized threshold $\tau_{m}^{*}$ for different $t^{\text{HD}}_{s}$ and we observe that it is slightly changing according to the value of $t^{\text{HD}}_{s}$. We also observe that the BER gain is relatively higher for changing $t^{\text{HD}}_{s}$ from $0.100\text{\,}\mathrm{s}$ to $0.150\text{\,}\mathrm{s}$ compared to from $0.150\text{\,}\mathrm{s}$ to $0.200\text{\,}\mathrm{s}$ due to the relative ISI difference.

Fig. 10(b) shows the simulation and theoretical BER of the full-duplex system with D-SIC. We can see that there is an optimal normalized detection threshold $\tau_{m}^{*}$ for different $t^{\text{FD}}_{s}$ values. Similar to the case of the half-duplex system, there is some similarity between the $\tau_{m}^{*}$ changes according to the $t^{\text{FD}}_{s}$ value and also the tendency of the BER gain with respect to

the $t^{\text{FD}}_{s}$ difference (see Fig. 10(a)).

Fig. 11 depicts the simulation and theoretical BERs of the full-duplex system while applying the different SIC techniques. The x-axis is the normalized detection threshold $\tau_{d}/N_{1}$. The performance of D-SIC is superior to the performance of A-SIC, because D-SIC deals with a wider time interval (i.e., when the discarding time is $T_{c}$, A-SIC ignores the molecules received in $[0,T_{c}]$, and D-SIC subtracts the expected number of the SI molecules received in $[T_{c},t_{s}]$). For the no-SIC case, we assume that $\text{Tx}_{i}$ is located on the surface of $\text{Rx}_{i}$ (i.e., $d_{1}$=$d_{2}$=$0\text{\,}\mathrm{\SIUnitSymbolMicro m}$), and we observe that the BER is nearly 0.5. The results indicate that the separation of transmitter and receiver is necessary in the two-way MCvD system. We observe that the BER of the full-duplex system with A-SIC is nearly 0.3. On the other hand, the BER of the full-duplex system with D-SIC becomes comparable with the half-duplex system in Fig.10(a). Moreover, we observe that when we apply both SIC techniques, the BER is slightly improved compared to the full-duplex system with only D-SIC. As we stated in (7) and (11), the number of received SI molecules from the current symbol can be modeled as a Gaussian random variable that has a variance $N_{1}P^{\varphi}_{ii}[0](1-P^{\varphi}_{ii}[0])x_{i}[n]$. Since the variance decreases as $P^{\varphi}_{ii}[0]$ decreases (when $P^{\varphi}_{ii}[0]\leq 0.5$), applying A-SIC can reduce the uncertainty of predicting the number of received SI molecules from the current symbol. For the full-duplex systems, since the BER of the A-SIC-only system and the no-SIC system are not reasonable, we consider those full-duplex systems with only D-SIC or D-SIC and A-SIC.

As derived in Section IV, the BER is a function of symbol duration ($t_{s}$), detection threshold ($\tau_{d}$), the number of molecules for encoding bit-1 ($N_{1}$), and the discarding time of molecular received signal for A-SIC ($T_{c}$). Since $N_{1}$ and $t_{s}$ are system parameters, we consider only $\tau_{d}$ and $T_{c}$ as variables to optimize. While we will show that optimal values for these parameters exist, we cannot derive them in closed-form and hence must resort to evaluating them numerically. To improve BER, we first need to see the structure of the BER of the full-duplex system with SIC in terms of $\tau_{d}$ and $T_{c}$ by the following analysis.

In Fig. 12, we depict a heatmap of the theoretical BER with respect to $T_{c}$ and $\tau_{m}$ for the full-duplex system with D-SIC and A-SIC. We observe that $\tau_{m}$ minimizing the BER with given $T_{c}$ decreases as $T_{c}$ increases. It is because that the molecules from the paired transmitter (i.e., desired signal) are also discarded by A-SIC. For each $t^{\text{FD}}_{s}$, we find global optimal normalized detection threshold and the discarding time in order to minimize the BER and denote them as red mark on the heatmap. We utilized these optimal values in the SIC algorithm to compare the throughput of the half-duplex system and the full-duplex system with SIC.

As was mentioned before, for comparison, we consider the following three cases: i) set $t^{\text{FD}}_{s}\!=\!t^{\text{HD}}_{s}/2$ (to observe a trade off between the throughput and BER); ii) set the same $t_{s}$ for both systems; iii) fix $t^{\text{HD}}_{s}$ and adjust $t^{\text{FD}}_{s}$ empirically to make the BER of both systems close. For the case iii), adjusted $t^{\text{FD}}_{s}$ values are in Table III.

For the case i), Fig. 13 shows that the throughputs of the full-duplex system with optimized SIC are almost double the throughput of the half-duplex system. Thus, we can achieve nearly double the transmission rate using the proposed SIC techniques without degrading the BER significantly. When $t^{\text{HD}}_{s}\!=\!t^{\text{FD}}_{s}$ is $0.100\text{\,}\mathrm{s}$ or $0.200\text{\,}\mathrm{s}$, the BER of the full-duplex system with optimized SIC is less than the BER of the half-duplex system. The converse is true when $t^{\text{HD}}_{s}\!=\!t^{\text{FD}}_{s}$ is $0.300\text{\,}\mathrm{s}$ or $0.400\text{\,}\mathrm{s}$. Adjusted $t^{\text{FD}}_{s}$ throughput values in Table III show the same tendency. The overall results imply that the full-duplex system with optimized SIC becomes better than the half-duplex system even in terms of BER when $t_{s}$ get decreased.

In Fig. 14 we depict the simulation and theoretical BER of the half-duplex system using QCSK and compare them to the full-duplex system with optimized SIC using BCSK where $t^{\text{FD}}_{s}\!=\!t^{\text{HD}}_{s}/2$. In this case, $M/t_{s}$ is the same for both systems. Hence, the BER determines the difference between the throughputs. For QCSK, we used an equally spaced number of molecules for encoding different bits (i.e., bit-0, 1, 2, 3) and three thresholds (i.e., $\tau_{d1}$, $\tau_{d2}$, $\tau_{d3}$) to detect the molecular received signal. Fig. 14 shows that the BER of the half-duplex system using QCSK is much higher than that of the full-duplex system using BCSK with optimized SIC. The theoretical BERs of the half-duplex system using QCSK are calculated by using a straightforward extension of (37). The throughput difference between the two systems can be seen in Table IV.

VI Conclusion

In this paper, we have investigated two different communication models of two-way MCvD—a half-duplex system and a full-duplex system. We derived the multi-receiver channel model. We also derived the BER formula and verified the formula by simulation. Theoretical analysis and simulations showed that severe SI occurs in the full-duplex system. Therefore, we proposed two SIC techniques to mitigate this interference: A-SIC and D-SIC. We analyzed the BER improvements in the full-duplex system with the proposed SIC techniques and numerically found the optimal values for the normalized detection threshold and the discarding time in order to minimize the system BER. To compare the half-duplex system with the full-duplex system, we evaluated the throughput in three different cases. The throughput of the full-duplex system with optimized SIC increased to more than that of the half-duplex system as $t_{s}$ decreased. With the proposed SIC techniques, we showed the possibility of full-duplex molecular communication using a single type of molecule. On the other hand, the BER analysis and simulation results revealed that using a concentration-based modulation technique of higher order significantly degrades the BER. Investigating a more effective modulation technique for the two-way MCvD will be a topic for the future work. We will also study the channel model for more general system configurations.

VII Appendix

VII-A Derivation of $d_{s^{\ell}_{1}}^{\text{Rx}_{2}}$ and $d_{s^{\ell}_{2}}^{\text{Rx}_{1}}$

Using (2), we rewrite (30) and combine with (13) as

[TABLE]

where $k_{1}$ and $k_{2}$ indicate $k_{1}(\mu,\eta,\phi)$ and $k_{1}(\mu,\eta,\phi)$, respectively. Note that (40) can be rewritten as quadratic equations of $d_{s^{\ell}_{1}}^{\text{Rx}_{2}}$ and $d_{s^{\ell}_{2}}^{\text{Rx}_{1}}$, and the solutions of the equations are derived as the following:

[TABLE]

VII-B Derivation of (33)

Through performing Laplace transform in (29) and solving simultaneous equations, we get

[TABLE]

where $F_{2}^{\text{Tx}_{1}}(s)$ and $F_{1}^{\text{Tx}_{1}}(s)$ are the Laplace transform of $F_{2}^{\text{Tx}_{1}}(t)$ and $F_{1}^{\text{Tx}_{1}}(t)$, respectively. To derive the inverse Laplace transforms of (42), we first derive the inverse Laplace transforms of $F_{2}^{\text{Tx}_{1}}(s^{2})$ and $F_{1}^{\text{Tx}_{1}}(s^{2})$, where

[TABLE]

Using [40, eqs. (5.1), (5.57), (1.18)], the inverse laplace transforms of (43) can be expressed as the following:

[TABLE]

where

[TABLE]

In (45), $u(t)$ is the unit step function. From (45), we derive the inverse Laplace transform of (42) as the following [40, eqs. (1.27)]:

[TABLE]

which give us (33).