Fazilet Gokbudak

Cevdet Bayram

Tuna Akdeniz

Ali Emre Tugcu

Pusane

Spatial Receptor Allocation for a Multiple Access Hub in Nanonetworks

Keywords: Molecular Communication via Diffusion, spatial transmitter allocation, nanonetworks

Following the same trends in other communication schemes, design of star network topologies will gain importance in order to communicate with a centralized hub in molecular communications as well. For a molecular communication multiple access network system, interlink interference (ILI) becomes a crucial issue in addition to intersymbol interference (ISI). This paper focuses on spatial allocation of the receptors of a centralized receiver among transmitters by proposing efficient solutions to reduce both ISI and ILI. The spatial allocation scheme has been designed using an optimization problem and verified with Monte Carlo simulations.

I. INTRODUCTION

In the last decade, molecular communication has increased its popularity as a possible solution of the communication of nanomachines [START_REF] Akyildiz | Nanonetworks: A new communication paradigm[END_REF], [START_REF] Farsad | A comprehensive survey of recent advancements in molecular communication[END_REF]. In recent molecular communication literature, different kinds of networks have been examined. These networks can be categorized into two network topologies as star networks and point-to-point networks. In pointto-point networks, main objective is to transmit an intended symbol from source to destination using relay nodes, while in star networks, multiple transmitters try to communicate with a centralized hub.

As point-to-point networks, one-way multi-hop diffusion nanonetworks are considered in [START_REF] Ahmadzadeh | Analysis and design of multi-hop diffusion-based molecular communication networks[END_REF], and the decode and forward method is proposed to increase the performance of the system. In [START_REF] Akdeniz | Novel network coding approaches for diffusion-based molecular nanonetworks[END_REF] and [START_REF] Aijaz | On error performance of network coding in diffusion-based molecular nanonetworks[END_REF], network coding is used for twoway communication to increase the performance of the system by increasing the data rate. In [START_REF] Tiwari | Error performance optimization using logarithmic barrier function in molecular nanonetworks[END_REF], an optimal barrier function is derived to minimize the error in point-to-point networks and a routing protocol for the network has been presented in [START_REF] Aijaz | Opportunistic routing in diffusion-based molecular nanonetworks[END_REF].

In the molecular communication literature, star network topologies have also been considered frequently. One of the possible applications of the star topology in nanonetworks can be healthcare monitoring by transmitting the sensed information to a centralized receiver, which can process these transmitted information as proposed in [START_REF] Moore | Diffusion-based multiple access by nano-transmitters to a micro-receiver[END_REF]. In [START_REF] Pierobon | A statistical-physical model of interference in diffusion-based molecular nanonetworks[END_REF], the statistical interference model of the channel is derived when multiple nanomachines release molecules to a centralized receiver. In [START_REF] Moore | Diffusion-based multiple access by nano-transmitters to a micro-receiver[END_REF], the bit error rate (BER) is optimized by scheduling the transmitters' release time of molecules and adjusting the concentration of these released molecules. In [START_REF] Zabini | Spatially distributed molecular communications: An asynchronous stochastic model[END_REF], a stochastic model of the network that involves a spherical receiver and randomly placed transmitters whose emission times are also random is derived. In [START_REF] Honary | Modeling molecular channel using queueing theory approach[END_REF], queueing network theory is used to model molecular nanonetworks.

Although there are various works related to nanonetworks with star network topology, to the best of our knowledge, none of them has examined the optimum spatial allocation of the receptors of the centralized receiver among transmitters to optimize communication performance. On the other hand, in relay network topology, optimum spatial position of the relay is determined to improve the performance [START_REF] Tavakkoli | Optimal positioning of relay node in cooperative molecular communication networks[END_REF]. In this paper, we propose a method to find the optimal spatial allocation of the receiver's receptors among the transmitters. As proposed in [START_REF] Akdeniz | Molecular signal modeling of a partially counting absorbing spherical receiver[END_REF], the transmitted molecules are more likely to be absorbed by the receptors that are closer to the transmitter. This is not surprising, since for these closer regions, molecules travel over shorter distances compared to the other regions. Hence, the molecules absorbed by the receptors of the receiver that are far from the transmitter require more time to reach these regions, and they are most likely to belong to the previous symbols, contributing to intersymbol interference (ISI). In [START_REF] Akdeniz | Molecular signal modeling of a partially counting absorbing spherical receiver[END_REF], it has been proved that the part of the receiver that is far from the transmitter should not be used for decoding. Inspired from this idea, we claim that there is plenty of room at the receiver for other transmitters. In this paper, we propose a method to allocate the parts of the receiver to different transmitters. The main motivation behind this allocation is determining the area of the receptors assigned to each transmitter and also determining the barriers between two consecutive receptor regions assigned to two neighboring transmitters in order to reduce both ISI and interlink interference (ILI).

II. SYSTEM MODEL

We consider N point nanotransmitters that are evenly placed around a centralized receiver. For the sake of simplicity, we consider N identical channels, but the proposed methods can be modified for different scenarios. Each transmitter aims to send its own message to the centralized spherical absorbing receiver with radius r r by releasing molecules in the diffusive fluidic channel. As shown in Fig. 1, different parts of the receiver surface are allocated for receiving messages from different transmitters. In Fig. 1a, the area of the allocated part of the receiver for each transmitter (also named as patch) is defined by the base of a spherical sector, whose apical angle is 2α r . Furthermore, as seen from the same figure, each pair of consecutive patches are separated by a neutral zone whose apical angle is α s , and the molecules absorbed at these neutral zones are not taken into account by the receiver for decision. In particular, these zones can be considered as barriers between consecutive patches in order to reduce ILI.

It is assumed that each transmitter sends its binary messages by on-off keying, which is a special case of concentration shift keying. In particular, each transmitter releases M molecules at the beginning of the transmission period t s for transmitting bit-1 and releases nothing to convey bit-0. The molecules in the channel move under the effects of free diffusion, which is explained in [START_REF] Yilmaz | Threedimensional channel characteristics for molecular communications with an absorbing receiver[END_REF]. At the end of the transmission period, the receiver makes a decision on the symbols transmitted from the transmitters by using the number of absorbed molecules at each patch.

In order to propose our approaches, a method for evaluating the number of the transmitters (hence, corresponding receivers) for a given α r and α s is required. As an illustrative example in Fig. 1, assuming that 2α r +α s = π/2, six transmitters exist for communications. Furthermore, we derive the analytical expression of the number of transmitter-receiver pairs for given α r and α s . Note that each receiver surface corresponds to a cone, whose apical angle is 2α r . Between all neighbour receivers, there is a α s gap to eliminate interference. Considering these, we can say that we require an angle of r + α s /2 + α s /2= 2α r + α s for each receiver. Therefore, in order to calculate the area required to place a receiver, we need to calculate the area of the spherical cap that coincides with the base of the cone whose apical angle is 2α r + α s , which is equal to 2πr2 r (1 -cos(α r + α s /2)). It is clear that the upper limit of the number of transmitters for given α r and α s can be obtained by dividing the area of the sphere by the area of the spherical cap that coincides with the base of the cone whose apical angle is 2α r + α s as

N = 4πr 2 r 2πr 2 r (1 -cos(α r + αs 2 )) = 2 (1 -cos(α r + αs 2 ))
.

(

) 1
Since each channel is identical (due to symmetry), it is enough to examine the communication performance of only one of them to determine the overall performance. The performance is evaluated using BER, which is affected by ISI due to the remaining molecules released by the transmitter for previous symbols and ILI caused by the molecules released by the neighboring transmitters. Considering the leftmost channel in Fig. 1a, ILI can most possibly occur due to the transmission from neighboring transmitters. Due to the nature of diffusion, amount of ILI coming from the rightmost channel to the leftmost channel is very limited compared to other channels. Therefore, the number of neighbors for any patch, N 1 , for a given α r and α s pair must be determined. In order to calculate N 1 , we need to define the transmitter orbit, which is a hypothetical circle, and all patches assigned to these neighboring transmitters are lined up on this orbit as shown in Fig. 1a and1b. Considering Fig. 1a, one can easily compute the radius of this orbit r or as r or = r r sin(2α r + α s ).

(

) 2
Therefore, the patches are lined up to a circle whose circumference is 2πr or . Accordingly, the number of patches placed in this orbit (N 1 ) can be obtained by dividing this circumference to the length of the arc of any patch that coincides with the orbit. Since the angle of the arc is (2α r + α s ), as can be seen in Fig. 1b, the arclength can be calculated as r (2α r + α s ). Therefore, N 1 can be written as

N 1 = 2πr r sin(2α r + α s ) r r (2α r + α s ) = 2π sin(2α r + α s ) (2α r + α s ) . (3)
It is possible that there can be many orbits satisfying 2kα r + kα s ≤ π/2 (for example, for the second closest neighbour receivers, the radius of the orbit would be r r sin(4α r + 2α s ) ).

Therefore, one can calculate the number of transmitters at the k -th orbit as

N k = 2πr r sin(2kα r + kα s ) r r (2α r + α s ) = 2π sin(2α r + α s ) (2α r + α ) . (4)
Hence, considering Fig. 1a, 2α r + α s = π/2, there is only one transmitter orbit, and the number of patches at this orbit can be obtained using (4) as 2π sin(π/2)/(π/2) = 4, which can also be verified by Fig. 1b. Once the number of neighboring transmitters (hence patches) is derived, the next step is to determine the fraction of absorbed molecules by the patch to the molecules released from the assigned transmitter to this patch at the n -th slot (p n ) and the fraction of absorbed molecules released from the neighbors of this transmitter to the molecules released from the assigned transmitter (l n ). In other words, p n and l n correspond to ISI and ILI, respectively, while p 1 is related to the intended signal. In [START_REF] Akdeniz | Molecular signal modeling of a partially counting absorbing spherical receiver[END_REF], the cumulative distribution of absorbed molecules by a spherical receiver until time t and angle α (F (α, t)) has already been derived. Using this function, the channel taps can be computed as pn(αr, ts) = F (αr, nts) -F (αr, (n-1)ts) , ln(αr, αs, ts) = (F (3αr + αs, nts) -F (3αr + αs, (n-1)ts)) -(F (αr + αs, nts) -F (αr + αs, (n-1)ts)) /N1.

(

Note that l n (α r , α s , t s ) corresponds to the ILI taps that correspond to the channel closest to the corresponding channel, whose channel taps are defined by p n (α r , t s ).

III. OPTIMUM SPATIAL ALLOCATION

In this paper, we investigate the optimum spatial allocation of the receiver to transmitters for two cases. The first one is the case where the number of transmitters, N , is fixed. For the second case, the total data rate transmitted to the receiver is fixed. In particular, the number of allowed bit transmissions per second is fixed and the number of transmitters is also an additional variable for this second problem. For both cases, our aim is to determine α r and α s angles that achieve the lowest BER. BER function in molecular communication channel can be obtained by considering all possible transmitted sequences, hence involves 2 L terms, where L is the channel memory. In multiple access channels, one also needs to take ILI into account by considering the sequences that the neighboring transmitters send. All in all, BER involves 2 LN1 terms, and it is not a tractable function. [START_REF] Akdeniz | Optimal reception delay in diffusion-based molecular communication[END_REF] proposes an objective function named signal-to-interference difference (SID), which is basically p 1 -L n=2 p n and provides approximately the same angles that BER gives. Inspired from this function, we propose the modified version of this function for multiple access channels. Since ILI is another factor that affects the performance of the system, signal-to-inter-symbol and -link interference difference function (SILD) is proposed to determine optimized α r and α s values as

SILD(α r , α s , t s ) =p 1 (α r , t s ) - L n=2 p n (α r , t s ) - L n=2 l n (α r , α s , t s ). (6)
A. Optimum Spatial Allocation for a Fixed Number of Transmitters

In the first case, it is assumed that there is a fixed number of transmitters (N ) that transmit their symbols during a fixed symbol interval t s , and the aim is to determine α r and α s to obtain the lowest BER or highest SILD. Since N is fixed, it can be concluded that 2α r + α s is also constrained by [START_REF] Akyildiz | Nanonetworks: A new communication paradigm[END_REF]. Considering this constraint, the optimization problem is reduced to finding either α r or α s by achieving maximum SILD as

α * r = argmax αr SILD(α r , α s , t s ) subject to: N < 2 (1 -cos(α r + α s /2)) < N + 1. (7)

B. Optimum Spatial Allocation for Fixed Data Rate

When the number of bits transmitted per second, namely data rate (R), is fixed, the number of transmitters becomes an additional variable for the optimization problem. In this case, for any α r and α s , the number of transmitters N can be obtained as given in [START_REF] Akyildiz | Nanonetworks: A new communication paradigm[END_REF]. Therefore, in order to achieve R bits per second, the required transmission symbol duration (t req s ) can be obtained as

t req s = R/N. (8)
Since the data rate (R) is fixed, the total number of transmitters N and the required symbol duration to achieve this data rate t req s are constraint variables. Therefore, considering [START_REF] Aijaz | Opportunistic routing in diffusion-based molecular nanonetworks[END_REF], the constraint turns out to be R = t req s 2/(1 -cos(α r + α s /2)) . This constraint implies that, unlike the optimization problem defined in [START_REF] Aijaz | Opportunistic routing in diffusion-based molecular nanonetworks[END_REF], there are two parameters that optimize the SILD, which are v = {α s , α r }. Therefore, the optimization problem for fixed data rate can be written as

v * = argmax
IV. PERFORMANCE ANALYSIS Although the analytical expression of SILD can be obtained by using the results in [START_REF] Akdeniz | Molecular signal modeling of a partially counting absorbing spherical receiver[END_REF], the corresponding SILD function is so complicated that it is not tractable to obtain an analytical solution. Therefore, performance analysis of the proposed systems is done by solving the optimization problems defined In order to evaluate BER, once the channel taps are obtained from [START_REF] Akdeniz | Molecular signal modeling of a partially counting absorbing spherical receiver[END_REF], normal approximations are done to model the number of arriving molecules to the receiver as done in [START_REF] Yilmaz | Arrival modelling for molecular communication via diffusion[END_REF]. In particular, to transmit bit-1, the transmitter releases M molecules and the arrival of these molecules can be modelled by using the channel taps, p n . Similarly, for modelling the arrival of the molecules from the neighbour channels, same methodology can be used by utilizing l n . For BER simulations, the threshold is always chosen to obtain the lowest BER in order to determine the ultimate performance of the systems. The performed simulations are conducted for 10 4 times and, for each simulation, 10 5 consecutive binary symbols are transmitted in a channel with memory L = 100, and M is chosen as 500 molecules.

It is observed from Fig. 2 that as the number of users N increases, the optimum angle α r decreases. This is expected since the allocated region to each transmitter should be narrower for more users. Furthermore, it can also be seen that the angle parameters of the minimum BER and maximum SILD are almost the same, which enables us to find the optimum receptor allocation scheme using SILD optimization without exhaustive BER simulations.

In Fig. 3, the optimization problem defined in ( 9) is examined for various data rates. In particular, as can be seen in ( 9), the data rate is considered as the amount of data transferred to the receiver. Therefore, both the number of transmitters N and corresponding transmission slot t req s are not constant, and they are easily determined by finding the optimum α r and α s . One can observe in Fig. 3 that, as the data rate increases, both α r and α s decrease. The decrease in α r and α s are expected since, as the amount of data transferred to the receiver increases, more transmitters are required and when the number of transmitters is increased, we need to compress them in smaller patches, which implies smaller α r and α s values. Furthermore, α r and α s values obtained by the proposed SILD function yield almost the same values with minimum BER obtained using exhaustive search.

V. CONCLUSION

In this paper, a star network communication scheme for nanonetworks is proposed for reducing the ISI and ILI. The proposed scheme is based on the allocation of the absorbing surface of the receiver (patch) to transmitters by assigning a neutral zone between neighboring receivers. In particular, the proposed approach consists of determining the surface area of the patch assigned to each transmitter as well as neutral zone. The optimization problems are constructed to determine the optimum surface area for each patch, and the results are verified with simulations that show that the solution of the proposed optimization problem yields minimum BER without performing exhaustive simulations for BER calculations.

Fig. 1 :

 1 Fig. 1: Left: Spatial allocation of the spherical receiver patches by point transmitters. The apical angle for each receiver is 2 α r and there are neutral zones with α s that separate neighboring receptor parts. Right: Expanded top view of the communication system.

2 ( 1 -

 21 αr,αs SILD(α r , α s , t req s )subject to: R = t req s cos(α r + α s /2)) .

Fig. 2 :Fig. 3 :

 23 Fig. 2: Optimum α to obtain lowest BER and maximum SILD

(a) Channel model for N = 6 Tx-Rx pair (b) Expanded view of the channel model


Fatih Dinc

Leander Thiele

Bayram Cevdet Akdeniz

The Effective Geometry Monte Carlo Algorithm: Applications to Molecular Communication

In this work, we address the systematic biases and random errors stemming from finite step sizes encountered in diffusion simulations. We introduce the Effective Geometry Monte Carlo (EG-MC) simulation algorithm which modifies the geometry of the receiver. We motivate our approach in a 1D toy model and then apply our findings to a spherical absorbing receiver in a 3D unbounded environment. We show that with minimal computational cost the impulse response of this receiver can be precisely simulated using EG-MC. Afterwards, we demonstrate the accuracy of our simulations and give tight constraints on the single free parameter in EG-MC. Finally, we comment on the range of applicability of our results. While we present the EG-MC algorithm for the specific case of molecular diffusion, we believe that analogous methods with effective geometry manipulations can be utilized to approach a variety of problems in other branches of physics such as condensed matter physics and cosmological large scale structure simulations.

Introduction

Understanding the communication between complex systems of diffusing molecules requires fast and reliable simulation methods. While many methods have been suggested in the literature to simulate diffusion processes inside such systems, they mainly focus on changing the path of the diffusing molecules to overcome systematic errors that arise from the finite step size. We first give a brief review of these approaches.

Molecular communication is growing in importance [START_REF] Farsad | A comprehensive survey of recent advancements in molecular communication[END_REF]. The understanding of communication between nanomachines is extremely important for further development in this area. In this pursuit, molecular communication via diffusion (MCvD) channels has been explored extensively in the literature [START_REF] Yilmaz | Arrival modelling for molecular communication via diffusion[END_REF][START_REF] Nakano | Channel model and capacity analysis of molecular communication with brownian motion[END_REF]. Although analytical results are derived for some simple diffusion channels, such as a point transmitter and a spherical absorbing receiver in an unbounded environment [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF] and in bounded environments [START_REF] Dinc | Impulse response of the channel with a spherical absorbing receiver and a spherical reflecting boundary[END_REF][START_REF] Al-Zubi | Modeling of ligand-receptor protein interaction in biodegradable spherical bounded biological micro-environments[END_REF], a deeper analysis including more complicated structures requires fast and reliable simulations of sophisticated diffusion channels. To meet this demand, many different simulation frameworks have been proposed in the literature, such as N3Sim [START_REF] Llatser | N3sim: Simulation framework for diffusion-based molecular communication nanonetworks[END_REF], NanoNS [START_REF] Gul | Nanons: A nanoscale network simulator framework for molecular communications[END_REF], BiNS2 [START_REF] Felicetti | Simulating an in vitro experiment on nanoscale communications by using bins2[END_REF] and AcCoRD [START_REF] Noel | Simulating with accord: Actor-based communication via reaction-diffusion[END_REF]. The most simple algorithm performing diffusion channel simulations is usually referred to as a Monte Carlo (MC) simulation algorithm. In this algorithm, the molecules iteratively diffuse inside the channel following a Gaussian probability distribution for their step sizes. In each iteration, the molecule position is calculated and molecules that are inside the volume of the receiver are absorbed [START_REF] Yilmaz | Simulation study of molecular communication systems with an absorbing receiver: Modulation and isi mitigation techniques[END_REF]. This algorithm performs quite accurately for small step sizes while requiring high computation power, i.e. iteration number, in order to simulate complex systems. For large step sizes, the accuracy of the simulation decreases due a combination of different effects. One such effect, called the intra-step absorption problem, is considered in [START_REF] Noel | Simulating with accord: Actor-based communication via reaction-diffusion[END_REF][START_REF] Arifler | Monte carlo analysis of molecule absorption probabilities in diffusion-based nanoscale communication systems with multiple receivers[END_REF][START_REF] Wang | A novel a priori simulation algorithm for absorbing receivers in diffusionbased molecular communication systems[END_REF]. This effect arises because some molecules close to the receiver are outside of the receiver boundary at two consecutive iteration steps, whereas in reality some portion of them should have been absorbed in-between these steps. In the literature, different absorption models approximating the behaviour of molecules near the receiver boundary (RMC in [START_REF] Arifler | Monte carlo analysis of molecule absorption probabilities in diffusion-based nanoscale communication systems with multiple receivers[END_REF] and APMC in [START_REF] Wang | A novel a priori simulation algorithm for absorbing receivers in diffusionbased molecular communication systems[END_REF]) are proposed to mitigate the intra-step absorption problem. In general, these methods perform random variable (RV) calculations to accurately describe the absorption probabilities of the molecules while decreasing the required iteration number. RMC requires that the surface of the receiver can be approximated by a flat surface, whereas APMC remains accurate only for large step sizes. Both methods work quite well in their respective regime of applicability at the expense of requiring an additional path-integral-like RV calculation which adds another layer of complexity to the simulation algorithm. However, intra-step absorption is only one of the problems which arise from the finite step size of MC simulations. In this work, we take a different perspective on the general problem of finite step size. As a first step, we invoke the principle of locality, which asserts that the state of an object can only be affected through its interactions with its nearby surroundings. While being utilized predominantly in physical contexts, such as field theory and relativity, the principle of locality is also pivotal to the Effective Geometry Monte Carlo simulation method (EG-MC) we develop in this article, where we exploit the locality of the absorption of molecules by a receiver. Building up on this intuition, instead of changing the diffusion of molecules near the receiver boundary we introduce an effective absorption surface such that the mean absorption position of the molecules becomes the receiver boundary. By doing so, we both solve the intra-step absorption problem described in the literature [START_REF] Noel | Simulating with accord: Actor-based communication via reaction-diffusion[END_REF][START_REF] Arifler | Monte carlo analysis of molecule absorption probabilities in diffusion-based nanoscale communication systems with multiple receivers[END_REF][START_REF] Wang | A novel a priori simulation algorithm for absorbing receivers in diffusionbased molecular communication systems[END_REF] and the negative absorption index problem (NAI) which we introduce below. The range of applicability of our method is constrained by the requirement of locality that the step size of the simulation is less than the distance between the receiver and the transmitter.

We organize this paper as follows. We first motivate the negative index problem and introduce the EG-MC algorithm in 1D while showing a way of fixing the free parameter that describes EG-MC. We then apply the EG-MC algorithm to a spherical absorbing receiver described in [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF]. Afterwards, we perform an error analysis to show the consistency of our results in different regimes, where we also illustrate the Poisson noise in the simulation. Finally, we conclude with the range of applicability and a quantitative measure of locality in the absorption process.

Problem Motivation: 1D Toy Model

We shall consider unbounded diffusion channels of arbitrary dimension first. In d-dimensional diffusion channels, the diffusion process is usually modelled as a Brownian motion with a step size

∆x i = N (0, 2D∆t), i = 1, 2, .., d, (1)
where D is the diffusion coefficient, ∆x i (i = 1, 2, .., d) are the incremental step sizes in the d-dimensions, ∆t is time step, and N (µ, σ 2 ) is the normal distribution with mean µ and variance σ 2 . The problems arising from the discrete time steps ∆t when diffusing molecules are encountering a boundary are illustrated in Fig. 1. Besides the well-known intra-step absorption problem, we introduce the novel concept of the negative absorption index problem: Definition (NAI Problem). We define the absorption index for 1D as where x absorbed is the absorbed position of the molecule, r x is the position of the receiver in 1D (r x = 0 unless otherwise stated) and denotes the average value. Finite step size simulations will always obtain AI 1D < 0, which we call the negative absorption index (NAI) problem.

AI 1D = x exp = x absorbed -r x (2)
Note that this definition of AI is meaningful only in 1D and an analogous definition in higher dimensions requires the surface to be approximately flat on scales comparable to the step size. For the scope of this paper, we deal only with this definition and the results that follow this intuition while leaving a more advanced treatment of the problem as future work.

Consider a conventional MC algorithm, as illustrated in Fig. 2. It is clear that the NAI problem will delay the reception of molecules at the receiver's boundary and therefore deteriorate the simulation's accuracy.

Conversely, a solution of the NAI problem would yield precise results for 1D simulations. This is where the idea of locality becomes significant. If the step size is smaller than the global geometry of the diffusion channel, the absorption process can be considered local. Consequently, an individual molecule's position at distances comparable to the step size from the receiver's boundary can be approximated as independent of the global geometry of the system, hence the local problem has as its only relevant length scale the simulation step size given by √ (D∆t). From this it follows immediately that formally shifting the receiver's boundary such that r x → r x +α √ (D∆t) as illustrated in Fig. 3 solves the NAI problem if the free parameter α is suitably chosen.

We conclude this section by illustrating the locality property of the receiver. To this aim, we perform simulations of 1D diffusion from a source of molecules 1 situated a distance L from the origin, where the receiver 1 All simulations in this work were performed with 10 5 particles. is located. We use 100 iteration steps. In Fig. 4, x exp / √ (D∆t) is plotted against the free parameter α. The shape of this plot is independent of the parameters D and L2 , and the interception of the linear graph with the line x exp = 0 yields the optimum choice for the free parameter α = 0.8235 ± 0.0005 in 1D. In Appendix B we show that an approximate calculation yields a similar value analytically.

Application: Spherical Absorbing Receiver

Having motivated the EG-MC algorithm from one-dimensional arguments, we now proceed to apply it to the practically more useful case of a 3D unbounded diffusion channel. As pointed out before, it must be borne in mind that although in 1D the NAI problem was the only systematic bias, in higher dimensions the intra-step absorption problem plays a role as well. We shall show that the EG-MC algorithm still gives results in excellent agreement with the analytic prediction for large step sizes.

We choose our geometry such that a spherical absorbing receiver of radius R is centred at the origin and receives particles from a point source situated a distance L away from the origin, as illustrated in Fig. 5(a).

In order to benchmark our simulation results, we use the analytic expressions available for this problem. They have already been derived in [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF] in quite a tedious fashion. We simplify the calculation considerably by introducing the generalized method of images, which we illustrate in Fig. 5(b). The complete derivation is given in Appendix A. It yields for the hitting rate, defined as the particle flux through the boundary

n hit (t) = R L L -R t √ 4πDt exp - (R -L) 2 4Dt , (3)
while the fraction of molecules absorbed by the receiver until time t is given by

N tot (t) = t 0 n hit (τ ) dτ = R L erfc L -R √ 4Dt . (4)
We now proceed by applying the EG-MC algorithm (Algorithm 1) to the 3D diffusion channel described above. The algorithm has the same structure as the Monte Carlo simulations, with the exception of line 7. Here, unlike usual Monte Carlo algorithms, we add a small thickness to the spherical receiver using the intuition we derived from the 1D toy model in the previous section. From now on, we denote the Monte Carlo simulations as the α = 0 limit of our algorithm. Moreover, as ∆t → 0, our algorithm converges to the Monte Carlo algorithm, as expected.

For this section, we pick α = 0.8235 for illustrative purposes (and shall explore the constraints on the parameter in Sec. 4). A comparison between simulation and analytical results is given in Figs. 6 and7. As can be seen from the figures, EG-MC outperforms the Monte Carlo algorithm significantly. Moreover, our algorithm gives robust results for far larger step sizes, hence decreasing the computational cost extremely. For example, in Fig. 7, the EG-MC algorithm converges in less than 300 iterations (we find that 20-30 iterations are the empirical cut-off for these cases where EG-MC performs quite well) while the Monte Carlo algorithm does not converge even for 10000 iterations. The comparison of n hit (t) obtained from both simulations with the analytical result. As can be seen, the EG-MC algorithm produces a result which follows the analytical curve closely. For both cases, the iteration number is taken as 300 and the step size can be calculated as ∆t = t f /iteration, where t f is the final time of the simulation.

Algorithm 1 Effective Geometry Monte Carlo (EG-MC)

1: procedure Absorption Detection(x i for i=1,2,3)

2:

Define iteration number, initial position and simulation parameters (D,L,∆t)

3: t f = ∆t • iteration 4:
for i = 1:iteration do 5:

Perform Diffusion Simulation

6: r = x 2 + y 2 + z 2 7: index=find(r < r x + α √ D∆t) 8:
N abs (i)= length(index) Before moving on to the error analysis, we shall first define some short-hand notations. From now on, we denote the discretized values by a square bracket, e.g. N [•]. Furthermore, we denote the number of particles absorbed between times t i and t i + ∆t as

∆N tot [t i ] = N p ti+∆t ti n hit (τ ) dτ. (5)
Here, N p denotes the number of particles initially released by the transmitter.

Error and Efficiency Analysis

In this section we discuss the uncertainty on the parameter α introduced above. To this aim, we compare the simulation result for a given α with the analytic result.

First, we define the integrated squared difference in the cumulative distribution (ISDCD) as

ISDCD = i (N sim tot [t i ] -N anl tot [t i ]) 2 , (6)
where the sum runs over all data points in a given simulation, and "sim" and "anl" denote the simulated and analytical results respectively. We choose this quantity as our primary measure for the accuracy of a given simulation for the following reasons: The cumulative distribution is a better measure than the bare probability density function (hit rate), since it mitigates the random noise in the simulation. Instead of comparing the cumulative distribution at a certain point in time, we use the sum, because it is less susceptible to random fluctuations, which are of particular importance in the steep increase of the cumulative distribution for small times.

In Fig. 8 we show the dependence of the ISDCD on the choice of the parameter α. It can clearly be seen that the ISDCD has a well defined minimum around α ∼ 0.8, and follows a quadratic dependence as confirmed by the fit.

Utilizing this behaviour, we can estimate the optimum value of α for different choices of simulation parameters (R, L, D) and different iteration numbers. This is easily done by retrieving the ISDCD, fitting it with a parabola and retrieving the value α that minimizes the ISDCD. A measure for the error on α is obtained by repeating this procedure (20 times) and taking the standard deviation of the results.

It may be argued that the better known reduced chi-squared χ 2 red of the simulated hitting rate with respect to the analytic hitting rate is a better estimator for the goodness of the simulation, since it directly describes the simulated data. However, as argued before, this quantity is very susceptible to random fluctuations. Nonetheless, we demonstrate that the optimal choice of α from the ISDCD corresponds to expected values of χ 2 red . We define We computed the error bars approximating the hitting rate as normally distributed. 4 The red line is the theoretical expectation for Poisson noise

χ 2 red = 1 iteration -1 i ∆N sim tot [t i ] -∆N anl tot [t i ] σ[t i ] 2 , (7)
√ (∆N anl tot [t i ]).
For clarity, only every fifth datapoint is shown.

where we take ∆N tot as the bare count rate in our simulation (i.e. it is integer valued and not normalized as in the previous sections).

We expect the random error σ[t i ] to be close to Poisson noise. 3 We confirm this by running the simulation repeatedly, retrieving for each datapoint the difference to the analytic result, and taking the standard deviation. The result is shown in Fig. 9, and confirms that we can set σ

2 [t i ] = ∆N anl tot [t i ]
. This is in agreement with the various models proposed in the literature [START_REF] Arjmandi | Diffusion-based nanonetworking: A new modulation technique and performance analysis[END_REF][START_REF] Etemadi | Compound poisson noise sources in diffusion-based molecular communication[END_REF].

In Fig. 10 we show the results for the optimum value of α as estimated from the ISDCD, together with the corresponding χ 2 red values from ∆N tot [t i ]. We first discuss the behaviour of the optimal parameter α O for different simulation inputs. It is clear that α O varies only very slightly, although the simulations were performed with very different input parameters. This confirms our expectation that to a good approximation we can treat α O as a constant. However, the deviations from the weighted mean are larger than expected from the errorbars. This is mainly due to a systematic effect at low iteration numbers, as for those the local problem at the receiver's boundary begins to loose its one-dimensional character which we invoked to motivate the EG-MC method. In general, we find the optimum α to be slightly lower in the 3D case in comparison to the 1D case, which we attribute to the fact that the α-correction increases the total number of particles received integrated over all time (according to R/L), while in 1D all particles are absorbed.

Next we discuss the behavior of χ 2 red . It is well known that an optimal fit is characterized by χ 2 red = 1, and it can be seen that this value is obtained -within statistical errors -for each simulation shown. We also list the values of χ 2 red for identical simulation parameters if α is set to 0. It is clear from these values that the EG-MC method outperforms a naive implementation by a large margin for iteration numbers of order 10 2 , while for larger iteration numbers the chi-squared is not a good measure for the goodness of fit due to the large standard deviation stemming from the Poisson noise.

After showing that the EG-MC method results in a better fit than the MC algorithm, we can now focus on comparing both algorithms in terms of efficiency. Before starting the comparison, we emphasize once more that the implementation of EG-MC algorithm is identical to the MC algorithm, with the only difference being the non-zero free parameter α in Algorithm 1. For the rest of this section, we will perform chi-squared and relative error analyses to assess the correctness of the simulations, which are performed with increasing iteration number. As a first order approximation, the simulation time linearly increases with the iteration number of the simulation. 10 5 particles are released from the transmitter for each simulation, whereas the time interval of the simulations is taken as t ∈ [0, t ch ]. Here, we define the characteristic time scale for each scenario as

t ch = (R -L) 2 D . ( 8
)
This definition is particularly important for comparison of different scenarios, as the time-scale should be The grey triangles indicate the value of χ 2 red obtained with this choice of α. We also list the values of χ 2 red for a simulation in which α = 0 at the bottom of the plot. The average α is weighted with the square of the individual standard deviations. The values (R, L, D, iter) are given in units of (µm, µm, µm 2 /s, 10 2 ). equivalent in all cases. We note that t ch corresponds to 6τ peak , where τ peak is the time n hit (t) reaches its maximum [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF]. Therefore, the simulations are space and time scale-invariant for the time interval [0, t ch ].

We compare the results of the chi-squared analysis for both algorithms with increasing iteration number in Fig. 11. Chi-squared analysis gives us the minimum iteration number after which the systematic error per time slot caused by the algorithm is less than the random error. Due to the scale invariance of the simulations, the χ 2 red values follow the same pattern with increasing iteration number (hence lower step size) as expected. The accuracy of the algorithm depends on the only scale-free parameter, which is the step-size adjusted with respect to the channel parameters (hence the iteration number). As can be seen from the figures, the iteration number required by MC algorithm to have similar accuracy as the EG-MC algorithm is around 1000. This shows that for the systematic error per time slot to be smaller than the random error, MC algorithm requires an order of magnitude larger simulation time with respect to the EGMC algorithm. While systematic error per time slot might be smaller than the random error, the algorithm can still make incorrect predictions for metrics where random error is integrated out. This is particularly clear from N tot (t) plotted in each subfigure, as the MC algorithm still cannot predict the correct values. The reason is as following: While the systematic error in n hit (t) is smaller than the random error, integrating out the random error to obtain N tot (t) accumulates the systematic error, resulting in the inaccurate prediction. However, the expected fraction of particles absorbed by the receiver is an important channel parameter which should be predicted accurately.

In order to test the algorithms in terms of their prediction of parameters where random errors are integrated out and systematic errors are accumulated, we propose the following metric, called the relative error:

ε rel = 1 iteration tmax t=∆t |N A tot (t) -N S tot (t)| N A tot (t) . ( 9
)
Here, S/A stands for simulation/analytic result. In implementing this metric, we occasionally use an offset, starting the sum from offset × ∆t instead of ∆t and adjusting the divison factor, as for the early times the random error still dominates for this metric. The result of the relative error analysis is shown in Fig. 12. As can be seen from the figure, the MC algorithm matches the accuracy of the EGMC algorithm around 30000 iterations, whereas the EGMC algorithm makes accurate predictions even for 300 iterations. This shows that the EGMC algorithm requires 2 orders of magnitude less iteration to converge to the analytic result. As to the first order approximation, simulation time linearly scales with the iteration number, the EGMC algorithm outperforms the MC algorithm while being more efficient in terms of simulation time by 2 orders of magnitude.

Range of Applicability and Application to Complex Structures

The last topic we shall address is the range of applicability of the EG-MC algorithm. As we have already discussed, the principle of locality dictates that √ (2D∆t) O(L). In other words, the step size should be small enough such that the molecules are not absorbed after only a few steps. To illustrate this concept, we define a comparison metric called relative inaccuracy as follows Relative inaccuracy = ISDCD(α = 0.8235) ISDCD(α = 0) .

Here ISDCD(α) is the integrated squared difference in the cumulative distribution defined in Eq. ( 6), which we evaluate for the EG-MC algorithm and compare it with a simple Monte Carlo implementation. In Fig. 13, we show the relative inaccuracy in dependence on step size for various choices of (R, L, D). It shows that our algorithm performs very well even for step sizes larger than the diameter of the receiver provided that a cutoff value related to the locality condition is not exceeded. In our calculations, we find that the cutoff condition can be given more precisely as

√ 2D∆t L -R. (11)
We can clearly see in Fig. 13 that this cutoff value is observed for the four different cases we present in this paper.

So far, we have focused on the case with a spherical absorbing receiver and a point transmitter. Now, we shall generalize our findings to more general systems with multiple transmitters and receivers having arbitrary shapes. We will start by showing that the EGMC algorithm converges to the MC algorithm for the time steps ∆t, for which MC algorithm is reliable. Let us re-consider Figs. 11 and12. The time step in these simulations can be given as

∆t = t ch iteration = (R -L) 2 D × iteration . ( 12
)
The incremental increase in the spherical receiver radius can be given as α √ D∆t = α(R -L)/ √ iteration. Then, we can find the fractional increase, ∆R, in the receiver radius as

∆R = α R -L R √ iteration ∼ 10 -2 , (13)
where we realize that α R-L R ∼ O(1) and iteration ∼ O(10 4 ) for the case where MC algorithm performs accurately. In the regime where MC algorithm is reliable, EGMC algorithm requires a single percent increase in the spherical radius, which can easily be neglected. This shows that as a lower bound, the EGMC algorithm performs as good as the MC algorithm in any complex structure.

Then, we argue using locality principle that the EGMC algorithm retains its superiority to the MC algorithm for a large group of complex structures. Let us assume that there are many receivers inside a channel. For the EGMC algorithm to be accurate, we require that the position step size, √ 2D∆t, in each iteration should be smaller than the distance smallest distance between the receivers. If this requirement is satisfied, changing the radius of one receiver does not affect the diffusion of molecules close to the other receivers5 , while correcting the systematic error caused by the MC algorithm around the original receiver. Since locality condition given in [START_REF] Yilmaz | Simulation study of molecular communication systems with an absorbing receiver: Modulation and isi mitigation techniques[END_REF] is less strict than the high iteration number condition of the MC algorithm, the EGMC algorithm is expected to outperform MC algorithm in any complex structure. This shows that using the EGMC algorithm is always preferred over using MC algorithm.

Finally, we conclude with specific examples where hte EGMC algorithm has been tested. In [START_REF] Dinc | Analytical estimation for the impulse response of an n-dimensional diffusion channel with an absorbing receiver[END_REF], the EGMC algorithm is applied for spherical absorbing receivers in an n-dimensional diffusion channel. Furthermore, in [Accepted article, in production], it is shown that MC algorithm leads to a systematic shift in the absorption angle of molecules. The EGMC algorithm fixes this problem, which shows that the EGMC algorithm truly solves the intra-step absorption problem 6 . Moreover, in this work, it is shown that the EGMC algorithm works accurately in the presence of constant flow as well, laying evidence for our claim that the EGMC algorithm can perform accurately in complex channels.

Conclusion

In this paper, we have shown that an effective geometry approach can lead to significantly faster and more reliable simulations of diffusion in molecular communication. As the complete method solely consists of including a small correction -α √ (D∆t) -to the receiver radius, it decreases the computational cost significantly. In fact, we have shown, in Fig. 7, that while EG-MC converges for iteration numbers ∼ 300, the Monte Carlo simulation does not converge to the analytical solution even for ∼ 10000 in the case of a spherical receiver.

From this perspective, our method outperforms Monte Carlo simulations by more than two orders of magnitude in terms of computation time.

We also note that the possibility of choosing a relatively large step size mitigates the problem of Poisson noise considerably, which implies that the simulation can be run with a smaller number of particles if identical noise levels for the hitting rate are to be obtained.

We believe that the EG-MC algorithm we present in this paper is just a first step towards a deeper understanding of effective geometries in the optimization of Monte Carlo type simulations. Similar algorithms can be implemented in different fields using Monte Carlo simulations where computational cost is of primary importance. As future work, we plan to generalize our findings to higher dimensional spaces while applying EG-MC algorithm for complex systems consisting of multiple receivers/transmitters with arbitrary shapes. Through the computational efficiency this method introduces, many systems deemed too complex to be simulated in reasonable time intervals should be considered for further research directions.

symmetry. This symmetry-motivated method is a more general version of the method of images which is used commonly in the literature. Hence, we obtain the probability distribution function for the molecule as

P (r, t) = exp -(r-L) 2 4Dt -exp -(r+L-2R) 2 4Dt 4πrL √ 4πDt . (18)
Then the hitting rate, defined as the particle flux through the boundary, is given by

n hit (t) = 4πR 2 D ∂ r P (r, t)| r=R = R L L -R t √ 4πDt exp - (R -L) 2 4Dt . ( 19
)
The fraction of molecules absorbed by the receiver until time t follows as

N tot (t) = t 0 n hit (τ ) dτ = R L erfc L -R √ 4Dt . (20)

B Analytical Approximation for α

In this appendix, we are deriving an approximate expression for α in the 1D case. We choose our geometry such that the receiver is located at x = 0, and the effective receiver at α . Let us consider a single time step in the simulation. At the beginning of the step, a particle is located at x 0 > α . We first find the expected absorption position x 1 (x 0 ) for this individual particle under the assumption that it is absorbed. We choose our units of length such that the probability density function for the location of this particle after the time step is given by P (x,

x 0 ) = N e -(x-x0) 2 , ( 21
)
where we omit normalization factors and lump them in N because they are of no importance for the calculation and would only serve to clutter the expressions. Hence we find that

x 1 (x 0 ) = α -∞ dx xP (x, x 0 ) α -∞ dx P (x, x 0 ) . ( 22
)
The probability of absorption is given by

P abs (x 0 ) = α -∞ dx P (x, x 0 ). ( 23
)
Denoting the instantaneous distribution of particles by ρ(x 0 ), the expected absorption position can then be given by using the probability distribution (P abs (x 0 )ρ(x 0 ) dx 0 ) as x abs = ∞ α dx 0 x 1 (x 0 ) P abs (x 0 )ρ(x 0 ) ∞ α dx 0 P abs (x 0 )ρ(x 0 ) ,

We now want to find α such that the expected absorption position (for all particles present at x > α ) equals zero. It follows that 0 = ∞ α dx 0 x 1 (x 0 ) P abs (x 0 )ρ(x 0 ),

We now invoke the approximation that the distribution ρ(x 0 ) is uniform. This approximation will yield a relatively precise result for α because only particles relatively close to the receiver have an appreciable chance of absorption, which justifies a zeroth order approximation.

Inserting the expressions given before, we find

We simplify the final integration using ξ ≡ x 0 -α , so that

0 = ∞ 0 dξ e -ξ 2 + ξ √ πerfc(ξ) + α √ πerfc(ξ) , (28)
from which it follows that α = Returning to the units of length used in the main body of the text, we find α = 2α = √ π/2 0.89. Hence, we see that using the approximation of a uniform ρ(x 0 ), a result within 7% of the empirically found value can be obtained.

Figure 1 :Figure 2 :

 12 Figure 1: An illustration of various problems resulting from finite step size. (i),(iii): Intra-step absorption problem, (ii): Negative absorption index (NAI) problem

Figure 3 :

 3 Figure 3: An Illustration of EG-MC Algorithm both in the general 3D case and in the 1D case. A new boundary is created such that the molecules are absorbed right at the boundary by the algorithm. The main motivation behind the shift of the boundary is the principle of locality, which results in linear relation between the free parameter α and the absorption position x exp (See Fig. 4).

Figure 4 :

 4 Figure 4: The linear fit describes the relationship between α and x exp (in 1/ √ D∆t) quite perfectly as expected. For this plot, a range of L (∈ [30, 200]µm) and D (∈ [80, 600]µm 2 /s) values are used and the mean values with error bars are plotted. The error bars are much smaller than the marker size. According to multiple trials with different parameters, the most optimum free parameter is α = 0.8235 ± 0.0005. The standard deviation is obtained from the linear fitting performed in this plot.

Figure 5 :

 5 Figure 5: (a) The 3D unbounded diffusion channel consisting of a point transmitter and a spherical receiver. (b) The 1D line, where the real space is projected onto the positive values and an imaginary space added to the negative values. In this pseudo-real (consisting of real and extended parts) space, we include a sink at r = 2R -L to use a more general version of method of images.

Figure 6 :

 6 Figure 6: (a) The proof of concept that EG-MC algorithm produces correct results for various transmitter distance L, where D = 80 µm 2 /s and R = 10µm. (b)The comparison of n hit (t) obtained from both simulations with the analytical result. As can be seen, the EG-MC algorithm produces a result which follows the analytical curve closely. For both cases, the iteration number is taken as 300 and the step size can be calculated as ∆t = t f /iteration, where t f is the final time of the simulation.

Figure 7 :

 7 Figure 7: The comparison of Monte Carlo simulation (a) and EG-MC simulation (b) with the analytical solution for different D values, where L = 30 µm and R = 10µm. For each individual case, the step size is different and therefore we note the value of total iteration. The step size can be calculated as ∆t = t f /iteration, where t f is the final time of the simulation.

9 :Figure 8 :

 98 Figure 8: Dependence of IS-DCD as defined in Eq. (6) on α for L = 35 µm, R = 10 µm, D = 80 µm 2 /s. The simulation was run with 100 iterations for each choice of α.Errorbars were obtained by running the simulation 20 times for each α and taking σ ISDCD = √ (Var(ISDCD)/20). The fit is a simple parabola.

Figure 9 :

 9 Figure 9: Standard deviation of the hitting rate taken from 30 simulations in dependence on time. The simulations were performed with L = 35 µm, R = 10 µm, D = 80 µm 2 /s, N = 10 3 iterations.We computed the error bars approximating the hitting rate as normally distributed.4 The red line is the theoretical expectation for Poisson noise√ (∆N anl tot [t i ]).For clarity, only every fifth datapoint is shown.

Figure 10 :

 10 Figure 10: Blue datapoints are the optimal value of α for different input parameters, estimated from the IS-DCD as described in the text. The grey triangles indicate the value of χ 2red obtained with this choice of α. We also list the values of χ 2 red for a simulation in which α = 0 at the bottom of the plot. The average α is weighted with the square of the individual standard deviations. The values (R, L, D, iter) are given in units of (µm, µm, µm 2 /s, 10 2 ).

Figure 11 :

 11 Figure11: The reduced chi-squared analysis to test the accuracy of the EGMC and MC algorithm. The final time is t ch for all cases and the simulation time step can be found by ∆t = t ch /iteration.

Figure 12 :

 12 Figure 12: The relative error analysis to test the accuracy of the EGMC and MC algorithm. The final time is t ch for all cases and the simulation time step can be found by ∆t = t ch /iteration.

Figure 13 :

 13 Figure 13: The relative inaccuracy between the Monte Carlo simulations and the EG-MC simulations for various step sizes. The iteration number is chosen such that either the final time is 6× the peak value of n hit (t) or at least 100 iterations to allow more accurate comparison between different cases. EG-MC algorithm clearly outperforms MC simulations for various step sizes. More importantly, EG-MC algorithm results in negligible error as long as √ (2D∆t) ≤ L -R, i.e. the locality assumption is valid.

  of variables χ ≡ x -x 0 in the second integral, dχ χe -χ 2 + x 0 e -χ 2 = ∞ α dx 0 1 2 e -(x0-α ) 2 + x 0 √ π 2 erfc(x 0 -α ) .

This has been found to hold for D ∈ [80, 600]µm 2 /s and L ∈ [30, 200]µm.

There is a small correction due to correlations between different time steps, as a high absorption rate in one time step in general leads to a lower absorption rate in the consecutive step. However, this effect is found to be very small.

According to the formula Var(σ) = σ 2 {1 -2Γ 2 (N/2)/[(N -1)Γ 2 ((N -1)/2)]}.

By this, we mean that the interference probability between the two receivers is suppressed by a Gaussian according to the distance between them and the time step ∆t.

This is because the EGMC algorithm predicts both the hitting rate and the absorption angle accurately, meaning that it mimics the diffusion process completely.

Acknowledgements

Research at the Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada, and by the Province of Ontario through the Ministry of Research and Innovation. LT acknowledges support by the Studienstiftung des Deutschen Volkes. FD and LT thank Matija Medvidovic for suggestions improving the code framework and Lauren Hayward Sierens for suggestions on the manuscript and error analysis. FD and BCA would like to thank Prof. Ali Emre Pusane and Prof. Tuna Tugcu for intellectually stimulating discussions on the molecular communication application of the EG-MC algorithm.

A Derivation of the Hitting Rate via Generalized Method of Images

In this appendix, we derive the analytic expression for the hitting rate n hit (t) and the integrated fraction of particles absorbed N tot (t) in the geometry described in Sec. [START_REF] Nakano | Channel model and capacity analysis of molecular communication with brownian motion[END_REF].

The impulse response for this diffusion channel has already been derived in [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF]. Here, we shall re-derive their findings while motivating the solution by symmetry arguments rather than lengthy calculations.

The diffusion of molecules in 3D is described by Fick's Law

where ∇ 2 is the Laplacian operator, P (r, t) is the probability distribution function of a diffusing molecule, D the diffusion coefficient defined in Eq. ( 1) and ∂ t is the partial derivative with respect to t.

The receiver absorbs all the molecules upon incident. Moreover, the molecules are initially localized at the position (0, 0, L). These two conditions specify a unique solution for the diffusion equation. Nonetheless, before proceeding, we shall exploit the SO(3) symmetry. The absorption angle of the molecules is dependent on the initial position of the transmitter, whereas the absorption rate is not. As we are only interested in the absorption rate, we can solve an equivalent but easier initial value problem in which the transmitter position is symmetrised over all angles. Then, the boundary conditions become

The calculation presented in [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF] derives its substantial length from the boundary condition in Eq. (15b). To avoid this, let us first write down the solution to the diffusion equation subject only to Eq. (15a), which is describing a freely diffusing particle [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF]

One important realization regarding this solution is that L is a free parameter, i.e. the diffusion equation is satisfied regardless of its choice, while L ensures that the boundary condition Eq. (15a) holds.

As the system is symmetric around all angles, we can project it onto a line stretching from 0 to ∞. This transformation maps 4πr 2 dr points onto dr points where r ∈ [0, ∞) is the parameter describing the real line. Usually this transformation is irreversible as it is not one-to-one. Nonetheless, by the SO(3) symmetry, no information is lost during the transformation. With this in mind, we can project the complete 3D space into an incomplete 1D real space r ∈ [0, ∞) with r = 0 describing the origin of the original space. Then, we complete this space by adding the negative values for r, as shown in Fig. 5(b), which do not have any physical meaning for the 3D space but are useful in solving the diffusion equation. As a final step, we include a sink at r = 2R -L with

while recalling that the diffusion space is r ≥ R. From the uniqueness theorem it follows that P (r, t) = P 0 (r, t) + P S (r, t) is the desired solution, as it satisfies the boundary conditions and the diffusion equation. Now, we can transform back to 3D space, where both P S (r, t) and P 0 (r, t) remain invariant due to the SO(3)


Student Memb, IEEE Fatih Dinc

Student Member, IEEE Bayram Cevdet Akdeniz

Member, IEEE Ali Emre Pusane

Member, IEEE Tuna Tugcu

Bayram Cevdet

B C Akdeniz

In this paper, the impulse response for a 3-D microfluidic channel in the presence of Poiseuille flow is obtained by solving the diffusion equation in radial coordinates. Using the radial distribution, the axial distribution is then approximated accordingly. Since Poiseuille flow velocity changes with radial position, molecules have different axial properties for different radial distributions. We, therefore, present a piecewise function for the axial distribution of the molecules in the channel considering this radial distribution. We lay evidence for our theoretical derivations for impulse response of the microfluidic channel and radial distribution of molecules through comparing them using various Monte Carlo simulations. Finally, the communication performance of the channel is examined. Index Terms-3-D microfluidic channel with flow, nonuniform diffusion, Poiseuille flow, impulse response of the microfluidic channel, molecular communication.

I. INTRODUCTION

M OLECULAR communication via diffusion (MCvD) is one of the most promising areas for nanonetworking due to its potential biocompatability. It is based on encoding information symbols by releasing messenger molecules into a fluidic environment. Released molecules diffuse through the environment under Brownian motion and the receiver makes a decision on the transmitted symbols by observing or absorbing the released molecules.

There are several different channel models proposed in the molecular communication literature. An extensive survey that involves the compilation of these channel models is presented

A General Analytical Approximation to Impulse Response of 3-D Microfluidic Channels in Molecular Communication

in [START_REF] Farsad | A comprehensive survey of recent advancements in molecular communication[END_REF]. Despite the large number of different diffusion channels in the literature, they all have one common problem due to the nature of the diffusion: inter symbol interference (ISI). Since the movement of molecules are slow and random in Brownian motion, some of the released molecules may not reach to the destination until the desired time. This possibly leads to an adverse effect on decoding. There are many modulation and equalization methods to eliminate the molecules that cause ISI [START_REF] Arjmandi | Diffusionbased nanonetworking: A new modulation technique and performance analysis[END_REF]- [START_REF] Akdeniz | Optimal reception delay in diffusion-based molecular communication[END_REF]. In addition to these methods, channel models that diminish ISI have been proposed by considering the reasons that cause ISI. It is clear that the dispersion and slow movement of the molecules can cause ISI in microfluidic channels. Especially in unbounded environments, the molecules are uniformly dispersed in the space and hence, the number of molecules reaching the receiver decreases. Therefore, using barriers that form a bounded channel can be a reasonable approach to keep the molecules closer to the destination and is a more realistic channel type considering biomedical applications. As proposed in [START_REF] Giné | Molecular communication options for long range nanonetworks[END_REF], vessel-like structures are good candidates for long range molecular communication since they preserve released molecules in a guided range. Another beneficial factor in molecular communication channel for reducing ISI is flow which increases the speed of the molecules [START_REF] Farsad | A comprehensive survey of recent advancements in molecular communication[END_REF]. Therefore, using microfluidic channels assisted by flow does not only diminish ISI but also increases the data rate.

There are various works related to microfluidic or vessellike channels. In [START_REF] Felicetti | Establishing digital molecular communications in blood vessels[END_REF], a decoding scheme for molecular communications in blood vessels has been proposed. In [START_REF] Bicen | Molecular transport in microfluidic channels for flow-induced molecular communication[END_REF], a rectangular microfluidic channel with flow has been modeled and analyzed. In [START_REF] Sun | Channel capacity modelling of blood capillary-based molecular communication with blood flow drift[END_REF], numerical capacity analysis of the vessel like molecular channel with flow is examined. In [START_REF] Turan | Channel model of molecular communication via diffusion in a vessellike environment considering a partially covering receiver[END_REF], a partially covering receiver in a vessel-like channel without flow is examined, and the channel characteristics are analyzed. Although the analytical channel impulse responses for unbounded channels are derived, the channel impulse response for microfluidic channel with Poiseuille flow has not been derived yet in the molecular communication literature except for some special cases. In [START_REF] Bicen | System-theoretic analysis and leastsquares design of microfluidic channels for flow-induced molecular communication[END_REF], the flow models of microfluidic channels with different cross-section area are presented and the impulse response is derived by solving a 1-D diffusionadvenction equation, which is only valid for some specific cases. In [START_REF] Wicke | Modeling duct flow for molecular communication[END_REF], channel impulse responses of point and planar transmitter in a microfluidic channel have been derived for dispersion and flow dominated cases. For the first case, radial distribution of the molecules is assumed to be uniform and the system is reduced to 1-D to solve the channel response. For the latter case, only the flow is considered by neglecting the effect of diffusion and the channel impulse response is obtained accordingly. Although that paper includes an elegant and extensive work for microfluidic channels with Poiseuille flow, it is only valid for the uniform radial distribution or flow dominated regions. This assumptions occurs, if Peclet number, a unitless number that compares the effects of flow and diffusion, is much higher or much lower than the ratio of the radius of the channel and the distance between the transmitter and receiver. Therefore, for other cases, derivation of the channel impulse response still remains as an open problem. In this paper, the analytical channel impulse response of microfluidic channel that involves Poiseuille flow is derived when an arbitrarily placed point transmitter and a planar observing receiver that fully covers the cross-section of the channel. In order to determine this function, firstly the analytical formula of the radial distribution of the molecules in time is derived. Using this distribution, the average velocity, displacement of a molecule and finally the probability of observation of a molecule by the receiver (which can also be regarded as channel impulse response) are obtained. All these functions enable us to examine a microfluidic channel without any real time simulation. Furthermore, using these functions, some channel properties like required time to reach uniform radial distribution in channel and required time for error-free communication is extracted.

The main contributions of this paper to the literature can be listed as follows:

• Derivation of the radial distribution of the released molecules in a vessel-like molecular communication channel under the effects of Poiseuille flow and diffusion.

• A piecewise approximation of the axial distribution of the released molecules for all cases. • Analytical derivation of the required communication interval for error-free communication

II. SYSTEM MODEL

The considered system model is depicted in Fig. 1. For the sake of simplicity, the coordinates in the channel are defined using cylindrical coordinates (r, x), where r = y 2 + z 2 ∈ [0, d] and x ∈ [-∞, ∞]. In this figure, a point transmitter placed at an arbitrary axial (x direction) and radial (r direction) position denoted by d 0 , and a circular planar observing receiver is placed x r away from the axial position of the point transmitter. The boundaries of the microfluidic channel are reflecting and there is a Poiseuille flow that changes with the axial position as

v(r ) = v m 1 - r 2 d 2 , ( 1
)
where v m is the maximum velocity that occurs at the center of the microfluidic channel. The planar observer receiver observes the number of molecules passing through its surface and makes a decision based on its observations without removing the messenger molecules from the environment. In addition to the flow, diffusion also exists in the channel and the diffusion process is modeled using the diffusion coefficient D, which is related to the variance of the Brownian motion. Since the flow is only available on the x axis, for each t seconds, the displacement of a molecule at a radial position r can be modeled in Cartesian coordinates as

x = tv(r ) + N (0, 2Dt) , y = N (0, 2Dt) , z = N (0, 2Dt) . ( 2
)
From (2), one can easily observe that the distribution of y and z (hence, radial displacement r ) are not dependent on x, but distribution of x is dependent on r . In other words, the movement of the molecules along the radial axis is purely diffusive while movement of the molecules along the axial axis is a combination of the diffusion and flow, whose velocity is determined by radial position. In order to compare the effect of flow and diffusion, Peclet number (Pe) is a useful dimensionless metric that can be obtained for the microfluidic channel as

Pe = v m d 2D . ( 3
)
Note that for pure diffusion Pe = 0, and for pure flow (i.e., advection) Pe approaches ∞.

III. CHANNEL IMPULSE RESPONSE

In order to derive the impulse response of the channel, we need to derive the joint radial, axial, and time distribution of the molecules, p (x, r, t|d 0 , d, x r ), in the microfluidic channel. As can be seen in ( 2), the radial distribution is independent of the axial distribution, but the axial distribution is dependent on the radial distribution. Considering this fact, we can rewrite p (x, r, t|d 0 , d, x r ) as

p (x, r, t|d 0 , d, x r ) = p (r, t|d 0 , d, x r ) • p (x, t|r, d 0 , d, x r ) . (4)
Therefore, our aim is finding the radial distribution p (r, t|d 0 , d, x r ) first, and then using this distribution to obtain the axial distribution p (x, t|r, d 0 , d, x r ). Once p (x, r, t|d 0 , d, x r ) and p (r, t|d 0 , d, x r ) are derived, the channel impulse response of the circular observer is obtained as

n hit (t|d 0 , d, x r ) = ∂ ∂t ∞ x r d 0 2π 0 dxdr dθ r p(x, r, t|d 0 , d, x r ).
(

) 5
Having found n hit (t|d 0 , d, x r ), one can easily find the fraction of the observed molecules by the receiver until time t, N hit (t|d 0 , d, x r ), by integrating n hit (t|d 0 , d, x r ) with respect to time as

N hit (t|d 0 , d, x r ) = t 0 n hit (τ |d 0 , d, x r )dτ. (6)
A. Derivation of the Radial Distribution p(r, t |d 0 , d, x r )

The radial distribution can be obtained by simply solving the diffusion equation or by drawing comparison to the heat flow, as discussed in [START_REF] Carslaw | Conduction of Heat in Solids[END_REF]. Here, we shall find the radial distribution of molecules using the former method. To describe the diffusion of the molecule inside the radial region, a solution to the Fick's Law, satisfying the necessary boundary conditions is needed. The equation is given as

D P (r, t|d 0 , d, x r ) = ∂ P (r, t|d 0 , d, x r ) ∂t , (7)
where P (r, t|d 0 , d, x r ) is the probability density of the molecule and is the Laplacian operator. The boundaries are reflecting, meaning that the probability current normal to the boundaries should be zero. Furthermore, the molecule is assumed to be situated at a distance d 0 away from the origin for t = 0, which would imply an initial condition for the molecules as P(r, t = 0|d 0 , d, x r ) = δ(rd 0 )δ(θ ), where θ is the angular coordinate in the y-z plane. Nonetheless, we note that we are ultimately interested in n hit (t|d 0 , d, x r ) and therefore transform our initial conditions in a way that would allow an easier solution with less computation. First, we note that the flow is independent of the angle θ at which the particle is situated at yz-plane, it only depends on r . Moreover, the circular observer observes all particles independent of the angle φ. Thus, finding the SO(2) invariant solution is sufficient enough 1 to find n hit (t|d 0 , d, x r ) and we shall focus on this 1 This simplification has been performed in the molecular communication literature several times, for example [START_REF] Yilmaz | Threedimensional channel characteristics for molecular communications with an absorbing receiver[END_REF], [START_REF] Al-Zu'bi | Modeling of ligand-receptor protein interaction in biodegradable spherical bounded biological micro-environments[END_REF].

problem from now on. The boundary conditions then can be given as

∂ P (r, t|d 0 , d, x r ) ∂r r = d = 0, (8a) P (r, t|d 0 , d, x r ) t = 0 = δ(r -d 0 ) 2πr , ( 8b
)
where we recall that, under Neumann boundary conditions, the Laplacian operator is guaranteed to have a unique solution up to an addition of a constant, which can be regarded as the normalization constant. In order to solve P (r, t|d 0 , d, x r ) separation of variable ansatz is used as

P (r, t|d 0 , d, x r ) = φ(r )T (t). ( 9
)
Here, the dependence on the initial and boundary conditions (related to d 0 , d and x r ) is dropped for brevity. This ansatz leads to the equation

D φ(r ) φ(r ) = Ṫ (t)
T (t) = -μ 2 ,
where Ṫ (t) denotes the derivative of T (t) with respect to t. We can deduce that

T (t) = Ce -μ 2 t ( 10
)
and arrive at the Neumann-eigenvalue problem for the Laplacian operator as

φ(r ) = - μ 2 D φ(r ). ( 11
)
At this point, it is important to recall some key properties of Laplacian operator [START_REF] Grebenkov | Geometrical structure of Laplacian eigenfunctions[END_REF]- [START_REF] Weinberger | An isoperimetric inequality for the N-dimensional free membrane problem[END_REF]. The eigenvalues μ 2 /D are nonnegative and real, as well as the eigenvectors corresponding to distinct eigenvalues being orthogonal and forming a basis for all possible solutions. The non-negativity of the eigenvalues ensures that T (t) does not tend to infinity as t → ∞, whereas the orthogonal basis guarantees a unique solution.

Rewriting the eigenvalue equation in polar coordinates, we obtain:

r 2 φ (r ) + r φ (r ) + μ 2 D r 2 φ(r ) = 0,
where φ (r ) denotes the derivative of φ(r ) with respect to r . The most general solution is

φ(r ) = J 0 μ √ D r + c 2 Y 0 μ √ D r ,
where J 0 and Y 0 are the Bessel's function of the first and second kind, respectively, and of the zeroth order, and c 2 is a constant to be determined by the boundary conditions.

Here, the coefficient of J 0 is chosen arbitrarily as the overall coefficient C of the solution is lumped into the time dependent part in [START_REF] Sun | Channel capacity modelling of blood capillary-based molecular communication with blood flow drift[END_REF]. The solution can now be shaped according to the boundary conditions given in [START_REF] Felicetti | Establishing digital molecular communications in blood vessels[END_REF]. One important observation is that, for t > 0, the probability density function P (r, t|d 0 , d, x r ) does not diverge at the origin of the y-z plane, resulting in c 2 = 0. The most general solution is then of the form:

P (r, t|d 0 , d, x r ) = ∞ n=0 C n J 0 μ n √ D r e -μ 2 n t ,
where C n and μ n are to be specified by the boundary conditions. We arbitrarily define the starting index as n = 0. Once the boundary condition given in (8a) is invoked, we arrive at

∂ P (r, t|d 0 , d, x r ) ∂r r=d = 0 ⇒ J 1 μ n √ D d = 0.
The first order Bessel function J 1 (r ) has infinitely many zeros. These zeros, defined as

β n = μ n √ D
d, then constitute infinitely many terms for the solution. With the notation β 0 = 0, the solution is a sum of infinitely many terms given as

P (r, t|d 0 , d, x r ) = ∞ n=0 C n J 0 β n r d e -μ 2 n t .
Before imposing the final condition given in (8b), it is useful to note the following normalization identity for Bessel functions [START_REF] Asmar | Partial Differential Equations and Boundary Value Problems[END_REF], which is given as

1 0 x J 0 (β n x)J 0 (β m x)dx = 0.5 J 0 (β n ) 2 δ nm ,
where δ nm is the Kronecker delta function. Using this identity for the initial condition given in (8b), we conclude that

C n = J 0 β n d 0 d πd 2 J 2 0 (β n )
, from which we find the final solution to be

P (r, t|d 0 , d, x r ) = ∞ n=0 J 0 β n d 0 d πd 2 J 2 0 (β n ) J 0 β n d r e -β 2 n Dt d 2 , ( 12
)
where the first term (β 0 = 0) corresponds to the uniform distribution P (r, t|d 0 , d, x r ) → 1/πd 2 , as t → ∞. Practically, it takes much shorter time to reach the uniform distribution.

Considering P (r, t|d 0 , d, x r ) in [START_REF] Bicen | System-theoretic analysis and leastsquares design of microfluidic channels for flow-induced molecular communication[END_REF], while n = 0 corresponds to the uniform distribution, n = 1 is the dominating term that makes P (r, t|d 0 , d, x r ) non-uniform. Therefore, if this term is arbitrarily small that can be assumed to be effectively zero, we can find the required time that radial distribution becomes uniform as

ke -β 2 1 Dt d 2 < , ( 13
)
where we define the bound

J 0 β 1 d 0 d J 2 0 (β 1 ) J 0 β 1 d r
≤ k and /πd 2 as the allowed deviation in the radial distribution from the uniformity. Here, we note that πd 2 is not lumped into the parameter to keep the measures unitless, since as t → ∞, P (r, t|d 0 , d, x r ) → 1/πd 2 and can be interpreted as the fractional error. Therefore, P (r, t|d 0 , d, x r ) becomes uniform if the following condition is satisfied:

t ≥ d 2 Dβ 2 1 log( k ). ( 14
)
We note that this bound is in line with the findings of [START_REF] Taylor | Dispersion of soluble matter in solvent flowing slowly through a tube[END_REF]. Taking = 10 -2 , we obtain

t ≥ d 2 Dβ 2 1 (log(100) + log(k)). ( 15
)
Since β 1 ≈ 3.83 and k is on the order of 1, a bound on time for uniform radial distribution can be approximated as

t ≥ t * ≈ d 2 3D . ( 16
)
Therefore, for t ≥ t * , we can expect radial homogeneity. Having found the probability density function P (r, t|d 0 , d, x r ), we can find the radial distribution as p (r, t|d 0 , d, x r ) = 2πr P (r, t|d 0 , d, x r ) .

(

B. Derivation of the Axial Distribution p(x, t|r, d 0 , d, x r )

The axial distribution of the molecules in the microfluidic channel has different characteristics for different p (r, t|d 0 , d, x r ). In particular, as discussed in [START_REF] Probstein | Physicochemical Hydrodynamics: An Introduction[END_REF], if p (r, t|d 0 , d, x r ) is uniformly distributed, p (x, t|r, d 0 , d, x r ) can be easily identified as

p (x, t|r, d 0 , d, x r ) = 1 √ 4π D e t exp - x -v m t 2 2 4D e t , (18)
which is equivalent to N t v m 2 , 2D e t and can be considered as a 1-D Brownian motion with effective diffusion coefficient D e = D 1 + Pe 2 48 shifted with average axial displacement for any time t.

The fact that ( 18) is valid for uniformly distributed p (r, t|d 0 , d, x r ) implies it is valid when following condition is satisfied,

Pe P c = 4x r d , ( 19
)
as explained in [START_REF] Probstein | Physicochemical Hydrodynamics: An Introduction[END_REF]. This implication makes sense since in order to obtain a uniform distribution in radial space, the channel should have either small radius (d) or high diffusion coefficient (D) so that the molecules disperse in the channel rapidly to reach the uniform state.

On the other hand, as Pe increases, it takes some time for p (r, t|d 0 , d, x r ) to become uniform. Hence, during this period, (18) cannot be used. We therefore propose two different regions, namely uniform and non-uniform radial regions, and using these regions a partial function for the derivation of p (x r , t|r, d 0 , d, x r ) is obtained: 1) p(x r , t|r, d 0 , d, x r ) for Uniform Radial Distribution Range:

As stated before, p (x, t|r, d 0 , d, x r ) can be written using [START_REF] Szegö | Inequalities for certain eigenvalues of a membrane of given area[END_REF], if radial distribution is uniform. We have already shown that for t > t * the radial distribution is uniform. Therefore for this period we can present the axial distribution using [START_REF] Szegö | Inequalities for certain eigenvalues of a membrane of given area[END_REF]. On the other hand since the average displacement is different for t > t * due to non-uniform radial distribution, we should evaluate the average displacement and plug it to (18) instead of v m t 2 . For this aim, the average axial displacement of a molecule x ex p (t) for t seconds can be evaluated using p (r, t|d 0 , d, x r ) as

x ex p (t) = t 0 d 0 p (r, τ |d 0 , d, x r ) v(r )drdτ . ( 20
)
Once x ex p (t) is obtained, p (x r , t|r, d 0 , d, x r ) is obtained using [START_REF] Szegö | Inequalities for certain eigenvalues of a membrane of given area[END_REF] as

p (x r , t|r, d 0 , d, x r ) = exp -(xr -x ex p (t )) 2 4D e t √ 4π D e t = f X t (x r ). ( 21
)
where f X t (x) ∼ N (x ex p (t), 2D e t).

2) p(x r , t|r, d 0 , d, x r ) for Non-Uniform Radial Distribution

Range: When Pe is comparable to or greater than P c , the released molecules need some time to reach a uniform radial distribution. Until that time, dispersion of molecules is limited; hence, a new axial distribution model should be defined. Let v(t) be the average velocity and v a (t) be the required average velocity to traverse a distance x r at time t, respectively. Then, one can easily define the following relation:

v a (t) = x r t . ( 22
)
Therefore, any molecule whose average velocity is higher than v a (t) until time t passes through the receiver. Considering this fact, the probability of exceeding x r distance until time t can be written as

Prob (x(t) ≥ x r ) = Prob (v(t) ≥ v a (t)) . ( 23
)
where x(t) is a random variable that defines the axial distribution of a molecule and this distribution is not known. Alternatively, one can calculate the probability of exceeding the average required velocity v a (t) for a given time t. Note that, the required average velocity involves two terms coming from the flow and diffusion as

v a (t) = v m 1 - r * (t) 2 d 2 + U t , ( 24
)
where r * (t) is the required radial position to achieve v a (t) with U , which is the axial displacement component coming from the diffusion distributed with f U (u) ∼ N (0, 2Dt). Using p (r, t|d 0 , d, x r ), we can obtain Prob (v(t) ≥ v a (t)) considering the velocity terms coming from the flow and diffusion using p (r, t|d 0 , d, x r ) as

Prob(v(t) ≥ v a (t)) = ∞ -∞ r * (t ) 0 1 t t 0 p (r, τ |d 0 , d, x r ) dτ dr f U (u)du, ( 25
)
where the result of the first two integral in (25) gives the average velocity distribution of a molecule due to flow and the outermost integral evaluates the contribution of the diffusion.

Once Prob (x(t) ≥ x r ) is obtained, the axial distribution p (x, t|r, d 0 , d, x r ) can also be obtained using this probability as

p (x, t|r, d 0 , d, x r ) = ∂ ∂ x k (1 -Prob (x(t) ≥ x k )) x k =x , ( 26
)

TABLE I PARAMETERS

where we find the probability density function p (x, t|r, d 0 , d, x r ) from the corresponding cumulative distribution function Prob (x(t) ≤ x).

Finally, using [START_REF] Arjmandi | ISI-avoiding modulation for diffusionbased molecular communication[END_REF], n hit (t|d 0 , d, x r ) can be partially represented as

n hit (t|d 0 , d, x r ) = ⎧ ⎪ ⎨ ⎪ ⎩ ∂ ∂t Prob (v(t) ≥ v a (t)) , t ≤ t * ∂ ∂t ∞ x r f X t (x)dx t > t * (27)
Accordingly, N hit (t|d 0 , d, x r ) can be obtained as

N hit (t|d 0 , d, x r ) = Prob (v(t) ≥ v a (t)) t ≤ t * 1 -F X t (x r ) t > t * (28)
where F X t (x r ) is the cumulative normal distribution whose mean and variance are x ex p (t) and 2D e t, respectively. Even though the derivations from ( 25)-( 28) are mainly heuristic and cumbersome, N hit (t|d 0 , d, x r ) is easy to evaluate and to compare with the simulation results, hence is of significant interest for the verification of our findings.

IV. SIMULATION RESULTS

We have verified the derived cumulative distribution of molecules observed by the receiver, N hit (t|d 0 , d, x r ) and the radial distribution of molecules at the channel p (r, t|d 0 , d, x r ) using Monte Carlo simulations. For both functions, the comparisons are made for three cases:

1) Pe P c 2) P c Pe 3) Pe ∼ P c , namely, Pe is much greater than P c , P c is much greater than Pe and Pe is comparable with P c , respectively. Corresponding simulation parameters are listed in Table I. For all simulations, 10 5 molecules are released and their positions are updated every t = 10 -3 s using [START_REF] Arjmandi | Diffusionbased nanonetworking: A new modulation technique and performance analysis[END_REF].

In Fig. 2, the derived analytical expression N hit (x r , t|d 0 , d, x r ) is verified using Monte Carlo simulations. As can be seen from the figure, the derived formula fits with simulations for all three regions. Especially for Fig. 2(b) and 2(c), the proposed piecewise function, which is separated by t * , can be easily identified. For the time duration before t * , the first equation in (28) is used while the second equation is used for the time duration after t * . On the other hand for other simulations t * is either too small or too high so that only one function is used. Furthermore, as d (hence Pe) increases, the time needed for the molecules to pass through the receiver decreases and converges to x r /v m . This is expected, since as the radius of the channel d increases, the radial position of the molecules does not disperse so much in the beginning. Hence, the flow speed affecting these molecules is around v m . On the other hand, as d decreases, the molecules can reach the boundary of the channel rapidly, which reduces the velocity of the flow, and thus, it takes more time to reach the receiver. Furthermore, as indicated in Fig. 2(e) and 2(f), the derived expression is verified for different initial radial positions d 0 . Another interesting observation from these figures can be obtained by comparing them with Fig. 2(b) and 2(c), respectively. As the initial radial position moves from center to boundary, in other words as d 0 increases, the initial speed of the molecules decreases. Hence, it takes more time to reach the receiver compared to the d 0 = 0 case. Nonetheless, for all cases, all molecules will be observed by the receiver in a short duration compared to the other channel models in the literature due to flow.

In Fig. 3, the radial distributions for different channel parameters and time are presented both theoretically and numerically. As can be seen from these figures, our proposed analytical formula is verified with the simulations. Furthermore, as indicated in Figs. 3(b) and 3(d), for higher times, p (r, t|d 0 , d, x r ) turns out to be linear; hence, there is a uniform radial distribution. In particular, for the parameters in Fig. 3(a) and 3(b), at t = 200ms, the radial distribution is nonuniform while at t = 2s, this distribution becomes uniform. This is expected since for these parameters, t * = d 2 3D ≈ 0.75s, resulting a non-uniform distribution for the time below t * . The similar situation can also be observed for Fig. 3(c) and 3(d) since for this case t * ≈ 5.3s.

V. COMMUNICATION ASPECT

Let us now investigate the communication performance of the 3-D microfluidic channels. As can be deduced from Fig. 2, all molecules will be observed if the receiver waits long enough. This leads to an ISI-free communication, if the symbol duration t s is chosen large enough. Therefore, our aim is to find the required transmission slot t f that leads to minimumerror communication in the channel. In order to achieve this, we are going to utilize [START_REF] Szegö | Inequalities for certain eigenvalues of a membrane of given area[END_REF] for all three cases (Pe P c , P c Pe and Pe ∼ P c ). Although [START_REF] Szegö | Inequalities for certain eigenvalues of a membrane of given area[END_REF] is only applicable for Pe P c , we can use it for the other cases to obtain t f . As Pe increases, the expected displacement of the molecules also increases, which implies that t f obtained using ( 18) is sufficient to achieve error free communication for other cases.

In particular, our aim is to determine t f using ( 18) that almost all molecules have passed from the receiver, which implies error-free communication as

∞ x r 1 4π D e t f exp ⎛ ⎜ ⎝- x - v m t f 2 2 4D e t f ⎞ ⎟ ⎠ dx = p sig , ( 29
)
where p sig denotes the ratio of the observed molecules by the receiver until t f and, in theory, it is desired to be close to 1 in order to eliminate ISI which can be obtained as 1-p sig . (29) can be rewritten as Consecutively, (30) can be rewritten as

Q x r - v m t f 2 4D e t f = p sig . (30)
x r - v m t f 2 = 4D e t f Q -1 ( p sig ), (31)
and the solution of (31) with respect to t f can be obtained as

t f = 2 2a 2 D e + v m x r + 2a a 2 D 2 e + D e v m x r v 2 m , (32)
where a = Q -1 p sig .

Based on the findings in (32), the communication performance of the channel is examined. In order to investigate the channel performance, bit error rate (BER) of the communication channel is determined by particle based simulations. The communication scheme is achieved by binary concentration shift keying (BCSK) by releasing M molecules at the beginning of the transmission slot t s to transmit bit-1 and releasing no molecules to transmit bit-0. For each simulation, 10 6 bits are transmitted and each simulation is repeated 100 times in order to obtain accurate results. Furthermore, the optimum threshold is chosen for each simulation in order to determine the ultimate performance of the channel and a counting noise has been added to the signal as presented in [START_REF] Tepekule | ISI mitigation techniques in molecular communication[END_REF].

As can be seen from Fig. 4, as D increases, the required t f decreases and BER performance increases. Furthermore, as seen in Fig. 4(a), when t s = 0.85s, the communication performance increases incredibly, since, as shown in Fig. 4(b), t f is lower than 0.85s for D values greater than 10 -10 m/s 2 . For t s = 0.75s, there are no any error free observations for any D, since t s = 0.75s is lower than all t f values and ISI leads to an adverse effect to the channel. Furthermore, as can be seen in Fig. 4(c), BER decreases as Pe/P c increases. This is expected, since for high Pe/P c ratio, the movement of the molecules is dominated by the flow instead of diffusion, which reduces randomness.

VI. CONCLUSION

Microfluidic channels with flow is a very good candidate for molecular communication since the boundaries and flow reduce the ISI and increase the data rate, which are still considered as open problems for many channel types. Therefore, the distribution of molecules in the channel and the impulse response are necessary to determine the characteristics of the channel. Although there are many different attempts for the derivation of the impulse response of microfluidic channels, all consider some specific cases by reducing the advectiondiffusion equation to the case for a 1-D channel, which is valid for uniform radial distribution. On the other hand, it may take some time to have uniform radial distribution, for smaller diffusion coefficient or higher channel radius. Until that time, the reduced 1-D model to solve axial distribution, hence that channel impulse response, is not applicable for these cases. We, therefore, firstly derive the radial distribution of the molecules inside the channel with respect to time, and accordingly determine the required time that the molecules reach the uniform radial state. Then, using the radial distribution and the time that molecules reach the uniform state, we derive the piecewise channel impulse response. Thereafter, we verify the derived formulas of the channel impulse response and radial distribution of molecules with Monte Carlo simulations and finally communication performance of the channel is examined.

Fig. 1 .

 1 Fig. 1. General channel model for a point transmitter and a planar (circular) observing receiver.

Fig. 2 .

 2 Fig. 2. Comparison of simulation and analytical results for cumulative axial molecule distribution, N hit (t|d 0 , d, x r ) for D = 100 × 10 -12 µm/s 2 , x r = 5mm, v m = 10mm/s and different d values which indicate different Pe regions. (a) d = 5µm, d 0 = 0, Pe = 250, P c = 4000. (b) d = 15µm, d 0 = 0, Pe = 750, P c = 1333, t * = 3/4s. (c) d = 20µm, d 0 = 0, Pe = 1000, P c = 1000, t * = 4/3s. (d) d = 100µm, d 0 = 0, Pe = 5000, P c = 200. (e) d = 15µm, d 0 = 0.5d Pe = 750, P c = 1333. (f) d = 40µm, d 0 = 0.25d Pe = 2000, P c = 500.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison of simulation and analytical results for radial molecule distribution, p(r, t|d 0 , d, x r ) for D = 100 × 10 -12 m/s 2 , x r = 5mm, v m = 10mm/s and different d values and time t. Note that linearity on p(r, t|d 0 , d, x r ) implies uniformity on P(r, t|d 0 , d, x r ). (a) d = 15 × 10 -6 m, t = 200ms. (b) d = 15 × 10 -6 m, t = 2s. (e) d = 40 × 10 -6 m, t = 2s. (f) d = 40 × 10 -6 m, t = 8s.

ACKNOWLEDGEMENTS

Fatih Dinç would like to thank Professor Achim Kempf for insightful discussions on the spectral theory of the Laplacian operator.

This work was supported in part by the Scientific and Technical Research Council of Turkey (TUBITAK) under Grant 116E916 and in part by the TETAM: Telecommunications and Informatics Technologies Research Center, Bo gaziçi University, under Grant DPT-2007K120610. Research at the Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada, and by the Province of Ontario through the Ministry of Research and Innovation. (