A Novel A Priori Simulation Algorithm for Absorbing Receivers in Diffusion-Based Molecular Communication Systems
aa r X i v : . [ c s . ET ] D ec SUBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 1
A Novel
A Priori
Simulation Algorithm forAbsorbing Receivers in Diffusion-Based MolecularCommunication Systems
Yiran Wang,
Student Member, IEEE,
Adam Noel,
Member, IEEE, and Nan Yang,
Member, IEEE
Abstract —A novel a priori
Monte Carlo (APMC) algorithmis proposed to accurately simulate the molecules absorbedat spherical receiver(s) with low computational complexity indiffusion-based molecular communication (MC) systems. It isdemonstrated that the APMC algorithm achieves high simulationefficiency since by using this algorithm, the fraction of moleculesabsorbed for a relatively large time step length precisely matchesthe analytical result. Therefore, the APMC algorithm overcomesthe shortcoming of the existing refined Monte Carlo (RMC)algorithm which enables accurate simulation for a relatively smalltime step length only. Moreover, for the RMC algorithm, anexpression is proposed to quickly predict the simulation accuracyas a function of the time step length and system parameters,which facilitates the choice of simulation time step for a givensystem. Furthermore, a likelihood threshold is proposed forboth the RMC and APMC algorithms to significantly savecomputational complexity while causing an extremely small lossin accuracy.
Index Terms —Diffusion-based molecular communication, ab-sorbing receivers, molecular communication simulation, MonteCarlo method.
I. I
NTRODUCTION M OLECULAR communication (MC) has emerged as anunderpinning paradigm of exchanging and conveyinginformation among nano-devices in very small dimensions orspecific environments, such as water, tunnels, and human bod-ies [2]. Unlike electromagnetic wave-enabled communication,MC delivers information based on chemical changes within theenvironment. In MC, transmitters first send out information-carrying molecules to propagate within the environment. Later,such molecules are captured by receivers to allow them toobtain the carried information. Biological examples of MCinclude chemotactic signaling, calcium signaling, and bacterialmigration [2]. Thus, MC offers the advantages of low energyconsumption and potential for biocompatibility in liquid andgaseous media. Notably, the development of MC is envisionedto support transformational nano-applications, e.g., intra-bodyhealth monitoring, target drug delivery, and food and waterquality monitoring [2].
The material in this paper was presented in part at the InternationalWorkshop on Molecular, Biological and Multiscale Communications, IEEEInternational Conference on Sensing, Communication and Networking (IEEESECON 2018), in Hong Kong, China in June 2018 [1].Y. Wang and N. Yang are with the Research School of Engineering,Australian National University, Canberra, ACT 2600, Australia (Emails:{yiran.wang, nan.yang}@anu.edu.au).A. Noel is with the School of Engineering, University of Warwick,Coventry, CV4 7AL, UK (Email: [email protected]).
Diffusion-based MC has been acknowledged as a sim-ple but commonly adopted MC system within the nano-communication research community [3]. In this system, infor-mation molecules propagate using kinetic energy only, whichpreserves a high energy efficiency. One of the major challengesin designing and analyzing a diffusion-based MC systemis receiver modeling. The majority of existing MC studieshave considered two types of receivers: passive receivers andactive receivers [3, 4]. Passive receivers do not impose anyimpact on molecule propagation, while active receivers havesome mechanism for molecules to react either within or atthe surface of the receiver. The active receiver model ismore general and is better representative of typical biologicalreceivers. This motivates us to investigate the properties ofabsorbing receivers, which are a common ideal approximationfor active receivers.The notion of diffusion with absorption is a long-existingbiological phenomenon that has been investigated in the liter-ature, e.g., [5, 6]. For example, [5] investigated the diffusioninto adsorbers, where the adsorbed molecules stop diffusionupon hitting the receiver, but later may desorb and resume dif-fusion. Specifically, [5] examined the diffusion into a sphericaladsorber, an ellipsoidal adsorber, and disk-like adsorber(s) inan infinite medium and obtained the diffusion currents for suchadsorbers by solving Fick’s equations. It is noted that someconclusions drawn for adsorbers can be applied to absorbers,due to the similarity between adsorption and absorption.Unlike [5], [6] discussed the absorption of molecules that caneither diffuse through or react chemically with the receiversurface. Based on such studies, [7, 8] considered a diffusion-based MC system with a single perfectly absorbing receiverwithin an unbounded three-dimensional (3D) environment.To be specific, [7] provided the numerical results of thehitting rate of molecules at different times and the fraction ofmolecules absorbed for a given time, the analysis of which wasperformed in [9], while [8] presented a simulation framework.Recently, [10, 11] evaluated the impact of a receiver withreversible adsorption on the performance of MC systems,where a molecule can be released back to the environmentat some time after being captured by the receiver.Apart from the theoretical analysis of MC systems, e.g.,[7, 10, 11], the simulation of MC systems also serves asan effective means for performance evaluation. Against thisbackground, a number of simulation frameworks, e.g., N3Sim[12], NanoNS [13], BiNS2 [14], and AcCoRD [15], havebeen developed to examine the behavior of information par-ticles in various MC environments. The development of such
UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 2 frameworks is based on two common simulation approaches,the microscopic approach and the mesoscopic approach [16].Considering a sub-region of the simulation environment, themicroscopic approach treats the molecules within the regionindividually, while the mesoscopic approach treats them inaggregate as uniformly distributed throughout the region. Awidely recognized example of the microscopic approach inMC is the particle-based simulation that uses Brownian motionto characterize particle propagation. A detailed comparisonbetween the microscopic and mesoscopic approaches andfurther information on how they have been implemented inMC frameworks were shown in [15, 17]. In this paper, wefocus on the particle-based microscopic approach due to itsprevalence for the physical simulation of an MC environment.It is worthwhile noting that the existing microscopic algo-rithms incur high computational complexity in the simulationof the fraction of molecules absorbed in MC systems withabsorbing receivers. This high complexity leads to a longsimulation run time since the existing algorithms require a verysmall time step length to accurately model the absorption. Totackle this issue, we develop a new particle-based microscopicalgorithm in this work. This algorithm aims to reduce thecomputational complexity for the simulation of an irreversibleperfectly absorbing receiver in a 3D MC system where everymolecule that arrives at the receiver is absorbed.In the simulation procedure of the particle-based micro-scopic algorithm for diffusion-based MC systems, moleculesare moved by adding Gaussian random variables (RVs) to their x -, y -, and z -coordinates at the end of every simulation timestep. This makes the coordinates discrete functions of time,although in reality the movements of molecules undergoingBrownian motion are continuous over both time and space.In practice, molecules may actually diffuse into an absorbingreceiver between two time steps. To determine this absorp-tion, some existing simulation algorithms simply compare theobserved coordinates of molecules after diffusing with thecoordinates of the receiver [8, 10, 12]. As a consequence,the molecules that hit a receiver between time steps cannotbe considered as “absorbed” by using these simulation algo-rithms. We name the probability that a molecule crosses theRX’s boundary between time steps as the intra-step absorptionprobability . Recently, [15] and [18] investigated the possibilityof intra-step absorption. Specifically, [15] declared a moleculeas “absorbed” if its straight-line trajectory within a time stepcrossed an absorbing surface. Alternatively, [18] approximatedthe intra-step absorption probability for spherical receiverboundaries using the equation for flat planar receiver bound-aries given by [19, Eq. (10)]. This approximation was referredto as the refined Monte Carlo (RMC) algorithm.Our paper focuses on improving the accuracy and reducingthe computational cost of simulating absorbing receivers. Wefirst run simulations using the RMC algorithm and measure thechange in accuracy as the time step length increases. Then, wefit the results to an expression to estimate the accuracy giventhe diffusion coefficient D , time step length ∆ t , RX’s radius r r , and transmitter-to-receiver distance r . After recognizingthe poor accuracy of the RMC algorithm when ∆ t increases,we propose a new method for simulations of absorbing re- ceiver(s) with a large time step length, namely, the a priori Monte Carlo (APMC) algorithm. The APMC algorithm usesthe a priori probability of a molecule being absorbed to decidewhether it is absorbed in the current simulation time step. Ifthe molecule is determined as “absorbed” using the a priori probability, then we omit the diffusion step. We show thatthis algorithm achieves very high accuracy when the diffusionstep length is large relative to the size of the receiver. Despitethat the RMC algorithm performs accurately for small timestep size, we observe from MATLAB that the simulation runtime of the RMC algorithm is higher than that of the APMCalgorithm for small time step size. Furthermore, we identifythat a major contributor to the computational complexity ofthe RMC and APMC algorithms is the generation of uniformRVs when assessing intra-step absorptions. Thus, we proposea likelihood threshold for both algorithms to reduce thecomputational cost caused by the generation of excessive RVs.Our contributions extend our preliminary work in [1] andare summarized as follows:1) As in [1], we present a new algorithm, i.e., the APMCalgorithm, for MC simulation with a relatively large timestep length. The advantages of the APMC algorithm are:a) For the case of a single perfectly absorbing receiver, weshow that by using the APMC algorithm, the fraction ofmolecules absorbed precisely matches the correspond-ing analytical result when √ D ∆ t/r r is relatively large.b) For the case of two perfectly absorbing receivers, weshow that by using the APMC algorithm, the fractionof molecules absorbed approaches the asymptotic an-alytical value when time grows large.We note that the advantages of the APMC algorithm inboth cases cannot be achieved by the existing algorithms.2) We propose polynomial fitting expressions to predictthe accuracy of the RMC algorithm introduced in [18]for the case of a single perfectly absorbing receiver.Aided by numerical results, we show that the third orderpolynomial fitting expression is the most accurate one foraccuracy prediction. This allows us to use this expressionto characterize the accuracy of the simulation withoutrunning it.3) We investigate the computational complexity of both theRMC algorithm and the APMC algorithm. Specifically,we compare their MATLAB run times and explore thetrade-off between simulation accuracy and computationalcomplexity for both algorithms, based on which we pro-pose a likelihood threshold for the absorption probabilityto save computation time. Molecule absorption for theRMC algorithm and the APMC algorithm is possible ifand only if the calculated absorption likelihood in thesimulation is higher than the likelihood threshold. Aidedby numerical results, we show that applying a likelihoodthreshold can save as much as 20% of the total number ofgenerated RVs with an extremely small loss in simulationaccuracy.Comparing to our preliminary work [1], which only pre-sented the received signals of the SMC, RMC and APMCalgorithms for the system with a single absorbing receiver, UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 3 y-axis z - a x i s x - a x i s Receiver(r ,0,0) MoleculesTransmitter(0,0,0) r r Ω r Fig. 1: Illustration of our system model. The TX is a point transmitter locatedat (0 , , , the RX is a spherical irreversible perfectly absorbing receiverlocated at ( r , , with r r being the radius and Ω r being the RX’s perfectlyabsorbing boundary. Molecules propagate in the environment according toBrownian motion. Absorbing ReceiverMolecule’s initial positionbefore propagationMolecule’s final locationafter propagation A c t u a l t r a j ec t o r y Fig. 2: Illustration of the intra-step molecule movement. There is a possibilitythat a molecule crosses an absorbing boundary within one time step, even if itsinitial and final positions during that time step are both outside the absorbingreceiver. this paper additionally presents the received signals for multi-receiver configurations, proposes prediction expressions forthe RMC algorithm, and investigates the computational com-plexity of the RMC and APMC algorithms. Furthermore, themajority of the numerical results and discussions shown in thispaper are not included in [1].The rest of the paper is organized as follows. The systemmodel of interests is described in Section II. Existing sim-ulation algorithms, our proposed APMC algorithm, and thelikelihood threshold are presented in Section III. Numericalresults and discussions are provided in Section IV. In SectionV, we present our conclusions.II. S
YSTEM M ODEL
We consider a diffusion-based MC system within a 3Dspace, as depicted in Fig. 1. In this system, a point transmitter(TX) is located at the origin of the space and a sphericalperfectly absorbing receiver (RX) is centered at location ( r , , . We denote r r as the RX’s radius and Ω r as theRX’s boundary. At the beginning of a transmission process, theTX instantaneously releases N molecules. We assume that themolecules are small enough to be considered as points. Oncereleased, molecules diffuse in the environment according toBrownian motion until hitting the RX’s boundary. We denote N hit (Ω r , t | r ) as the number of molecules released from theTX at time t = 0 s and absorbed by the RX by time t . As per [9, Eq. (3.116)], we express N hit (Ω r , t | r ) as N hit (Ω r , t | r ) = N r r r erfc (cid:18) r − r r √ Dt (cid:19) , (1)where D is the diffusion coefficient and erfc ( · ) is the com-plementary error function. D describes the proportionalityconstant between the flux due to molecular diffusion and thegradient in the concentration of molecules.The majority of the existing simulation algorithms forabsorbing RXs did not consider the possibility of intra-stepmolecule absorption [18]. This absorption is depicted in Fig. 2,which shows that the actual trajectory of a molecule may crossthe RX’s boundary during one simulation time step, even if itsinitial position at the beginning of the time step and its finalposition at the end of the same time step are both outsidethe absorbing RX. If this crossing occurs, the molecule isabsorbed by the RX in practice. When the possibility of thisabsorption is ignored, the number of absorbed molecules isunderestimated, thus deteriorating the accuracy of simulation.In this work, we refer to the probability that a moleculeis absorbed during a simulation time step as the intra-stepabsorption probability . For MC systems, the discussion of thisprobability was limited. Two of the papers that consideredthe intra-step absorption probability are [15, 18]. In [15],the absorption of a molecule was determined by whetherits straight-line trajectory has crossed the absorbing surface.In [18], the intra-step absorption probability of a perfectlyabsorbing RX with a spherical boundary was approximatedas that of a perfectly absorbing RX with an infinite planar boundary, given by [19, Eq. (10)]Pr RMC = exp (cid:18) − l i l f D ∆ t (cid:19) , (2)where l i is the initial distance of a molecule from the absorbingboundary at the beginning of a time step, l f is the final distanceof a molecule from the absorbing boundary at the end of thesame time step, and ∆ t is the time step length.Apart from the aforementioned MC system which containsonly one absorbing spherical RX, we also consider a systemcontaining two identical absorbing spherical RXs. In thissystem, TX stays at the origin of the space and two absorbingRXs are located on opposite sides of and equidistant fromthe TX. Focusing on this two-RX system, [20] evaluated thefraction of molecules absorbed, which is defined as the ratiobetween the number of molecules absorbed at the RX and thetotal number of molecules released. As the time elapsed goesto infinity, the asymptotic fraction of molecules absorbed byeach RX is given by [20]Pr t →∞ = p µ − cos η ) × ∞ X n =0 e − ( n + ) µ sinh (cid:0) n + (cid:1) ( µ − µ )sinh (cid:0) n + (cid:1) ( µ − µ ) P n (cos( η )) , (3)where ( µ, η, φ ) are the bispherical coordinates correspond-ing to (0 , , in the natural coordinate system, µ =cosh − ( r /r r ) , µ = − cosh − ( r /r r ) , r r is the RX’s radius, r is the distance from the center of one RX to the TX, and P n is the n th-degree Legendre polynomial. The two-RX system UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 4
Absorbing Receiver Molecule’s initial positionbefore propagationMolecule’s final positionafter propagation A c t u a l t r a j ec t o r y O b s e r v e d t r a j ec t o r y Fig. 3: An example which shows the inaccuracy of the determination criterionin [15]. In this example, the line segment from the molecule’s initial positionto its final position crosses the RX’s surface, indicating that this molecule isabsorbed by the RX as per the determination criterion in [15]. However, themolecule’s actual trajectory indicates that this molecule is not absorbed.
Molecule’s initial locationbefore propagationMolecule’s final locationafter propagation l f l i ReceiverAbsorbing A c t u a l t r a j ec t o r y Fig. 4: Illustration of the distances used in the intra-step absorption probabilitycalculation in the RMC algorithm [18]. allows us to show that our proposed simulation algorithm canbe applied to not only a single-RX system but a multi-RXsystem, which will be demonstrated in Section IV.III. S
IMULATION A LGORITHMS
In this section, we first describe the common structure ofthe existing simulation algorithms for absorbing RXs. Then wepropose polynomial functions to predict the accuracy of theRMC algorithm. After that, we present our proposed APMCalgorithm and determine the likelihood threshold to be appliedto the APMC and RMC algorithms. The algorithms discussedin this paper are summarized in Table I.
A. Existing Simulation Algorithms
By observing the existing simulation algorithms for micro-scopic molecule absorption (such as those in [8, 15, 18]), weidentify that they follow a common structure. This structureis presented in
Algorithm 1 . For simulating molecules thatfollow Brownian motion, ∆ t needs to be carefully chosensuch that the diffusion step length is relatively small comparedto the distances that define reflective or absorbing boundaries[19]. We note that the ratio of the step length to the distance TABLE I: Comparison of Simulation Algorithms
Algorithm
Diffusion first? Intra-step absorption?SMC [6, 8] Yes NoRMC [18] Yes Yes. Absorption probability is givenby (2).AcCoRD [15] Yes Yes. Absorption occurs when themolecule trajectory crosses boundary.APMC No No. Molecules are absorbed beforebeing diffused.
Algorithm 1
A common structure of simulation algorithmsfor a single absorbing RX Determine the end time of simulation. for all simulation time steps do if t = 0 then Add N molecules to environment. end if Scan all not-yet-absorbed molecules, i.e., moleculeswhich are not absorbed by the RX. for all not-yet-absorbed molecules do Propagate each molecule for one step according toBrownian motion. Determine if the molecule is absorbed. end for end for affects the performance of the simulation algorithms, as willbe demonstrated in Section IV.We clarify that each algorithm has its own criterion fordetermining whether or not a molecule is absorbed by the RX,as given in Line 9 of
Algorithm 1 . Such determination criteriafor algorithms in [8, 15, 18] are summarized as follows: • As per the determination criterion in [8], the moleculesbeing observed inside the RX at the end of a time stepare absorbed. This criterion is referred to by [18] asthe simplistic Monte Carlo (SMC) algorithm. We notethat the performance of the SMC algorithm is inaccurate,unless the time step length is very small. This is becausethe SMC algorithm ignores the possibility of intra-stepabsorption and thus, underestimates the number of ab-sorbed molecules. • As per the determination criterion in [15], a molecule isabsorbed if the line segment from its initial position toits final position crosses the RX’s boundary. However, wenote that the line segment crossing the RX’s surface isneither sufficient nor necessary to correctly detect intra-step absorption. For example, in Fig. 2 the moleculeabsorption that actually occurs cannot be detected bythe criterion in [15]. In another case, shown in Fig. 3,molecule absorption is determined by the criterion in [15]but does not actually occur. The accuracy of the algorithmin [15] is better than but still comparable to the SMCalgorithm. • As per the determination criterion in [18], referred to asthe RMC algorithm and described in
Algorithm 2 , (2) isused to calculate the intra-step absorption probability ofan RX with a perfectly absorbing spherical boundary. As
UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 5
Algorithm 2
Molecule absorption determination in [18] if The molecule’s distance to ( r , , is smaller than r r , then The molecule is absorbed. else Calculate the intra-step absorption probability, Pr
RMC ,using (2). Generate a uniform RV u . if Pr RMC ≥ u then The molecule is absorbed. end if end if depicted in Fig. 4, in the RMC algorithm, l i denotes theshortest distance between the molecule’s initial positionand the RX’s boundary and l f denotes the shortest dis-tance between the molecule’s final position and the RX’sboundary. As shown in [18], the accuracy of the SMCalgorithm is comparable to that of the RMC algorithmwhen ∆ t is relatively small. When ∆ t increases whileother parameters remain the same, the accuracy of theRMC algorithm becomes measurably higher than that ofthe SMC algorithm. However, in addition to the resultspresented in [18], we run simulations of the SMC andRMC algorithms and observe that the simulated fractionof absorbed molecules of the RMC algorithm deviatesfrom the analytical one given by (2) when the root meansquare (RMS) of the diffusion step length, √ D ∆ t , isrelatively larger than r r . We will show the impact of ∆ t on the accuracy of the RMC algorithm in Section IV. B. Performance Prediction of RMC Algorithm
In this subsection, we propose a rule-of-thumb expressionto predict the accuracy of the RMC algorithm when theparameters of the MC system are within specified ranges.Specifically, our aim is to make a rough but fast predictionof the algorithm’s accuracy for given parameters, withoutresorting to simulations. To collect necessary data for accuracyprediction, we adopt the following procedure (the intermediateplots mentioned in the procedure are not presented in thispaper due to page limit):1) For the parameters listed in Table II, we run simulationsusing the RMC algorithm for one time step. We notethat the fraction of molecules absorbed at every currenttime step is based on the results generated in the previoustime step. In order to eliminate the impact of resultsfrom previous time steps, we choose to test for the firstsimulation time step only, which means that we calculatewith only two samples, one at the very beginning of thetransmission, and another at the end of the first time step.Simulation is repeated for N molecules and the resultobtained from the simulation of each molecule is calleda realization. For each set of D , ∆ t , and r , we plot R versus r r , where the measure of accuracy, R , is calculated according to R = 1 − P i =1 ( Pr hit ( i − − Pr sim ( i − P i =1 (cid:0) Pr sim ( i − − Pr sim (cid:1) , (4)where Pr hit ( i −
1) = N hit (Ω r , ( i − t | r ) /N, (5)is the analytical fraction of absorbed molecules at time ( i − t , as obtained from (1). We define Pr sim ( i − as the simulated fraction of absorbed molecules by theRMC algorithm at time ( i − t , which is obtained byaveraging the total number of absorbed molecules out ofall realizations at ( i − t over N , and refer to it as the i th simulated sample. Additionally, we define Pr sim as themean of all simulated samples.2) We plot reference lines that indicate different values of R in the same figure.3) We observe the values of r , r r , D , and ∆ t that achievethe selected values of R .4) For each selected value of R , we plot a 3D figurewhere the z -coordinate is the recorded value of r r , the x -coordinate is the recorded value of D ∆ t , and the y -coordinate is the recorded value of r .By observing the obtained results, we find that the accuracyof the RMC algorithm decreases when D ∆ t increases, r increases or r r decreases. To predict the accuracy, we define adimensionless variable κ that is calculated from the parametersof the MC system. We assume linear dependency between thevariables, i.e., κ and r r , for the sake of a quick and simplecalculation. Using the curve fitting tool in MATLAB, we fitthe 3D plots to equations of κ . Here, we define κ as κ = r r ( r D ∆ t ) − , (6)where the exponent of − is chosen such that the right handside is dimensionless. We next fit R to a polynomial functionof κ . We test the first, second, and third order polynomial fits ,which are given by R ≈ κ + 47100 , (7) R ≈ − κ + 392 κ − , (8)and R ≈ κ − κ + 813 κ − , (9)respectively. Given that R is the accuracy of the RMCalgorithm, which implies ≤ R ≤ , we impose . ≤ κ ≤ . for (9). If κ > . , then the RMC algorithm ispredicted to have asymptotically high accuracy. If κ < . ,then the RMC is predicted to have asymptotically low accu-racy. We then arrive at R ≈ , if κ < . , κ − κ +813 κ − , if . ≤ κ ≤ . , , otherwise . (10) We do not test higher order fits in this paper because we seek an easy-to-compute expression.
UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 6
TABLE II: Range of MC System Parameters for Performance Prediction
Parameter Notation and Range
TX-RX distance µ m ≤ r ≤ µ m Radius of RX µ m < r r ≤ r µ m RMS of diffusion step length µ m ≤ √ D ∆ t ≤ √ µ m Accuracy R ∈ { . , . , . , . , . , . , . , . , . } We will demonstrate in Section IV that the third order poly-nomial fit is the best match among (7), (8), and (10). Theaccuracy of the three polynomial fits will be shown in Fig. 8and the corresponding measured RMS error (RMSE) will beshown in Table III.
C. New A Priori Monte Carlo Algorithm
In this subsection, we propose a new simulation algorithmfor approximating the fraction of molecules absorbed at aperfectly absorbing RX when √ D ∆ t/r r is larger than thatconsidered for the RMC algorithm in [18]. We refer to thenewly proposed algorithm as the APMC algorithm.The procedure of the APMC algorithm is to first calculate, before the j th molecule diffuses, the probability that thismolecule will be absorbed in the current time step. Thisprobability depends on the distance between this moleculeand the center of the RX, d j , and the time step length, ∆ t .Specifically, the probability that this molecule will be absorbedin the current time step is calculated asPr APMC = r r d j erfc (cid:18) d j − r r √ D ∆ t (cid:19) , (11)which is obtained by scaling (1) by N , replacing the totalsimulation time t with ∆ t , and replacing r with d j . In (11),Pr APMC denotes the fraction of absorbed molecules releasedfrom a location d j away from the RX at time t = 0 s and absorbed by the RX by time ∆ t . We note that (11) iscalculated repeatedly in every time step for each of the freemolecules with the molecule’s current updated location whenit diffuses. Then the molecule absorption is determined bygenerating a uniform RV u , where ≤ u ≤ , and comparingits value with the probability obtained by (11). A moleculeis marked as “absorbed” if u ≤ Pr APMC . After determiningthe molecules absorbed, each not-yet-absorbed molecule ispropagated according to Brownian motion. If any propagatedmolecule is inside the RX’s boundary at the end of the currenttime step, we revert the diffusion of this molecule and letit propagate again, until this molecule diffuses to a locationoutside the RX. This is because if a molecule propagates to alocation inside the RX, it contradicts the preconditioning thatthe molecule is not absorbed. The APMC algorithm is detailedin
Algorithm 3 . D. Likelihood Threshold for Simulation Complexity Reduction
In this subsection, we propose a likelihood threshold to re-duce the computational complexity for the RMC algorithm and
Algorithm 3
The APMC algorithm for molecule absorption Determine the end time of simulation. for all simulation time steps do if t = 0 then Release N molecules into environment. end if Scan all not-yet-absorbed molecules. for all not-yet-absorbed molecules do Calculate the distance between the j th molecule to ( r , , , denote by d j . Calculate the absorbed probability Pr
APMC for each not-yet-absorbed molecule using (11) with r r , d j , D , and ∆ t . if Pr APMC ≥ u then The molecule is absorbed. end if end for for all not-yet-absorbed molecules do Propagate the molecule for one step. while the molecule’s distance to ( r , , ≤ r r do Revert the movement of this molecule to thelocation before propagation.
Propagate this molecule again. end while end for end for the APMC algorithm. For both algorithms, the computationalcomplexity consists of three parts and can be written as C sim = O (cid:18) t end ∆ t ( c diffuse + c absorb + c locate ) (cid:19) , (12)where c diffuse denotes the complexity of diffusing all molecules, c absorb denotes the complexity of determining the molecules tobe absorbed and absorbing these molecules, and c locate denotesthe complexity of determining whether a molecule is inside thespherical RX. We observe from MATLAB run time profilesthat the time for generating uniform RVs alone can take uparound of the total simulation run time of both algorithmsand that an increase in the number of generated RVs is closelyrelated to an increase in the total simulation run time for bothalgorithms. Therefore, we use the number of generated RVs tocharacterize the computational complexity of both algorithms.There are two types of RVs for the APMC algorithmand the RMC algorithm, namely, the uniform RVs used fortesting the molecule absorption and the Gaussian RVs used forpropagating molecules according to Brownian motion. Thus,we denote N u and N g as the number of uniform RVs and thenumber of Gaussian RVs, respectively. For each determinationof molecule absorption, as shown in Line 6 in Algorithm2 and Line 10 in
Algorithm 3 , we add to N u . For eachmolecule propagation, as shown in Line 8 in Algorithm 1 and Lines 15 and 18 in
Algorithm 3 , we add to N g sincethere are three Gaussian RVs added to the x -, y -, and z -coordinates of a molecule. Here, we adopt the commonly usedBox-Muller transform [21] to convert the number of generatedGaussian RVs to an equivalent number of generated uniform UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 7
RVs. According to [21], the generation of a pair of GaussianRVs requires /π pairs of uniform RVs. Thus, we use /π asa conversion factor and then the equivalent number of totallygenerated uniform RVs, which characterizes the computationalcomplexity of simulation, is given by N total = N u + (4 /π ) N g .Based on the aforementioned computational complexitycharacterization, we now propose a likelihood threshold, ξ ,to reduce the computational complexity. With this likelihoodthreshold, molecule absorption for the RMC algorithm ispossible if and only if Pr RMC ≥ ξ , while molecule absorptionfor the APMC algorithm is possible if and only if Pr APMC ≥ ξ .This can significantly reduce N u in simulation. Although theaccuracy of the RMC and APMC algorithms may slightlydecrease when ξ applies, an appropriate value of ξ canprovide a good trade-off between simulation accuracy andcomputational complexity. Such trade-off of the RMC andAPMC algorithms will be illustrated in Section IV.IV. N UMERICAL R ESULTS
In this section, we compare the fraction of molecules ab-sorbed determined by the SMC algorithm, the RMC algorithm,and the APMC algorithm, with the aid of particle-based simu-lations, to show the benefits of our APMC algorithm. We alsoinvestigate the accuracy of the prediction expression obtainedfrom empirical simulations in Section III-B by examining theperformance of the RMC algorithm. Furthermore, we comparesimulation run times of the RMC and APMC algorithmsand explore the trade-off between simulation accuracy andcomputational complexity of algorithms. The algorithms areimplemented in MATLAB.Throughout this section, we denote M as the number oftime-varying samples, and hence M − as the number oftime steps. Unless otherwise stated, we adopt the diffusioncoefficient of D = 10 − m / s , the TX-RX distance of r =50 µ m , and the number of molecules released of N = 10 . A. A Single Absorbing Receiver
In this subsection, we focus on the MC system with a singleabsorbing RX. We first examine the impact of r r and ∆ t on the fraction of molecules absorbed produced by the threealgorithms in Fig. 5 and Fig. 6, respectively.Fig. 5 plots the simulated fraction of molecules absorbedversus time t for M = 100 , ∆ t = 0 . , and r r = 20 µ m or . µ m . In this figure, the analytical result obtained using(1) is also plotted for examining the accuracy of the algo-rithms. From this figure, we observe that our APMC algorithmachieves a higher accuracy when r r decreases while ∆ t , D ,and r remain unchanged. Specifically, Fig. 5(a) shows that theRMC algorithm matches the fraction of molecules absorbedas predicted by the analytical result when √ D ∆ t/r r is small,which meets our expectation. Indeed, the performance of theRMC algorithm depends on the value of √ D ∆ t/r r . When √ D ∆ t/r r approaches 0 and r is larger than r r , the surfacearea of the RX can be approximated by an infinite plane.Therefore, the probability of a molecule entering the RX be-tween time steps is comparable to the probability of a moleculecrossing a flat planar boundary between time steps. We also Time [s] F r ac ti on o f M o l ec u l e s A b s o r b e d RMCSMCAPMCAnalytical0 42 6 8 (a) r r = 20 µ m Time [s] F r ac ti on o f M o l ec u l e s A b s o r b e d RMCSMCAPMCAnalytical (b) r r = 0 . µ m Fig. 5: Comparison of the fraction of molecules absorbed results producedby the SMC algorithm, the RMC algorithm, and the APMC algorithm versustime for different r r with M = 100 and ∆ t = 0 . . Time [s] F r ac ti on o f M o l ec u l e s A b s o r b e d (a) ∆ t = 0 . Time [s] F r ac ti on o f M o l ec u l e s A b s o r b e d (b) ∆ t = 5 s Fig. 6: Comparison of the fraction of molecules absorbed results producedby the SMC algorithm, the RMC algorithm, and the APMC algorithm versustime for different ∆ t with M = 10 and r r = 10 µ m . observe from Fig. 5(a) that when √ D ∆ t/r r is small, theSMC algorithm and the APMC algorithm underestimates andoverestimates the fraction of molecules absorbed, respectively.Unlike Fig. 5(a), Fig. 5(b) shows that the APMC algorithmmatches the fraction of molecules absorbed predicted bythe analytical result while the other two algorithms do not.Particularly, the fractions of molecules absorbed produced bythe RMC and SMC algorithms are very far away from theanalytical result when t > . Indeed, when √ D ∆ t/r r islarge, the spherical RX’s boundary cannot be approximatedby a flat planar boundary. Therefore, the RMC algorithmoverestimates the fraction of molecules absorbed.Fig. 6 plots the simulated fraction of molecules absorbed,together with the analytical result obtained using (1), versustime t for M = 10 , r r = 10 µ m , and ∆ t = 0 . or . We observe from this figure that when ∆ t increasesfrom . to , our APMC algorithm achieves a higheraccuracy. Specifically, Fig. 6(a) shows that the gap betweenthe fraction of molecules absorbed produced by the APMCalgorithm and that produced by the RMC algorithm is verysmall. Also, Fig. 6(a) shows that both the APMC and RMCalgorithms overestimate the fraction of molecules absorbed. UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 8
Time [s] N u m b e r o f N e w l y - A b s o r b e d M o l ec u l e s (a) RMC Algorithm Time [s] N u m b e r o f N e w l y - a b s o r b e d M o l ec u l e s (b) APMC Algorithm Fig. 7: Distribution of the number of newly-absorbed molecules during eachtime step for the RMC and APMC algorithms versus time for ∆ t = 0 . , M = 10 , and r r = 10 µ m . Simulations are repeated times and N = 10 molecules are released each time. The spectrum bars represent observationprobabilities.TABLE III: RMSE Measurements for Polynomial Fits in Fig. 8 Fig. 8(a) 0.0638 0.1487 0.0440Fig. 8(b) 0.0255 0.4378 0.0218
Unlike Fig. 6(a), Fig. 6(b) shows that the APMC algorithmmatches the fraction of molecules absorbed predicted by theanalytical result. The fraction of molecules absorbed producedby the RMC algorithm is approximately twice that producedby the APMC algorithm when t ≥
10 s in Fig. 6(b). Thisdemonstrates the accuracy of our APMC algorithm when √ D ∆ t/r r is large.Fig. 7 depicts the distribution of newly-absorbed moleculesduring each time step (e.g., from t = 1 s to t = 1 . )for both the RMC and APMC algorithms versus time t for M = 10 , r r = 10 µ m , and ∆ t = 0 . . The analytical resultin this figure is obtained using (1). Comparing Fig. 7(a) withFig. 7(b), we observe that the RMC algorithm significantlyoverestimates the number of newly-absorbed molecules when t = 1 s , while our APMC algorithm gives an improvedaccuracy when t = 1 s by absorbing molecules accordingto (11). When t increases, the overestimation of the RMCalgorithm becomes slightly less severe than that of the APMCalgorithm, this is in accordance with the observation fromFig. 6(a). More importantly, we observe that the distribution ofnewly-absorbed molecules in Fig. 7(b) is very similar to thatin Fig. 7(a), which demonstrates that our APMC algorithmdoes not noticeably disrupt the statistical distribution.We now examine the predicted accuracy of the RMC algo-rithm by comparing it with the measured one. The expressionsfor the first, second, and third order polynomial fits to predictaccuracy are given by (7), (8), and (10), respectively. We notethat the calculation of (7), (8), and (10) requires κ only. Asgiven by (6), κ is a function of D , ∆ t , r , and r r .Fig. 8 plots the predicted accuracy given by (7), (8), and(10) as well as the measured accuracy of the APMC and RMCalgorithms versus D ∆ t for r = 40 µ m , M = 2 , and r r =15 µ m or µ m . We observe from Fig. 8 that when D ∆ t increases, for a fixed r r , the accuracy of the RMC algorithm D t [ m ] A cc u r ac y , R (a) r r = 15 µ m D t [ m ] A cc u r ac y , R (b) r r = 20 µ m Fig. 8: Measured and predicted accuracy of the RMC algorithm and themeasured accuracy of the APMC algorithm versus D ∆ t for different r r with r = 40 µ m and M = 2 . D t [ m ] r r = 20 m A cc u r ac y , R r r = 25 m r r = 30 m Fig. 9: Measured and predicted accuracy of the RMC algorithm versus D ∆ t for different r r with M = 2 . is at first higher than that of the APMC algorithm, but soonbecomes much lower than that of the APMC algorithm. Wealso observe that the accuracy of the APMC stays close to 1for most of the range of D ∆ t considered. This demonstratesthe accuracy and robustness of our APMC algorithm when √ D ∆ t/r r is large.To quantitatively assess the accuracy of the first, second, andthird order polynomial fits in Fig. 8, we calculate the RMSEmeasurements of these fits viaRMSE = s P i =1 ( Pr hit ( i − − Pr sim ( i − , (13)where Pr hit ( i − , Pr sim ( i − , and i are defined the sameas in (4). The RMSE measurements for both subfigures aregiven in Table III. Based on both subfigures and Table III,we find that the third order polynomial fit for R is moreaccurate than the first and second order polynomial fits for r r = 15 µ m and µ m , since the third order polynomial fitis, on average, closer to the measured accuracy than the firstand second order polynomial fits, as shown in both subfigures,and achieves the lowest RMSE fitting measurement, as shownin Table III. Therefore, in Figs. 9 and 10 we only consider thethird order polynomial fit for the RMC algorithm. UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 9 Receiver radius r r [ m] A cc u r ac y , R RMC - MeasurementRMC - 3rd order predictionAPMC - Measurement t = 3s t = 1s t = 1s t = 3s
10 20 30 40 50
Fig. 10: Measured and predicted accuracy of the RMC algorithm and themeasured accuracy of the APMC algorithm versus r r for different ∆ t with M = 2 and r = 80 µ m . Fig. 9 plots the measured accuracy and the predicted ac-curacy of the RMC algorithm versus D ∆ t for r r = 20 µ m , µ m , and µ m . For the measured accuracy at a given D ∆ t , we run simulations for different combinations of D and ∆ t which achieve the same D ∆ t , e.g., the combination of D = 2 × − m / s and ∆ t = 0 . or the combination of D = 1 × − m / s and ∆ t = 1 s . The simulated resultfor each combination is displayed by a point (i.e., circle)in a scatter plot. We observe that the points overlap witheach other for a given D ∆ t , which demonstrates that theperformance of the RMC algorithm depends on D ∆ t , but not D or ∆ t separately. We further find that the average absolutevalue of the difference between the measured accuracy andthe predicted accuracy is . for r r = 20 µ m , . for r r = 25 µ m , and . for r r = 30 µ m , all of which arelower than .Fig. 10 plots the measured and predicted accuracy of theRMC algorithm together with the measured accuracy of theAPMC algorithm versus r r for large time step lengths ∆ t = 1 s and . As we observe from this figure, the accuracy of theAPMC algorithm stays close to 1 for most of the range of r r considered. This demonstrates that for large time step lengths,the APMC algorithm preserves a very high accuracy, whichagrees with the discussion on Fig. 6(b), and the variance inaccuracy is small. In addition, we observe that the measuredaccuracy of the RMC algorithm agrees well with the thirdorder polynomial fit given in (10). This is in accordancewith the observation made from Fig. 9 that the third orderpolynomial fit well approximates the measured accuracy when D ∆ t ≤ .Finally, we compare the computational complexity of theAPMC algorithm with that of the RMC algorithm. Fig. 11plots the average run time per realization of MATLAB sim-ulation for the APMC and RMC algorithms for ∆ t = 0 .
01 s , .
05 s , . , . , , , and . We clarify that each runtime is averaged over at least 20 realizations on an Intel i7desktop PC. We first observe from the figure that the averagedrun time per realization decreases when ∆ t becomes higher.This is due to the fact that when ∆ t increases, the number oftime steps in the simulation decreases. This leads to a decreasein the number of propagation and absorption operations, which -2 -1 Run Time per Realization [s] t = 5s t = 2s t = 1s t = 0.5s t = 0.1s t = 0.05s t = 0.01s APMCRMC (5.7294, 5.1504)s(1.1976, 1.0466)s(0.6234, 0.5340)s(0.1338, 0.1184)s(0.0718, 0.0632)s(0.0378, 0.0348)s(0.0174, 0.0180)s
Fig. 11: Average run time per realization of MATLAB simulation for theAPMC and RMC algorithms for r r = 20 µ m , N = 10 , and ∆ t = 0 .
01 s , .
05 s , . , . , , , and . -2 -1 A cc u r ac y , R N u m b e r o f R V s AccuracyRVs = 0.0404 (a) r r = 10 µ m -2 -1 A cc u r ac y , R N u m b e r o f R V s AccuracyRVs = 0.05455 (b) r r = 20 µ m Fig. 12: Number of RVs generated and measured accuracy of the APMCalgorithm versus ξ for r = 50 µ m , M = 10 , N = 10 , D = 10 − m / s , ∆ t = 10 s , and r r = 10 µ m or r r = 20 µ m . affects the computational complexity. Second, we observethat the APMC algorithm requires a slightly shorter run timeto simulate a system transmission process than the RMCalgorithm, except for ∆ t = 5 s . This is due to the fact thatthe APMC algorithm overestimates the number of absorbedmolecules when √ D ∆ t/r r is small. An overestimation meansthat more molecules are absorbed and removed permanentlyfrom the environment than expected, which leads to fewermolecules to be propagated and tracked by the system andthus reducing the computational complexity of the simulation.When the time step length increases to ∆ t = 5 s , √ D ∆ t/r r becomes large and the RMC algorithm leads to a more severeoverestimation than the APMC algorithm. Therefore, the RMCalgorithm requires a shorter time to simulate than the APMCalgorithm.Figs. 12 and 13 plot the computational complexity and themeasured accuracy of the APMC algorithm and the RMCalgorithm, respectively, versus the likelihood threshold ξ for r = 50 µ m , M = 10 , N = 10 , D = 10 − m / s , and r r = 10 µ m or r r = 20 µ m . We observe in all subfiguresthat when ξ increases, the measured accuracy decreases. Thisis because when a higher likelihood threshold is applied,the number of molecules that have absorption probabilities UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 10 -2 -1 A cc u r ac y , R N u m b e r o f R V s AccuracyRVs = 0.0545 (a) r r = 10 µ m -2 -1 A cc u r ac y , R N u m b e r o f R V s AccuracyRVs = 0.08283 (b) r r = 20 µ m Fig. 13: Number of RVs generated and measured accuracy of the RMCalgorithm versus ξ for r = 50 µ m , M = 10 , N = 10 , D = 10 − m / s , ∆ t = 0 . , and r r = 10 µ m or r r = 20 µ m . lower than the threshold increases. We note that ignoring lowabsorption probabilities leads to an underestimation of thenumber of absorbed molecules, resulting in lower accuracy.It also leads to a larger number of N g since N g dependson the number of molecules that need to be propagated. Wealso observe that the growth in N g is negligible compared tothe decrease in N u when ξ is low. We further observe thatwhen ξ increases beyond a certain value, such as ξ = 0 . in Fig. 12(a), ξ = 0 . in Fig. 12(b), ξ = 0 . inFig. 13(a), and ξ = 0 . in Fig. 13(b), the number ofgenerated RVs begins to increase, since the growth in N g eventually outnumbers the reduction in N u .In addition, we observe from Fig. 12 and Fig. 13 that when r r increases from µ m to µ m , the change in measuredaccuracy becomes less severe when applying the same ξ . Thisis because a larger RX leads to a higher a priori absorptionprobability in the APMC algorithm and a higher intra-stepabsorption probability in the RMC algorithm. Therefore, theimpact of the likelihood threshold on the number of absorp-tions is smaller for a larger RX. We clarify that Figs. 12 and 13reveal the impact of the likelihood threshold on the numberof required RVs and the measured accuracy. With the aid ofthese figures and similar results, a suitable likelihood thresholdcan be identified to achieve an acceptable loss of accuracy forsimulating other system parameters. B. Multiple Absorbing Receivers
In this subsection, we focus on the MC system with multipleabsorbing RXs. We run simulations for the system with twoabsorbing RXs, placed symmetrically as described in the lastparagraph of Section II, and the system with four absorbingRX, where the TX is a point located at the origin, and the RXsare located at ( r , , , ( − r , , , (0 , r , , and (0 , − r , in a Cartesian coordinate system with the radius as r r . Forboth systems, we examine the fraction of molecules absorbedas produced by the RMC and APMC algorithms. Due to thesymmetry of the systems, the probabilities that a molecule isabsorbed by each of the RXs are identical. Fig. 14: Comparison of the fraction of molecules absorbed for a system withtwo perfectly absorbing RXs and a system with four perfectly absorbing RXs.Results are produced by the RMC algorithm and the APMC algorithm versustime for r = 100 µ m , M = 5 × , D = 1 . × − m / s , ∆ t = 2 s ,with r r = 10 µ m and r r = 40 µ m . The asymptotic fraction of absorbedmolecules of one of the RXs in a system with two perfectly absorbing RXsas time goes to infinity is also plotted. Fig. 14 plots the simulated fraction of molecules ab-sorbed by one of the RXs versus time t for r = 100 µ m , M = 5 × , ∆ t = 2 s , and r r = { , } µ m for thesystem with two absorbing RXs and the system with fourabsorbing RXs. For the system with two absorbing RXs, wecompare the simulated results with the asymptotic fractionsof molecules absorbed calculated using (3) to examine theaccuracy of simulations. This is due to a lack of time-varyinganalytical results for the fraction of molecules absorbed bytwo absorbing RXs in the literature. We observe that when r r = 10 µ m , the fraction of molecules absorbed as producedby the APMC algorithm approaches the asymptotic resultswhen t is large, while the fraction of molecules absorbedas produced by the RMC algorithm is much higher than theasymptotic results. This is because when r r is small, thespherical RX’s boundary cannot be accurately approximatedby a flat planar boundary. We also observe that when r r increases to µ m , the fractions of molecules absorbed asproduced by both algorithms approach the asymptotic resultswhen t is large. Indeed, when r r increases, the approximationof the spherical RX’s boundary by a flat planar boundarybecomes more accurate. For the system with four absorbingRXs, we observe that when r r = 10 µ m , the fraction ofmolecules absorbed as produced by the APMC algorithm ismuch lower than that as produced by the RMC algorithm.Specifically, the RMC algorithm approaches a value whichis even higher than the asymptotic probability of a moleculebeing absorbed by one of the RXs in the two-RX systemwhen t is large, which is an intuitive overestimation. When r r increases to µ m , the fractions of molecules absorbed asproduced by both algorithms approach the same value, whichis similar to the system with two absorbing RXs.V. C ONCLUSION
In this paper, we proposed a new algorithm named the a priori
Monte Carlo (APMC) algorithm to use the a priori
UBMITTED TO IEEE TRANSACTIONS ON NANOBIOSCIENCE 11 approach for simulating the fraction of molecules being ab-sorbed by spherical receiver(s). Based on numerical results, wedemonstrated that the APMC algorithm outperforms the RMCalgorithm for large √ D ∆ t/r r . Thus, the APMC algorithm issuitable to be used in simulations for both the single absorbingreceiver case and the two absorbing receivers case when √ D ∆ t/r r is large. Moreover, for the single absorbing receivercase, we proposed an expression to predict the simulationaccuracy of an existing algorithm, the refined Monte Carlo(RMC) algorithm. We further investigated the MATLAB runtime of the APMC and RMC algorithms and applied likeli-hood thresholds to reduce the computational complexity ofboth algorithms. This investigation revealed that the APMCalgorithm significantly reduces the computational complexityof simulating absorbing receivers in MC systems withoutcompromising accuracy. Thus, it is worthwhile to extend theAPMC algorithm to other reactive surfaces where the a priori probability of surface interaction is known.R EFERENCES[1] Y. Wang, A. Noel, and N. Yang, “A new simulation algorithm forabsorbing receiver in molecular communication,” in
Proc. IEEE SECON2018 Int. Workshop Molecular, Biol. Multi-Scale Commun. , Hong Kong,China, June 2018, pp. 1–4.[2] T. Nakano, A. W. Eckford, and T. Haraguchi,
Molecular Communica-tion . Cambridge University Press, 2013.[3] N. Farsad, H. B. Yilmaz, A. Eckford, C. B. Chae, and W. Guo, “Acomprehensive survey of recent advancements in molecular communi-cation,”
IEEE Commun. Surv. Tut. , vol. 18, no. 3, pp. 1887–1919, Thirdquarter 2016.[4] T. Nakano, M. J. Moore, F. Wei, A. V. Vasilakos, and J. Shuai, “Molecu-lar communication and networking: Opportunities and challenges,”
IEEETrans. Nanobiosci. , vol. 11, no. 2, pp. 135–148, May 2012.[5] H. C. Berg,
Random Walks in Biology . Princeton University Press,1993.[6] J. Crank,
The Mathematics of Diffusion . Oxford University Press, 1979.[7] H. B. Yilmaz, A. C. Heren, T. Tugcu, and C. B. Chae, “Three-dimensional channel characteristics for molecular communications withan absorbing receiver,”
IEEE Commun. Lett. , vol. 18, no. 6, pp. 929–932,June 2014.[8] H. B. Yilmaz and C.-B. Chae, “Simulation study of molecular com-munication systems with an absorbing receiver: Modulation and ISImitigation techniques,”
Simulation Model. Practice. Th. , vol. 49, pp.136–150, Dec. 2014.[9] K. Schulten and I. Kosztin, “Lectures in theoretical biophysics,”
Uni-versity of Illinois , vol. 117, 2000.[10] Y. Deng, A. Noel, M. Elkashlan, A. Nallanathan, and K. C. Cheung,“Modeling and simulation of molecular communication systems witha reversible adsorption receiver,”
IEEE Trans. Mol. Biol. Multi-ScaleCommun. , vol. 1, no. 4, pp. 347–362, Dec. 2015.[11] A. Ahmadzadeh, H. Arjmandi, A. Burkovski, and R. Schober, “Com-prehensive reactive receiver modeling for diffusive molecular commu-nication systems: Reversible binding, molecule degradation, and finitenumber of receptors,”
IEEE Trans. Nanobiosci. , vol. 15, no. 7, pp. 713–727, Dec. 2016.[12] I. Llatser, D. Demiray, A. Cabellos-Aparicio, D. T. Altilar, and E. Alar-cón, “N3Sim: Simulation framework for diffusion-based molecular com-munication nanonetworks,”
Simulation Model. Practice. Th. , vol. 42, pp.210–222, Mar. 2014.[13] E. Gul, B. Atakan, and O. B. Akan, “Nanons: A nanoscale networksimulator framework for molecular communications,”
Nano Commun.Netw. , vol. 1, no. 2, pp. 138–156, June 2010.[14] L. Felicetti, M. Femminella, G. Reali, P. Gresele, and M. Malvestiti,“Simulating an in vitro experiment on nanoscale communications byusing BiNS2,”
Nano Commun. Netw. , vol. 4, no. 4, pp. 172–180, Dec.2013.[15] A. Noel, K. C. Cheung, R. Schober, D. Makrakis, and A. Hafid,“Simulating with AcCoRD: Actor-based communication via reaction–diffusion,”
Nano Commun. Netw. , vol. 11, pp. 44–75, Mar. 2017.[16] A. Hellander, S. Hellander, and P. LoÌ´Ltstedt, “Coupled mesoscopicand microscopic simulation of stochastic reaction-diffusion processes in mixed dimensions,”
Multi-scale Model. Sim. , vol. 10, no. 2, pp. 585–611, May 2012.[17] N. Farsad, H. B. Yilmaz, A. Eckford, C.-B. Chae, and W. Guo, “A com-prehensive survey of recent advancements in molecular communication,”
IEEE Commun. Surv. Tut. , vol. 18, no. 3, pp. 1887–1919, Third Quarter2016.[18] D. Arifler and D. Arifler, “Monte Carlo analysis of molecule absorptionprobabilities in diffusion-based nanoscale communication systems withmultiple receivers,”
IEEE Trans. Nanobiosci. , vol. 16, no. 3, pp. 157–165, Apr. 2017.[19] S. S. Andrews and D. Bray, “Stochastic simulation of chemical reactionswith spatial resolution and single molecule detail,”
Phys. Biol. , vol. 1,no. 3, p. 137, Aug. 2004.[20] H. Sano, “Solutions to the smoluchowski equation for problems involv-ing the anisotropic diffusion or absorption of a particle,”
J. Chem. Phys. ,vol. 74, no. 2, pp. 1394–1400, Aug. 1981.[21] G. Marsaglia and T. A. Bray, “A convenient method for generatingnormal variables,”