Diffusive Mobile MC with Absorbing Receivers: Stochastic Analysis and Applications
Trang Ngoc Cao, Arman Ahmadzadeh, Vahid Jamali, Wayan Wicke, Phee Lep Yeoh, Jamie Evans, Robert Schober
11 Diffusive Mobile MC with Absorbing Receivers:Stochastic Analysis and Applications
Trang Ngoc Cao, Arman Ahmadzadeh, Vahid Jamali, Wayan Wicke,Phee Lep Yeoh, Jamie Evans, and Robert Schober
Abstract
This paper presents a stochastic analysis of the time-variant channel impulse response (CIR) of a threedimensional diffusive mobile molecular communication (MC) system where the transmitter, the absorbingreceiver, and the molecules can freely diffuse. In our analysis, we derive the mean, variance, probabilitydensity function (PDF), and cumulative distribution function (CDF) of the CIR. We also derive the PDF andCDF 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 imperfectchannel state information. For the first application, we exploit the mean and variance of the CIR to optimize acontrolled-release drug delivery system employing a mobile drug carrier. We evaluate the performance of theproposed release design based on the PDF and CDF of the CIR. We demonstrate significant savings in theamount of released drugs compared to a constant-release scheme and reveal the necessity of accounting forthe drug-carrier’s mobility to ensure reliable drug delivery. For the second application, we exploit the PDF ofthe distance between the mobile transceivers and the CDF of p to optimize three design parameters of an MCsystem employing on-off keying modulation and threshold detection. Specifically, we optimize the detectionthreshold at the receiver, the release profile at the transmitter, and the time duration of a bit frame. We showthat the proposed optimal designs can significantly improve the system performance in terms of the bit errorrate and the efficiency of molecule usage. I. I
NTRODUCTION
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 referencestherein. For example, the simple diffusive channel model of an unbounded three-dimensional (3D) MCsystem with impulsive point release of information carrying molecules [2] has been widely used forsystem analysis and design, see [3], [4], and references therein. Diffusion channel models with drift
This paper has been presented in part at IEEE ICC 2019. a r X i v : . [ c s . I T ] A ug [5] and chemical reactions [6] have also been considered. However, most of the previously studiedMC 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 designof mobile MC systems have gained considerable attention, e.g., see [1], [8]–[15], and references therein.In [8], a mobile ad hoc nanonetwork was considered where mobile nanomachines collect environmentalinformation and deliver it to a mobile central control unit. The mobility of the nanomachines wasdescribed 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 detectionwith non-diffusive information molecules was proposed. The authors in [10] considered adaptivedetection and inter-symbol interference (ISI) mitigation in mobile MC systems, while [11] analyzedthe 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 systemnumerically. In [12], a comprehensive framework for modeling the time-variant channels of diffusivemobile MC systems with diffusive transceivers was developed. However, all of the works mentionedabove assumed a passive receiver.On the other hand, for many MC applications, a fully absorbing receiver is considered to be a morerealistic model compared to a passive receiver as it captures the interaction between the receiver andthe information molecules, e.g., the conversion of the information molecules to a new type of moleculeor the absorption and removal of the information molecules from the environment [2], [3]. Since themolecules are removed from the environment after being absorbed by the receiver, the channel impulseresponse (CIR) for absorbing receivers is a more complicated function of the distance between thetransceivers and the receiver’s radius compared to passive receivers. Therefore, the stochastic analysisof mobile MC systems with absorbing receiver is very challenging. For the fully absorbing receiver indiffusive mobile MC systems, theoretical expressions for the average distribution of the first hittingtime, i.e., the mean of the CIR, were derived for a one-dimensional (1D) environment without driftin [13] and with drift in [14]. Based on the 1D model in [13], the error rate and channel capacityof the system were examined in [15]. However, none of these works provides a statistical analysisof 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 thetime-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 beinginjected or extravasated from the cardiovascular system into the tissue surrounding a targeted diseasedcell site, the drug carriers may not be anchored at the targeted site but may move, mostly via diffusion[16]–[19]. The diffusion of the drug carriers results in a time-variant absorption rate of the drugs evenif the drug release rate is constant. Furthermore, experimental and theoretical studies have indicatedthat the total drug dosage as well as the rate and time period of drug absorption by the receptors ofthe diseased cells are critical factors in the healing process [18], [20]. Therefore, to satisfy reducedrug 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 amountof released drugs is minimized while a desired rate of drug absorption at the diseased site during aprescribed time period is achieved. To this end, the mobility of the drug carriers and the absorptionrate of the drugs have to be accurately taken into account. This can be accomplished by exploitingthe MC paradigm where the drug carriers, diseased cells, and drug molecules are modeled as mobiletransmitters, absorbing receivers, and signaling molecules, respectively [3]. Release profile designs fordrug delivery systems based on an MC framework were proposed in [4], [21]–[23]. However, in theseworks, the transceivers were fixed and only the movement of the drug molecules was considered. Inthis paper, we exploit the analytical results obtained for the stochastic parameters of the time-variantMC channel with absorbing receiver for the optimization of the release profile of drug delivery systemswith 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 randommovements 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 communicationsystem design challenging. In this paper, we consider a mobile MC system employing on-off keyingand threshold detection and optimize three design parameters to improve the system performanceunder imperfect CSI. First, we optimize the detection threshold at the receiver for minimization of themaximum bit error rate (BER) in a frame when the number of molecules available for transmission isuniformly 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 ofthe maximum BER in a frame given a fixed number of molecules available for transmission of theentire 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 designensures that molecules are used efficiently as a molecule release occurs only if the released molecule isobserved at the receiver with sufficiently high probability. For the proposed design tasks, the results ofthe stochastic analysis of the transceivers’ positions and of the probability that a molecule is absorbedduring 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 systememploying an absorbing receiver. In particular, we derive the mean, variance, PDF, and CDF ofthe corresponding CIR. Moreover, we derive the PDF and CDF of the probability that a moleculeis absorbed during a given time period. The stochastic channel analysis is exploited for the designof drug delivery and MC systems. ‚ For drug delivery systems, the release profile is optimized for the minimization of the amount ofreleased drugs while ensuring that the absorption rate at the diseased cells does not fall below aprescribed threshold for a given period of time. We show that the proposed design requires asignificantly 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 imperfectCSI, we optimize three design parameters, namely the detection threshold at the receiver, therelease profile at the transmitter, and the time duration of a bit frame. Simulation results showsignificant performance gains for the proposed designs in terms of BER and the efficiencyof molecule usage compared to baseline systems with uniform molecule release and withoutlimitation on time duration of a bit frame, respectively. ‚ Our results reveal that the transceivers’ mobility has a significant impact on the system performanceand 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 withabsorbing receiver are expected to be useful not only for the design of the drug delivery and MCsystems considered in this paper but also for the design of detection schemes and the evaluation ofthe 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 amolecule is absorbed during a given time period, the MC system design for imperfect CSI, and thecorresponding simulation results are not included in [24].The remainder of this paper is organized as follows. In Section II, we introduce the considereddiffusive mobile MC system with absorbing receiver and the time-variant channel model. In Section
Drug molecule Diseased cellsSpherical RxSpherical Tx Absorbed drug molecule
Fig. 1. System model for drug delivery. The drug carrier and the diseased cells of a tumor are modeled as diffusive spherical transmitter( Tx ) and spherical absorbing receiver ( Rx ), respectively. The drug molecules are absorbed by Rx , when they hit its surface. Differentdistances between Tx and Rx over time are due to Tx ’s diffusion. 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. G
ENERAL S YSTEM AND C HANNEL M ODEL
In this section, we first introduce the model for a general diffusive mobile MC system with absorbingreceiver. Subsequently, we specialize the model to drug delivery and communication systems withimperfect CSI. Finally, we define the time-variant CIR and the received signal.
A. System Model
We consider a linear diffusive mobile MC system in an unbounded 3D environment with constanttemperature and viscosity. The system comprises one mobile spherical transparent transmitter, denotedby Tx , with radius a tx , one mobile spherical absorbing receiver, denoted by Rx , with radius a rx , andthe signaling molecules of type X .The movements of Tx , Rx , and X molecules are assumed to be mutually independent and followBrownian motion with diffusion coefficients D Tx , D Rx , and D X , respectively. This assumption, whichwas also made in [12] and [13], is motivated by the fact that the mobility of small objects is governedby Brownian motion.We assume that Tx releases molecules at its center instantaneously and discretely during theconsidered period of time denoted by T . Let t i and T b denote the time instant of the i -th releaseand the duration of the interval between the i -th and the p i ` q -th release, respectively. We have t i “ p i ´ q T b and i P t , . . . , I u , where I is the total number of releases during T . We denote thetime-varying distance between the centers of Tx and Rx at time t by r p t q . Furthermore, let α i and A “ ř Ii “ α i denote the number of molecules released at time t i and the total number of moleculesreleased during T , respectively. For concreteness, we specialize the considered general model to twoapplication scenarios.
1) Drug Delivery Systems:
A drug delivery system comprises a drug carrier releasing drug moleculesand diseased cells absorbing them. We model the drug carrier and diseased cells as Tx and Rx of thegeneral MC system, respectively, see Fig. 1. The drug carriers in drug delivery systems are typicallynanoparticles, such as spherical polymers or polymer chains, having a size not smaller than
100 nm [17]. Moreover, drug carriers are designed to carry drug molecules and interaction with the drug orthe receiver is not intended. Hence, the drug carriers can be modeled as mobile spherical transparenttransmitters, Tx . When the drug molecules hit the tumor, they are absorbed by receptors on thesurface of the diseased cells [18], [20]. For convenience, we model the tumor as a spherical absorbingreceiver, 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 receiverwith radius a 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 bloodinto the interstitial tissue near the diseased cells, where they start to move. We assume that theinjection position can be estimated and thus r p t “ q is known. The movement of the drug carryingnanoparticles in the tissue is caused by diffusion and convection mechanisms but diffusion is expectedto be dominant in most cases [16]–[19]. At the tumor site, the drug carrier releases drug moleculesof type X , which also diffuse in the tissue [18]. Hence, we can adopt Brownian motion to modelthe diffusion of Tx and X molecules with diffusion coefficients D Tx , and D X , respectively [7]. Weconsider a rooted tumor and thus D Rx “ , which is a special case of the considered general systemmodel.We assume the instantaneous and discrete release of drugs. After releasing for a period, the drugcarrier may be removed by blood circulation or run out of drugs. Thus, for drug delivery systems, T , T 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 continuousrelease can be approximated by letting T b Ñ , i.e., I Ñ 8 . Moreover, A and α i denote the totalnumber of drug molecules released during T and the number of drug molecules released at time t i ,respectively.
2) Molecular Communication System:
For the considered MC system, we assume Brownian motionof the transceivers and signaling molecules. We assume multi-frame communication between mobile Tx and Rx with instantaneous molecule release for each bit transmission. Hence, T , T 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 P t , . . . , I u , denotethe i -th bit in the bit frame. We assume that symbols and are transmitted independently andwith equal probability. Thus, the probability of transmitting ˜ b i is Pr p ˜ b i q “ { , where Pr p¨q denotesprobability and ˜ b i P t , u is a realization of b i . We assume that on-off keying modulation is employed.At time t i , Tx releases α i molecules to transmit bit 1 and no molecules for bit 0. Then, A “ ř Ii “ α i is the total number of molecules available for transmission in a given bit frame. B. Time-variant CIR and Received Signal
Considering again the general system model, we now model the channel between Tx and Rx aswell as the received signals at Rx for drug delivery and MC systems, respectively.
1) Time-variant CIR:
Let h p t, τ q denote the hitting rate, i.e., the absorption rate of a given molecule,at time τ after its release at time t at the center of Tx . Then, for an infinitesimally small observationwindow ∆ τ , i.e., ∆ τ Ñ , we can interpret h p t, τ q ∆ τ as the probability of absorption of a moleculeby Rx between times τ and τ ` ∆ τ after its release at time t . The hitting rate h p t, τ q is also referredto as the CIR since it completely characterizes the time-variant channel, which is assumed to be linear.For a given distance between Tx and Rx , r p t q , the CIR h p t, τ q of a diffusive mobile MC system attime τ is given by [1], [2] h p t, τ q “ a rx ? πD τ ˆ ´ a rx r p t q ˙ exp ˜ ´ p r p t q ´ a rx q D τ ¸ , τ ą , (1)where h p t, τ q “ , for τ ď . Here, D is the effective diffusion coefficient capturing the relativemotion of the signaling molecules and Rx , i.e., D “ D X ` D Rx , see [25, Eq. (8)]. In the consideredMC system, due to the motion of the transceivers, the distance r p t q is a random variable, and thus,the CIR h p t, τ q is time-variant and should be modeled as a stochastic process [12].
2) Received Signal for Drug Delivery System:
In drug delivery, the absorption rate ultimatelydetermines the therapeutic impact of the drug [18], [20]. Thus, we formally define the absorptionrate as the desired received signal, and make achieving a desired absorption rate the objective forsystem design. Recall that h p t, τ q ∆ τ , ∆ τ Ñ , is the probability of absorption of a molecule by Rx between times τ and τ ` ∆ τ after the release at time t . If α i molecules are released at Tx at time t i ,the expected number of molecules absorbed at Rx between times t and t ` ∆ t , for ∆ t Ñ , due to this release is α i h p t i , t ´ t i q ∆ t . During the period r , t s , the total number of released drug moleculesis A t “ ř i α i , @ i | t i ă t, and the expected number of drug molecules absorbed between times t and t ` ∆ t , for ∆ t Ñ is given by s p t q “ ř i α i h p t i , t ´ t i q ∆ t, @ i | t i ă t . Let g p t q denote the absorptionrate of drug molecules X at Rx at time t , i.e., g p t q “ s p t q{ ∆ t , ∆ t Ñ . Then, we have g p t q “ ÿ @ i | t i ă t α i h p t i , t ´ t i q . (2)As mentioned before, the absorption rate g p t q , i.e., the received signal, of the tumor cells directlyaffects the healing efficacy of the drug. Hence, we will design the drug delivery system such that g p t q does not fall below a prescribed value. Since g p t q is a function of h p t i , t ´ t i q , it is random dueto the diffusion of Tx . Therefore, the design of the drug delivery system has to take into accountthe statistical properties of g p t q , which can be obtained from the results of the statistical analysis of h p t i , t ´ t i q .
3) Received Signal for MC System:
For the MC system design, the received signal, denoted by q i ,is defined as the number of X molecules absorbed at Rx during bit interval T b after the transmissionof the i -th bit at t i by Tx as the received signal, denoted by q i . We detect the transmitted informationbased on the received signal, q i . It has been shown in [6] that q i follows a Binomial distribution thatcan be accurately approximated by a Gaussian distribution when α i is large, which we assume here.We focus on the effect of the transceivers’ movements on the MC system performance and designthe optimal release profile of Tx to account for these movements. We assume the bit interval to besufficiently long such that most of the molecules have been captured by or have moved far away from Rx before the following bit is transmitted, i.e., ISI is negligible. We note that enzymes [6] and reactiveinformation molecules, such as acid/base molecules [26], [27], may be used to speed up the moleculeremoval process and to increase the accuracy of the ISI-free assumption. Moreover, we model externalnoise sources in the environment as Gaussian background noise with mean and variance equal to η [1]. Thus, we have q i „ N ´ µ i, ˜ b i , σ i, ˜ b i ¯ for b i “ ˜ b i , (3)where µ i, “ σ i, “ η , µ i, “ α i p p t i , T b q ` η , σ i, “ α i p p t i , T b qp ´ p p t i , T b qq ` η . Here, N p µ, σ q denotes a Gaussian distribution with mean µ and variance σ . p p t, T b q denotes the probability that asignaling molecule is absorbed during bit interval T b after its release at time t at the center of Tx .For a given distance r p t q , p p t, T b q is given by [2] p p t, T b q “ ż T b h p t, τ q d τ “ a rx r p t q erfc ˆ r p t q ´ a rx ? D T b ˙ , (4) where erfc p¨q is the complementary error function. Since r p t q is a random variable and p p t, T b q is afunction of r p t q , p p t, T b q and any function of p p t, T b q , e.g., the received signal q i , are random processes.Moreover, p p t, T b q is also a function of h p t, τ q . Hence, for MC system design, we have to take intoaccount the statistical properties of p p t, T b q , which can be obtained based on the proposed statisticalanalysis of r p t q and h p t, τ q .In summary, the design of both drug delivery and MC systems depends on the statistical propertiesof the CIR, h p t, τ q , and r p t q including their means, variances, PDFs, and CDFs, which will be analyzedin the next section. III. S TOCHASTIC C HANNEL A NALYSIS
In this section, we first analyze the distribution of the distance between the transceivers, r p t q , andthen use it to derive the statistics of the time-variant CIR, h p t, τ q , and p p t, T b q as a function h p t, τ q .In particular, we develop analytical expressions for the mean, variance, PDF, and CDF of h p t, τ q andthe PDF and CDF of p p t, T b q . A. Distribution of the Tx - Rx Distance for a Diffusive System
In the 3D space, r p t q is given by r p t q “ bř d Pt x,y,z u p r d, Rx p t q ´ r d, Tx p t qq , where r d, Tx p t q and r d, Rx p t q , d P t x, y, z u , are the Cartesian coordinates representing the positions of Tx and Rx attime t , respectively. Let us assume, without loss of generality, that the diffusion of Tx and Rx starts at t “ . Then, given the Brownian motion model for the mobility of Tx and Rx , we have r d, Tx p t q „ N p r d, Tx p t “ q , D Tx t q and r d, Rx p t q „ N p r d, Rx p t “ q , D Rx t q , where we assume that r d, Tx p t “ q and r d, Rx p t “ q are known. Let us define r d p t q “ r d, Rx p t q ´ r d, Tx p t q . Then, we have r d p t q „ N p r d p t “ q , D t q , where D “ D Tx ` D Rx is the effective diffusion coefficient capturingthe relative motion of Tx and Rx , see [25, Eq. (10)]. Given the Gaussian distribution of r d p t q , weknow that [28] γ “ r p t q? D t “ d ř d Pt x,y,z u r d p t q D t , (5)follows a noncentral chi-distribution, i.e., γ „ X k p λ q , with k “ degrees of freedom and parameter λ “ b ř d Pt x,y,z u r d p t “ q D t “ r ? D t , where r denotes r p t “ q . The statistical properties of randomvariable r p t q are provided in the following lemma. Lemma 1:
The mean, variance, PDF, and CDF of random variable r p t q , which represents the distancebetween the centers of the diffusive mobile Tx and Rx , are given by, respectively, E t r p t qu “ c D tπ exp ˆ ´ r D t ˙ ` ˆ r ` D tr ˙ erf ˆ r ? D t ˙ , (6) Var t r p t qu “ r ` D t ´ E t r p t qu , (7) f r p t q p r q “ rr ? πD t exp ˆ ´ r ` r D t ˙ sinh ˆ r r D t ˙ , (8)and F r p t q p r q “ ´ Q ˆ λ, r ? D Tx t ˙ . (9)where erf p¨q is the error function, Q M p a, b q is the Marcum Q-function [29], E t¨u denotes statisticalexpectation, Var t¨u denotes variance, and f t¨u p¨q and F t¨u p¨q denote the PDF and CDF of the randomvariable in the subscript, respectively. Proof:
Please refer to Appendix A.
Remark 1:
From (6) and (7), we can observe that when t Ñ 8 , we have exp ´ ´ r D t ¯ Ñ and erf ´ r ? D t ¯ Ñ and, as a result, E t r p t qu Ñ 8 . Intuitively, because of diffusion, the transceiverseventually move far away from each other on average. Remark 2:
We note that (8) was derived under the assumption that Tx can diffuse in the entire 3Denvironment. However, in reality, Tx cannot move inside Rx , i.e., it does not interact with Rx , andthus will be reflected when it hits Rx ’s boundary. Hence, the actual f r p t q p r q , derived in [25], differsfrom (8), e.g., f r p t q p r q “ for r ă a tx ` a rx . However, for very small r , i.e., r « , (8) approacheszero. Hence, (8) is a valid approximation for the actual f r p t q p r q . The validity of this approximation isevaluated in Section VI via simulations, where, in our particle-based simulation, Tx is reflected uponcollision with Rx [30]. B. Statistical Moments of Time-variant CIR
In this subsection, we derive the statistical moments of the time-variant CIR, i.e., mean m p t, τ q andvariance σ p t, τ q . In particular, the mean of the time-variant CIR, m p t, τ q , can be written as m p t, τ q “ ż h p t, τ q ˇˇ r p t q“ r f r p t q p r q d r. (10)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 moleculesreleased by a diffusive transparent transmitter and captured by a diffusive absorbing receiver is givenby m p t, τ q “ a rx a π p D τ ` D t q r τ exp ˆ ´ a D τ ´ r D t ˙ „ ´ e v p t,τ q u p t,τ q ˆ v p t, τ q u p t, τ q ` a rx ˙ (11) ˆ erfc ˜ v p t, τ q a u p t, τ q ¸ ` e w p t,τ q u p t,τ q ˆ w p t, τ q u p t, τ q ` a rx ˙ erfc ˜ w p t, τ q a u p t, τ q ¸ff , where u p t, τ q , v p t, τ q , and w p t, τ q are defined, for compactness, as follows u p t, τ q “ D τ ` D t , v p t, τ q “ ´ a rx D τ ´ r D t , w p t, τ q “ ´ a rx D τ ` r D t . (12) Proof:
Substituting (1) and (8) into (10) and using the integrals given by [31, Eq. (2.3.15.4) andEq. (2.3.15.7)], we obtain the expression for m p t, τ q in (11). Remark 3: m p t, τ q is a function of time t . Hence, h p t, τ q is a non-stationary stochastic process. Ingeneral, at large t , m p t, τ q decreases when t increases and eventually approaches zero when t Ñ 8 .This means that as t increases, the molecules released by Tx , on average, have a decreasing chance ofbeing absorbed by Rx since the transceivers move away from each other as mentioned in Remark 1.In order to obtain the variance of h p t, τ q , σ p t, τ q “ φ p t, τ q ´ m p t, τ q , (13)we first need to find an expression for the second moment φ p t, τ q , defined as φ p t, τ q “ E t h p t, τ qu .The following corollary provides an analytical expression for φ p t, τ q . Corollary 1: φ p t, τ q is given by φ p t, τ q “ c p t, τ q ż ` exp ` ´ ˆ u p t, τ q r ´ ˆ v p t, τ q r ˘ ´ exp ` ´ ˆ u p t, τ q r ´ ˆ w p t, τ q r ˘˘ (14) ˆ ˆ r ´ a rx ` a r ˙ d r , where c p t, τ q “ a e ´ a D τ ´ r D t D πτ r ? πD t , ˆ u p t, τ q “ D τ ` D t , (15) ˆ v p t, τ q “ ´ a rx D τ ´ r D t , ˆ w p t, τ q “ ´ a rx D τ ` r D t . Proof:
From the definition, we have φ p t, τ q “ E (cid:32) h p t, τ q ( “ ż h p t, τ q ˇˇ r p t q“ r f r p t q p r q d r . (16)Substituting (1) and (8) into (16) and simplifying the expression, we obtain (14). Remark 4:
The expression in (14) comprises integrals of the form ş exp p ax ` bx q { x d x , where a and b are constants. Such integrals cannot be obtained in closed form. However, the integrals canbe evaluated numerically in a straightforward manner. C. Distribution Functions of the Time-variant CIR
In this subsection, we derive analytical expressions for the PDF and CDF of h p t, τ q . The PDF of h p t, τ q is given in the following theorem. Theorem 2:
The PDF of the impulse response of a time-variant channel with diffusive moleculesreleased by a diffusive transparent transmitter and captured by a diffusive absorbing receiver is givenby $’’’’&’’’’% f h p t,τ q p h q “ f r p t q p r p h qq ˆ h p r p h q ,τ q ´ f r p t q p r p h qq ˆ h p r p h q ,τ q , for ď h ă h ‹ ,f h p t,τ q p h q Ñ 8 , for h “ h ‹ ,f h p t,τ q p h q “ , otherwise, (17)where ˆ h p r, τ q denotes h p t, τ q , given by (1), as a function of r p t q and τ , f r p t q p r q is given by (8), r p h q and r p h q , r p h q ă r p h q , are the solutions of the equation ˆ h p r, τ q “ h , h ‹ is the maximum value of ˆ h p r, τ q for all values of r p t q , and ˆ h p r, τ q is given by ˆ h p r, τ q “ a rx ? πD τ exp ˜ ´ p r ´ a rx q D τ ¸ ˆ a rx r ´ p r ´ a rx q D τ ´ ´ a rx r ¯˙ . (18) Proof:
Please refer to Appendix B.As stated in the proof of Theorem 2, there are two different values of r p t q , r and r , leading to thesame value of ˆ h p r, τ q , i.e., h p t, τ q , when ď h ă h ‹ . Hence, the PDF of h p t, τ q is a function of the PDFsof these two values of r p t q . However, when h p t, τ q reaches its maximum, f h p t,τ q p h q approaches infinityand does not depend on f r p t q p r p h qq since the probability of h “ h ‹ , i.e., Pr p h “ h ‹ q “ f h p t,τ q p h q d h ,is finite and d h approaches at h “ h ‹ .The CDF of h p t, τ q is given in the following corollary. Corollary 2:
The CDF of the impulse response of a time-variant channel with diffusive moleculesreleased by a diffusive transparent transmitter and captured by a diffusive absorbing receiver is givenby $’’’’&’’’’% F h p t,τ q p h q “ F r p t q p r p h qq ` ´ F r p t q p r p h qq , for ď h ď h ‹ ,F h p t,τ q p h q “ , for h ă ,F h p t,τ q p h q “ , for h ą h ‹ , (19) where F r p t q p r q is given by (9). Proof:
From the definition of the CDF and (17), we have F h p t,τ q p h q “ ż h f h p t,τ q p ˇ h q dˇ h “ ż h f r p t q p ˇ r p ˇ h qqB ˆ h p ˇ r , τ q{B ˇ r ´ f r p t q p ˇ r p ˇ h qqB ˆ h p ˇ r , τ q{B ˇ r dˇ h (20) “ ż r p h q f r p t q p ˇ r q dˇ r ´ ż r p h q8 f r p t q p ˇ r q dˇ r “ F r p t q p r p h qq ` ´ F r p t q p r p h qq , where ˇ r and ˇ r , ˇ r ă ˇ r , are the solutions of the equation ˆ h p ˇ r, τ q “ ˇ h . This completes the proof.Similar to the PDF, the CDF of h p t, τ q also depends on the CDFs of two values of r p t q , i.e., r p h q and r p h q . D. Distribution Functions of p p t, T b q Calculating the mean of p p t, T b q involves an integral of the form ş erfc p ax q exp p´ b x ` cx q d x ,with appropriate constants a, b, c ą , 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 p t, T b q in thefollowing corollary. Corollary 3:
The PDF and CDF of the probability that a diffusive molecule is absorbed by adiffusive absorbing receiver during an interval T b after its release at time t by a diffusive transparenttransmitter are, respectively, given by f p p t,T b q p p q “ ´ f r p t q p ˜ r p p qq p p ˜ r q , (21) F p p t,T b q p p q “ ´ F r p t q p ˜ r p p qq , (22)where f r p t q p r q and F r p t q p r q are given by (8) and (9), respectively. Here, ˜ r p p q is the solution of theequation p p t, T b q “ p and p p ˜ r q is given by p p ˜ r q “ ´ a rx ˜ r erfc ˆ ˜ r ´ a rx ? D T b ˙ ´ a rx ˜ r ? πD T b exp ˆ ´ p ˜ r ´ a rx q D T b ˙ . (23) Proof:
The proof of Corollary 3 follows the same steps as the proof of Theorem 2 and Corollary 2and exploits that p p t, T b q is a function of r p t q as shown in (4). From (23), we observe that p p ˜ r q ă so the equation p p t, T b q “ p has only one solution. Then, we apply the relations for the PDFs andCDFs of functions of random variables [32] to obtain (21) and (22).The mean, variance, PDF, and CDF of h p t, τ q and p p t, T b q can be exploited to design efficient andreliable synthetic MC systems. As examples, we consider the design and analysis of drug deliveryand MC systems in the following two sections. IV. D
RUG D ELIVERY S YSTEM D ESIGN
In this section, we apply the derived stochastic parameters of the time-variant CIR for absorbingreceivers for the design and performance evaluation of drug delivery systems.
A. Controlled-Release Design
The treatment of many diseases requires the diseased cells to absorb a minimum rate of drugsduring a prescribed time period at minimum cost [20]. To design an efficient drug delivery systemsatisfying this requirement, we minimize the total number of released drug molecules, A “ ř Ii “ α i ,subject to the constraint that the absorption rate g p t q is equal to or larger than a target rate, θ p t q ,for a period of time, denoted by T Rx . We allow θ p t q to be a function of time so that the designedsystem can satisfy different treatment requirements over time. Since g p t q is random, we cannot alwaysguarantee g p t q ě θ p t q . Hence, we will design the system based on the first and second order momentsof g p t q and use the PDF and CDF of g p t q to evaluate the system performance. In particular, wereformulate the constraint such that the mean of g p t q minus a certain deviation is equal to or abovethe threshold θ p t q during T Rx , i.e., E t g p t qu ´ β V t g p t qu ě θ p t q , ď t ď T Rx , where V t¨u denotesstandard deviation and β is a coefficient determining how much deviation from the mean is taken intoaccount. Based on (2), the constraint can be written as a function of α i as follows E t g p t qu ´ β V t g p t qu p a q ě ÿ @ i | t i ă t α i p E t h p t i , t ´ t i qu ´ β V t h p t i , t ´ t i quq ě θ p t q , (24)for ď t ď T Rx . Inequality p a q in (24) is due to E t g p t qu “ E " ř @ i | t i ă t α i h p t i , t ´ t i q * “ ř @ i | t i ă t α i ˆ E t h p t i , t ´ t i qu and Minkowski’s inequality [33]: V t g p t qu “ »– E $&%¨˝ ÿ @ i | t i ă t p α i h p t i , t ´ t i q ´ E t α i h p t i , t ´ t i quq ˛‚ ,.-fifl { (25) ď ÿ @ i | t i ă t α i “ E (cid:32) p h p t i , t ´ t i q ´ E t h p t i , t ´ t i quq (‰ { “ ÿ @ i | t i ă t α i V t h p t i , t ´ t i qu . Note that we may not be able to find α i such that (24) holds for all values of β and θ p t q . Forexample, when β is too large, E t g p t qu ´ β V t g p t qu can be negative and hence (24) cannot be satisfiedfor θ p t q ą . However, when E t h p t i , t ´ t i qu ą β V t h p t i , t ´ t i qu , i.e., either β or V t h p t i , t ´ t i qu are small, such that β V t h p t i , t ´ t i qu is sufficiently small, we can always find α i so that (24) holdsfor arbitrary θ p t q . Since time t is a continuous variable, the constraint in (24) has to be satisfiedfor all values of t , ď t ď T 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 holdonly for a finite number of time instants t “ t n “ n ∆ t n , where n “ , . . . , N and ∆ t n “ T Rx { N .Then, the proposed optimization problem for the design of α i is formulated as follows min α i ě , @ i A “ I ÿ i “ α i (26)s.t. ÿ i,t i ă t α i p m p t i , n ∆ t n ´ t i q ´ βσ p t i , n ∆ t n ´ t i qq ě θ p n ∆ t n q , for n “ , . . . , N, where m p t, τ q and σ p t, τ q are the mean (11) and the standard deviation (13) of h p t, τ q , respectively.Since m p t, τ q and σ p t, τ q do not oscillate but are well-behaved and smooth functions of t as shown inSection VI, a small value of N (e.g., N “ ) is usually enough to meet the continuous constraint (24)for all t . Having m p t, τ q in (11) and σ p t, τ q in (13) and treating the α i as real numbers, (26) can bereadily solved numerically as a linear program. We note that although the numbers of drug molecules α i are integers, for tractability, we solve (26) for real α i and quantize the results to the nearest integervalues.We note that the problem in (26) is statistical in nature and provides guidance for the offline designof the drug delivery system. B. System Performance
Since g p t q ě θ p t q is required for proper operation of the system, we evaluate the system performancein terms of the probability that the drug absorption rate satisfies the target rate g p t q ě θ p t q , denotedby P θ p t q “ Pr t g p t q ě θ p t qu . P θ p t q is given in the following theorem. Theorem 3:
The system performance metric P θ p t q can be expressed as P θ p t q “ ´ f α h p t ,t ´ t q p θ p t qq ˚ ¨ ¨ ¨ ˚ f α i ´ h p t i ´ ,t ´ t i ´ q p θ p t qq ˚ F α i h p t i ,t ´ t i q p θ p t qq , (27)where ˚ denotes convolution, and i “ , , . . . satisfies t i ď t . In (27), we define f α i h p t i ,t ´ t i q p θ p t qq “ { α i ˆ f h p t i ,t ´ t i q p θ p t q{ α i q and F α i h p t i ,t ´ t i q p θ p t qq “ F h p t i ,t ´ t i q p θ p t q{ α i q . Proof:
From the definition of the CDF, we have P θ p t q “ ´ F g p t q t θ p t qu “ ´ ż θ p t q f g p t q p g q d g. (28)Due to the summation of independent random variables in (2), i.e., independent releases at t i , we have f g p t q p g q “ ` f α h p t ,t ´ t q ˚ ¨ ¨ ¨ ˚ f α i h p t i ,t ´ t i q ˘ p g q . (29) Substituting (29) into (28), then using the integration property of the convolution, i.e., ż θ p t q ` f α h p t ,t ´ t q ˚ ¨ ¨ ¨ ˚ f α i h p t i ,t ´ t i q ˘ p g q d g “ f α h p t ,t ´ t q p θ p t qq ˚ ¨ ¨ ¨ ˚ ż θ p t q f α i h p t i ,t ´ t i q p g q d g, (30)and using the definition of the CDF, we obtain (27).We note that the analytical expressions for the PDF and CDF of h p t, τ q 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 θ p t q can be guaranteed based on the statisticalmoments 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 θ p t q “ Pr t g p t q ě θ p t qu is given asfollows P θ p t q ě ´ β . (31) Proof:
By using (24) and the Chebyshev inequality [34], we obtain P θ p t q p a q ě Pr ! | g p t q ´ E t g p t qu| ď E t g p t qu ´ θ p t q ) p b q ě ´ V t g p t qup E t g p t qu ´ θ p t qq p c q ě ´ β , (32)where p a q can be obtained by expanding the absolute value on the right-hand side, p b q is due to theChebyshev inequality, and p c q is due to (24). This completes the proof. Remark 5:
Proposition 1 is not only useful for evaluating the system performance, but also providesa guideline for the design of the release profile of drugs in (26). For example, to ensure a highabsorption rate probability of P θ p t q ě . , from (31), we need to set the β coefficient in (26) as β “ . Note that a useful bound can only be obtained based on (31) when β ą and (24) is satisfied.V. MC S YSTEM D ESIGN FOR I MPERFECT
CSIIn this section, we apply the stochastic analysis presented in Section III for the design of MC systemswith imperfect CSI. The CSI is imperfect due to the movement of the transceivers and assumed to beknown only at the beginning of a bit frame. In particular, we optimize three design parameters of adiffusive mobile MC system employing on-off keying modulation and threshold detection, namely thedetection threshold at Rx , the release profile at Tx , and the time duration of a bit frame. By choosingthe optimal values of those three parameters, we can improve the system performance while keepingthe overall system relatively simple. First, we optimize the detection threshold for minimization ofthe maximum BER in a frame assuming a uniform release profile. This approach can be employedin very simple MC systems where Tx is not capable of adjusting the number of released molecules. Second, we optimize the release profile at 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 theframe. This second approach to MC optimization improves the system performance in terms of BERbut requires a mechanism to control the number of molecules released at Tx . Third, we design theoptimal duration of the bit frame satisfying a constraint on the efficiency of molecule usage. Thus, thisthird approach improves the system performance in terms of the efficiency of molecule usage. Thethree proposed designs can be performed offline. Furthermore, they can be combined with each otheror carried out separately depending on the capabilities and requirements of the system. For all threedesigns, as a first step, we need to derive the BER as a function of the number of released molecules. A. Detection and BER
We consider a simple threshold detector at Rx , where the received signal q i is compared with adetection threshold, denoted by ξ , in order to determine the detected bit ˆ b i as follows ˆ b i “ $’&’% if q i ą ξ, if q i ď ξ. (33)Given the assumption of no ISI and Pr p ˜ b i q “ { , from (3) and (33), the error probability of the i -th bit, denoted by P b p b i q , can be simplified as [25, Eq. (12)] P b p b i q “ ´
14 erf ˆ ξ ´ η ? η ˙ ` ż f r p t q p r i q erf p ζ i p ξ, α i qq d r i , (34)where f r p t q p r i q is given in (8), r i is r p t i q for brevity, and ζ i p ξ, α i q “ ξ ´ µ i, σ i, ? “ ξ ´p α i p p t i ,T b q` η q ? p α i p p t i ,T b qp ´ p p t i ,T b qq` η q .Note that P b p b i q depends on i since the distance r p t i q between the transceivers is a function of releasetime t i . B. Optimal Detection Threshold for Uniform Release
We first consider system design for uniform release, where the number of available molecules isuniformly allocated across all bits of a frame. To facilitate reliable communication, our objective is tooptimize the detection threshold, ξ , such that the maximum value of the error rate of the bits in aframe is minimized, given the total number of available molecules in a frame, A , i.e., min ξ max i t P b p b i qu s . t . α i “ A { I. (35)From (34), the problem is equivalent to min ξ max i "ż f r p t q p r i q erf p ζ i p ξ, α i qq d r i ´ erf ˆ ξ ´ η ? η ˙* s . t . α i “ A { I. (36) The following lemma reveals the convexity of the problem in (36).
Lemma 2:
For η ă ξ ă µ i, , the objective function in (36) is convex in ξ . Proof:
Please refer to Appendix C.Note that η ă ξ ă µ i, is intuitively satisfied for typical system parameters since the decision thresholdshould be higher than the average noise level when bit ” ” is sent but should not exceed the meanof the received signal when bit ” ” is sent. Otherwise, a high error rate would result. Due to theconvexity of problem (36), the global optimum ξ can be easily obtained by numerical methods suchas the interior-point method [32]. C. Optimal Release with Fixed Detection Threshold
For the second proposed design, we aim to optimize the release profile, i.e., the number of moleculesavailable for release for each bit, α i , such that max i t P b p b i qu is minimized given a total number ofmolecules A that are available for release in a frame min α max i t P b p b i qu s . t . I ÿ i “ α i “ A, (37)where α “ r α i , α , . . . , α I s .For a given threshold ξ , we can re-express (37) based on (34) as min α max i "ż f r p t q p r i q erf p ζ i p ξ, α i qq d r i * s . t . I ÿ i “ α i “ A. (38)The following lemma states the convexity of the optimization problem in (38). Lemma 3:
For η ă ξ ă µ i, , the objective function in (38) is convex in α . Proof:
Please refer to Appendix D.Hence, the global optimum of (38) can be readily obtained by numerical methods such as theinterior-point method.Note that, for tractability, similar to the proposed drug delivery design, we solve (36) and (38) forreal α i and quantize the results to the nearest integer values. D. Optimal Time Duration of a Bit Frame
In the third proposed design, we consider the molecule usage efficiency for communication. Weevaluate the efficiency based on p p t, T b q , i.e., the probability that a signaling molecule is absorbedduring bit interval T b after its release at time t . If p p t, T b q is too small, none of the released moleculesmay actually reach the receiver and thus the molecules are wasted, i.e., the system has low efficiency.Hence, we want to keep p p t, T b q above a certain value, denoted by ψ . Intuitively, as on average h p t, τ q decreases over time t , p p t, T b q , which is the integral over h p t, τ q with respect to τ , also on averagedecreases over time. Therefore, our objective is to choose the maximum duration of a bit frame,denoted by T ‹ , such that p p t, T b q ą ψ for t ď T ‹ ´ T b , where t is the release time.Since p p t, T b q is a random process, we cannot enforce p p t, T b q ą ψ but can only bound theprobability that p p t, T b q ą ψ is satisfied, i.e., Pr p p p t, T b q ą ψ q ě P , where P is a design parameter.Moreover, we have Pr p p p t, T b q ą ψ q “ ´ F p p t,T b q p ψ q p a q “ F r p t q p ˜ r p ψ qq , (39)where equality p a q is due to (22). As such, we can re-express the problem as maximizing the durationof a bit frame such that F r p t q p ˜ r p ψ qq ě P holds. To this end, in the following lemma, we analyze F r p t q p ˜ r p ψ qq as a function of time t . Lemma 4: F r p t q p ˜ r p ψ qq is a decreasing function of time t . Proof:
Please refer to Appendix E.Since Lemma 4 shows that F r p t q p ˜ r p ψ qq is a decreasing function of time, the maximum durationof a bit frame satisfying F r p t q p ˜ r p ψ qq ě P can be found by solving F r p T ‹ ´ T b q p ˜ r p ψ qq “ P , where F r p T ‹ ´ T b q p ˜ r p ψ qq is given in (9). Remark 6:
If multiple frames are transmitted, the proposed design framework can be applied toeach frame, respectively. However, the optimal designs may be different for different frames due to themoving 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 optimaldesign framework can also be applied to transparent receivers. For a transparent receiver, p p t, T b q isthe probability that a molecule is observed inside the volume of the transparent receiver at time T b after its release at time t at the center of Tx .VI. N UMERICAL R ESULTS
In this section, we provide numerical results to evaluate the accuracy of the derived expressions andanalyze the performance of the MC systems in the considered application scenarios. We use the setof simulation parameters summarized in Table I, unless stated otherwise. The parameters are chosento match the actual system parameters in drug delivery systems, as will be explained in detail inSubsection VI-B.
A. Time-variant Channel Analysis
In this subsection, we numerically analyze the time-variant MC channel. For verification of theaccuracy of the expressions derived in Section III, we employ a hybrid particle-based simulation TABLE IS
YSTEM PARAMETERS USED FOR NUMERICAL RESULTS
Parameter Value Parameter Value D Tx [ m { s ] ˆ ´ D Rx [ m { s ] D X [ m { s ] ˆ ´ r [ m ] ˆ ´ a tx [ m ] ˆ ´ a rx [ m ] ˆ ´ T [ h ] T Rx [ h ] I N T b [ s ] . θ p t qr s ´ s . . . . . . . . . . t = { , , , , } s τ [s] m ( t , τ ) AnalysisSimulation
Fig. 2. Mean of the CIR h p r p t q , τ q as a function of time τ . .
05 0 . .
15 0 . .
25 0 . t = { , , } s τ = 0 . h ( t, τ ) f h ( t , τ ) ( h ) AnalysisSimulation
Fig. 3. PDF of the CIR f h p t,τ q p h q for τ “ .
17 s and t “t , , u s . approach. In particular, we use particle-based simulation of the Brownian motion of the transceiversto generate realizations of the random distance between Tx and Rx , r p t q . Then, we use Monte Carlosimulation to obtain the desired statistical results by suitably averaging the CIRs (1) obtained forthe realizations of r p t q . For particle-based simulation of the Brownian motion of Tx , Tx performs arandom walk with a random step size in space in every discrete time step of length ∆ t st “ . Thelength of each step in space is modeled as a Gaussian random variable with zero mean and standarddeviation ? D Tx ∆ t st . Furthermore, we also take into account the reflection of Tx upon collision with Rx . When Tx hits Rx , we assume that it bounces back to the position it had at the beginning of theconsidered simulation step [30].Fig. 2 shows the mean of the CIR, m p t, τ q , as a function of τ . In general, for large τ , m p t, τ q decreases when t increases as expected since the transceivers move away from each other on average.For large τ , m p t, τ q also decreases when τ increases as would be the case in a static system. Notethat in the simulations, unlike the analysis, we have taken into account the reflection of Tx when ithits Rx . Therefore, the good agreement between simulation and analytical results in Fig. 2 suggeststhat the reflection of Tx does not have a significant impact on the statistical properties of h p t, τ q andthe 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 “ t , , u s and τ “ .
17 s . Fig. 3shows that when t increases, a smaller value of h p t, τ q is more likely to occur since, on average, thetransceivers move away from each other. When t is very large, it is likely that the molecules cannotreach Rx , and hence, cannot be absorbed, consequently h p t, τ q Ñ . We also observe that h p t, τ q hasa maximum value and f h p t,τ q p h p t, τ qq Ñ 8 when h p t, τ q approaches the maximum value as statedin (17). For example, h p t “
36 s , τ “ .
17 s q is random but its maximum possible value is . and f h p t “
36 s ,τ “ .
17 s q p . q Ñ 8 .Fig. 2 and 3 show a perfect match between simulation and analytical results. This confirms theaccuracy of our analysis of the time-variant CIR in Sections III. Since particle-based simulation iscostly, in the following subsections, we adopt Monte-Carlo simulation by averaging our results over independent realizations of both the distance r p t q and the CIR. The distance r p t q is calculated fromthe locations of the transceivers, which are generated from Gaussian distributions, see Subsection III-A.In particular, r p t q “ bř d Pt x,y,z u p r d, Rx p t q ´ r d, Tx p t qq , where r d, Tx p t q „ N p r d, Tx p t “ q , D Tx t q , r d, Rx p t q „ N p r d, Rx p t “ q , D Rx t q , r d, Tx p t “ q “ , and r d, Rx p t “ q “ {? r . The CIR is givenby (1) for each realization of r p t q . B. Drug Delivery System Design
In this section, we provide numerical results for the considered drug delivery system. As mentionedabove, the parameters in Table I are chosen to match real system parameters, e.g., the diffusioncoefficient D X of drug molecules vary from ´ to ´ m { s [19], drug carriers have sizes a tx ě
100 nm [17], the size of tumor cells is on the order of µ m , and drug carriers can be injected orextravasated from the cardiovascular system into the tissue surrounding the targeted diseased cell site[18], i.e., close to the tumor cells. The dosing periods in drug delivery systems are on the order ofdays [35], i.e.,
24 h . For simplicity, we set N “ and the value of the required absorption rate is setto θ p t q “ ´ . We choose I relatively large to obtain small intervals T b .In Fig. 4, we plot the number of released molecules α i versus the corresponding release time t i [ h ]for different system parameters. The coefficients are obtained by solving the optimization problem , , , , , D Tx = 10 − m / s D Tx = 5 × − m / s D Tx = 10 − m / s D Tx = 0m / s BenchmarkRelease time t i [ h ] α i β = 0 β = 0 . β = 1 β = 2 Fig. 4. Optimal number of released molecules α i as a function ofrelease time t i [ h ] for different system parameters and T “
24 h .The black horizontal dotted line is the benchmark when the α i arenot optimized. .
992 7 .
994 7 .
996 7 .
998 8 8 .
002 8 .
004 8 .
006 8 . . . . . design 1design 2 design 3Time t [ h ] E { g ( t ) } a nd V { g ( t ) } E { g ( t ) } analysis E { g ( t ) } simulation V { g ( t ) } analysis V { g ( t ) } simulation θ ( t ) Fig. 5. E t g p t qu and V t g p t qu between the 1000-th release and the1002-th release, i.e., at about , for three different designs. Design1 (green line): naive design without considering Tx ’s movementwith D Tx “ ´ m { s and β “ ; design 2 (blue line) and 3(red line): optimal design for ` D Tx r m { s s , β ˘ “ ` ´ , ˘ , and ` ´ , ˘ , respectively. in (26) with β “ t , . , , u for D Tx “ ´ m { s and β “ t , . u for D Tx “ t , u ˆ ´ m { s . As mentioned in the discussion of (24), we cannot choose large values of β when the diffusioncoefficient is large, i.e., the standard deviation is large, as the problem in (26) may become infeasible.Fig. 4 shows that for all considered parameter settings, we first have to release a large number ofmolecules for the absorption rate to exceed the threshold. Then, in the static system with D Tx “ { s ,the optimal coefficient decreases with increasing time, since a fraction of the molecules previouslyreleased from Tx linger around Rx and are absorbed later. However, for the time-variant channel, Tx eventually diffuses away from Rx as time t increases and hence, molecules released at later times by Tx will be far away from Rx and may not reach it. Therefore, at later times, the amount of drugsreleased has to be increased for the absorption rate to not fall below the threshold. For larger D Tx , Tx diffuses away from Rx faster and thus, the number of released molecules α i have to increase faster.This type of drug release, i.e., first releasing a large amount of drugs, then reducing and eventuallyincreasing the amount of released drugs again, is called a tri-phasic release [36]. Once we havedesigned the release profile, we can implement it by choosing a suitable drug carrier as shown in [36].Moreover, as expected, for larger β , we need to release more drugs to ensure that (26) is feasible. Theblack horizontal dotted line in Fig. 4 is a benchmark where the α i , @ i, are not optimized but naively .
992 7 .
994 7 .
996 7 .
998 8 8 .
002 8 .
004 8 .
006 8 . . . . . D Tx = 10 − m / s D Tx = 5 × − m / s D Tx = 10 − m / s D Tx = 10 − m / s (naive design)Lower bound D Tx = 10 − m / s , β = 2 Time t [ h ] P θ β = 0 β = 0 . β = 1 β = 2 Analysis
Fig. 6. P θ p t q as a function of time t [ h ] between the 1000-threlease and the 1002-th release, i.e., at about . A = 10 A = 10 A = 10 Bit index i N u m b e r o fr e l ea s e d m o l ec u l e s α i Uniform releaseOptimal release
Fig. 7. The number of molecules available for release for each bitfor uniform and optimal release when A “ (cid:32) , , ( and T “
300 s . set to α i “ α “ . For this naive design, A “ α I « . ˆ , whereas with the optimal α i ,for β “ and D Tx “ ´ m { s , A “ . ˆ , i.e., less than the A required for the naivedesign, and for β “ and D Tx “ ´ m { s , A “ . ˆ , i.e., less than the A required forthe naive design. This highlights that applying the optimal release profile can save significant amountsof drugs and still satisfy the therapeutic requirements. Moreover, as observed in Fig. 4, at later times,e.g., t i ą
15 h for D Tx “ ´ m { s , the values of α i required to satisfy the desired absorption rateare higher than the fixed α i used in the naive design, i.e., the benchmark, which means that the naivedesign cannot provide the required absorption rate.In Fig. 5, we plot the mean and standard deviation of the absorption rate, E t g p t qu and V t g p t qu ,between the -th release and the -th release for three designs. For designs 1 and 2, we assumed D Tx “ ´ m { s and β “ , and for design 3, we adopted D Tx “ ´ m { s and β “ . Notethat the considered time window, e.g., between the -th release and the -th release, is chosenarbitrarily in the middle of T to analyze the system behavior between individual releases. For design1, Tx diffuses with D Tx “ ´ m { s but the release profile is designed without accounting for Tx ’smobility, i.e., the adopted α i are given by the green line in Fig. 4 obtained under the assumptionof D Tx “ { s . For designs 2 and 3, the mobility of Tx is taken into account. The black dashedline marks the threshold θ p t q that g p t q should not fall below. It is observed from Fig. 5 that when Tx diffuses but the design does not take into account the mobility, the requirement that the expected absorption rate, E t g p t qu , exceeds θ p t q , is not satisfied for most of the time. For design 2 with β “ ,we observe that E t g p t qu ą θ p t q always holds but E t g p t qu ´ V t g p t qu ą θ p t q does not always hold. Fordesign 3 with β “ , we observe that E t g p t qu ´ V t g p t qu ą θ p t q always holds since β ą enforces agap between E t g p t qu and θ p t q . In other words, even if g p t q deviates from the mean, it can still exceed θ p t q .In Fig. 6, we present the system performance in terms of the probability that g p t q ě θ p t q , P θ p t q , forthe time period between the -th and -th releases, i.e., at about . The lines and markersdenote simulation and analytical results, respectively. Fig. 6 shows a good agreement between analyticaland simulation results. In Fig. 6, we observe that P θ p t q increases with increasing β because the designfor larger β enforces a larger gap between E t g p t qu and θ p t q , as can be seen in Fig. 5. Moreover, fora given β , P θ p t q will be different for different D Tx . In particular, for larger D Tx , P θ p t q is smallerdue to the faster diffusion and increasing randomness of the CIR. Moreover, in Fig. 6, the green lineshows that the naive design, i.e., design 1 in Fig. 5, has very poor performance. In Fig. 6, we alsoobserve that between two releases, P θ p t q first increases due to the released drugs and then decreasesdue to drug diffusion. Furthermore, in Fig. 6, we also show the lower bound on P θ p t q derived inProposition 1 for D Tx “ ´ m { s and β “ , where (31) yields P θ p t q ě . . Fig. 6 shows thatthe red dash-dotted line, i.e., P θ p t q for D Tx “ ´ m { s and β “ , is indeed above the horizontalblack dashed line, i.e., P θ p t q “ . . C. Molecular Communication System Design
In this subsection, we show numerical results for the second application scenario, i.e., an MCsystem with imperfect CSI. We apply again the system parameters in Table I except that here we set D Rx “ ´ m { s , I “ , η “ , T “
300 s , and T b “ T { I “
10 s to also allow Rx to move and toreduce 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 availablefor transmission of each bit in a frame, for an MC system with fixed detection thresholds and fixed A , A “ t , , u . The fixed detection thresholds ξ are obtained from Subsection V-B by assuminguniform release. Fig. 7 reveals that in order to minimize the maximum BER in a frame, fewer moleculesshould be released at the beginning of the frame and the number of released molecules graduallyincreases with time. This is expected since, on average, for later release times, more molecules areneeded to compensate for the increasing distance between the transceivers.Fig. 8 shows the maximum BER within a frame for uniform release and the proposed releasedesign obtained from (38), with a fixed detection threshold obtained from (36), as a function of A for − − − − − − − T = 3000 s T = 300 s T = 100 s − − − A = 10 Bit index i B E R Uniform release, analyticalOptimal release, analyticalSimulation
Total number of molecules A M a x i m u m B E R Uniform release, analyticalOptimal release, analyticalSimulation
Fig. 8. Maximum BER in a frame as a function of A with uniformand optimal release. The inset shows the BER for each bit in aframe for uniform and optimal release for A “ and T “
300 s . , . . . . ψ = { . , . , . , . , . } Release time t [s] P r ( p ( t , T b ) > ψ ) AnalysisSimulation
Fig. 9. The probability that p p t, T b q is larger than a given value ψ as a function of t for T b “
10 s . T “ t , , u s . As can be observed, the proposed optimal release profile leads to significantperformance improvements compared to uniform release, especially for large A . For example, for A “ and T “ , the maximum BER is reduced by a factor of for optimal release comparedto 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. 8, we showthe BER as a function of bit index i in one frame for uniform and optimal release for A “ and T “
300 s . We observe that the optimal release achieves a lower maximum BER compared to theuniform release. We also observe that the optimal release leads to approximately the same BER foreach bit in a frame which highlights the benefits of the proposed design.Fig. 9 shows the probability that p p t, T b q is larger than a given value ψ , Pr p p p t, T b q ą ψ q , asa function of time t . Pr p p p t, T b q ą ψ q provides information about the probability that a releasedmolecule is absorbed at the Rx, i.e., the efficiency of molecule usage. We observe from Fig. 9 that Pr p p p t, T b q ą ψ q is a decreasing function of time as expected from the analysis in Subsection V-D.Moreover, for a given t , Pr p p p t, T b q ą ψ q is smaller for larger ψ . Furthermore, we can deduce themaximum time duration of a bit frame, T ‹ , satisfying a required molecule usage efficiency from Fig. 9.For example, for t ď
300 s , Pr p p p t, T b q ą . q ě . holds. Thus, T ‹ “ t ` T b “
310 s guarantees amolecule usage efficiency of Pr p p p t, T b q ą . q ě . . VII. C
ONCLUSIONS
In this paper, we considered a diffusive mobile MC system with an absorbing receiver, in whichboth the transceivers and the molecules diffuse. We provided a statistical analysis of the time-variantCIR and its integral, i.e., the probability that a molecule is absorbed by Rx during a given time period.We applied this statistical analysis to two system design problems, namely drug delivery and on-offkeying based MC with imperfect CSI. For the drug delivery system, we proposed an optimal releaseprofile which minimizes the number of released drug molecules while ensuring a target absorptionrate for the drugs at the diseased site during a prescribed time period. The probability of satisfyingthe constraint on the absorption rate was adopted as a system performance criterion and evaluated.We observed that ignoring the reality of the Tx ’s mobility for designing the release profile leadsto unsatisfactory performance. For the MC system with imperfect CSI, we optimized three designparameters, i.e., the detection threshold at Rx , the release profile at Tx , and the time duration of abit frame. Our simulation results revealed that the proposed MC system designs achieved a betterperformance in terms of BER and molecule usage efficiency compared to a uniform-release systemand a system without limitation on molecule usage, respectively. Overall, our results showed thattaking into account the time-variance of the channel of mobile MC systems is crucial for achievinghigh performance. A PPENDIX
A: P
ROOF OF L EMMA E t r p t qu p a q “ a D t E t γ u p b q “ a D t ? ´ λ ÿ n “ Γ p n ` q n !Γ p n ` { q ˆ λ ˙ n (40) “ c D tπ ´ λ ÿ n “ p n ` q ! p n ` q ! ` λ ˘ n “ c D tπ ´ λ ÿ n “ ˆ p n q ! p n q ! ` p n q ! p n ` q ! ˙ ` λ ˘ n p c q “ c D tπ ´ λ ˆ ` ? π ? λ e λ erf ˆ λ ? ˙ ` ? π ? λ e λ erf ˆ λ ? ˙˙ , where Γ p¨q denotes the Gamma function, equality p a q is due to (5), equality p b q is due to [37, Eq. (1.5)],and equality p c q is obtained by applying [31, Eq. (5.2.11.6) and Eq. (5.2.11.7)]. Simplifying the finalexpression, we obtain (6). To prove (7), we have E (cid:32) r p t q ( p a q “ D t E (cid:32) γ ( p b q “ D t ´ λ ÿ n “ Γ p n ` { q n !Γ p n ` { q ˆ λ ˙ n (41) “ D te ´ λ ˜ λ ÿ n “ p n ´ q ! ˆ λ ˙ n ´ ` ÿ n “ n ! ˆ λ ˙ n ¸ p c q “ D te ´ λ ˆ λ e λ ` e λ ˙ “ r ` D t, where equality p a q is obtained due to (5), equality p b q is due to [37, Eq. (1.5)], and equality p c q isobtained by applying the Maclaurin series of the exponential function. From (40) and (41), we obtain(7) since Var t r p t qu “ E t r p t qu ´ E t r p t qu .To prove (8), we have f r p t q p r q p a q “ ? D t f γ p γ q p b q “ rr ? πD t exp ˆ ´ r ` r D t ˙ sinh ˆ r r D t ˙ , (42)where f γ p γ q is the PDF of γ . Equality p a q in (42) exploits the fact that γ is a function of r p t q [34,Eq. (5-16)]. Equality p b q in (42) is obtained from the expression for PDF f γ p γ q [28, Eq. (1.6)] and therelation I { p x q “ b πx sinh p x q , where I { p x q is the Bessel function of the first kind and order { .Moreover, since r p t q? D Tx t follows a noncentral chi distribution, we obtain (9) as [29, Eq. (1)] F r p t q p r q “ F r p t q ? D t ˆ r ? D t ˙ “ ´ Q ˆ λ, r ? D t ˙ . (43)A PPENDIX
B: P
ROOF OF T HEOREM ˆ h p r, τ q and h p t, τ q are two functions of different variablesbut 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 ˆ h p r, τ q “ is equivalent to a cubic equation in r , given by ar ` br ` cr ` d “ , with properly defined coefficients a , b , c , d and discriminant ∆ “ abcd ´ b d ` b c ´ ac ´ a d . From (18), we have ∆ ă and thus ˆ h p r, τ q “ hasonly one real valued solution, denoted by r ‹ , which corresponds to the maximum value of ˆ h p r, τ q ,denoted by h ‹ . Then, from (18), we observe that ˆ h p r, τ q ą for r ă r ‹ and ˆ h p r, τ q ă for r ą r ‹ .Therefore, the equation ˆ h p r, τ q “ h has two solutions r p h q and r p h q , r p h q ă r p h q , when h ă h ‹ ,and has only one solution r ‹ when h “ h ‹ “ ˆ h p r ‹ , τ q . Finally, we derive (17) by exploiting [34,Eq. (5-16)] for the PDF of functions of random variables. Moreover, for h “ h ‹ , ˆ h p r, τ q “ so f h p t,τ q p h q “ f r p t q p r ‹ q ˆ h p r ‹ ,τ q Ñ 8 . A PPENDIX
C: P
ROOF OF L EMMA erf p ζ i p ξ, α i qq ´ erf ´ ξ ´ η ? η ¯ is convex in ξ . erf p¨q is aconvex and non-decreasing function for negative arguments and a concave and non-decreasing functionfor positive arguments. For η ă ξ ă µ i, , we have ζ i p ξ, α i q ă and ξ ´ η ? η ą . Moreover, ζ i p ξ, α i q and ξ ´ η ? η are affine functions of ξ . Then, erf p ζ i p ξ, α i qq is a convex and non-decreasing function of an affinefunction and thus is convex in ξ [32, Eq. (3.10)]. erf ´ ξ ´ η ? η ¯ is a concave and non-decreasing functionof an affine function and thus is concave in ξ [32, Eq. (3.10)]. Therefore, erf p ζ i p ξ, α i qq ´ erf ´ ξ ´ η ? η ¯ is convex, which concludes the proof.A PPENDIX
D: P
ROOF OF L EMMA erf p ζ i p ξ, α i qq is convex in α . First, erf p¨q is a convexand non-decreasing function for negative arguments and ζ i p ξ, α i q ă for ξ ă µ i, . Second, when η ă ξ , ζ i p ξ, α i q is a convex function in α since its Hessian matrix is positive semi-definite. Then, erf p ζ i p ξ, α i qq is a convex and non-decreasing function of a convex function, and thus, it is convex in α for η ă ξ ă µ i, [32, Eq. (3.10)], which concludes the proof.A PPENDIX
E: P
ROOF OF L EMMA B F r p t q p r qB t ă . Using the Taylor series expansion of the sinh p¨q function in (8), we obtain B F r p t q p r qB t “ BB t ˜ż r xr ? πD t exp ˆ ´ x ` r D t ˙ ÿ n “ p n ` q ! ˆ r x D t ˙ n ` d x ¸ (44) “ ż r xr ? πD ÿ n “ « p n ` q ! ˆ r x D ˙ n ` BB t ˆ exp ˆ ´ x ` r D t ˙ t ´ n ´ { ˙ff d x “ ż r xr ? πD ÿ n “ « p n ` q ! ˆ r x D ˙ n ` exp ˆ ´ x ` r D t ˙ t ´ n ´ { ˆ ˆ ´ n ´ { ` x ` r D t ˙ d x ď . R EFERENCES [1] V. Jamali, A. Ahmadzadeh, W. Wicke, A. Noel, and R. Schober, “Channel modeling for diffusive molecular communication–Atutorial review,”
Proceedings of the IEEE,
Early Access, pp. 1–46, 2019.[2] H. B. Yilmaz, A. C. Heren, T. Tugcu, and C. Chae, “Three-dimensional channel characteristics for molecular communicationswith an absorbing receiver,”
IEEE Commun. Lett. , vol. 18, no. 6, pp. 929–932, Jun. 2014. [3] N. Farsad, H. B. Yilmaz, A. Eckford, C. B. Chae, and W. Guo, “A comprehensive survey of recent advancements in molecularcommunication,” IEEE Commun. Surveys Tuts. , vol. 18, no. 3, pp. 1887–1919, Feb. 2016.[4] M. Femminella, G. Reali, and A. V. Vasilakos, “A molecular communications model for drug delivery,”
IEEE Trans. Nanobiosci. ,vol. 14, no. 8, pp. 935–945, Dec. 2015.[5] S. Kadloor, R. S. Adve, and A. W. Eckford, “Molecular communication using brownian motion with drift,” vol. 3, no. 1, pp.89–99, Jun. 2012.[6] A. Noel, K. Cheung, and R. Schober, “Improving receiver performance of diffusive molecular communication with enzymes,”
IEEE Trans. Nanobiosci. , vol. 13, no. 1, pp. 31–43, Mar. 2014.[7] U. A. K. Chude-Okonkwo, R. Malekian, B. T. Maharaj, and A. V. Vasilakos, “Molecular communication and nanonetwork fortargeted drug delivery: A survey,”
IEEE Commun. Surveys Tuts. , vol. 19, no. 4, pp. 3046–3096, May 2017.[8] A. Guney, B. Atakan, and O. B. Akan, “Mobile ad hoc nanonetworks with collision-based molecular communication,”
IEEETrans. Mobile Comput. , vol. 11, no. 3, pp. 353–366, Mar. 2012.[9] T. Nakano, Y. Okaie, S. Kobayashi, T. Koujin, C. Chan, Y. Hsu, T. Obuchi, T. Hara, Y. Hiraoka, and T. Haraguchi, “Performanceevaluation of leader-follower-based mobile molecular communication networks for target detection applications,”
IEEE Trans.Commun. , vol. 65, no. 2, pp. 663–676, Feb. 2017.[10] G. Chang, L. Lin, and H. Yan, “Adaptive detection and ISI mitigation for mobile molecular communication,”
IEEE Trans.Nanobiosci. , vol. 17, no. 1, pp. 21–35, Jan. 2018.[11] L. Lin, Q. Wu, F. Liu, and H. Yan, “Mutual information and maximum achievable rate for mobile molecular communicationsystems,”
IEEE Trans. Nanobiosci. , vol. 17, no. 4, pp. 507–517, Oct. 2018.[12] A. Ahmadzadeh, V. Jamali, and R. Schober, “Stochastic channel modeling for diffusive mobile molecular communication systems,”
IEEE Trans. Commun. , vol. 66, no. 12, pp. 6205–6220, Dec. 2018.[13] W. Haselmayr, S. M. H. Aejaz, A. T. Asyhari, A. Springer, and W. Guo, “Transposition errors in diffusion-based mobile molecularcommunication,”
IEEE Commun. Lett. , vol. 21, no. 9, pp. 1973–1976, Sep. 2017.[14] N. Varshney, W. Haselmayr, and W. Guo, “On flow-induced diffusive mobile molecular communication: First hitting time andperformance analysis,”
IEEE Trans. Mol. Biol. Multi-Scale Commun.,
Early Access, pp. 1–1, Jul. 2019.[15] N. Varshney, A. K. Jagannatham, and P. K. Varshney, “On diffusive molecular communication with mobile nanomachines,” in , Mar. 2018, pp. 1–6.[16] M. Sefidgar, M. Soltani, K. Raahemifar, H. Bazmara, S. Nayinian, and M. Bazargan, “Effect of tumor shape, size, and tissuetransport properties on drug delivery to solid tumors,”
Journal of Biological Engineering , vol. 8, no. 12, Jun. 2014.[17] A. Pluen, P. A. Netti, R. K. Jain, and D. A. Berk, “Diffusion of macromolecules in agarose gels: Comparison of linear andglobular configurations,”
Biophysical Journal , vol. 77, no. 1, pp. 542–552, 1999.[18] B. K. Lee, Y. H. Yun, and K. Park, “Smart nanoparticles for drug delivery: Boundaries and opportunities,”
Chemical EngineeringScience , vol. 125, pp. 158–164, Apr. 2015.[19] X. Wang, Y. Chen, L. Xue, N. Pothayee, R. Zhang, J. S. Riffle, T. M. Reineke, and L. A. Madsen, “Diffusion of drug deliverynanoparticles into biogels using time-resolved microMRI,”
The Journal of Physical Chemistry Letters , vol. 5, no. 21, pp. 3825–3830,Oct. 2014.[20] K. B. Sutradhar and C. D. Sumi, “Implantable microchip: the futuristic controlled drug delivery system,”
Drug Delivery , vol. 23,no. 1, pp. 1–11, Apr. 2014.[21] Y. Chahibi, M. Pierobon, and I. F. Akyildiz, “Pharmacokinetic modeling and biodistribution estimation through the molecularcommunication paradigm,”
IEEE Trans. Biomed. Eng. , vol. 62, no. 10, pp. 2410–2420, Oct. 2015.[22] S. Salehi, N. S. Moayedian, S. S. Assaf, R. G. Cid-Fuentes, J. Sol-Pareta, and E. Alarcn, “Releasing rate optimization in a singleand multiple transmitter local drug delivery system with limited resources,”
Nano Commun. Netw. , vol. 11, pp. 114–122, Mar.2017. [23] S. Salehi, N. S. Moayedian, S. H. Javanmard, and E. Alarcn, “Lifetime improvement of a multiple transmitter local drug deliverysystem based on diffusive molecular communication,” IEEE Trans. Nanobiosci. , vol. 17, no. 3, pp. 352–360, Jul. 2018.[24] T. N. Cao, A. Ahmadzadeh, V. Jamali, W. Wicke, P. L. Yeoh, J. Evans, and R. Schober, “Diffusive mobile MC for controlled-releasedrug delivery with absorbing receiver,” in , May 2019.[25] A. Ahmadzadeh, V. Jamali, A. Noel, and R. Schober, “Diffusive mobile molecular communications over time-variant channels,”
IEEE Commun. Lett. , vol. 21, no. 6, pp. 1265–1268, Jun 2017.[26] N. Farsad and A. Goldsmith, “A molecular communication system using acids, bases and hydrogen ions,” in , July 2016, pp. 1–6.[27] V. Jamali, N. Farsad, R. Schober, and A. Goldsmith, “Diffusive molecular communications with reactive signaling,” in , May 2018, pp. 1–7.[28] K. S. Miller, R. I. Bernstein, and L. Blumenson, “Generalized Rayleigh processes,”
Quarterly of Applied Mathematics , vol. 16,no. 2, pp. 137–145, Jul. 1958.[29] G. H. Robertson, “Computation of the noncentral chi-square distribution,”
The Bell System Technical Journal , vol. 48, no. 1, pp.201–207, Jan 1969.[30] Y. Deng, A. Noel, M. Elkashlan, A. Nallanathan, and K. C. Cheung, “Modeling and simulation of molecular communicationsystems with a reversible adsorption receiver,”
IEEE Trans. Mol. Biol. Multi-Scale Commun. , vol. 1, no. 4, pp. 347–362, Dec 2015.[31] A. P. Prudnikov, Y. A. Brychkov, and O. I. Marichev,
Integrals and Series . New York: Gordon and Breach Science, 1986, vol. 1.[32] S. Boyd and L. Vandenberghe,
Convex Optimization . USA: Cambridge University Press, 2004.[33] D. A. Stephens, “Math 556 Mathematical Statistics I - Some Inequalities,” Lecture Notes, Fall 2008.[34] A. Papoulis and S. U. Pillai,
Probability, Random Variables and Stochastic Processes . USA: McGraw-Hill, 2002.[35] D. Y. Arifin, L. Y. Lee, and C. Wang, “Mathematical modeling and simulation of drug release from microspheres: Implications todrug delivery systems,”
Advanced Drug Delivery Reviews , vol. 58, no. 12, pp. 1274–1325, Sep 2006.[36] S. Fredenberg, M. Wahlgren, M. Reslow, and A. Axelsson, “The mechanisms of drug release in poly(lactic-co-glycolic acid)-baseddrug delivery systems - A review,”
International Journal of Pharmaceutics , vol. 415, no. 1, pp. 34–52, May 2011.[37] J. H. Park, “Moments of the generalized Rayleigh distribution,”