Optimizations for Mobile MIMO Relay Molecular Communication via Diffusion with Network Coding

Abstract

We investigate mobile multiple-input multiple-output (MIMO) molecular communication via diffusion (MCvD) system which is consisted of two source nodes, two destination nodes and one relay node in the mobile three-dimensional channel. First, the combinations of decode-and-forward (DF) relaying protocol and network coding (NC) scheme are implemented at relay node. The adaptive thresholds at relay node and destination nodes can be obtained by maximum a posteriori (MAP) probability detection method. Then the mathematical expressions of the average bit error probability (BEP) of this mobile MIMO MCvD system based on DF and NC scheme are derived. Furthermore, in order to minimize the average BEP, we establish the optimization problem with optimization variables which include the ratio of the number of emitted molecules at two source nodes and the initial position of relay node. We put forward an iterative scheme based on block coordinate descent algorithm which can be used to solve the optimization problem and get optimal values of the optimization variables simultaneously. Finally, the numerical results reveal that the proposed iterative method has good convergence behavior. The average BEP performance of this system can be improved by performing the joint optimizations.

1. Introduction

In molecular communication via diffusion (MCvD), the information transmission is finished via molecules based on diffusion in the fluidic environment [1]. MCvD system is considered as a promising communication paradigm in the field of microscale communication among nanomachines due to the biological compatibility and energy-efficient technique [2]. One of the most important applications of MCvD in the area of biomedical domain, is targeted drug delivery, which has been laying the foundation for the advancement of internet of bio-nano things (IoBNT) [3-5].

Recently, many researchers pay close attentions to the mobile MCvD because of nodes mobilities. In 2018, Ahmadzadeh et al. [6] established a three-dimensional (3D) time-varying random channel model of point-to-point mobile MCvD, and proposed a mathematical framework model to study the channel impulse response (CIR). Lin et al. [7] presented the capacity analysis of mobile one-hop MCvD system under two factors including inter-symbol interference (ISI) and noise. In 2019, Cao et al. [8] studied a 3D mobile MCvD system with active absorption nanomachine and deduced the characteristic of CIR. Varshney et al. [9] considered a 1D multi-hop mobile MCvD model that can be applied to human intravascular communication, and proposed the optimal threshold detection rules. Chouhan et al. [10] studied the optimization problem of obtaining the optimal detection threshold and the optimal number of emitted molecules for a 1D multi-hop mobile MCvD model. In 2020, Wang et al. [11] used depleted molecule shift keying (D-MoSK) modulation to suppress the ISI in a mobile diffusive-drift MCvD model. Huang et al. [12] used the estimated initial distance to realize signal detection in point-to-point mobile MCvD system.

In this paper, we consider mobile multiple-input multiple-output (MIMO) MCvD with one relay and propose an iterative algorithm based on block coordinate descent algorithm (BCDA) to jointly optimize the number of emitted molecules at two source nodes and the initial position of relay node for this system. The motivation of this paper is based on the following considerations:

(1) To the best of our knowledge, the mobile MIMO MCvD with one relay has not been yet investigated. In order to improve the transmission reliability, a network coding (NC) scheme [13-14] is used at relay node. The positioning of relay node in the dimension of micrometers is realistic and can be used in the area of drug delivery. In addition, how to get the optimal number of emitted molecules at two source nodes and the position of relay node for minimizing the average bit error probability (BEP) is a challenge work.

(2) Compared with the recent work in [15], there is only one optimization variable to get the position of relay node. For the mobile MIMO MCvD system, there are three optimization variables which represent the ratio relationship of initial distance from different axes to compute the coordinate of relay node. Then there are four optimization variables including the ratio of number of emitted molecules at two source nodes and the coordinates from three axes of relay node to be optimized to achieve the minimum value of average BEP. It is still an open problem in this mobile MIMO MCvD system.

In our paper, we make the following contributions:

(1) The decode-and-forward (DF) relay protocol and NC scheme are implemented at relay node R which can use fewer time slots to complete one bit information transmission compared to the scheme without NC. The adaptive thresholds at relay node and destination nodes can be obtained by maximum a posteriori (MAP) probability detection scheme.

(2) The mathematical expression of the average BEP of the mobile MIMO MCvD system with one relay is derived. Then we put forward an optimization problem for minimizing the average BEP and solve it by using an iterative algorithm based on BCDA method in order to obtain the values of these optimization variables which include the optimal ratio of the number of emitted molecules at two source nodes and the optimal initial position of relay node.

(3)The numerical results reveal that this iterative algorithm has good convergence behavior. In particular, when the system parameters are given, the optimal ratio of the number of emitted molecules at two source nodes and the optimal initial position of relay node can be determined simultaneously. Then the performance of this system is improved by solving the optimization problem.

We organize the rest of the paper as follows. Section 2 introduces the overview of related works. Section 3 describes the MIMO MCvD system with one relay using DF relay protocol and NC scheme in 3D environment. The optimization problem is implemented in Section 4. The average BEP performance of mobile MIMO MCvD system is evaluated in Section 5. We conclude the paper in Section 6.

2. Related Works

In this section, existing works on the optimization and network coding in MCvD literature are introduced.

Many researchers discussed how to optimize decision thresholds at each nanomachine to achieve better performance of MCvD system. In [16], Tavakkoli et al. used Bisection methods to solve the problem of optimizing threshold at destination of two-hop MCvD system in 1D channel with drift. On this basis, they also proposed the bisection algorithm for joint optimizations of relay location and destination threshold for two-hop MCvD system in 3D channel [15]. Fang et al. [17] optimized the error performance of one-hop MCvD system with multiple receiver nanomachines to search the optimal decision thresholds. Varshney et al. [18] proposed the optimal decision rules at each cooperative nanomachine and destination nanomachine for mobile multi-hop MCvD system to get the decision threshold. The authors in [19]-[20] focused on different MCvD models and studied the optimization problems to obtain optimal decision threshold and optimal number of transmitted molecules.

The existing works about NC focused on MCvD system which is consisted of two nodes communicating via a relay node. Unluturk et al. [21] studied the tradeoff between rate and delay of static MCvD system with three nodes by using NC mechanism. Akdeniz et al. [13] proposed two novel NC schemes and obtained the analytical error probability expressions of the MCvD system. Farahnak-Ghazani et al. [14] realized a novel molecular physical-layer NC by using chemical reactions which were explored to mitigate noise and ISI. Kwak et al. [22] considered a two-way MCvD system which was composed of two transmitters and two receivers without relay and then proposed analog self-interference (SI) and digital-SI cancellation strategies for this system.

In order to improve the performance of mobile MIMO MCvD system, the optimization method and network coding at the relay node are jointly used.

3. The Mobile MIMO MCvD System With One Relay

The mobile MIMO MCvD system with one relay by using DF protocol and NC scheme is considered in our work. Fig. 1 shows the system model. Two source nodes Tx₁ and Tx₂ want to transmit information to two destination nodes Rx₁ and Rx₂, respectively. We suppose the node Tx₁ is not in the transmission range of node Rx₁ and node Tx₂ is not in the transmission range of node Rx₂. We assume that the overall fluid environment is large enough and it is assumed that five nodes are passive observers. Each node is a sphere with fixed radius r_q and volume V_q, here q represents nodes Tx_p (p=1, 2), Rx_p (p=1, 2) and relay node R. Another assumption is that the source nodes, relay node and destination nodes can be perfectly synchronized.

Fig. 1. The mobile MIMO MCvD system with one relay by using network coding.

In the process of transmission, we adopt half-duplex communication protocol. The relay node R does not transmit and receive molecules in the same time slot. When the scheme without NC is adopted in this MIMO MCvD system with one relay, it needs four time slots to complete one bit information transmission (Tx₁ → R , R → Rx₁ , Tx₂ → R , R → Rx₂). For the DF-NC scheme, one bit transmission can be completed in two time slots. In this MIMO MCvD system, node Tx₁ transmits one bit W_Tx1^j to Rx₁ and node Tx₂ transmits one bit W_Tx2^j to Rx₂. We give the following transmission process.

Step 1, the source node Tx₁ transmits information W_Tx1^j with type A₁ molecules and node Tx₂ transmits information W_Tx2^j with type A₂ molecules to node R in time slot j, respectively. Then the molecules transmitted by nodes Tx₁ and Tx₂ can arrive at node R in time slot j. At the same time, the nodes Rx₁ and Rx₂ can receive the information from nodes Tx₂ and Tx₁, respectively.

Step 2, the received information from Tx₁ is decoded by node R and Tx₂ into Ŵ_Tx1^j and Ŵ_Tx2^j , respectively. Node R broadcasts W_R^j by releasing molecules with type A₃ which is the XOR result of Ŵ_Tx1^j and Ŵ_Tx2^j to node Rx₁ and Rx₂ at the same time in time slot (j+1).

Step 3, node Rx₁ decodes the information from the relay node as W_Tx1^j ⊕ W_R^j and finds the information transmitted by node Tx₁. By using the same method, node Rx₂ can obtain the information transmitted by node Tx₂.

We consider the optimizations of the initial position of node R and the optimal number of emitted molecules at node Tx₁ and node Tx₂ in this mobile MIMO MCvD system. On one hand, when node Tx₁ and node Tx₂ both send the information bit 1, the number of molecules from node Tx₁ and node Tx₂ which arrive at relay node R depends on the initial position of node R. In such a case, how to get the initial position of node R in order to minimize the average BEP is studied. On the other hand, the optimal ratio of the number of emitted molecules at node Tx₁ and node Tx₂ used to minimize the average BEP is optimized. Therefore, the ratio of the number of emitted molecules at Tx₁ and Tx₂ and the initial position of node R can be optimized simultaneously by using the iteration algorithm based on BCDA.

4. MIMO MCvD System with One Relay

Let r_R and V_R be the radius and volume of node R, respectively. d_{(Tx1, R)}(t) represents the distance between node Tx₁ and node R at time t. For the case that node Tx₁ and node R are mobile and the molecules are released by node Tx₁ at time t=0, the probability of receiving one molecule at time t>0 for the link Tx₁→R is derived by [11]

\(\begin{aligned}h\left(t, \tau_{s}\right)=\frac{V_{R}}{\left(4 \pi D_{1} \tau_{s}\right)^{\frac{3}{2}}} \exp \left(-\frac{d_{\left(\mathrm{Tx}_{1}, R\right)}^{2}(t)}{4 D_{1} \tau_{s}}\right), \end{aligned}\)       (1)

where D₁ describes the relative motion between the molecules and node R. D₁ = D_A1 + D_R, here D_A1 and D_R are diffusion coefficients of type A₁ molecules and node R, respectively. τ_s is the relative time with respect to t. Let 1 d_(Tx1,R)⁰ represent the initial distance between two nodes Tx₁ and R at t=0. Based on [26], µ_h(t,τ_s) is the mean of h(t,τ_s) in (1) and it is obtained by

\(\begin{aligned}\mu_{h}\left(t, \tau_{s}\right)=\frac{V_{R} \exp \left(-\frac{\left(d_{\left(\mathrm{Tx}_{1}, R\right)}^{0}\right)^{2}}{4\left(D_{1} \tau_{s}+D_{2} t\right)}\right)}{\left(4 \pi\left(D_{1} \tau_{s}+D_{2} t\right)\right)^{3 / 2}},\end{aligned}\)       (2)

where D₂ = D_Tx1 + D_R. Let w = D₁τ_s + D₂t and n = D₁τ_s + 2D₂t. We use 2 σ_h²(t , τ_s) to represent the variance of h(t , τ_s ) , which is derived by [23]

\(\begin{aligned}\sigma_{h}^{2}\left(t, \tau_{s}\right)=\frac{\left(V_{R}\right)^{2}}{(4 \pi)^{3}}\left(\frac{\exp \left(-\frac{\left(d_{\left(\mathrm{Tx}_{1}, R\right)}^{0}\right)^{2}}{2 n}\right)}{\left(2 w n-n^{2}\right)^{3 / 2}}-\frac{\exp \left(-\frac{\left(d_{\left(\mathrm{Tx}_{1}, R\right)}^{0}\right)^{2}}{2 w}\right)}{w^{3}}\right).\end{aligned}\)       (3)

4.1 Derivation of the mathematical expression of the average BEP

W_Tx1^j represents the transmitted bit 1 or 0 at node Tx₁ in time slot j. N_(Tx1,R)^C[j] is defined as the number of molecules arriving at node R which are released by node Tx₁ both in time slot j. Then we have

N_(Tx1,R)^C[j] = N_A1[j]W_Tx1^jh(jT_s,τ_s),       (4)

where N_A1[j] represents the transmission of bit 1 at node Tx₁ by using type A₁ molecules. The node Tx₁ transmits the information bit W_Tx1^j in time slot j. N_(Tx1,R)^ISI[j] is defined as the number of ISI molecules arriving at node R. It is expressed by

\(\begin{aligned}N_{\left(\mathrm{Tx}_{1}, R\right)}^{I S I}[j]=\sum_{i=1}^{j-1} N_{A_{1}}[i] W_{\mathrm{Tx}_{1}}^{i} h\left(i T_{s},(j-i) T_{s}+\tau_{s}\right).\end{aligned}\)       (5)

N_(Tx1,R)^A₁[j] represents the number of emitted molecules at node Tx₁ which arrive at node R in time slot j. It is composed of three parts including N_(Tx1,R)^C[j], N_(Tx1,R)^ISI[j] and the noise for the link Tx₁ → R . Then we have

N_(Tx1,R)^A₁[j] = N_(Tx1,R)^C[j] + N_(Tx1,R)^ISI[j] + N_(Tx1,R)^Noise,       (6)

where N_(Tx1,R)^Noise is the noise generated for the link Tx₁ → R which follows Normal distribution 𝒩(µ_(Tx1,R)^Noise, (σ_(Tx1,R)^Noise)²) with mean µ_(Tx1,R)^Noise = 0 and variance (σ_(Tx1,R)^Noise)² [24]. Let W_Tx1ⁱ be the transmitted bit by node Tx₁ in time slot i. According to (4)-(6), we can get

\(\begin{aligned} N_{\left(\mathrm{Tx}_{1}, R\right)}^{A_{1}}[j] & =N_{A_{1}}[j] W_{\mathrm{Tx}_{1}}^{j} h\left(j T_{s}, \tau_{s}\right)+\sum_{i=1}^{j-1} N_{A_{1}}[i] W_{\mathrm{Tx}_{1}}^{i} h\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)+N_{\left(\mathrm{T}_{1}, R\right)}^{\text {Noise }} \\ & =\sum_{i=1}^{j} N_{A_{1}}[i] W_{\mathrm{Tx}_{1}}^{i} h\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)+N_{\left(\mathrm{Tx}_{1}, R\right)}^{\text {Noise }} .\end{aligned}\)       (7)

Similarly, we can get the number of molecules arriving at node R with type A₂ in time slot j for the link Tx₂ → R which is defined as N_(Tx2,R)^A₂[j]. It is computed by

\(\begin{aligned}N_{\left(\mathrm{Tx}_{2}, R\right)}^{A_{2}}[j]=\sum_{i=1}^{j} N_{A_{2}}[i] W_{\mathrm{Tx}_{2}}^{i} h\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)+N_{\left(\mathrm{Tx}_{2}, R\right)}^{\text {Noise }},\end{aligned}\)      (8)

where W_Tx2ⁱ is the transmitted bit at node Tx₂ in the i-th time slot. N_(Tx2,R)^Nosce is defined as the noise for the link Tx₂ → R . N_Rxp^A₃[j + 1] is the number of type A₃ molecules arriving at node Rx_p in time slot (j+1), which is represented by

\(\begin{aligned}N_{\mathrm{Rx}_{a}}^{A_{3}}[j+1]=\sum_{i=2}^{j+1} N_{A_{3}}[i] W_{R}^{i} h\left(i T_{s},(j+1-i) T_{s}+\tau_{s}\right),\end{aligned}\)       (9)

where N_A3[i] is the number of emitted type A₃ molecules at relay R to represent the transmission of bit 1. W_Rⁱ is the XOR result of Ŵ_Tx1ⁱ and Ŵ_Tx2ⁱ . Here, Ŵ_Tx1ⁱ and Ŵ_Tx2ⁱ are the corresponding decoded information of W_Tx1ⁱ and W_Tx2ⁱ , respectively. Based on the decision threshold at node R corresponding to received type A_p molecules which is θ_(Txp,R)^A_p[j](p ∈ {1, 2} ), we set up the signal detection rule at relay R as follows:

\(\begin{aligned}\hat{W}_{\mathrm{Tx}_{p}}^{j}=\left\{\begin{array}{l}1, \text { if } N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j]>\theta_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j], \\ 0, \text { otherwise },\end{array}\right.\end{aligned}\)       (10)

where θ_(Txp,R)^A_p[j] (p ∈ {1, 2} ) is used to decode bit W_Txp^j and Ŵ_Txp^j is the decoded information of bit W_Txp^j .

Let N_(Txp,R)^A_p[j] represent the number of molecules arriving at node R. It follows Normal distributions. Then the hypothesis testing corresponding to the received signal at relay node R is given by

H₀ : N_(Txp,R)^A_p[j] ~ N(μ_R,0^A_p, (σ_R,0^A_p)²),

H₁ : N_(Txp,R)^A_p[j] ~ N(μ_R,1^A_p, (σ_R,1^A_p)²),      (11)

where μ_R,0^A_p, (σ_R,0^A_p)² ,μ_R,1^A_pand (σ_R,1^A_p)² are obtained by

\(\begin{aligned}\begin{array}{l}\mu_{R, 0}^{A_{p}}=\sum_{i=1}^{j-1} N_{A_{p}}[i] \beta_{\mathrm{Tx}_{p}}^{i} \mu_{h}\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)+\mu_{\left(\mathrm{Tx}_{p}, R\right)}^{\text {Noise }}, \\ \left(\sigma_{R, 0}^{A_{p}}\right)^{2}=\sum_{i=1}^{j-1} N_{A_{p}}[i] \beta_{\mathrm{Tx}_{p}}^{i} \mu_{h}\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)\left(1-\mu_{h}\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)\right) \\ +N_{A_{p}}^{2}[i] \beta_{\mathrm{Tx}_{p}}^{i}\left(1-\beta_{\mathrm{Tx}_{p}}^{i}\right) \mu_{h}^{2}\left(m T_{s},(i-m) T_{s}+\tau_{s}\right)+\left(\sigma_{\left(\mathrm{Tx}_{p}, R\right)}^{\text {Noise }}\right)^{2}, \\ \mu_{R, 1}^{A_{p}}=\sum_{i=1}^{j} N_{A_{a}}[i] \beta_{\mathrm{Tx}_{p}}^{i} \mu_{h}\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)+\mu_{\left(\mathrm{Tx}_{p}, R\right)}^{\text {Noise }} \text {, } \\ \left(\sigma_{R, 1}^{A_{p}}\right)^{2}=\sum_{i=1}^{j} N_{A_{p}}[i] \beta_{\mathrm{Tx}_{p}}^{i} \mu_{h}\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)\left(1-\mu_{h}\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)\right) \\ +N_{A_{p}}^{2}[i] \beta_{\mathrm{Tx}_{p}}^{i}\left(1-\beta_{\mathrm{Tx}_{p}}^{i}\right) \mu_{h}^{2}\left(m T_{s},(i-m) T_{s}+\tau_{s}\right)+\left(\sigma_{\left(\mathrm{T}_{p}, R\right)}^{\text {Noise }}\right)^{2}, \\\end{array}\end{aligned}\)      (12)

where β_Txpⁱ represents the probability that Tx_p transmits bit 1 in time slot i. We use the maximum a posteriori (MAP) method to decrease the BEP. Λ(N_(Txp,R)^A_p[j]) is the likelihood-ratio and can be obtained by

\(\begin{aligned}\Lambda\left(N_{\left(\operatorname{Tx}_{p}, R\right)}^{A_{p}}[j]\right)=\frac{P\left(N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j] \mid H_{1}\right)}{P\left(N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j] \mid H_{0}\right)}=\frac{f^{1}\left(N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j]\right)}{f^{0}\left(N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j]\right)} \underset{H_{0}}{\gtrless} \frac{1-\beta_{\mathrm{Tx}_{p}}[j]}{\beta_{\mathrm{Tx}_{p}}[j]}\end{aligned}\)       (13)

where P(H₁) = β_Txp[j] and P(H₀) = 1-β_Txp[j] represent the probability that node Tx₁ transmits bit 1 and bit 0 in time slot j, respectively. f⁰(N_(Txp,R)^A_p[j]) and f¹ (N_(Txp,R)^A_p[j]) are the probability density function of (N_(Txp,R)^A_p[j] under H₀ and H₁, respectively, which are given by

\(\begin{aligned}\begin{array}{l}f^{0}\left(N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j]\right)=\frac{1}{\sqrt{2 \pi\left(\sigma_{R, 0}^{A_{p}}\right)^{2}}} \times e^{-\frac{\left(N_{\left(\mathrm{T}_{p}, R\right)}^{A_{p}}[j]-\mu_{R, 0}^{A_{p}}\right)^{2}}{2\left(\sigma_{R, 0}^{A_{p}}\right)^{2}}}, \\ f^{1}\left(N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j]\right)=\frac{1}{\sqrt{2 \pi\left(\sigma_{R, 1}^{A_{p}}\right)^{2}}} \times e^{-\frac{\left(N_{\left(\mathrm{T}_{p}, R\right)}^{A_{p}}[j]-\mu_{R, 1}\right)^{A_{p}}}{2\left(\sigma_{R, 1}^{A_{p}}\right)^{2}}} .\end{array}\end{aligned}\)       (14)

Therefore, θ_(Txp,R)^A_p[j]) j at node R for the link Tx_p → R is computed as follows:

\(\begin{aligned}N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j] \underset{H_{0}}{\gtrless} \operatorname{round}\left(\frac{v+\sqrt{v^{2}-u q}}{u}\right) \equiv \theta_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j],\end{aligned}\)      (15)

where the parameters u, v and q in (15) can be computed by

\(\begin{aligned}\begin{array}{l}u=\left(\sigma_{R, 1}^{A_{p}}\right)^{2}-\left(\sigma_{R, 0}^{A_{p}}\right)^{2}, \\ v=\mu_{R, 0}^{A_{p}}\left(\sigma_{R, 1}^{A_{p}}\right)^{2}-\mu_{R, 1}^{A_{p}}\left(\sigma_{R, 0}^{A_{p}}\right)^{2}, \\ q=\left(\mu_{R, 0}^{A_{p}}\right)^{2}\left(\sigma_{R, 1}^{A_{p}}\right)^{2}-\left(\mu_{R, 1}^{A_{p}}\right)^{2}\left(\sigma_{R, 0}^{A_{p}}\right)^{2}-2\left(\sigma_{R, 0}^{A_{p}}\right)^{2}\left(\sigma_{R, 1}^{A_{p}}\right)^{2}\left(\ln \frac{1-\beta_{\mathrm{Tx}_{p}}[j]}{\beta_{\mathrm{Tx}_{p}}[j]}-\ln \frac{\sigma_{R, 0}^{A_{p}}}{\sigma_{R, 1}^{A_{p}}}\right) .\end{array}\end{aligned}\)       (16)

When W_Tx1 ^j and W_Tx2^j are known, the probability of error occurred at node R is defined by P(E_R|W_Tx1^j, W_Tx2^j). Then we have

P(E_R|W_Tx1^j, W_Tx2^j) = P(E_R^A₁|W_Tx1^j)(1-P(E_R^A₂|W_Tx2^j)) + (1-P(E_R^A₁|W_Tx1^j))P(E_R^A₂|W_Tx2^j),       (17)

where P(E_R^A_p|W_Txp^j) represents the error probability that W_Txp^j is transmitted by node Tx_p in time slot j. Therefore, we can get

\(\begin{aligned}\begin{array}{l}P\left(E_{R}^{A_{p}} \mid W_{\mathrm{Tx}_{p}}^{j}=0\right)=P\left\{N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j]>\theta_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j] \mid W_{\mathrm{Tx}_{1}}^{j}=0\right\}=Q\left(\frac{\theta_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j]-\mu_{R, 0}^{A_{p}}}{\sigma_{R, 0}^{A_{p}}}\right), \\ P\left(E_{R}^{A_{p}} \mid W_{\mathrm{Tx}_{p}}^{j}=1\right)=P\left\{N_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j] \leq \theta_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{p}}[j] \mid W_{\mathrm{Tx}_{p}}^{j}=1\right\}=1-Q\left(\frac{\theta_{\left(\mathrm{Tx}_{p}, R\right)}^{A_{1}}[j]-\mu_{R, 1}^{A_{p}}}{\sigma_{R, 1}^{A_{p}}}\right) .\end{array}\end{aligned}\)       (18)

According to (17) and (18), the error probabilities at node R in (17) when W_Tx1^j and W_Tx2^j are known can be derived by Appendix A.

Let 3 θ_Rxp^A₃ be the decision threshold at node Rx_p (p∈ {1, 2}) corresponding to the number of molecules with type A₃. The signal detection at node Rx_p (p∈ {1, 2}) is obtained by

\(\begin{aligned}\hat{W}_{\mathrm{Rx}_{p}}^{j+1}=\left\{\begin{array}{l}1, \text { if } N_{\mathrm{Rx}_{p}}^{A_{3}}[j+1]>\theta_{\mathrm{Rx}_{p}}^{A_{3}}, \\ 0, \text { otherwise }\end{array}\right.\end{aligned}\)       (19)

where θ_Rxp^A₃ is used to get the decoded information ŷ_Rxp^j+1 in time slot (j+1) at node Rx_p which is transmitted to node Rx_p (p ∈ {1, 2}) . The value of ŷ_Rxp^j+1 is 0 or 1, then we can establish a binary hypothesis test model at node Rx_p (p ∈ {1, 2}) as follows:

H₀ : N_Rxp^A₃[j+1] ~ N(µ_Rxp,0^A₃, (σ_Rxp,0^A₃)²),

H₁ : N_Rxp^A₃[j+1] ~ N(µ_Rxp,1^A₃, (σ_Rxp,1^A₃)²),       (20)

where µ_Rxp,0^A₃ ,(σ_Rxp,0^A₃)², µ_Rxp,1^A₃ and (σ_Rxp,1^A₃)² are the corresponding mean and variance of the Normal distributions which N_Rxp^A₃[j+1] follows under H₀ and H₁, respectively, which are obtained by

\(\begin{aligned}\begin{array}{l}\mu_{\mathrm{Rx}_{p}, 0}^{A_{3}}=\sum_{i=2}^{j} N_{A_{3}}[i] \beta_{R}^{i} \mu_{h}\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)+\mu_{\left(R, \mathrm{Rx}_{p}\right)}^{\text {Noise }}, \\ \left(\sigma_{\mathrm{Rx}_{p}, 0}^{A_{3}}\right)^{2}=\sum_{i=2}^{j} N_{A_{3}}[i] \beta_{R}^{i}\left(\sigma_{h}\left(i T_{s},(j-i) T_{s}+\tau_{s}\right)\right)^{2}+\left(\sigma_{\left(R, \mathrm{Rx}_{p}\right)}^{\text {Noise }}\right)^{2}, \\ \mu_{\mathrm{Rx}_{p}, 1}^{A_{3}}=\sum_{i=2}^{j+1} N_{A_{3}}[i] \beta_{R}^{i} \mu_{h}\left(i T_{s},(j+1-i) T_{s}+\tau_{s}\right)+\mu_{\left(R, \mathrm{Rx}_{p}\right)}^{\text {Noise }}, \\ \left(\sigma_{\mathrm{Rx}_{p}, 1}^{A_{3}}\right)^{2}=\sum_{i=2}^{j+1} N_{A_{3}}[i] \beta_{R}^{i}\left(\sigma_{h}\left(i T_{s},(j+1-i) T_{s}+\tau_{s}\right)\right)^{2}+\left(\sigma_{\left(R, \mathrm{Rx}_{p}\right)}^{\text {Noose }}\right)^{2},\end{array}\end{aligned}\)       (21)

where β_Rⁱ is the probability that node R transmits bit 1 in time slot i. Let P(E_Rxp^A₃ |W_R^j = 0) and P(E_Rxp^A₃|W_R^j = 1) represent error probability at node Rx_p (p ∈ {1, 2}) corresponding to the molecules with type A₃ when W_R^j = 0 and W_R^j = 1, respectively. Then we have

\(\begin{aligned}\begin{array}{l}P\left(E_{\mathrm{Rx}_{p}} \mid W_{R}^{j+1}=0\right)=P\left\{N_{\mathrm{Rx}_{p}}^{A_{3}}[j+1]>\theta_{\mathrm{Rx}_{p}}^{A_{3}}[j+1] \mid W_{R}^{j+1}=0\right\}=Q\left(\frac{\theta_{\mathrm{Rx}_{p}}^{A_{3}}[j+1]-\mu_{\mathrm{Rx}_{p}, 0}^{A_{3}}}{\sigma_{\mathrm{Rx}_{p}, 0}^{A_{3}}}\right), \\ P\left(E_{\mathrm{Rx}_{p}} \mid W_{R}^{j+1}=1\right)=P\left\{N_{\mathrm{Rx}_{p}}^{A_{3}}[j+1] \leq \theta_{\mathrm{Rx}_{p}}^{A_{3}}[j+1] \mid W_{R}^{j+1}=1\right\}=1-Q\left(\frac{\theta_{\mathrm{Rx}_{p}}^{A_{3}}[j+1]-\mu_{\mathrm{Rx}_{p}, 1}^{A_{3}}}{\sigma_{\mathrm{Rx}_{p}, 1}^{A_{3}}}\right) .\end{array}\end{aligned}\)       (22)

The relay node R decodes information W_Tx1^j and W_Tx2^j with error when W_R^j ≠ Ŵ_tx1^j ⊕ Ŵ_tx2^j. Under this case, there is one decoding error of W_Tx1^j or W_Tx2^j , and the other is decoded correctly. When W_Tx1^j and W_Tx2^j are known, the error probability at node Rx_p (p ∈ {1, 2}) is denoted by P(E_Rxp|W_Tx1^j,W_Tx2^j ). Four combinations of values of W_Tx1^j and W_Tx2^j are generated whenW_Tx1^j ,W_Tx2^j ∈ {0,1}. Then the results of Pe_(Tx1,Rx1)[j + 1] can be obtained by Appendix B. Similarly, the error probability Pe_(Tx2,Rx2)[j + 1] can be computed. The average BEP of transmitting one bit for mobile MIMO MCvD system with one relay is formulated by

\(\begin{aligned}P e_{A v g}[j+1]=\frac{1}{2} \sum_{p=1}^{2} P e_{\left(\mathrm{Tx}_{p}, \mathrm{Rx}_{p}\right)}[j+1]\end{aligned}\)       (23)

4.2 Optimizing the ratio of the number of emitted molecules at two source nodes and the initial position of relay node

The goal of the optimization problem is to minimize the average BEP of transmitting one bit in mobile MIMO MCvD system with one relay and the optimization variables are the ratio of the number of emitted molecules at two source nodes and the initial position of relay node. Then this optimization problem is formulated as follows:

\(\begin{aligned}\min _{m, x_{R}^{0}, y_{R}^{0}, z_{R}^{0}} P e_{A v g}[j+1]=\min _{m, x_{R}^{0}, y_{R}^{0}, z_{R}^{0}} \frac{1}{2} \sum_{p=1}^{2} P e_{\left(\mathrm{Tx}_{p}, \mathrm{Rx}_{p}\right)}[j+1]\end{aligned}\)       (24)

where m = N_A1 / (N_A1 + N_A2) . x_R⁰, y_R⁰ and z_R⁰ are the coordinates of initial position of node R from three different axes, respectively. Then the initial position of relay node is computed by the optimal values of x_R⁰, y_R⁰ and z_R⁰ .

Algorithm 1: Bisection method

Choose ε which satisfies 0 < ε < 1. Then we set δ = 0, γ = 1. In addition, the optimal value of m which is represented by m*, can satisfy the relationships δ ≤ Pe_Avg[j + 1]|_m=m* and γ ≥ Pe_Avg[j + 1]|_m=m*.

Iterations:

Step 1, α = (δ + γ)/2 .

Step 2, when x_R⁰, y_R⁰ and z_R⁰ are fixed, check whether the problem in (24) is convex feasibility problem.

Step 3, if the convex problem in (24) is feasible, we can update δ = α ; else γ = α .

Until ||δ − γ|| ≤ ε .

Algorithm 2: The iterative algorithm for joint optimizations of m, x_R⁰, y_R⁰ and z_R⁰

The initialization is expressed as follows: Set 0 < ε < 1 and iteration number k=0, x_R⁰ =0.15μm, y_R⁰ =0.15μm and z_R⁰ =0. Obtain m⁰ on the basis of bisection method in Algorithm 1.

The process of iterations is composed of three steps which is described as follows:

Step 1, when m^k is fixed, search (x_R⁰)^k+1 according to the bisection method in Algorithm 1, and then find (y_R⁰)^k+1 and (z_R⁰)^k+1 by Algorithm 1, respectively.

Step 2, for fixed (x_R⁰)^k+1, (y_R⁰)^k+1 and (z_R⁰)^k+1 , find m^k+1 by using the bisection method shown in Algorithm 1.

Step 3, set k=k+1.

Until ||m^k - m^k-1||≤ ε, ||(x_R⁰)^k - (x_R⁰)^k-1||≤ ε , ||(y_R⁰)^k - (y_R⁰)^k-1||≤ ε and ||(z_R⁰)^k - (z_R⁰)^k-1||≤ ε .

5. Numerical Results

The numerical results are used to evaluate the average BEP of MIMO MCvD system with one relay based on DF and NC scheme by optimizing the ratio of the number of emitted molecules at node Tx₁ and Tx₂ and the initial distance between Tx₁ and node R. The numerical parameters for the evaluations are given in Table 1. The coordinates of the initial position of four nodes are set in Table 2 as follows:

Table 1. The numerical simulation parameters

Table 2. The initial positions of nodes Tx_p and Rx_p (p=1, 2)

In Fig. 2, it is obvious that the average BEP is a function of m and the initial distance between node R and Tx₁. We can see that the average BEP with respect to m and the initial distance between node R and Tx₁ is convex since it has minimum value [25]. The value of m is N_A1 / (N_A1 + N_A2) which represents the ratio of the number of emitted molecules at node Tx₁. The initial distance between Tx₁ and R can be computed by the positions of node Tx₁ and R. In such a case, the local optimum is the global optimum. Then we use the iteration algorithm to find the optimal value of m and the coordinates of the initial position of node R from three different axes x_R⁰ , y_R⁰ and z_R⁰ . Then the initial position of node R can be obtained. The parameters in Fig. 2 are set as N_A1 + N_A2 = 3 × 10⁴ , N_A3 = 2 × 10⁴. D_Txp = D_Rxp = D_R = 10^-13 m² /s (p=1, 2).

Fig. 2. The average BEP is varying with m and the initial distance between Tx₁ and R.

In Fig. 3, the convergence result of the average BEP with respect to number of iterations by using the iteration algorithm is shown. It is obvious that the average BEP is decreasing with the increase of the number of iterations, and finally it converges a stable value. The result shows that the error probability performance exhibits fast convergence under different values of ratios between N_A1 and N_A2 in Fig. 3(a) and the signal-to-noise ratio (SNR) in Fig. 3(b). However, the convergence speeds of average BEP by implementing the iteration algorithm with different numbers of molecules released by Tx₁ and Tx₂ are different. We set ε = 0.0001 for the iterative algorithm for different values of N_A1 and N_A2 . In Fig. 3(a), for the cases N_A1 = N_A2 , N_A1 =2 N_A2 and N_A1 =4 N_A2 , the average BEP requires 15, 13 and 10 iterations to achieve convergence, respectively. In addition, when N_A1 = N_A2 , the average BEP has the minimum value. In Fig. 3(b), the larger value of SNR, the average BEP decreases faster and can converge at the minimum value. The parameters in Fig. 3 are set as N_{A1 +} N_A2 = 3 ×10⁴, N_A3 = 2 ×10⁴.

Fig. 3. Convergence result of average BEP vs number of iterations with different values of (a) ratios between N_A1 and N_A2 ; (b) SNR.

Fig. 4 explores the variation of average BEP versus the total number of emitted molecules at node Tx₁ and node Tx₂ which is N_A1 + N_A2 with different values of noise variance of each link. It can be concluded that the average BEP decreases with an increase of N_A1 + N_A2 . After some value of N_A1 + N_A2, the average BEP reaches its minimum value. The average BEP declines slowly when the noise variance is larger. Therefore, we can see that the average BEP with (σ^Noise)²=500 reaches its minimum value more slowly than that with (σ^Noise)²=1500. When the value of N_A1 + N_A2 is fixed, the larger value of (σ^Noise)² , the larger value of the average BEP. The parameters in Fig. 4 are set as N_A3 = 2 × 10⁴.

Fig. 4. The average BEP is varying with the total number of released molecules at Tx₁ and Tx₂ with different value of variance of noise.

Fig. 5 shows that the average BEP decreases with the total number of emitted molecules at sources nodes Tx₁ and Tx₂. In addition, T_s has an influence on the average BEP. In Fig. 5, the average BEP is decreasing with . T_s When other parameters are set the same, we observe that the average BEP with T_s=0.2s reduces the fastest than those with T_s=0.1s and T_s=0.01s. It can be explained by the following results: the receiving probability of one molecule for the transmission of each link will increase with T_s. Under this case, more molecules are received by the receiver node with a larger value of T_s which leads to a decrease of the average BEP.

Fig. 5. The average BEP vs the total number of released molecules at Tx₁ and Tx₂ with different values of T_s.

N_A3 is the number of molecules released by R. Fig. 6 plots the average BEP decreases with N_A3 under different values of N_A1 + N_A2. When the other parameters are fixed, the receiving probability of one molecule for the transmission links R → Rx₁ and R → Rx₂ will increase with N_A3 . Then the average BEP will decrease. For the same value of N_A3 , the average BEP is decreasing with the values of N_A1 + N_A2, which is consistent with the results in Fig. 5. For each value of N_A1 + N_A2 , the values of optimal m are around 0.5 and the optimal initial distance between Tx₁ and R under different values of N_A3 are different.

Fig. 6. The average BEP vs number of molecules released by R with different values of N_A1 + N_A2 .

In Fig. 7, we consider three different cases with N_A1 + N_A2 =2 × 10⁴, N_A1 + N_A2 = 3 × 10⁴ and N_A1 + N_A2 =4 × 10⁴ to evaluate the average BEP of the mobile MIMO MCvD system with one relay. According to Fig. 7, the average BEP decreases first with thresholds at Rx₁ and Rx₂ and then arrives at its minimum value. When the thresholds at Rx₁ and Rx₂ reach some value, the average BEP begins to increase and then reaches the maximum value. Especially, the minimum value of average BEP under the cases N_A1 + N_A2 = 2 × 10⁴ and N_A1 + N_A2 = 3 × 10⁴ are smaller than that under the case N_A1 + N_A2 = 4 × 10⁴ . Furthermore, we also find that when the value of N_A1 + N_A2 is larger, the optimal decision threshold at relay node R is also larger. This result is based on the facts that the decision threshold at relay node R is directly related to the value of N_A1 or N_A2 . Similarly, when the value of N_A3 is larger, the optimal decision threshold at node Rx₁ or node Rx₂ is also larger. Furthermore, when N_A3 = 2 × 10⁴ , the optimal value of N_A1 / (N_A1 + N_A2) is 0.5, 0.5 and 0.47 under three cases when N_A1 + N_A2 = {2, 3, 4} × 10⁴ , respectively. The optimal initial distance between Tx₁ and R under the corresponding three cases is 9.89μm, 10.29μm and 10.50μm, respectively.

Fig. 7. The average BEP is varying with decision thresholds at Rx₁ and Rx₂ with different values of total number of molecules released by Tx₁ and Tx₂.

The result in Fig. 8 reveals that the average BEP is varying with m under different optimal detection threshold schemes between MAP and fixed detection threshold scheme. Under the four cases, with the increasing of m, it is obviously that the average BEP decreases first and reaches its minimum value and then increases. When the value of m continues to increase, the average BEP gets its maximum value. It is easily to obtain that the average BEP is smaller under MAP scheme than the case under fixed detection thresholds which are 100 and 120. Therefore, we use MAP scheme to obtain the optimal detection threshold at each node. The minimum value of m is 0.5 under the three cases. In addition, when the variance between the detection threshold obtained by MAP scheme and the value of fixed detection threshold is larger, the average BEP is also larger. The parameters in Fig. 8 are set as N_A1 + N_A2 = 3 × 10⁴, N_A3 = 2 × 10⁴, d_(Tx1,R)⁰ = 10μm.

Fig. 8. The average BEP vs m with different detection threshold schemes.​​​​​​​

6. Conclusions

In this work, the mobile MIMO MCvD system with one relay node by using DF relaying protocol and NC scheme in 3D environment is investigated. In particular, we minimize the average BEP of this system by using the proposed iteration algorithm based on BCDA. First, we derive the mathematical expressions of the average BEP of this system. Second, we establish the optimization problem whose objective is to achieve the minimum value of average BEP. Third, four optimization variables which are the number of emitted molecules at nodes Tx₁ and Tx₂ and the coordinates of relay node R from three axes can be jointly optimized by using the proposed iteration algorithm. Furthermore, the optimization variables can be simultaneously obtained by using the iteration algorithm to obtain the minimum average BEP when the other system parameters are determined.

According to the numerical results, we demonstrate that the iteration algorithm based on BCDA has good convergence behaviors with different values of m and SNR. By using the proposed iteration algorithm, the optimal ratio of the number of emitted molecules at two source nodes and the optimal initial position of relay node can be obtained. Then the number of emitted molecules at two source nodes can be determined. Moreover, it is observed that the impacts of some parameters on the average BEP of this system are different. Specifically, increasing the number of emitted molecules by source nodes or by relay node in each time slot and the size of each time slot can improve the performance of this system. MAP scheme is superior to fixed detection threshold scheme with different values of m. We plan to explore the optimization problems of MIMO MCvD system with multiple relays to minimize the average BEP of this system.

Appendix A

Appendix B