Diffusive Mobile MC with Absorbing Receivers: Stochastic Analysis and Applications

TL;DR

This paper provides a comprehensive stochastic analysis of a diffusive mobile molecular communication system with absorbing receivers, deriving key statistical properties of the channel and applying them to optimize drug delivery and communication system performance.

Contribution

It introduces the first analytical derivation of the CIR's distribution and absorption probability in a mobile MC system, enabling optimized system designs.

Findings

Significant reduction in drug release compared to constant schemes.

Improved bit error rate and molecule efficiency with optimized parameters.

Demonstrated importance of accounting for mobility in system design.

Abstract

This paper presents a stochastic analysis of the time-variant channel impulse response (CIR) of a three dimensional diffusive mobile molecular communication (MC) system where the transmitter, the absorbing receiver, and the molecules can freely diffuse. In our analysis, we derive the mean, variance, probability density function (PDF), and cumulative distribution function (CDF) of the CIR. We also derive the PDF and CDF of the probability p that a released molecule is absorbed at the receiver during a given time period. The obtained analytical results are employed for the design of drug delivery and MC systems with imperfect channel state information. For the first application, we exploit the mean and variance of the CIR to optimize a controlled-release drug delivery system employing a mobile drug carrier. We evaluate the performance of the proposed release design based on the PDF and CDF of the CIR. We demonstrate significant savings in the amount of released drugs compared to a constant-release scheme and reveal the necessity of accounting for the drug-carrier's mobility to ensure reliable drug delivery. For the second application, we exploit the PDF of the distance between the mobile transceivers and the CDF of p to optimize three design parameters of an MC system employing on-off keying modulation and threshold detection. Specifically, we optimize the detection threshold at the receiver, the release profile at the transmitter, and the time duration of a bit frame. We show that the proposed optimal designs can significantly improve the system performance in terms of the bit error rate and the efficiency of molecule usage.

Tables

Table 1. TABLE I: System parameters used for numerical results

Parameter	Value	Parameter	Value
$D_{\mathrm{Tx}}$ [$m^2/s$]	$1\times {10}^{-14}$	$D_{\mathrm{Rx}}$ [$m^2/s$]	$0$
$D_X$ [$m^2/s$]	$8\times {10}^{-11}$	$r_0$ [$m$]	$10\times {10}^{-6}$
$a_{\mathrm{tx}}$ [$m$]	$1\times {10}^{-7}$	$a_{\mathrm{rx}}$ [$m$]	$1\times {10}^{-6}$
$T$ [$h$]	$24$	$T_{\mathrm{Rx}}$ [$h$]	$24$
$I$	$3000$	$N$	$5$
$T_b$[$s$]	$28.8$	$\theta \left(t\right)\left[s^{-1}\right]$	$1$

Equations

$$h(t,\tau )=\frac{a_{\mathrm{rx}}}{\sqrt{4\pi D_1{\tau }^3}}\left(1-\frac{a_{\mathrm{rx}}}{r(t)}\right)\exp \left(-\frac{{\left(r(t)-a_{\mathrm{rx}}\right)}^2}{4D_1\tau }\right),\tau >0,$$

$$g\left(t\right)=\sum _{\forall i|t_i<t}{\alpha }_{i}h\left(t_i,t-t_i\right).$$

$$q_i\sim \mathcal{N}\left({\mu }_{i,{\tilde{b}}_i},{\sigma }_{i,{\tilde{b}}_i}^2\right)\text{ for }{b}_i={\tilde{b}}_i,$$

$$p(t,T_{\mathrm{b}})={\int }_0^{T_{\mathrm{b}}}h(t,\tau )\mathrm{d}\tau =\frac{a_{\mathrm{rx}}}{r(t)}\mathrm{erfc}\left(\frac{r(t)-a_{\mathrm{rx}}}{2\sqrt{D_1{T}_{\mathrm{b}}}}\right),$$

$$\gamma =\frac{r(t)}{\sqrt{2D_{2}t}}=\sqrt{\frac{{\sum }_{d\in \{x,y,z\}}{r}_d^2(t)}{2D_{2}t}},$$

$$\mathrm{E}\left\{r(t)\right\}$$

$$\mathrm{Var}\left\{r(t)\right\}$$

$$f_{r(t)}(r)$$

$$\text{and}{F}_{r(t)}(r)$$

$$m(t,\tau )={\int }_0^{\infty }h(t,\tau )|{}_{r(t)=r}{f}_{r(t)}(r)\mathrm{d}r.$$

$$m(t,\tau )$$

$$\times \mathrm{erfc}\left(\frac{v(t,\tau )}{2\sqrt{u(t,\tau )}}\right)+{\mathrm{e}}^{\frac{w(t,\tau {)}^2}{4u(t,\tau )}}\left(\frac{w(t,\tau )}{2u(t,\tau )}+a_{\mathrm{rx}}\right)\mathrm{erfc}\left(\frac{w(t,\tau )}{2\sqrt{u(t,\tau )}}\right)],$$

$$u(t,\tau )=\frac{1}{4D_1\tau }+\frac{1}{4D_{2}t},v(t,\tau )=-\frac{a_{\mathrm{rx}}}{2D_1\tau }-\frac{r_0}{2D_{2}t},w(t,\tau )=-\frac{a_{\mathrm{rx}}}{2D_1\tau }+\frac{r_0}{2D_{2}t}.$$

$${\sigma }^2(t,\tau )=\varphi (t,\tau )-m^2(t,\tau ),$$

$$\varphi (t,\tau )=$$

$$\times \left(r_1-2a_{\mathrm{rx}}+\frac{a_{\mathrm{rx}}^2}{r_1}\right)\mathrm{d}{r}_1,$$

$$c(t,\tau )=\frac{a_{\mathrm{rx}}^2{\mathrm{e}}^{-\frac{a_{\mathrm{rx}}^2}{2D_1\tau }-\frac{r_0^2}{4D_{2}t}}}{8D_1\pi {\tau }^3{r}_0\sqrt{\pi D_{2}t}},$$

$$\hat{v}(t,\tau )=-\frac{a_{\mathrm{rx}}}{D_1\tau }-\frac{r_0}{2D_{2}t},$$

$$\varphi (t,\tau )$$

$$\{\begin{array}{ll}{f}_{h\left(t,\tau \right)}(h)=\frac{f_{r(t)}(r_1(h))}{{\hat{h}}^{\prime }(r_1(h),\tau )}-\frac{f_{r(t)}(r_2(h))}{{\hat{h}}^{\prime }(r_2(h),\tau )},& \text{for }0\le h<h^{\star },\\ f_{h\left(t,\tau \right)}(h)\to \infty ,& \text{for }h=h^{\star },\\ f_{h\left(t,\tau \right)}(h)=0,& \text{otherwise,}\end{array}$$

$${\hat{h}}^{\prime }(r,\tau )=$$

$$\{\begin{array}{ll}{F}_{h\left(t,\tau \right)}(h)=F_{r(t)}(r_1(h))+1-F_{r(t)}(r_2(h)),& \text{for }0\le h\le h^{\star },\\ F_{h\left(t,\tau \right)}(h)=0,& \text{for }h<0,\\ F_{h\left(t,\tau \right)}(h)=1,& \text{for }h>h^{\star },\end{array}$$

$$F_{h\left(t,\tau \right)}(h)=$$

$$=$$

$$f_{p(t,T_{\mathrm{b}})}(p)=-\frac{f_{r(t)}(\tilde{r}(p))}{p^{\prime }(\tilde{r})},$$

$$F_{p(t,T_{\mathrm{b}})}(p)=1-F_{r(t)}(\tilde{r}(p)),$$

$$p^{\prime }(\tilde{r})=$$

$$\mathrm{E}\left\{g(t)\right\}-\beta \mathrm{V}\left\{g(t)\right\}\stackrel{(a)}{\ge }\sum _{\forall i|t_i<t}{\alpha }_i\left(\mathrm{E}\left\{h\left(t_i,t-t_i\right)\right\}-\beta \mathrm{V}\left\{h\left(t_i,t-t_i\right)\right\}\right)\ge \theta (t),$$

$$\mathrm{V}\left\{g(t)\right\}$$

$$\le \sum _{\forall i|t_i<t}{\alpha }_i{\left[\mathrm{E}\left\{{\left(h\left(t_i,t-t_i\right)-\mathrm{E}\left\{h\left(t_i,t-t_i\right)\right\}\right)}^2\right\}\right]}^{1/2}=\sum _{\forall i|t_i<t}{\alpha }_i\mathrm{V}\left\{h\left(t_i,t-t_i\right)\right\}.$$

$$\underset{{\alpha }_i\ge 0,\forall i}{\min }$$

$$\sum _{i,t_i<t}{\alpha }_i\left(m\left(t_i,n\Delta t_n-t_i\right)-\beta \sigma \left(t_i,n\Delta t_n-t_i\right)\right)\ge \theta (n\Delta t_n),\text{ for }n=1,\dots ,N,$$

$$P_{\theta }(t)=$$

$$P_{\theta }(t)=1-F_{g(t)}\left\{\theta (t)\right\}=1-{\int }_0^{\theta (t)}{f}_{g(t)}(g)\mathrm{d}g.$$

$$f_{g(t)}(g)=\left(f_{{\alpha }_{1}h\left(t_1,t-t_1\right)}\ast \cdots \ast f_{{\alpha }_{i}h\left(t_i,t-t_i\right)}\right)(g).$$

Full text

I Introduction

As appropriate channel models are essential for the analysis and design of molecular communication (MC) systems, MC channel modeling has been extensively studied in the literature, see [1] and references therein. For example, the simple diffusive channel model of an unbounded three-dimensional (3D) MC system with impulsive point release of information carrying molecules [2] has been widely used for system analysis and design, see [3, 4], and references therein. Diffusion channel models with drift [5] and chemical reactions [6] have also been considered. However, most of the previously studied MC channel models assume static communication systems where the transceivers do not move.

Recently, many applications have emerged where the transceivers are mobile, including drug delivery [7], mobile ad hoc networks [8], and detection of mobile targets [9]. Hence, the modeling and design of mobile MC systems have gained considerable attention, e.g., see [8, 9, 10, 11, 12, 1, 13, 14, 15], and references therein. In [8], a mobile ad hoc nanonetwork was considered where mobile nanomachines collect environmental information and deliver it to a mobile central control unit. The mobility of the nanomachines was described by a 3D model but information was only exchanged when two nanomachines collided. In [9], a leader-follower-based model for two-dimensional mobile MC networks for target detection with non-diffusive information molecules was proposed. The authors in [10] considered adaptive detection and inter-symbol interference (ISI) mitigation in mobile MC systems, while [11] analyzed the mutual information and maximum achievable rate in such systems. However, the authors of [10] and [11] did not provide a stochastic analysis of the time-variant channel but analyzed the system numerically. In [12], a comprehensive framework for modeling the time-variant channels of diffusive mobile MC systems with diffusive transceivers was developed. However, all of the works mentioned above assumed a passive receiver.

On the other hand, for many MC applications, a fully absorbing receiver is considered to be a more realistic model compared to a passive receiver as it captures the interaction between the receiver and the information molecules, e.g., the conversion of the information molecules to a new type of molecule or the absorption and removal of the information molecules from the environment [3, 2]. Since the molecules are removed from the environment after being absorbed by the receiver, the channel impulse response (CIR) for absorbing receivers is a more complicated function of the distance between the transceivers and the receiver’s radius compared to passive receivers. Therefore, the stochastic analysis of mobile MC systems with absorbing receiver is very challenging. For the fully absorbing receiver in diffusive mobile MC systems, theoretical expressions for the average distribution of the first hitting time, i.e., the mean of the CIR, were derived for a one-dimensional (1D) environment without drift in [13] and with drift in [14]. Based on the 1D model in [13], the error rate and channel capacity of the system were examined in [15]. However, none of these works provides a statistical analysis of the time-variant CIR of a 3D diffusive mobile MC system with absorbing receiver. In this paper, we address this issue and exploit the obtained analytical results for the stochastic parameters of the time-variant MC channel for the design of drug delivery and MC systems.

In drug delivery systems, drug molecules are carried to diseased cell sites by nanoparticle drug carriers, so that the drug is delivered to the targeted site without affecting healthy cells [7]. After being injected or extravasated from the cardiovascular system into the tissue surrounding a targeted diseased cell site, the drug carriers may not be anchored at the targeted site but may move, mostly via diffusion [16, 17, 18, 19]. The diffusion of the drug carriers results in a time-variant absorption rate of the drugs even if the drug release rate is constant. Furthermore, experimental and theoretical studies have indicated that the total drug dosage as well as the rate and time period of drug absorption by the receptors of the diseased cells are critical factors in the healing process [18, 20]. Therefore, to satisfy reduce drug cost, over-dosing, and negative side effects to healthy cells yet satisfy the treatment requirements, it is important to optimize the release profile of drug delivery systems such that the total amount of released drugs is minimized while a desired rate of drug absorption at the diseased site during a prescribed time period is achieved. To this end, the mobility of the drug carriers and the absorption rate of the drugs have to be accurately taken into account. This can be accomplished by exploiting the MC paradigm where the drug carriers, diseased cells, and drug molecules are modeled as mobile transmitters, absorbing receivers, and signaling molecules, respectively [3]. Release profile designs for drug delivery systems based on an MC framework were proposed in [21, 4, 22, 23]. However, in these works, the transceivers were fixed and only the movement of the drug molecules was considered. In this paper, we exploit the analytical results obtained for the stochastic parameters of the time-variant MC channel with absorbing receiver for the optimization of the release profile of drug delivery systems with mobile drug carriers.

In diffusive mobile MC systems, knowledge of the CIR is needed for reliable communication design. However, the CIR may not always be available in a diffusive mobile MC system due to the random movements of the transceivers. In particular, the distance between the transceivers at the time of release, on which the CIR depends, may only be known at the start of a transmission frame. In other words, the movement of the transceivers causes the CSI to become outdated, which makes communication system design challenging. In this paper, we consider a mobile MC system employing on-off keying and threshold detection and optimize three design parameters to improve the system performance under imperfect CSI. First, we optimize the detection threshold at the receiver for minimization of the maximum bit error rate (BER) in a frame when the number of molecules available for transmission is uniformly allocated to each bit of the frame. Second, we optimize the release profile at the transmitter, i.e., the optimal number of molecules available for the transmission of each bit, for minimization of the maximum BER in a frame given a fixed number of molecules available for transmission of the entire frame. Third, we maximize the frame duration under the constraint that the probability that a released molecule is absorbed by the receiver does not fall below a prescribed value. Such a design ensures that molecules are used efficiently as a molecule release occurs only if the released molecule is observed at the receiver with sufficiently high probability. For the proposed design tasks, the results of the stochastic analysis of the transceivers’ positions and of the probability that a molecule is absorbed during a given time period are exploited.

In summary, the main contributions of this paper are as follows:

• •

We provide a statistical analysis of the time-variant channel of a 3D diffusive mobile MC system employing an absorbing receiver. In particular, we derive the mean, variance, PDF, and CDF of the corresponding CIR. Moreover, we derive the PDF and CDF of the probability that a molecule is absorbed during a given time period. The stochastic channel analysis is exploited for the design of drug delivery and MC systems.

• •

For drug delivery systems, the release profile is optimized for the minimization of the amount of released drugs while ensuring that the absorption rate at the diseased cells does not fall below a prescribed threshold for a given period of time. We show that the proposed design requires a significantly lower amount of released drugs compared to a design with constant-release rate.

• •

For MC systems employing on-off keying modulation and threshold detection based on imperfect CSI, we optimize three design parameters, namely the detection threshold at the receiver, the release profile at the transmitter, and the time duration of a bit frame. Simulation results show significant performance gains for the proposed designs in terms of BER and the efficiency of molecule usage compared to baseline systems with uniform molecule release and without limitation on time duration of a bit frame, respectively.

• •

Our results reveal that the transceivers’ mobility has a significant impact on the system performance and should be taken into account for MC system design.

We note that the derived analytical results for the time-variant CIR of mobile MC systems with absorbing receiver are expected to be useful not only for the design of the drug delivery and MC systems considered in this paper but also for the design of detection schemes and the evaluation of the performance (e.g., the capacity and throughput) of such systems.

This paper expands its conference version [24]. In particular, the analysis of the probability that a molecule is absorbed during a given time period, the MC system design for imperfect CSI, and the corresponding simulation results are not included in [24].

The remainder of this paper is organized as follows. In Section II, we introduce the considered diffusive mobile MC system with absorbing receiver and the time-variant channel model. In Section III, we provide the proposed statistical analysis of the time-variant channel. In Sections IV and V, we apply the derived results for optimization of drug delivery and MC systems with imperfect CSI, respectively. Numerical results are presented in Section VI, and Section VII concludes the paper.

II General System and Channel Model

In this section, we first introduce the model for a general diffusive mobile MC system with absorbing receiver. Subsequently, we specialize the model to drug delivery and communication systems with imperfect CSI. Finally, we define the time-variant CIR and the received signal.

II-A System Model

We consider a linear diffusive mobile MC system in an unbounded 3D environment with constant temperature and viscosity. The system comprises one mobile spherical transparent transmitter, denoted by $\mathrm{Tx}$, with radius $a_{\mathrm{tx}}$, one mobile spherical absorbing receiver, denoted by $\mathrm{Rx}$, with radius $a_{\mathrm{rx}}$, and the signaling molecules of type $\mathrm{X}$.

The movements of $\mathrm{Tx}$, $\mathrm{Rx}$, and $\mathrm{X}$ molecules are assumed to be mutually independent and follow Brownian motion with diffusion coefficients $D_{\mathrm{Tx}}$, $D_{\mathrm{Rx}}$, and $D_{\mathrm{X}}$, respectively. This assumption, which was also made in [12] and [13], is motivated by the fact that the mobility of small objects is governed by Brownian motion.

We assume that $\mathrm{Tx}$ releases molecules at its center instantaneously and discretely during the considered period of time denoted by $T$. Let $t_{i}$ and $T_{\mathrm{b}}$ denote the time instant of the $i$-th release and the duration of the interval between the $i$-th and the $(i+1)$-th release, respectively. We have $t_{i}=(i-1)T_{\mathrm{b}}$ and $i\in\{1,\dots,I\}$, where $I$ is the total number of releases during $T$. We denote the time-varying distance between the centers of $\mathrm{Tx}$ and $\mathrm{Rx}$ at time $t$ by $r(t)$. Furthermore, let $\alpha_{i}$ and $A=\sum_{i=1}^{I}\alpha_{i}$ denote the number of molecules released at time $t_{i}$ and the total number of molecules released during $T$, respectively. For concreteness, we specialize the considered general model to two application scenarios.

II-A1 Drug Delivery Systems

A drug delivery system comprises a drug carrier releasing drug molecules and diseased cells absorbing them. We model the drug carrier and diseased cells as $\mathrm{Tx}$ and $\mathrm{Rx}$ of the general MC system, respectively, see Fig. 1. The drug carriers in drug delivery systems are typically nanoparticles, such as spherical polymers or polymer chains, having a size not smaller than $100\text{\,}\mathrm{nm}$ [17]. Moreover, drug carriers are designed to carry drug molecules and interaction with the drug or the receiver is not intended. Hence, the drug carriers can be modeled as mobile spherical transparent transmitters, $\mathrm{Tx}$. When the drug molecules hit the tumor, they are absorbed by receptors on the surface of the diseased cells [18, 20]. For convenience, we model the tumor as a spherical absorbing receiver, $\mathrm{Rx}$. In reality, the colony of cancer cells may potentially have a different geometry, of course. However, as an abstract approximation, we model the cancer cells as one effective spherical receiver with radius $a_{\mathrm{rx}}$ and with a surface area equivalent to the total surface area of the tumor (see Fig. 1). Hence, the absorption on the actual and the modeled surfaces is expected to be comparable [16].

In a drug delivery system, the drug carriers can be directly injected or extravasated from the blood into the interstitial tissue near the diseased cells, where they start to move. We assume that the injection position can be estimated and thus $r(t=0)$ is known. The movement of the drug carrying nanoparticles in the tissue is caused by diffusion and convection mechanisms but diffusion is expected to be dominant in most cases [16, 17, 18, 19]. At the tumor site, the drug carrier releases drug molecules of type $\mathrm{X}$, which also diffuse in the tissue [18]. Hence, we can adopt Brownian motion to model the diffusion of $\mathrm{Tx}$ and $\mathrm{X}$ molecules with diffusion coefficients $D_{\mathrm{Tx}}$, and $D_{\mathrm{X}}$, respectively [7]. We consider a rooted tumor and thus $D_{\mathrm{Rx}}=0$, which is a special case of the considered general system model.

We assume the instantaneous and discrete release of drugs. After releasing for a period, the drug carrier may be removed by blood circulation or run out of drugs. Thus, for drug delivery systems, $T$, $T_{\mathrm{b}}$, $t_{i}$, and $I$ denote the release period of the drug, the duration of the interval between two releases, the release instants of the drug molecules, and the number of releases, respectively. A continuous release can be approximated by letting $T_{\mathrm{b}}\rightarrow 0$, i.e., $I\rightarrow\infty$. Moreover, $A$ and $\alpha_{i}$ denote the total number of drug molecules released during $T$ and the number of drug molecules released at time $t_{i}$, respectively.

II-A2 Molecular Communication System

For the considered MC system, we assume Brownian motion of the transceivers and signaling molecules. We assume multi-frame communication between mobile $\mathrm{Tx}$ and $\mathrm{Rx}$ with instantaneous molecule release for each bit transmission. Hence, $T$, $T_{\mathrm{b}}$, $t_{i}$, and $I$ denote the duration of a bit frame, the duration of one bit interval, the beginning of the $i$-th bit interval, and the number of bits in a frame, respectively. For an arbitrary bit frame, let $b_{i}$, $i\in\{1,\dots,I\}$, denote the $i$-th bit in the bit frame. We assume that symbols [math] and $1$ are transmitted independently and with equal probability. Thus, the probability of transmitting $\tilde{b}_{i}$ is $\Pr(\tilde{b}_{i})=1/2$, where $\Pr(\cdot)$ denotes probability and $\tilde{b}_{i}\in\left\{0,1\right\}$ is a realization of $b_{i}$. We assume that on-off keying modulation is employed. At time $t_{i}$, $\mathrm{Tx}$ releases $\alpha_{i}$ molecules to transmit bit 1 and no molecules for bit 0. Then, $A=\sum_{i=1}^{I}\alpha_{i}$ is the total number of molecules available for transmission in a given bit frame.

II-B Time-variant CIR and Received Signal

Considering again the general system model, we now model the channel between $\mathrm{Tx}$ and $\mathrm{Rx}$ as well as the received signals at $\mathrm{Rx}$ for drug delivery and MC systems, respectively.

II-B1 Time-variant CIR

Let $h(t,\tau)$ denote the hitting rate, i.e., the absorption rate of a given molecule, at time $\tau$ after its release at time $t$ at the center of $\mathrm{Tx}$. Then, for an infinitesimally small observation window $\Delta\tau$, i.e., $\Delta\tau\rightarrow 0$, we can interpret $h(t,\tau)\Delta\tau$ as the probability of absorption of a molecule by $\mathrm{Rx}$ between times $\tau$ and $\tau+\Delta\tau$ after its release at time $t$. The hitting rate $h(t,\tau)$ is also referred to as the CIR since it completely characterizes the time-variant channel, which is assumed to be linear.

For a given distance between $\mathrm{Tx}$ and $\mathrm{Rx}$, $r(t)$, the CIR $h(t,\tau)$ of a diffusive mobile MC system at time $\tau$ is given by [2, 1]

[TABLE]

where $h(t,\tau)=0$, for $\tau\leq 0$. Here, $D_{1}$ is the effective diffusion coefficient capturing the relative motion of the signaling molecules and $\mathrm{Rx}$, i.e., $D_{1}=D_{X}+D_{\mathrm{Rx}}$, see [25, Eq. (8)]. In the considered MC system, due to the motion of the transceivers, the distance $r(t)$ is a random variable, and thus, the CIR $h(t,\tau)$ is time-variant and should be modeled as a stochastic process [12].

II-B2 Received Signal for Drug Delivery System

In drug delivery, the absorption rate ultimately determines the therapeutic impact of the drug [18, 20]. Thus, we formally define the absorption rate as the desired received signal, and make achieving a desired absorption rate the objective for system design. Recall that $h(t,\tau)\Delta\tau$, $\Delta\tau\rightarrow 0$, is the probability of absorption of a molecule by $\mathrm{Rx}$ between times $\tau$ and $\tau+\Delta\tau$ after the release at time $t$. If $\alpha_{i}$ molecules are released at $\mathrm{Tx}$ at time $t_{i}$, the expected number of molecules absorbed at $\mathrm{Rx}$ between times $t$ and $t+\Delta t$, for $\Delta t\rightarrow 0$, due to this release is $\alpha_{i}h(t_{i},t-t_{i})\Delta t$. During the period $[0,t]$, the total number of released drug molecules is $A_{t}=\sum_{i}\alpha_{i},\forall i|t_{i}<t,$ and the expected number of drug molecules absorbed between times $t$ and $t+\Delta t$, for $\Delta t\rightarrow 0$ is given by $s(t)=\sum_{i}\alpha_{i}h(t_{i},t-t_{i})\Delta t,\forall i|t_{i}<t$. Let $g(t)$ denote the absorption rate of drug molecules $\mathrm{X}$ at $\mathrm{Rx}$ at time $t$, i.e., $g(t)=s(t)/\Delta t$, $\Delta t\rightarrow 0$. Then, we have

[TABLE]

As mentioned before, the absorption rate $g(t)$, i.e., the received signal, of the tumor cells directly affects the healing efficacy of the drug. Hence, we will design the drug delivery system such that $g(t)$ does not fall below a prescribed value. Since $g(t)$ is a function of $h(t_{i},t-t_{i})$, it is random due to the diffusion of $\mathrm{Tx}$. Therefore, the design of the drug delivery system has to take into account the statistical properties of $g(t)$, which can be obtained from the results of the statistical analysis of $h(t_{i},t-t_{i})$.

II-B3 Received Signal for MC System

For the MC system design, the received signal, denoted by $q_{i}$, is defined as the number of $\mathrm{X}$ molecules absorbed at $\mathrm{Rx}$ during bit interval $T_{\mathrm{b}}$ after the transmission of the $i$-th bit at $t_{i}$ by $\mathrm{Tx}$ as the received signal, denoted by $q_{i}$. We detect the transmitted information based on the received signal, $q_{i}$. It has been shown in [6] that $q_{i}$ follows a Binomial distribution that can be accurately approximated by a Gaussian distribution when $\alpha_{i}$ is large, which we assume here. We focus on the effect of the transceivers’ movements on the MC system performance and design the optimal release profile of $\mathrm{Tx}$ to account for these movements. We assume the bit interval to be sufficiently long such that most of the molecules have been captured by or have moved far away from $\mathrm{Rx}$ before the following bit is transmitted, i.e., ISI is negligible. We note that enzymes [6] and reactive information molecules, such as acid/base molecules [26, 27], may be used to speed up the molecule removal process and to increase the accuracy of the ISI-free assumption. Moreover, we model external noise sources in the environment as Gaussian background noise with mean and variance equal to $\eta$ [1]. Thus, we have

[TABLE]

where $\mu_{i,0}=\sigma_{i,0}^{2}=\eta$, $\mu_{i,1}=\alpha_{i}p(t_{i},T_{\mathrm{b}})+\eta$, $\sigma_{i,1}^{2}=\alpha_{i}p(t_{i},T_{\mathrm{b}})(1-p(t_{i},T_{\mathrm{b}}))+\eta$. Here, $\mathcal{N}\left(\mu,\sigma^{2}\right)$ denotes a Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$. $p(t,T_{\mathrm{b}})$ denotes the probability that a signaling molecule is absorbed during bit interval $T_{\mathrm{b}}$ after its release at time $t$ at the center of $\mathrm{Tx}$. For a given distance $r(t)$, $p(t,T_{\mathrm{b}})$ is given by [2]

[TABLE]

where $\mathrm{erfc}(\cdot)$ is the complementary error function. Since $r(t)$ is a random variable and $p(t,T_{\mathrm{b}})$ is a function of $r(t)$, $p(t,T_{\mathrm{b}})$ and any function of $p(t,T_{\mathrm{b}})$, e.g., the received signal $q_{i}$, are random processes. Moreover, $p(t,T_{\mathrm{b}})$ is also a function of $h(t,\tau)$. Hence, for MC system design, we have to take into account the statistical properties of $p(t,T_{\mathrm{b}})$, which can be obtained based on the proposed statistical analysis of $r(t)$ and $h(t,\tau)$.

In summary, the design of both drug delivery and MC systems depends on the statistical properties of the CIR, $h(t,\tau)$, and $r(t)$ including their means, variances, PDFs, and CDFs, which will be analyzed in the next section.

III Stochastic Channel Analysis

In this section, we first analyze the distribution of the distance between the transceivers, $r(t)$, and then use it to derive the statistics of the time-variant CIR, $h(t,\tau)$, and $p(t,T_{\mathrm{b}})$ as a function $h(t,\tau)$. In particular, we develop analytical expressions for the mean, variance, PDF, and CDF of $h(t,\tau)$ and the PDF and CDF of $p(t,T_{\mathrm{b}})$.

III-A Distribution of the $\mathrm{Tx}$-$\mathrm{Rx}$ Distance for a Diffusive System

In the 3D space, $r(t)$ is given by $r(t)=\sqrt{\sum_{d\in\{x,y,z\}}(r_{d,\mathrm{Rx}}(t)-r_{d,\mathrm{Tx}}(t))^{2}}$, where $r_{d,\mathrm{Tx}}(t)$ and $r_{d,\mathrm{Rx}}(t)$, $d\in\{x,y,z\}$, are the Cartesian coordinates representing the positions of $\mathrm{Tx}$ and $\mathrm{Rx}$ at time $t$, respectively. Let us assume, without loss of generality, that the diffusion of $\mathrm{Tx}$ and $\mathrm{Rx}$ starts at $t=0$. Then, given the Brownian motion model for the mobility of $\mathrm{Tx}$ and $\mathrm{Rx}$, we have $r_{d,\mathrm{Tx}}(t)\sim\mathcal{N}\left(r_{d,\mathrm{Tx}}(t=0),2D_{\mathrm{Tx}}t\right)$ and $r_{d,\mathrm{Rx}}(t)\sim\mathcal{N}\left(r_{d,\mathrm{Rx}}(t=0),2D_{\mathrm{Rx}}t\right)$, where we assume that $r_{d,\mathrm{Tx}}(t=0)$ and $r_{d,\mathrm{Rx}}(t=0)$ are known. Let us define $r_{d}(t)=r_{d,\mathrm{Rx}}(t)-r_{d,\mathrm{Tx}}(t)$. Then, we have $r_{d}(t)\sim\mathcal{N}\left(r_{d}(t=0),2D_{2}t\right)$, where $D_{2}=D_{\mathrm{Tx}}+D_{\mathrm{Rx}}$ is the effective diffusion coefficient capturing the relative motion of $\mathrm{Tx}$ and $\mathrm{Rx}$, see [25, Eq. (10)]. Given the Gaussian distribution of $r_{d}(t)$, we know that [28]

[TABLE]

follows a noncentral chi-distribution, i.e., $\gamma\sim\mathcal{X}_{k}(\lambda)$, with $k=3$ degrees of freedom and parameter $\lambda=\sqrt{\frac{\sum_{d\in\{x,y,z\}}r_{d}^{2}(t=0)}{2D_{2}t}}=\frac{r_{0}}{\sqrt{2D_{2}t}}$, where $r_{0}$ denotes $r(t=0)$. The statistical properties of random variable $r(t)$ are provided in the following lemma.

Lemma 1

The mean, variance, PDF, and CDF of random variable $r(t)$, which represents the distance between the centers of the diffusive mobile $\mathrm{Tx}$ and $\mathrm{Rx}$, are given by, respectively,

[TABLE]

where $\mathrm{erf}(\cdot)$ is the error function, $\mathbf{Q}_{M}\left(a,b\right)$ is the Marcum Q-function [29], $\mathrm{E}\left\{\cdot\right\}$ denotes statistical expectation, $\mathrm{Var}\{\cdot\}$ denotes variance, and $f_{\left\{\cdot\right\}}(\cdot)$ and $F_{\left\{\cdot\right\}}(\cdot)$ denote the PDF and CDF of the random variable in the subscript, respectively.

Proof:

Please refer to Appendix A. ∎

Remark 1

From (6) and (7), we can observe that when $t\rightarrow\infty$, we have $\exp\left(-\frac{r_{0}^{2}}{4D_{2}t}\right)\rightarrow 1$ and $\mathrm{erf}\left(\frac{r_{0}}{\sqrt{4D_{2}t}}\right)$ $\rightarrow~{}0$ and, as a result, $\mathrm{E}\left\{r(t)\right\}\rightarrow\infty$. Intuitively, because of diffusion, the transceivers eventually move far away from each other on average.

Remark 2

We note that (8) was derived under the assumption that $\mathrm{Tx}$ can diffuse in the entire 3D environment. However, in reality, $\mathrm{Tx}$ cannot move inside $\mathrm{Rx}$, i.e., it does not interact with $\mathrm{Rx}$, and thus will be reflected when it hits $\mathrm{Rx}$’s boundary. Hence, the actual $f_{r(t)}(r)$, derived in [25], differs from (8), e.g., $f_{r(t)}(r)=0$ for $r<a_{\mathrm{tx}}+a_{\mathrm{rx}}$. However, for very small $r$, i.e., $r\approx 0$, (8) approaches zero. Hence, (8) is a valid approximation for the actual $f_{r(t)}(r)$. The validity of this approximation is evaluated in Section VI via simulations, where, in our particle-based simulation, $\mathrm{Tx}$ is reflected upon collision with $\mathrm{Rx}$ [30].

III-B *Statistical Moments of Time-variant CIR *

In this subsection, we derive the statistical moments of the time-variant CIR, i.e., mean $m(t,\tau)$ and variance $\sigma^{2}(t,\tau)$. In particular, the mean of the time-variant CIR, $m(t,\tau)$, can be written as

[TABLE]

A closed-form expression for (10) is provided in the following theorem.

Theorem 1

The mean of the impulse response of a time-variant channel with diffusive molecules released by a diffusive transparent transmitter and captured by a diffusive absorbing receiver is given by

[TABLE]

where $u(t,\tau)$, $v(t,\tau)$, and $w(t,\tau)$ are defined, for compactness, as follows

[TABLE]

Proof:

Substituting (1) and (8) into (10) and using the integrals given by [31, Eq. (2.3.15.4) and Eq. (2.3.15.7)], we obtain the expression for $m(t,\tau)$ in (11). ∎

Remark 3

$m(t,\tau)$* is a function of time $t$. Hence, $h(t,\tau)$ is a non-stationary stochastic process. In general, at large $t$, $m(t,\tau)$ decreases when $t$ increases and eventually approaches zero when $t~{}\rightarrow~{}\infty$. This means that as $t$ increases, the molecules released by $\mathrm{Tx}$, on average, have a decreasing chance of being absorbed by $\mathrm{Rx}$ since the transceivers move away from each other as mentioned in Remark 1.*

In order to obtain the variance of $h(t,\tau)$,

[TABLE]

we first need to find an expression for the second moment $\phi(t,\tau)$, defined as $\phi(t,\tau)=\mathrm{E}\left\{h^{2}(t,\tau)\right\}$. The following corollary provides an analytical expression for $\phi(t,\tau)$.

Corollary 1

$\phi(t,\tau)$* is given by*

[TABLE]

where

[TABLE]

Proof:

From the definition, we have

[TABLE]

Substituting (1) and (8) into (16) and simplifying the expression, we obtain (14). ∎

Remark 4

The expression in (14) comprises integrals of the form $\int_{0}^{\infty}\exp\left(ax^{2}+bx\right)/x\ \mathrm{d}x$, where $a$ and $b$ are constants. Such integrals cannot be obtained in closed form. However, the integrals can be evaluated numerically in a straightforward manner.

III-C *Distribution Functions of the Time-variant CIR *

In this subsection, we derive analytical expressions for the PDF and CDF of $h\left(t,\tau\right)$. The PDF of $h\left(t,\tau\right)$ is given in the following theorem.

Theorem 2

The PDF of the impulse response of a time-variant channel with diffusive molecules released by a diffusive transparent transmitter and captured by a diffusive absorbing receiver is given by

[TABLE]

where $\hat{h}\left(r,\tau\right)$ denotes $h\left(t,\tau\right)$, given by (1), as a function of $r(t)$ and $\tau$, $f_{r(t)}(r)$ is given by (8), $r_{1}(h)$ and $r_{2}(h)$, $r_{1}(h)<r_{2}(h)$, are the solutions of the equation $\hat{h}\left(r,\tau\right)=h$, $h^{\star}$ is the maximum value of $\hat{h}\left(r,\tau\right)$ for all values of $r(t)$, and $\hat{h}^{\prime}(r,\tau)$ is given by

[TABLE]

Proof:

Please refer to Appendix B.

∎

As stated in the proof of Theorem 2, there are two different values of $r(t)$, $r_{1}$ and $r_{2}$, leading to the same value of $\hat{h}(r,\tau)$, i.e., $h(t,\tau)$, when $0\leq h<h^{\star}$. Hence, the PDF of $h(t,\tau)$ is a function of the PDFs of these two values of $r(t)$. However, when $h(t,\tau)$ reaches its maximum, $f_{h\left(t,\tau\right)}(h)$ approaches infinity and does not depend on $f_{r(t)}(r(h))$ since the probability of $h=h^{\star}$, i.e., $\Pr(h=h^{\star})={f_{h\left(t,\tau\right)}(h)}\ \mathrm{d}h$, is finite and $\mathrm{d}h$ approaches [math] at $h=h^{\star}$.

The CDF of $h(t,\tau)$ is given in the following corollary.

Corollary 2

The CDF of the impulse response of a time-variant channel with diffusive molecules released by a diffusive transparent transmitter and captured by a diffusive absorbing receiver is given by

[TABLE]

where $F_{r(t)}(r)$ is given by (9).

Proof:

From the definition of the CDF and (17), we have

[TABLE]

where $\check{r}_{1}$ and $\check{r}_{2}$, $\check{r}_{1}<\check{r}_{2}$, are the solutions of the equation $\hat{h}\left(\check{r},\tau\right)=\check{h}$. This completes the proof. ∎

Similar to the PDF, the CDF of $h\left(t,\tau\right)$ also depends on the CDFs of two values of $r(t)$, i.e., $r_{1}(h)$ and $r_{2}(h)$.

III-D Distribution Functions of $p(t,T_{\mathrm{b}})$

Calculating the mean of $p(t,T_{\mathrm{b}})$ involves an integral of the form $\int_{0}^{\infty}\mathrm{erfc}\left(ax\right)\exp(-b^{2}x^{2}+cx)\mathrm{d}x$, with appropriate constants $a,b,c>0$, for which a closed-form expression is not known. However, based on the results in Subsections III-A and III-C, we obtain the PDF and CDF of $p(t,T_{\mathrm{b}})$ in the following corollary.

Corollary 3

The PDF and CDF of the probability that a diffusive molecule is absorbed by a diffusive absorbing receiver during an interval $T_{\mathrm{b}}$ after its release at time $t$ by a diffusive transparent transmitter are, respectively, given by

[TABLE]

[TABLE]

where $f_{r(t)}(r)$ and $F_{r(t)}(r)$ are given by (8) and (9), respectively. Here, $\tilde{r}(p)$ is the solution of the equation $p(t,T_{\mathrm{b}})=p$ and $p^{\prime}(\tilde{r})$ is given by

[TABLE]

Proof:

The proof of Corollary 3 follows the same steps as the proof of Theorem 2 and Corollary 2 and exploits that $p(t,T_{\mathrm{b}})$ is a function of $r(t)$ as shown in (4). From (23), we observe that $p^{\prime}(\tilde{r})<0$ so the equation $p(t,T_{\mathrm{b}})=p$ has only one solution. Then, we apply the relations for the PDFs and CDFs of functions of random variables [32] to obtain (21) and (22). ∎

The mean, variance, PDF, and CDF of $h\left(t,\tau\right)$ and $p(t,T_{\mathrm{b}})$ can be exploited to design efficient and reliable synthetic MC systems. As examples, we consider the design and analysis of drug delivery and MC systems in the following two sections.

IV Drug Delivery System Design

In this section, we apply the derived stochastic parameters of the time-variant CIR for absorbing receivers for the design and performance evaluation of drug delivery systems.

IV-A Controlled-Release Design

The treatment of many diseases requires the diseased cells to absorb a minimum rate of drugs during a prescribed time period at minimum cost [20]. To design an efficient drug delivery system satisfying this requirement, we minimize the total number of released drug molecules, $A=\sum_{i=1}^{I}\alpha_{i}$, subject to the constraint that the absorption rate $g(t)$ is equal to or larger than a target rate, $\theta(t)$, for a period of time, denoted by $T_{\mathrm{Rx}}$. We allow $\theta(t)$ to be a function of time so that the designed system can satisfy different treatment requirements over time. Since $g(t)$ is random, we cannot always guarantee $g(t)\geq\theta(t)$. Hence, we will design the system based on the first and second order moments of $g(t)$ and use the PDF and CDF of $g(t)$ to evaluate the system performance. In particular, we reformulate the constraint such that the mean of $g(t)$ minus a certain deviation is equal to or above the threshold $\theta(t)$ during $T_{\mathrm{Rx}}$, i.e., $\mathrm{E}\left\{g(t)\right\}-\beta\mathrm{V}\left\{g(t)\right\}\geq\theta(t)$, $0\leq t\leq T_{\mathrm{Rx}}$, where $\mathrm{V}\{\cdot\}$ denotes standard deviation and $\beta$ is a coefficient determining how much deviation from the mean is taken into account. Based on (2), the constraint can be written as a function of $\alpha_{i}$ as follows

[TABLE]

for $0\leq t\leq T_{\mathrm{Rx}}$. Inequality $(a)$ in (24) is due to \mathrm{E}\left\{g(t)\right\}=\mathrm{E}\bigg{\{}\underset{\forall i|t_{i}<t}{\sum}\alpha_{i}h\left(t_{i},t-t_{i}\right)\bigg{\}}=\sum_{\forall i|t_{i}<t}\alpha_{i}  $\times\mathrm{E}\left\{h\left(t_{i},t-t_{i}\right)\right\}$ and Minkowski’s inequality [33]:

[TABLE]

Note that we may not be able to find $\alpha_{i}$ such that (24) holds for all values of $\beta$ and $\theta(t)$. For example, when $\beta$ is too large, $\mathrm{E}\left\{g(t)\right\}-\beta\mathrm{V}\left\{g(t)\right\}$ can be negative and hence (24) cannot be satisfied for $\theta(t)>0$. However, when $\mathrm{E}\left\{h\left(t_{i},t-t_{i}\right)\right\}>\beta\mathrm{V}\left\{h\left(t_{i},t-t_{i}\right)\right\}$, i.e., either $\beta$ or $\mathrm{V}\left\{h\left(t_{i},t-t_{i}\right)\right\}$ are small, such that $\beta\mathrm{V}\left\{h\left(t_{i},t-t_{i}\right)\right\}$ is sufficiently small, we can always find $\alpha_{i}$ so that (24) holds for arbitrary $\theta(t)$. Since time $t$ is a continuous variable, the constraint in (24) has to be satisfied for all values of $t$, $0\leq t\leq T_{\mathrm{Rx}}$, and thus there is an infinite number of constraints, each of which corresponds to one value of $t$. Therefore, we simplify the problem by relaxing the constraints to hold only for a finite number of time instants $t=t_{n}=n\Delta t_{n}$, where $n=1,\dots,N$ and $\Delta t_{n}=T_{\mathrm{Rx}}/N$. Then, the proposed optimization problem for the design of $\alpha_{i}$ is formulated as follows

[TABLE]

where $m\left(t,\tau\right)$ and $\sigma\left(t,\tau\right)$ are the mean (11) and the standard deviation (13) of $h\left(t,\tau\right)$, respectively. Since $m(t,\tau)$ and $\sigma(t,\tau)$ do not oscillate but are well-behaved and smooth functions of $t$ as shown in Section VI, a small value of $N$ (e.g., $N=5$) is usually enough to meet the continuous constraint (24) for all $t$. Having $m\left(t,\tau\right)$ in (11) and $\sigma\left(t,\tau\right)$ in (13) and treating the $\alpha_{i}$ as real numbers, (26) can be readily solved numerically as a linear program. We note that although the numbers of drug molecules $\alpha_{i}$ are integers, for tractability, we solve (26) for real $\alpha_{i}$ and quantize the results to the nearest integer values.

We note that the problem in (26) is statistical in nature and provides guidance for the offline design of the drug delivery system.

IV-B System Performance

Since $g(t)\geq\theta(t)$ is required for proper operation of the system, we evaluate the system performance in terms of the probability that the drug absorption rate satisfies the target rate $g(t)\geq\theta(t)$, denoted by $P_{\theta}(t)=\mathrm{Pr}\left\{g(t)\geq\theta(t)\right\}$. $P_{\theta}(t)$ is given in the following theorem.

Theorem 3

The system performance metric $P_{\theta}(t)$ can be expressed as

[TABLE]

where $\ast$ denotes convolution, and $i=1,2,\dots$ satisfies $t_{i}\leq t$. In (27), we define $f_{\alpha_{i}h\left(t_{i},t-t_{i}\right)}\left(\theta(t)\right)=1/\alpha_{i}\times f_{h\left(t_{i},t-t_{i}\right)}\left(\theta(t)/\alpha_{i}\right)$ and $F_{\alpha_{i}h\left(t_{i},t-t_{i}\right)}\left(\theta(t)\right)=F_{h\left(t_{i},t-t_{i}\right)}\left(\theta(t)/\alpha_{i}\right)$.

Proof:

From the definition of the CDF, we have

[TABLE]

Due to the summation of independent random variables in (2), i.e., independent releases at $t_{i}$, we have

[TABLE]

Substituting (29) into (28), then using the integration property of the convolution, i.e.,

[TABLE]

and using the definition of the CDF, we obtain (27). ∎

We note that the analytical expressions for the PDF and CDF of $h(t,\tau)$ in Theorem 2 and Corollary 2, respectively, are not in closed form. Nevertheless, the evaluation of the system performance in (27) can be approximated by a discrete convolution which can be easily evaluated numerically.

Furthermore, we note that a minimum value of $P_{\theta}(t)$ can be guaranteed based on the statistical moments of the CIR without knowledge of the PDF and the CDF as shown in the following proposition.

Proposition 1

For a given solution of (24), a lower bound on $P_{\theta}(t)=\mathrm{Pr}\left\{g(t)\geq\theta(t)\right\}$ is given as follows

[TABLE]

Proof:

By using (24) and the Chebyshev inequality [34], we obtain

[TABLE]

where $(a)$ can be obtained by expanding the absolute value on the right-hand side, $(b)$ is due to the Chebyshev inequality, and $(c)$ is due to (24). This completes the proof. ∎

Remark 5

Proposition 1 is not only useful for evaluating the system performance, but also provides a guideline for the design of the release profile of drugs in (26). For example, to ensure a high absorption rate probability of $P_{\theta}(t)\geq 0.75$, from (31), we need to set the $\beta$ coefficient in (26) as $\beta=2$. Note that a useful bound can only be obtained based on (31) when $\beta>1$ and (24) is satisfied.

V MC System Design for Imperfect CSI

In this section, we apply the stochastic analysis presented in Section III for the design of MC systems with imperfect CSI. The CSI is imperfect due to the movement of the transceivers and assumed to be known only at the beginning of a bit frame. In particular, we optimize three design parameters of a diffusive mobile MC system employing on-off keying modulation and threshold detection, namely the detection threshold at $\mathrm{Rx}$, the release profile at $\mathrm{Tx}$, and the time duration of a bit frame. By choosing the optimal values of those three parameters, we can improve the system performance while keeping the overall system relatively simple. First, we optimize the detection threshold for minimization of the maximum BER in a frame assuming a uniform release profile. This approach can be employed in very simple MC systems where $\mathrm{Tx}$ is not capable of adjusting the number of released molecules. Second, we optimize the release profile at $\mathrm{Tx}$ for minimization of the maximum BER in a frame, assuming a fixed detection threshold and a fixed number of molecules available for transmission in the frame. This second approach to MC optimization improves the system performance in terms of BER but requires a mechanism to control the number of molecules released at $\mathrm{Tx}$. Third, we design the optimal duration of the bit frame satisfying a constraint on the efficiency of molecule usage. Thus, this third approach improves the system performance in terms of the efficiency of molecule usage. The three proposed designs can be performed offline. Furthermore, they can be combined with each other or carried out separately depending on the capabilities and requirements of the system. For all three designs, as a first step, we need to derive the BER as a function of the number of released molecules.

V-A Detection and BER

We consider a simple threshold detector at $\mathrm{Rx}$, where the received signal $q_{i}$ is compared with a detection threshold, denoted by $\xi$, in order to determine the detected bit $\hat{b}_{i}$ as follows

[TABLE]

Given the assumption of no ISI and $\Pr(\tilde{b}_{i})=1/2$, from (3) and (33), the error probability of the $i$-th bit, denoted by $P_{\mathrm{b}}(b_{i})$, can be simplified as [25, Eq. (12)]

[TABLE]

where $f_{r(t)}(r_{i})$ is given in (8), $r_{i}$ is $r(t_{i})$ for brevity, and $\zeta_{i}(\xi,\alpha_{i})=\frac{\xi-\mu_{i,1}}{\sigma_{i,1}\sqrt{2}}=\frac{\xi-\left(\alpha_{i}p(t_{i},T_{\mathrm{b}})+\eta\right)}{\sqrt{2\left(\alpha_{i}p(t_{i},T_{\mathrm{b}})(1-p(t_{i},T_{\mathrm{b}}))+\eta\right)}}$. Note that $P_{\mathrm{b}}(b_{i})$ depends on $i$ since the distance $r(t_{i})$ between the transceivers is a function of release time $t_{i}$.

V-B Optimal Detection Threshold for Uniform Release

We first consider system design for uniform release, where the number of available molecules is uniformly allocated across all bits of a frame. To facilitate reliable communication, our objective is to optimize the detection threshold, $\xi$, such that the maximum value of the error rate of the bits in a frame is minimized, given the total number of available molecules in a frame, $A$, i.e.,

[TABLE]

From (34), the problem is equivalent to

[TABLE]

The following lemma reveals the convexity of the problem in (36).

Lemma 2

For $\eta<\xi<\mu_{i,1}$, the objective function in (36) is convex in $\xi$.

Proof:

Please refer to Appendix C. ∎

Note that $\eta<\xi<\mu_{i,1}$ is intuitively satisfied for typical system parameters since the decision threshold should be higher than the average noise level when bit ”[math]” is sent but should not exceed the mean of the received signal when bit ”$1$” is sent. Otherwise, a high error rate would result. Due to the convexity of problem (36), the global optimum $\xi$ can be easily obtained by numerical methods such as the interior-point method [32].

V-C Optimal Release with Fixed Detection Threshold

For the second proposed design, we aim to optimize the release profile, i.e., the number of molecules available for release for each bit, $\alpha_{i}$, such that $\max_{i}\left\{P_{\mathrm{b}}(b_{i})\right\}$ is minimized given a total number of molecules $A$ that are available for release in a frame

[TABLE]

where $\bm{\alpha}=\left[\alpha_{i},\alpha_{2},\dots,\alpha_{I}\right]$.

For a given threshold $\xi$, we can re-express (37) based on (34) as

[TABLE]

The following lemma states the convexity of the optimization problem in (38).

Lemma 3

For $\eta<\xi<\mu_{i,1}$, the objective function in (38) is convex in $\bm{\alpha}$.

Proof:

Please refer to Appendix D. ∎

Hence, the global optimum of (38) can be readily obtained by numerical methods such as the interior-point method.

Note that, for tractability, similar to the proposed drug delivery design, we solve (36) and (38) for real $\alpha_{i}$ and quantize the results to the nearest integer values.

V-D Optimal Time Duration of a Bit Frame

In the third proposed design, we consider the molecule usage efficiency for communication. We evaluate the efficiency based on $p(t,T_{\mathrm{b}})$, i.e., the probability that a signaling molecule is absorbed during bit interval $T_{\mathrm{b}}$ after its release at time $t$. If $p(t,T_{\mathrm{b}})$ is too small, none of the released molecules may actually reach the receiver and thus the molecules are wasted, i.e., the system has low efficiency. Hence, we want to keep $p(t,T_{\mathrm{b}})$ above a certain value, denoted by $\psi$. Intuitively, as on average $h(t,\tau)$ decreases over time $t$, $p(t,T_{\mathrm{b}})$, which is the integral over $h(t,\tau)$ with respect to $\tau$, also on average decreases over time. Therefore, our objective is to choose the maximum duration of a bit frame, denoted by $T^{\star}$, such that $p(t,T_{\mathrm{b}})>\psi$ for $t\leq T^{\star}-T_{\mathrm{b}}$, where $t$ is the release time.

Since $p(t,T_{\mathrm{b}})$ is a random process, we cannot enforce $p(t,T_{\mathrm{b}})>\psi$ but can only bound the probability that $p(t,T_{\mathrm{b}})>\psi$ is satisfied, i.e., $\Pr\left(p(t,T_{\mathrm{b}})>\psi\right)\geq P$, where $P$ is a design parameter. Moreover, we have

[TABLE]

where equality $(a)$ is due to (22). As such, we can re-express the problem as maximizing the duration of a bit frame such that $F_{r(t)}(\tilde{r}(\psi))\geq P$ holds. To this end, in the following lemma, we analyze $F_{r(t)}(\tilde{r}(\psi))$ as a function of time $t$.

Lemma 4

$F_{r(t)}(\tilde{r}(\psi))$* is a decreasing function of time $t$.*

Proof:

Please refer to Appendix E. ∎

Since Lemma 4 shows that $F_{r(t)}(\tilde{r}(\psi))$ is a decreasing function of time, the maximum duration of a bit frame satisfying $F_{r(t)}(\tilde{r}(\psi))\geq P$ can be found by solving $F_{r(T^{\star}-T_{\mathrm{b}})}(\tilde{r}(\psi))=P$, where $F_{r(T^{\star}-T_{\mathrm{b}})}(\tilde{r}(\psi))$ is given in (9).

Remark 6

If multiple frames are transmitted, the proposed design framework can be applied to each frame, respectively. However, the optimal designs may be different for different frames due to the moving transceivers, whose distances are assumed to be perfectly estimated at the start of each frame.

Remark 7

Here, we discuss a system with an absorbing receiver. Nevertheless, the proposed optimal design framework can also be applied to transparent receivers. For a transparent receiver, $p(t,T_{\mathrm{b}})$ is the probability that a molecule is observed inside the volume of the transparent receiver at time $T_{\mathrm{b}}$ after its release at time $t$ at the center of $\mathrm{Tx}$.

VI Numerical Results

In this section, we provide numerical results to evaluate the accuracy of the derived expressions and analyze the performance of the MC systems in the considered application scenarios. We use the set of simulation parameters summarized in Table I, unless stated otherwise. The parameters are chosen to match the actual system parameters in drug delivery systems, as will be explained in detail in Subsection VI-B.

VI-A Time-variant Channel Analysis

In this subsection, we numerically analyze the time-variant MC channel. For verification of the accuracy of the expressions derived in Section III, we employ a hybrid particle-based simulation approach. In particular, we use particle-based simulation of the Brownian motion of the transceivers to generate realizations of the random distance between $\mathrm{Tx}$ and $\mathrm{Rx}$, $r(t)$. Then, we use Monte Carlo simulation to obtain the desired statistical results by suitably averaging the CIRs (1) obtained for the realizations of $r(t)$. For particle-based simulation of the Brownian motion of $\mathrm{Tx}$, $\mathrm{Tx}$ performs a random walk with a random step size in space in every discrete time step of length $\Delta t^{\mathrm{st}}=$1\text{,}\mathrm{ms}$$. The length of each step in space is modeled as a Gaussian random variable with zero mean and standard deviation $\sqrt{2D_{\mathrm{Tx}}\Delta t^{\mathrm{st}}}$. Furthermore, we also take into account the reflection of $\mathrm{Tx}$ upon collision with $\mathrm{Rx}$. When $\mathrm{Tx}$ hits $\mathrm{Rx}$, we assume that it bounces back to the position it had at the beginning of the considered simulation step [30].

Fig. 3 shows the mean of the CIR, $m(t,\tau)$, as a function of $\tau$. In general, for large $\tau$, $m(t,\tau)$ decreases when $t$ increases as expected since the transceivers move away from each other on average. For large $\tau$, $m(t,\tau)$ also decreases when $\tau$ increases as would be the case in a static system. Note that in the simulations, unlike the analysis, we have taken into account the reflection of $\mathrm{Tx}$ when it hits $\mathrm{Rx}$. Therefore, the good agreement between simulation and analytical results in Fig. 3 suggests that the reflection of $\mathrm{Tx}$ does not have a significant impact on the statistical properties of $h(t,\tau)$ and the approximation in (8) and the analytical results obtained based on it are valid.

In Fig. 3, we plot the PDF of the CIR for time instances $t=\left\{36,360,3600\right\}$\mathrm{s} and $\tau=$0.17\text{\,}\mathrm{s}. Fig. 3 shows that when $t$ increases, a smaller value of $h(t,\tau)$ is more likely to occur since, on average, the transceivers move away from each other. When $t$ is very large, it is likely that the molecules cannot reach $\mathrm{Rx}$, and hence, cannot be absorbed, consequently $h(t,\tau)\rightarrow 0$. We also observe that $h(t,\tau)$ has a maximum value and $f_{h(t,\tau)}\left(h(t,\tau)\right)\rightarrow\infty$ when $h(t,\tau)$ approaches the maximum value as stated in (17). For example, $h(t=$36\text{,}\mathrm{s}$,\tau=$0.17\text{,}\mathrm{s}$)$ is random but its maximum possible value is $0.29$ and f_{h(t=36\text{,}\mathrm{s}$,\tau=$0.17\text{,}\mathrm{s})}\left(0.29\right)\rightarrow\infty.

Fig. 3 and 3 show a perfect match between simulation and analytical results. This confirms the accuracy of our analysis of the time-variant CIR in Sections III. Since particle-based simulation is costly, in the following subsections, we adopt Monte-Carlo simulation by averaging our results over $10^{5}$ independent realizations of both the distance $r(t)$ and the CIR. The distance $r(t)$ is calculated from the locations of the transceivers, which are generated from Gaussian distributions, see Subsection III-A. In particular, $r(t)=\sqrt{\sum_{d\in\{x,y,z\}}(r_{d,\mathrm{Rx}}(t)-r_{d,\mathrm{Tx}}(t))^{2}}$, where $r_{d,\mathrm{Tx}}(t)\sim\mathcal{N}\left(r_{d,\mathrm{Tx}}(t=0),2D_{\mathrm{Tx}}t\right)$, $r_{d,\mathrm{Rx}}(t)\sim\mathcal{N}\left(r_{d,\mathrm{Rx}}(t=0),2D_{\mathrm{Rx}}t\right)$, $r_{d,\mathrm{Tx}}(t=0)=0$, and $r_{d,\mathrm{Rx}}(t=0)=1/\sqrt{3}r_{0}$. The CIR is given by (1) for each realization of $r(t)$.

VI-B Drug Delivery System Design

In this section, we provide numerical results for the considered drug delivery system. As mentioned above, the parameters in Table I are chosen to match real system parameters, e.g., the diffusion coefficient $D_{\mathrm{X}}$ of drug molecules vary from $10^{-9}$ to $10^{-14}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ [19], drug carriers have sizes $a_{\mathrm{tx}}\geq$100\text{,}\mathrm{nm} [[17](#bib.bib17)], the size of tumor cells is on the order of $\mathrm{\SIUnitSymbolMicro m}$, and drug carriers can be injected or extravasated from the cardiovascular system into the tissue surrounding the targeted diseased cell site [[18](#bib.bib18)], i.e., close to the tumor cells. The dosing periods in drug delivery systems are on the order of days [[35](#bib.bib35)], i.e., $24\text{\,}\mathrm{h}$. For simplicity, we set $N=5$ and the value of the required absorption rate is set to $\theta(t)=$1\text{\,}\mathrm{s}^{-1}. We choose $I$ relatively large to obtain small intervals $T_{\mathrm{b}}$.

In Fig. 5, we plot the number of released molecules $\alpha_{i}$ versus the corresponding release time $t_{i}$ [$\mathrm{h}$] for different system parameters. The coefficients are obtained by solving the optimization problem in (26) with $\beta=\{0,0.4,1,2\}$ for $D_{\mathrm{Tx}}=10^{-14}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ and $\beta=\{0,0.4\}$ for $D_{\mathrm{Tx}}=\{5,10\}\times 10^{-14}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$. As mentioned in the discussion of (24), we cannot choose large values of $\beta$ when the diffusion coefficient is large, i.e., the standard deviation is large, as the problem in (26) may become infeasible. Fig. 5 shows that for all considered parameter settings, we first have to release a large number of molecules for the absorption rate to exceed the threshold. Then, in the static system with $D_{\mathrm{Tx}}=$0\text{,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}, the optimal coefficient decreases with increasing time, since a fraction of the molecules previously released from $\mathrm{Tx}$ linger around $\mathrm{Rx}$ and are absorbed later. However, for the time-variant channel, $\mathrm{Tx}$ eventually diffuses away from $\mathrm{Rx}$ as time $t$ increases and hence, molecules released at later times by $\mathrm{Tx}$ will be far away from $\mathrm{Rx}$ and may not reach it. Therefore, at later times, the amount of drugs released has to be increased for the absorption rate to not fall below the threshold. For larger $D_{\mathrm{Tx}}$, $\mathrm{Tx}$ diffuses away from $\mathrm{Rx}$ faster and thus, the number of released molecules $\alpha_{i}$ have to increase faster. This type of drug release, i.e., first releasing a large amount of drugs, then reducing and eventually increasing the amount of released drugs again, is called a tri-phasic release [[36](#bib.bib36)]. Once we have designed the release profile, we can implement it by choosing a suitable drug carrier as shown in [[36](#bib.bib36)]. Moreover, as expected, for larger $\beta$, we need to release more drugs to ensure that ([26](#S4.E26)) is feasible. The black horizontal dotted line in Fig. [5](#S6.F5) is a benchmark where the $\alpha_{i},\forall i,$ are not optimized but naively set to $\alpha_{i}=\alpha_{1}=5493$. For this naive design, $A=\alpha_{1}I\approx 1.65\times 10^{7}$, whereas with the optimal $\alpha_{i}$, for $\beta=0$ and $D_{\mathrm{Tx}}=10^{-13}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$, $A=1.2\times 10^{7}$, i.e., $27\%$ less than the $A$ required for the naive design, and for $\beta=1$ and $D_{\mathrm{Tx}}=10^{-14}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$, $A=7.6\times 10^{6}$, i.e., $54\%$ less than the $A$ required for the naive design. This highlights that applying the optimal release profile can save significant amounts of drugs and still satisfy the therapeutic requirements. Moreover, as observed in Fig. [5](#S6.F5), at later times, e.g., $t_{i}>$15\text{\,}\mathrm{h} for $D_{\mathrm{Tx}}={10^{-13}}$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$$, the values of $\alpha_{i}$ required to satisfy the desired absorption rate are higher than the fixed $\alpha_{i}$ used in the naive design, i.e., the benchmark, which means that the naive design cannot provide the required absorption rate.

In Fig. 5, we plot the mean and standard deviation of the absorption rate, $\mathrm{E}\left\{g(t)\right\}$ and $\mathrm{V}\left\{g(t)\right\}$, between the $1000$-th release and the $1002$-th release for three designs. For designs 1 and 2, we assumed $D_{\mathrm{Tx}}=10^{-13}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ and $\beta=0$, and for design 3, we adopted $D_{\mathrm{Tx}}=10^{-14}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ and $\beta=1$. Note that the considered time window, e.g., between the $1000$-th release and the $1002$-th release, is chosen arbitrarily in the middle of $T$ to analyze the system behavior between individual releases. For design 1, $\mathrm{Tx}$ diffuses with $D_{\mathrm{Tx}}=10^{-13}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ but the release profile is designed without accounting for $\mathrm{Tx}$’s mobility, i.e., the adopted $\alpha_{i}$ are given by the green line in Fig. 5 obtained under the assumption of $D_{\mathrm{Tx}}=$0\text{,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$$. For designs 2 and 3, the mobility of $\mathrm{Tx}$ is taken into account. The black dashed line marks the threshold $\theta(t)$ that $g(t)$ should not fall below. It is observed from Fig. 5 that when $\mathrm{Tx}$ diffuses but the design does not take into account the mobility, the requirement that the expected absorption rate, $\mathrm{E}\left\{g(t)\right\}$, exceeds $\theta(t)$, is not satisfied for most of the time. For design 2 with $\beta=0$, we observe that $\mathrm{E}\left\{g(t)\right\}>\theta(t)$ always holds but $\mathrm{E}\left\{g(t)\right\}-\mathrm{V}\left\{g(t)\right\}>\theta(t)$ does not always hold. For design 3 with $\beta=1$, we observe that $\mathrm{E}\left\{g(t)\right\}-\mathrm{V}\left\{g(t)\right\}>\theta(t)$ always holds since $\beta>0$ enforces a gap between $\mathrm{E}\left\{g(t)\right\}$ and $\theta(t)$. In other words, even if $g(t)$ deviates from the mean, it can still exceed $\theta(t)$.

In Fig. 7, we present the system performance in terms of the probability that $g(t)\geq\theta(t)$, $P_{\theta}(t)$, for the time period between the $1000$-th and $1002$-th releases, i.e., at about $8\text{\,}\mathrm{h}$. The lines and markers denote simulation and analytical results, respectively. Fig. 7 shows a good agreement between analytical and simulation results. In Fig. 7, we observe that $P_{\theta}(t)$ increases with increasing $\beta$ because the design for larger $\beta$ enforces a larger gap between $\mathrm{E}\left\{g(t)\right\}$ and $\theta(t)$, as can be seen in Fig. 5. Moreover, for a given $\beta$, $P_{\theta}(t)$ will be different for different $D_{\mathrm{Tx}}$. In particular, for larger $D_{\mathrm{Tx}}$, $P_{\theta}(t)$ is smaller due to the faster diffusion and increasing randomness of the CIR. Moreover, in Fig. 7, the green line shows that the naive design, i.e., design 1 in Fig. 5, has very poor performance. In Fig. 7, we also observe that between two releases, $P_{\theta}(t)$ first increases due to the released drugs and then decreases due to drug diffusion. Furthermore, in Fig. 7, we also show the lower bound on $P_{\theta}(t)$ derived in Proposition 1 for $D_{\mathrm{Tx}}=10^{-14}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ and $\beta=2$, where (31) yields $P_{\theta}(t)\geq 0.75$. Fig. 7 shows that the red dash-dotted line, i.e., $P_{\theta}(t)$ for $D_{\mathrm{Tx}}=10^{-14}$ $\mathrm{m}^{2}\mathrm{/}\mathrm{s}$ and $\beta=2$, is indeed above the horizontal black dashed line, i.e., $P_{\theta}(t)=0.75$.

VI-C Molecular Communication System Design

In this subsection, we show numerical results for the second application scenario, i.e., an MC system with imperfect CSI. We apply again the system parameters in Table I except that here we set $D_{\mathrm{Rx}}=10^{-11}~{}$\mathrm{m}^{2}\mathrm{/}\mathrm{s}, $I=30$, $\eta=1$, $T=$300\text{\,}\mathrm{s}, and $T_{\mathrm{b}}=T/I=$10\text{,}\mathrm{s}$$ to also allow $\mathrm{Rx}$ to move and to reduce the transmission window compared to the drug delivery system.

In Fig. 7, we consider the optimal release design, i.e., the optimal number of molecules available for transmission of each bit in a frame, for an MC system with fixed detection thresholds and fixed $A$, $A=\left\{10^{3},10^{4},10^{5}\right\}$. The fixed detection thresholds $\xi$ are obtained from Subsection V-B by assuming uniform release. Fig. 7 reveals that in order to minimize the maximum BER in a frame, fewer molecules should be released at the beginning of the frame and the number of released molecules gradually increases with time. This is expected since, on average, for later release times, more molecules are needed to compensate for the increasing distance between the transceivers.

Fig. 9 shows the maximum BER within a frame for uniform release and the proposed release design obtained from (38), with a fixed detection threshold obtained from (36), as a function of $A$ for $T=\left\{100,300,3000\right\}$\mathrm{s}. As can be observed, the proposed optimal release profile leads to significant performance improvements compared to uniform release, especially for large $A$. For example, for $A=10^{5}$ and $T=$3000\text{\,}\mathrm{s}, the maximum BER is reduced by a factor of $8$ for optimal release compared to uniform release. On the other hand, to achieve a given desired BER, the total number of molecules $A$ required for optimal release is lower than that for uniform release. In the inset of Fig. 9, we show the BER as a function of bit index $i$ in one frame for uniform and optimal release for $A=10^{4}$ and $T=$300\text{,}\mathrm{s}$$. We observe that the optimal release achieves a lower maximum BER compared to the uniform release. We also observe that the optimal release leads to approximately the same BER for each bit in a frame which highlights the benefits of the proposed design.

Fig. 9 shows the probability that $p(t,T_{\mathrm{b}})$ is larger than a given value $\psi$, $\Pr\left(p(t,T_{\mathrm{b}})>\psi\right)$, as a function of time $t$. $\Pr\left(p(t,T_{\mathrm{b}})>\psi\right)$ provides information about the probability that a released molecule is absorbed at the Rx, i.e., the efficiency of molecule usage. We observe from Fig. 9 that $\Pr\left(p(t,T_{\mathrm{b}})>\psi\right)$ is a decreasing function of time as expected from the analysis in Subsection V-D. Moreover, for a given $t$, $\Pr\left(p(t,T_{\mathrm{b}})>\psi\right)$ is smaller for larger $\psi$. Furthermore, we can deduce the maximum time duration of a bit frame, $T^{\star}$, satisfying a required molecule usage efficiency from Fig. 9. For example, for $t\leq$300\text{,}\mathrm{s}, $\Pr\left(p(t,T_{\mathrm{b}})>0.02\right)\geq 0.8$ holds. Thus, $T^{\star}=t+T_{\mathrm{b}}=$310\text{\,}\mathrm{s} guarantees a molecule usage efficiency of $\Pr\left(p(t,T_{\mathrm{b}})>0.02\right)\geq 0.8$.

VII Conclusions

In this paper, we considered a diffusive mobile MC system with an absorbing receiver, in which both the transceivers and the molecules diffuse. We provided a statistical analysis of the time-variant CIR and its integral, i.e., the probability that a molecule is absorbed by $\mathrm{Rx}$ during a given time period. We applied this statistical analysis to two system design problems, namely drug delivery and on-off keying based MC with imperfect CSI. For the drug delivery system, we proposed an optimal release profile which minimizes the number of released drug molecules while ensuring a target absorption rate for the drugs at the diseased site during a prescribed time period. The probability of satisfying the constraint on the absorption rate was adopted as a system performance criterion and evaluated. We observed that ignoring the reality of the $\mathrm{Tx}$’s mobility for designing the release profile leads to unsatisfactory performance. For the MC system with imperfect CSI, we optimized three design parameters, i.e., the detection threshold at $\mathrm{Rx}$, the release profile at $\mathrm{Tx}$, and the time duration of a bit frame. Our simulation results revealed that the proposed MC system designs achieved a better performance in terms of BER and molecule usage efficiency compared to a uniform-release system and a system without limitation on molecule usage, respectively. Overall, our results showed that taking into account the time-variance of the channel of mobile MC systems is crucial for achieving high performance.

Appendix A Proof of Lemma 1

To prove (6), we have

[TABLE]

where $\Gamma(\cdot)$ denotes the Gamma function, equality $(a)$ is due to (5), equality $(b)$ is due to [37, Eq. (1.5)], and equality $(c)$ is obtained by applying [31, Eq. (5.2.11.6) and Eq. (5.2.11.7)]. Simplifying the final expression, we obtain (6).

To prove (7), we have

[TABLE]

where equality $(a)$ is obtained due to (5), equality $(b)$ is due to [37, Eq. (1.5)], and equality $(c)$ is obtained by applying the Maclaurin series of the exponential function. From (40) and (41), we obtain (7) since $\mathrm{Var}\left\{r(t)\right\}=\mathrm{E}\left\{r^{2}(t)\right\}-\mathrm{E}^{2}\left\{r(t)\right\}$.

To prove (8), we have

[TABLE]

where $f_{\gamma}\left(\gamma\right)$ is the PDF of $\gamma$. Equality $(a)$ in (42) exploits the fact that $\gamma$ is a function of $r(t)$ [34, Eq. (5-16)]. Equality $(b)$ in (42) is obtained from the expression for PDF $f_{\gamma}\left(\gamma\right)$ [28, Eq. (1.6)] and the relation $I_{1/2}(x)=\sqrt{\frac{2}{\pi x}}\sinh(x)$, where $I_{1/2}(x)$ is the Bessel function of the first kind and order $1/2$.

Moreover, since $\frac{r(t)}{\sqrt{2D_{\mathrm{Tx}}t}}$ follows a noncentral chi distribution, we obtain (9) as [29, Eq. (1)]

[TABLE]

Appendix B Proof of Theorem 2

For this proof, we keep in mind that $\hat{h}\left(r,\tau\right)$ and $h\left(t,\tau\right)$ are two functions of different variables but give the same value $h$ since $r$ is a function of $t$. Taking the derivative of (1) with respect to $r$, we obtain (18). From (18), we observe that $\hat{h}^{\prime}(r,\tau)=0$ is equivalent to a cubic equation in $r$, given by $ar^{3}+br^{2}+cr+d=0$, with properly defined coefficients $a$, $b$, $c$, $d$ and discriminant $\Delta=18abcd-4b^{3}d+b^{2}c^{2}-4ac^{3}-27a^{2}d^{2}$. From (18), we have $\Delta<0$ and thus $\hat{h}^{\prime}(r,\tau)=0$ has only one real valued solution, denoted by $r^{\star}$, which corresponds to the maximum value of $\hat{h}(r,\tau)$, denoted by $h^{\star}$. Then, from (18), we observe that $\hat{h}^{\prime}(r,\tau)>0$ for $r<r^{\star}$ and $\hat{h}^{\prime}(r,\tau)<0$ for $r>r^{\star}$. Therefore, the equation $\hat{h}\left(r,\tau\right)=h$ has two solutions $r_{1}(h)$ and $r_{2}(h)$, $r_{1}(h)<r_{2}(h)$, when $h<h^{\star}$, and has only one solution $r^{\star}$ when $h=h^{\star}=\hat{h}\left(r^{\star},\tau\right)$. Finally, we derive (17) by exploiting [34, Eq. (5-16)] for the PDF of functions of random variables. Moreover, for $h=h^{\star}$, $\hat{h}^{\prime}(r,\tau)=0$ so $f_{h\left(t,\tau\right)}(h)=\frac{f_{r(t)}(r^{\star})}{\hat{h}^{\prime}(r^{\star},\tau)}\rightarrow\infty$.

Appendix C Proof of Lemma 2

To prove Lemma 2, we need to prove that $\mathrm{erf}\left(\zeta_{i}(\xi,\alpha_{i})\right)-\mathrm{erf}\left(\frac{\xi-\eta}{\sqrt{2\eta}}\right)$ is convex in $\xi$. $\mathrm{erf}(\cdot)$ is a convex and non-decreasing function for negative arguments and a concave and non-decreasing function for positive arguments. For $\eta<\xi<\mu_{i,1}$, we have $\zeta_{i}(\xi,\alpha_{i})<0$ and $\frac{\xi-\eta}{\sqrt{2\eta}}>0$. Moreover, $\zeta_{i}(\xi,\alpha_{i})$ and $\frac{\xi-\eta}{\sqrt{2\eta}}$ are affine functions of $\xi$. Then, $\mathrm{erf}\left(\zeta_{i}(\xi,\alpha_{i})\right)$ is a convex and non-decreasing function of an affine function and thus is convex in $\xi$ [32, Eq. (3.10)]. $\mathrm{erf}\left(\frac{\xi-\eta}{\sqrt{2\eta}}\right)$ is a concave and non-decreasing function of an affine function and thus is concave in $\xi$ [32, Eq. (3.10)]. Therefore, $\mathrm{erf}\left(\zeta_{i}(\xi,\alpha_{i})\right)-\mathrm{erf}\left(\frac{\xi-\eta}{\sqrt{2\eta}}\right)$ is convex, which concludes the proof.

Appendix D Proof of Lemma 3

To prove Lemma 3, we need to prove that $\mathrm{erf}\left(\zeta_{i}(\xi,\alpha_{i})\right)$ is convex in $\bm{\alpha}$. First, $\mathrm{erf}(\cdot)$ is a convex and non-decreasing function for negative arguments and $\zeta_{i}(\xi,\alpha_{i})<0$ for $\xi<\mu_{i,1}$. Second, when $\eta<\xi$, $\zeta_{i}(\xi,\alpha_{i})$ is a convex function in $\bm{\alpha}$ since its Hessian matrix is positive semi-definite. Then, $\mathrm{erf}\left(\zeta_{i}(\xi,\alpha_{i})\right)$ is a convex and non-decreasing function of a convex function, and thus, it is convex in $\bm{\alpha}$ for $\eta<\xi<\mu_{i,1}$ [32, Eq. (3.10)], which concludes the proof.

Appendix E Proof of Lemma 4

To prove Lemma 4, we need to prove $\frac{\partial F_{r(t)}(r)}{\partial t}<0$. Using the Taylor series expansion of the $\sinh(\cdot)$ function in (8), we obtain

[TABLE]