Diffusive Molecular Communications with Reactive Molecules: Channel Modeling and Signal Design

TL;DR

This paper models and designs signal detection for molecular communication systems using reactive molecules, introducing new modulation schemes, detectors, and an algorithm for channel response computation to improve ISI mitigation.

Contribution

It presents a comprehensive framework for reactive molecular communication, including novel modulation schemes, detectors, and an efficient algorithm for channel response calculation.

Findings

Proposed multiple modulation schemes for reactive MC systems.

Derived ML detectors for different receiver types.

Developed an algorithm for channel response computation validated by simulations.

Abstract

This paper focuses on molecular communication (MC) systems using two types of signaling molecules which may participate in a reversible bimolecular reaction in the channel. The motivation for studying these MC systems is that they can realize the concept of constructive and destructive signal superposition, which leads to favorable properties such as inter-symbol interference (ISI) reduction and avoiding environmental contamination due to continuous release of signaling molecules into the channel. This work first presents a general formulation for binary modulation schemes that employ two types of signaling molecules and proposes several modulation schemes as special cases. Moreover, two types of receivers are considered: a receiver that is able to observe both types of signaling molecules (2TM), and a simpler receiver that can observe only one type of signaling molecules (1TM). For both of these receivers, the maximum likelihood (ML) detector for general binary modulation is derived under the assumption that the detector has perfect knowledge of the ISI-causing sequence. In addition, two suboptimal detectors of different complexity are proposed, namely an ML-based detector that employs an estimate of the ISI-causing sequence and a detector that neglects the effect of ISI. The proposed detectors, except for the detector that neglects ISI for the 2TM receiver, require knowledge of the channel response (CR) of the considered MC system. Moreover, the CR is needed for performance evaluation of all proposed detectors. However, deriving the CR of MC systems with reactive signaling molecules is challenging since the underlying partial differential equations that describe the reaction-diffusion mechanism are coupled and non-linear. Therefore, we develop an algorithm for efficient computation of the CR and validate its accuracy via particle-based simulation.

Figures

Tables

Table 1. TABLE I: Default Values of the System Parameters [ 11 , 12 ] .

Variable	Definition	Value
$N_A^{\mathrm{tx}},N_B^{\mathrm{tx}}$	Number of released molecules	$5\times {10}^3$ molecules
$d$	Distance between transmitter and receiver	$250$ nm
$D_A,D_B$	Diffusion coefficient	${10}^{-10}$ ${\text{m}}^2\cdot {\text{s}}^{-1}$
${\kappa }_f$	Forward reaction rate	${10}^{-17}$ molecule${}_{}{}^{-1}\cdot$m${}_{}{}^3\cdot$s-1
${\kappa }_b$	Backward reaction rate	${10}^{25}$ molecule$\cdot$m${}_{}{}^{-3}\cdot$s-1
$r$	Receiver radius	$50$ nm
$T^{\mathrm{symb}}$	Symbol duration	$200$ $\mu$s

Equations

$${\mathcal{T}}_A=\left\{t|t=(k-1)T^{\mathrm{symb}},\forall k\mathrm{and}s[k]=1\right\}$$

$${\mathcal{T}}_B=\left\{t|t=(k-1)T^{\mathrm{symb}},\forall k\mathrm{and}s[k]=0\right\},$$

$${\mathcal{T}}_A=\left\{t|t=(k-1)T^{\mathrm{symb}},\forall k\mathrm{and}s[k]=1\right\}$$

$${\mathcal{T}}_B=\left\{t|t=(k-1)T^{\mathrm{symb}}+{\tau }_1,\forall k\mathrm{and}s[k]=1\right\}.$$

$${\mathcal{T}}_A=\left\{t|t=(k-1)T^{\mathrm{symb}}+(1-s[k]){\tau }_0,\forall k\right\}$$

$${\mathcal{T}}_B=\left\{t|t=(k-1)T^{\mathrm{symb}}+s[k]{\tau }_1,\forall k\right\}.$$

$$A+B\underset{{\kappa }_b}{\stackrel{{\kappa }_f}{\rightleftharpoons }}\varnothing ,$$

$$y_i\sim \mathcal{P}\left(s{\bar{y}}_i^{(1)}(\mathbf{s})+(1-s){\bar{y}}_i^{(0)}(\mathbf{s})\right),i\in \{A,B\},$$

$${\hat{s}}^{\mathrm{ml}}$$

$${\hat{s}}^{\mathrm{ml}}=\{\begin{array}{ll}1,& \mathrm{if}{y}_A\ge \alpha (\mathbf{s})y_B+\beta (\mathbf{s})\\ 0,& \mathrm{otherwise},\end{array}$$

$${\hat{s}}^{\mathrm{ml}}$$

$${\hat{s}}^{\mathrm{ml}}=\{\begin{array}{ll}1,& \mathrm{if}{y}_A\ge \gamma (\mathbf{s})\\ 0,& \mathrm{otherwise},\end{array}$$

$$\hat{s}=\{\begin{array}{ll}1,& \mathrm{if}{y}_A\ge y_B\\ 0,& \mathrm{otherwise},\end{array}$$

$$\hat{s}=\{\begin{array}{ll}1,& \mathrm{if}{y}_A\ge \gamma \\ 0,& \mathrm{otherwise},\end{array}$$

$${\bar{y}}_i\left(t\right)=\underset{\mathbf{r}\in {\mathcal{V}}^{\mathrm{rx}}}{\int }{C}_i(\mathbf{r},t)\mathrm{d}\mathbf{r},i\in \{A,B\},$$

$$\frac{\partial C_A(\mathbf{r},t)}{\partial t}=$$

$$\frac{\partial C_B(\mathbf{r},t)}{\partial t}=$$

$$C_A^{\mathrm{eq}}{C}_B^{\mathrm{eq}}=\frac{{\kappa }_b}{{\kappa }_f}\triangleq {\kappa }_w.$$

$$C_i(\mathbf{r},t+\Delta t)=$$

$$\underset{\Delta t\to 0}{\lim }\Delta C_i^{\mathrm{x}}(\mathbf{r},t+\Delta t)=C_i^{\mathrm{x}}(\mathbf{r},t+\Delta t)-C_i(\mathbf{r},t),\mathrm{x}\in \{\mathrm{df},\mathrm{rc}\},$$

$$C_i(\mathbf{r},t+\Delta t)$$

$$={\bar{G}}_i(\mathbf{r})+C_i^{\mathrm{df}}(\mathbf{r},t+\Delta t)+C_i^{\mathrm{rc}}(\mathbf{r},t+\Delta t)-C_i(\mathbf{r},t),$$

$$C_i^{\mathrm{df}}(\mathbf{r},t+\Delta t)$$

$$C_A^{\mathrm{rc}}(\mathbf{r},t+\Delta t)$$

$$C_B^{\mathrm{rc}}(\mathbf{r},t+\Delta t)$$

$$C_i^{\mathrm{df}}(\mathbf{r},t+\Delta t)$$

$$C_i^{\mathrm{rc}}(\mathbf{r},t+\Delta t)$$

$$C_i^{\mathrm{df}}(\mathbf{r},t+\Delta t)$$

$$W_i(r,\tilde{r})=$$

$$C_A^{\mathrm{fr}}(\mathbf{r},t+\Delta t)$$

$$C_B^{\mathrm{fr}}(\mathbf{r},t+\Delta t)$$

$$C_i^{\mathrm{br}}(\mathbf{r},t+\Delta t)$$

$$\text{LLR}=$$

$$=$$

$$=\log \left(\frac{f_{\mathcal{P}}(y_A|s=1,\mathbf{s})}{f_{\mathcal{P}}(y_A|s=0,\mathbf{s})}\right)=\log \left(\frac{({\bar{y}}_A^{(1)}(\mathbf{s}){)}^{y_A}{e}^{-{\bar{y}}_A^{(1)}(\mathbf{s})}}{({\bar{y}}_A^{(0)}(\mathbf{s}){)}^{y_A}{e}^{-{\bar{y}}_A^{(0)}(\mathbf{s})}}\right)$$

Full text

MC molecular communication ISI inter-symbol interference ML maximum likelihood CR channel response PDE partial differential equation FDM finite difference method BER bit error rate RV random variable PMF probability mass function LLR log likelihood ratio

I Introduction

Recent advances in biology, nanotechnology, and medicine have given rise to the need for communication in nano/micrometer scale dimensions [2, 3]. In nature, a common strategy for communication between nano/microscale entities such as bacteria, cells, and organelles (i.e., components of cells) is diffusive MC [3, 4]. In contrast to conventional wireless communication systems that encode information into electromagnetic waves, MC systems embed information in the characteristics of signaling molecules such as their type and concentration. Therefore, diffusive MC has been considered as a bio-inspired approach for communication between small-scale nodes for applications where conventional wireless communication may be inefficient or infeasible. In fact, over the last few years, several testbeds have been developed as proof-of-concept for biological [5, 6, 7] and non-biological [8, 9, 10] diffusive MC systems.

One characteristic of MC is that the receiver always observes a constructive superposition of the molecules released in previous symbol intervals or by different transmitters. This feature leads to several undesirable effects. First, many concepts in conventional communications that rely on both constructive and destructive superposition of signals, such as precoding, beamforming, and orthogonal sequences, are not applicable in MC. Second, the release of signaling molecules in consecutive symbol intervals introduces significant ISI as the channel impulse response of MC channels is heavy-tailed. Third, if molecules are continuously released, particularly into a bounded environment, the concentration of the signaling molecules increases over time and contaminates the environment.

One solution to cope with these challenges is to use enzymes to degrade the signaling molecules in the environment [11]. It has been shown in [11] that ISI is significantly reduced if enzymes are uniformly present in the channel. However, having uniformly distributed enzymes in the environment has two main drawbacks. First, degradation of the signaling molecules via enzymes cannot be controlled, which may hurt performance. Second, the ISI reduction comes at the expense of reducing the peak concentration of the signaling molecules observed at the receiver. In [12], the authors proposed to employ acids and bases for signaling. This MC system has the advantage that the release of the molecules can be controlled by the transmitter and acids and bases can react to cancel each other out. Hence, the contamination of the environment by signaling molecules is avoided. Moreover, the use of acids and bases implies the destructive and constructive superposition of signaling molecules in the channel (not at the receiver) which can be exploited to reduce ISI. In fact, the effectiveness of this reactive signaling for ISI reduction was experimentally verified in [8]. These advantages of the MC system in [12, 8] motivate us to consider MC systems with reactive signaling molecules in this paper.

Knowledge of the CR is typically needed for receiver design and performance evaluation. However, for MC with reactive signaling molecules, deriving the CR is challenging since the underlying partial differential equations that describe the reaction-diffusion mechanism are coupled and non-linear. A closed-form solution to these equations has not yet been found, which has led to various approximations [11, 13, 14, 15, 16], see [17] for a recent overview of various channel models that have been developed so far for MC systems. For instance, in [11], the spatial and temporal distribution of the enzyme concentration was assumed to be constant to obtain an approximate solution. However, for MC systems in which the transmitter releases reactive signaling molecules into the channel, the concentrations of the molecules are temporally and spatially non-uniform and hence the constant distribution assumption does not hold. In the absence of closed-form solutions, numerical methods are commonly used to solve reaction-diffusion equations in the chemistry and physics literature [14]. This approach was applied to MCs in [12] where the authors employed a finite difference method (FDM) to solve the reaction-diffusion equation for a one-dimensional environment. Another approach to compute the expected concentrations of molecules is to average many realizations of concentrations obtained via a stochastic simulation [15, 18, 19]. However, the computational complexity of these numerical and simulation methods is very high. In [16], data is encoded in the concentration difference of two types of molecules and it is shown that if identical diffusion coefficients for both types of signaling molecules are assumed then the resulting PDE for the concentration difference is linear. However, the statistical model for the difference of the observed molecules is still a function of the concentrations of both types of molecules.

In this paper, we consider an MC system that employs two types of molecules for signaling where the signaling molecules may participate in a reversible bimolecular reaction, such as the acid and base reaction in [12]. The considered reversible bimolecular reaction involves two reactions with different rates, namely the reaction of two reactant molecules that yields a product molecule and the decomposition of this product molecule into the two reactant molecules. Moreover, we assume an unbounded (one-, two-, or three-dimensional111One- and two-dimensional environments are used as first-order approximations of three-dimensional environments in special cases to facilitate the analysis. For instance, propagation along long tubes can be approximated as one-dimensional [20] and propagation in a low-depth petri dish or thin tissues can be approximated as two-dimensional [21].) environment and a passive receiver for simplicity. We first present a general formulation for binary modulation schemes that employ two types of signaling molecules and also propose several modulation schemes as special cases. Moreover, we consider two types of receivers, namely a receiver that is able to observe both types of signaling molecules (2TM) and a simpler receiver that can observe only one type of signaling molecules (1TM). For both types of receivers, we derive a genie-aided ML detector that assumes perfect knowledge of the previously transmitted symbols, i.e., the ISI-causing symbols. The performance of this genie-aided ML detector yields an upper bound on the performance of any practical detector. In addition, we propose two suboptimal detectors of different complexity, namely an ML-based detector that employs an estimate of ISI-causing symbols and an ISI-neglecting detector that simply ignores the effect of ISI. All of the proposed detectors, except the ISI-neglecting detector for the 2TM receiver, require knowledge of the CR of the considered MC system. Moreover, the CR is needed for performance evaluation of the proposed detectors. As discussed earlier, the computation of the CR is complicated by the non-linearity that arises due to the bimolecular reaction and must be characterized for all possible sequences of numbers of molecules released by the transmitter into the channel. To address this issue, we develop an algorithm for efficient computation of the CR of the considered MC system for any arbitrary sequence of released numbers of molecules. This algorithm discretizes only the time variable and solves the resulting PDEs analytically in the space variables by decoupling the underlying reaction-diffusion equations and efficiently exploiting the simplifying characteristics inherent to the considered unbounded environment and passive receiver. This numerical method is much faster than algorithms that discretize both space and time, e.g., FDM used in [12], and avoids their inherent numerical stability issues [14]. The accuracy of the proposed algorithm for CR computation is validated using particle-based simulation. Simulation results obtained using this CR computation algorithm show that MC systems with reactive signaling have superior performance relative to those with non-reactive signaling due to the reduction of ISI enabled by the chemical reactions. Moreover, we show that unlike the MC system in [11], ISI is reduced in the considered MC system without reducing the peak of the CR. Furthermore, due to the considerable ISI reduction by chemical reactions, the proposed suboptimal ML-based detector performs very close to the genie-aided bound, and the proposed ISI-neglecting detector causes only a small performance degradation with respect to the genie-aided bound.

We note that this paper expands its conference version [1] in several directions. First, only 2TM receivers were considered in [1], whereas in this paper, we also consider 1TM receivers which are less complex. Moreover, in this paper, we derive a genie-aided detector and propose suboptimal detectors not only for 2TM receivers but also for 1TM receivers. Furthermore, we derive CR computation algorithms for one-, two-, and three-dimensional environments for completeness, whereas only a three-dimensional environment was considered in [1]. Finally, unlike [1], in this paper, the concentration update rules of the proposed CR computation algorithm are further simplified for several special cases.

The remainder of this paper is organized as follows. In Section II, the considered system is presented. The proposed optimal and suboptimal detectors are derived in Section III and the proposed algorithm for computing the CR is introduced in Section IV. Simulation results are reported in Section V, and conclusions are drawn in Section VI.

II System Model

The considered MC system consists of a transmitter, a receiver, and a channel which are introduced in detail in the following. Fig. 1 illustrates the considered system model in a three-dimensional environment.

II-A Transmitter

We assume a point-source transmitter located at the origin of the $n$-dimensional Cartesian coordinate system, i.e., $\mathbf{r}=(r_{1},\dots,r_{n})=(0,\dots,0),\,\,n\in\{1,2,3\}$. The transmitter employs two types of molecules for signaling, namely type-$A$ and type-$B$ molecules. In particular, the transmitter releases $N_{i}^{\mathrm{tx}}$ type-$i$ molecules into the channel at time instances $t\in\mathcal{T}_{i},\,i\in\{A,B\}$. By properly defining $\mathcal{T}_{i}$, different modulation schemes can be realized. Let $s[k]\in\{0,1\}$ denote the binary symbol at the $k$-th symbol interval. In this paper, we focus on the following three binary modulation schemes: conventional molecule shift keying (MoSK) [22, 23, 24], a proposed variant of on-off keying (OOK), and a new modulation scheme named order shift keying (OSK), see Fig. 2.

Conventional MoSK: For binary one, $s[k]=1$, the transmitter releases $N_{A}^{\mathrm{tx}}$ type-$A$ molecules at the beginning of the symbol interval and no type-$B$ molecules whereas for binary zero, $s[k]=0$, the transmitter releases $N_{B}^{\mathrm{tx}}$ type-$B$ molecules at the beginning of the symbol interval and no type-$A$ molecules. For this modulation scheme, we obtain

[TABLE]

where $T^{\mathrm{symb}}$ denotes the length of a symbol interval.

Proposed OOK: For binary one, $s[k]=1$, the transmitter releases $N_{A}^{\mathrm{tx}}$ type-$A$ molecules at the beginning of the symbol interval and $N_{B}^{\mathrm{tx}}$ type-$B$ molecules $\tau_{1}$ seconds after the start of the symbol interval. For binary zero, $s[k]=0$, the transmitter releases no molecules. Here, we choose $\tau_{1}$ as the peak of the CR for the case where only $N_{A}^{\mathrm{tx}}$ type-$A$ molecules are instantaneously released at $t=0$. For this modulation scheme, we obtain

[TABLE]

Note that in contrast to conventional OOK modulation which employs only one type of molecule, here two types of molecules are used for $s[k]=1$. The proposed modulation scheme effectively removes type-$A$ molecules that remain in the environment. As will be shown in Section V, this modulation scheme considerably reduces ISI compared to conventional OOK modulation. A similar idea was introduced in [25] where the release patterns of type-$A$ and type-$B$ molecules were studied for efficient pulse shaping. However, the analysis in [25] relies on the linearity of the system and is not applicable to the non-linear system considered in this paper.

Proposed OSK: For binary one, $s[k]=1$, the transmitter releases $N_{A}^{\mathrm{tx}}$ type-$A$ molecules at the beginning of the symbol interval and $N_{B}^{\mathrm{tx}}$ type-$B$ molecules $\tau_{1}$ seconds after the start of the symbol interval. In a similar manner, for binary zero, $s[k]=0$, the transmitter releases $N_{B}^{\mathrm{tx}}$ type-$B$ molecules at the beginning of the symbol interval and $N_{A}^{\mathrm{tx}}$ type-$A$ molecules $\tau_{0}$ seconds after the start of the symbol interval. In particular, we choose $\tau_{1}$ ($\tau_{0}$) as the peak of the CR assuming instantaneous release of only $N_{A}^{\mathrm{tx}}$ ($N_{B}^{\mathrm{tx}}$) type-$A$ (type-$B$) molecules at $t=0$. For this modulation scheme, we obtain

[TABLE]

In other words, for OSK modulation, the information bit is encoded in the order in which the type-$A$ and type-$B$ molecules are released in each symbol interval.

To summarize, the basic idea behind the proposed OOK and OSK modulation schemes is that after the release of a given type of signaling molecules (referred to as primary molecules), their reactive counterpart molecules (referred to as secondary/cleaning molecules) are released such that both types of molecules cancel each other out. This effectively shortens the CR and reduces ISI. The further the release times of both types of molecules are apart, the less efficient the cancellation and, as a result, the ISI reduction become. However, if the two release times are too close together, the peak concentration of the primary molecules at the receiver will be severely reduced. To account for this trade-off while maintaining a simple design, we propose to release the secondary/cleaning molecules at the transmitter when the peak concentration of the primary molecules has already been observed at the receiver. To formalize this, let $\bar{c}_{i}(t),\,\,i\in\{A,B\}$, denote the expected concentration of type-$i$ molecules at the receiver when only $N_{i}^{\mathrm{tx}}$ type-$i$ molecules are released at time $t=0$. Then, we choose $\tau_{1}=\mathrm{argmax}_{t}\,\bar{c}_{A}(t)$ and $\tau_{0}=\mathrm{argmax}_{t}\,\bar{c}_{B}(t)$.

II-B Channel

We assume an unbounded $n$-dimensional environment for $n\in\{1,2,3\}$. The type-$A$ and type-$B$ molecules released by the transmitter diffuse in the environment with diffusion coefficients $D_{A}$ and $D_{B}$ in m${}^{2}\cdot$s*-1*, respectively, and may participate in the following bimolecular reaction

[TABLE]

where $\kappa_{f}$ and $\kappa_{b}$ denote the forward reaction rate constant in molecule${}^{-1}\cdot$m${}^{n}\cdot$s${}^{-1},\,\,n\in\{1,2,3\}$, and the backward reaction rate constant in molecule$\cdot$m${}^{-n}\cdot$s${}^{-1},\,\,n\in\{1,2,3\}$, respectively. Moreover, symbol $\varnothing$ denotes chemical species which are of no interest for communication. Note that (7) includes the reactions considered in [12, 16]. Moreover, if type-$B$ molecules represent enzymes and only type-$A$ molecules are used for signaling, (7) includes the degradation reaction in [11] when the enzyme concentration is constant everywhere and is much larger than the concentration of the type-$A$ molecules such that the reaction in (7) does not change the enzyme concentration.

II-C Receiver

Let us assume that the receiver is located at distance $d$ from the transmitter. We consider two types of passive receivers: i) a receiver that is able to count the numbers of both type-$A$ and type-$B$ molecules within its volume, denoted by 2TM receiver; and ii) a receiver that is able to count only type-$A$ molecules (1TM) within its volume, denoted by 1TM receiver. 1TM receivers are simpler to implement whereas 2TM receivers offer better performance since they exploit the diversity gain that observing two types of molecules provides. However, as is shown in Section V, regardless of whether a 1TM or a 2TM receiver is employed, reactive MC systems benefit from the ISI reduction enabled by employing two types of reactive signaling molecules.

Let $\bar{y}_{A}\left(t\right)$ and $\bar{y}_{B}\left(t\right)$ denote the expected numbers of type-$A$ and type-$B$ molecules observed at the receiver at time $t$, respectively, due to the release of a known sequence of numbers of molecules by the transmitter. We refer to $\bar{y}_{i}\left(t\right),\,\,i\in\{A,B\}$, as the CR of the considered MC system. Note that depending on the length of the symbol interval, at a given time, the receiver may observe molecules released by the transmitter in multiple previous symbol intervals, i.e., ISI may exist. Let $L$ be the length of the channel memory222Theoretically, the memory length of the considered MC channel is infinite; however, from a practical point-of-view, the effect of the previous symbols becomes negligible after several symbol intervals. and $\mathbf{s}\in\{0,1\}^{L-1}$ denote the vector of the $L-1$ previously transmitted symbols. We also refer to $\mathbf{s}$ as the ISI-causing sequence333For notational simplicity, we drop the symbol index $k$ in the remainder of the paper.. Therefore, given $\mathbf{s}$ and $s$, the number of type-$i$ molecules counted at the receiver at sample time $t_{s}$ is modelled as [26, 27]

[TABLE]

where $\mathcal{P}(\lambda)$ denotes a Poisson random variable (RV) with mean $\lambda$. Moreover, $\bar{y}_{i}^{(s)}(\mathbf{s})$ is $\bar{y}_{i}\left(t_{s}\right)$ under the condition that the symbol in the current symbol interval is $s\in\{0,1\}$ and the ISI-causing sequence is $\mathbf{s}\in\{0,1\}^{L-1}$. We note that the Poisson distribution assumed for $y_{i}$ in (8) is an approximation which has been shown to be accurate for reaction-diffusion processes in the chemistry and physics literature [28, 29]. In Section V, we will validate the Poisson model in (8) using the particle-based simulator developed in Appendix F.

Note that due to the reaction process, the CR of the considered MC system $\bar{y}_{i}^{(s)}(\mathbf{s})$ is a non-linear function of the transmitted data symbols $s$ and $\mathbf{s}$. Therefore, we cannot simply compute the CR for one shot transmission and use convolution to capture the effect of ISI [27]. In particular, to fully characterize the average behavior of the system, one has to compute the CR for both symbol hypotheses $s\in\{0,1\}$ and all $2^{L-1}$ possible ISI-causing sequences. In Section III, we derive optimal and suboptimal detectors assuming $\bar{y}_{i}^{(s)}(\mathbf{s}),\,\,\forall s,\mathbf{s}$, is known. Then, in Section IV, we derive an efficient numerical algorithm for computation of CR $\bar{y}_{i}^{(s)}(\mathbf{s}),\,\,\forall s,\mathbf{s}$.

III Detection Methods for Binary Modulation

In this section, we derive the genie-aided ML detector for binary modulation assuming the ISI-causing sequence is known. This provides an upper bound on performance for any practical detector. Subsequently, we propose two suboptimal practical detectors having different complexities.

III-A Optimal Genie-Aided ML Detector

In the following, we focus on symbol-by-symbol detection. We consider a genie-aided ML detector that has perfect knowledge of the ISI-causing sequence.

2TM Receivers: For 2TM receivers, the genie-aided ML detection problem for the considered binary modulation schemes is given by

[TABLE]

where $\Pr(\cdot)$ denotes probability and $f_{\mathcal{P}}(x)=\frac{\lambda^{x}e^{-\lambda}}{x!}$ is the probability mass function (PMF) of a Poisson RV with mean $\lambda$. Equality $(a)$ follows from the fact that conditioned on $\bar{y}_{i}^{(s)}(\mathbf{s}),\,\,i\in\{A,B\}$, and $(s,\mathbf{s})$, RVs $y_{A}$ and $y_{B}$ are independent. The optimal detector for 2TM receivers is given in the following proposition.

Proposition 1

The genie-aided ML detector for 2TM receivers as a solution of (9) is given by

[TABLE]

where \alpha(\mathbf{s})=\frac{1}{\chi(\mathbf{s})}\log\Big{(}\frac{\bar{y}_{B}^{(0)}(\mathbf{s})}{\bar{y}_{B}^{(1)}(\mathbf{s})}\Big{)}, \beta(\mathbf{s})=\frac{1}{\chi(\mathbf{s})}\big{(}\bar{y}_{A}^{(1)}(\mathbf{s})+\bar{y}_{B}^{(1)}(\mathbf{s})-\bar{y}_{A}^{(0)}(\mathbf{s})-\bar{y}_{B}^{(0)}(\mathbf{s})\big{)}, and \chi(\mathbf{s})=\log\Big{(}\frac{\bar{y}_{A}^{(1)}(\mathbf{s})}{\bar{y}_{A}^{(0)}(\mathbf{s})}\Big{)}.

Proof:

The proof is given in Appendix A. ∎

Proposition 1 reveals that the optimal detection rule for 2TM receivers involves only linear processing of the received signals $y_{A}$ and $y_{B}$ and a threshold comparison. The coefficients of the corresponding linear operation, i.e., $\alpha(\mathbf{s})$ and $\beta(\mathbf{s})$, depend on the CR $\bar{y}_{i}^{(s)}(\mathbf{s})$ and the ISI-causing sequence $\mathbf{s}$.

1TM Receivers: For 1TM receivers, the ML problem is formulated as

[TABLE]

The solution to the above problem is given in the following corollary.

Corollary 1

The genie-aided ML detector as a solution of (11) is given by

[TABLE]

where \gamma(\mathbf{s})=\left(\bar{y}_{A}^{(1)}(\mathbf{s})-\bar{y}_{A}^{(0)}(\mathbf{s})\right)\big{/}\log\Big{(}\frac{\bar{y}_{A}^{(1)}(\mathbf{s})}{\bar{y}_{A}^{(0)}(\mathbf{s})}\Big{)}.

Proof:

The proof is given in Appendix B. ∎

Corollary 1 shows that the optimal detection rule for 1TM receivers is a simple threshold detector. The detection threshold, $\gamma(\mathbf{s})$, depends on the CR $\bar{y}_{A}^{(s)}(\mathbf{s})$ and the ISI-causing sequence $\mathbf{s}$.

Remark 1

We note that the detectors in Proposition 1 and Corollary 1 are not limited to the specific modulation schemes introduced in Section II-A and can be applied in any (reactive or non-reactive) MC system that employs binary modulation based on two types and one type of signaling molecules, respectively. The adopted modulation scheme and whether or not the signaling molecules are reactive are reflected in the value of the CR $\bar{y}_{i}^{(s)}(\mathbf{s})$, of course.

The optimal detectors in Proposition 1 and Corollary 1 provide performance upper bounds since they require knowledge of the ISI-causing sequence $\mathbf{s}$ which is not available in practice. Therefore, in the following, we propose some practical suboptimal detectors.

III-B Suboptimal Detectors

In the following, we propose two categories of suboptimal detectors.

III-B1 ML-Based Detectors

Here, we assume that the receiver employs the detectors in Proposition 1 or Corollary 1 but uses its own estimates of the previous symbols as the ISI-causing sequence. This leads to suboptimal detectors which we refer to as “ML-based detectors with estimated ISI”. We show in Section V that the performance of these ML-based detectors with estimated ISI is very close to the performance upper bound provided by the corresponding genie-aided ML detector.

III-B2 ISI-Neglecting Detectors

Recall that the main motivation for considering MC with reactive signaling as well as the proposed OOK and OSK modulation schemes was to reduce ISI. In the following, we simplify the detectors given in Proposition 1 and Corollary 1 assuming an ISI-free channel and propose to use them when the ISI is sufficiently small. In particular, for an ISI-free channel, parameters $\alpha(\mathbf{s})$, $\beta(\mathbf{s})$, and $\gamma(\mathbf{s})$ do not depend on previous symbols and are denoted by $\alpha$, $\beta$, and $\gamma$, respectively. This simplifies the detector in Proposition 1 for 2TM receivers to

[TABLE]

and the detector in Corollary 1 for 1TM receivers to

[TABLE]

where in (13), we further assume symmetry with respect to the two molecule types, i.e., $D_{A}=D_{B}$ and $N_{A}=N_{B}$ hold444Note that if $D_{A}=D_{B}$ and $N_{A}^{\mathrm{tx}}=N_{B}^{\mathrm{tx}}$ hold, we obtain $\bar{y}_{A}^{(1)}(\mathbf{s})=\bar{y}_{B}^{(0)}(\mathbf{s})\triangleq a$ and $\bar{y}_{A}^{(0)}(\mathbf{s})=\bar{y}_{B}^{(1)}(\mathbf{s})\triangleq b$. This leads to $\alpha=\frac{1}{\log(a/b)}\log(a/b)=1$ and $\beta=\frac{1}{\log(a/b)}(a+b-b-a)=0$ in Proposition 1..

Remark 2

Note that the detector in (13) does not need knowledge of the CR which makes it suitable for MC systems with limited computational capabilities. However, for the detector in (14), the computation of $\gamma$ requires knowledge of the CR, see e.g. [30, Eq. (5)] for the optimal value of $\gamma$ for an ISI-free MC channel. On the other hand, the advantage of the detector in (14) over that in (13) is that the corresponding receiver is simpler because it needs to be able to count only one type of molecule.

IV Computation of Channel Response

for MC Systems with Reactive Signaling Molecules

In this section, we first formally present the problem statement for CR computation. Next, we derive a numerical algorithm for computing the CR and discuss its complexity with respect to the available numerical methods for CR computation.

IV-A Problem Statement

Let $C_{A}(\mathbf{r},t)$ and $C_{B}(\mathbf{r},t)$ denote the concentrations of type-$A$ and type-$B$ molecules in point $\mathbf{r}$ and at time $t$. Considering a passive receiver, $\bar{y}_{i}\left(t\right)$ is obtained as

[TABLE]

where $\mathcal{V}^{\mathrm{rx}}$ is the set of points within the receiver space. Concentrations $C_{A}(\mathbf{r},t)$ and $C_{B}(\mathbf{r},t)$ can be found using the following reaction-diffusion equations [12, 15]

[TABLE]

where $\nabla^{2}=\frac{\partial^{2}}{\partial r_{1}^{2}}+\cdots+\frac{\partial^{2}}{\partial r_{n}^{2}},\,\,n\in\{1,2,3\}$, is the Laplace operator. Moreover, $G_{i}(\mathbf{r},t)=\sum_{t_{i}\in\mathcal{T}_{i}}N_{i}^{\mathrm{tx}}\delta_{\mathrm{d}}(\mathbf{r},t-t_{i}),\,\,i\in\{A,B\}$, represents the concentration of type-$i$ molecules that are released by the transmitter into the channel where $\delta_{\mathrm{d}}(\mathbf{r},t)=\delta_{\mathrm{d}}(r_{1})\cdots\delta_{\mathrm{d}}(r_{n})\delta_{\mathrm{d}}(t),\,\,n\in\{1,2,3\}$, and $\delta_{\mathrm{d}}(\cdot)$ is the Dirac delta function. In addition, without loss of generality, we assume initial condition $\underset{\epsilon\to 0}{\lim}\,\,C_{i}(\mathbf{r},-\epsilon)=C_{i}^{\mathrm{init}}(\mathbf{r})$ where $C_{i}^{\mathrm{init}}(\mathbf{r})$ is the spatial concentration of type-$i$ molecules before transmission starts at $t=0$.

Remark 3

Suppose that the MC environment is in steady state equilibrium before transmission starts at $t=0$. In this case, the concentrations of type-$A$ and type-$B$ molecules are temporally and spatially uniform, i.e., $C_{i}^{\mathrm{init}}(\mathbf{r})\triangleq C_{i}^{\mathrm{eq}}$. Therefore, substituting $\frac{\partial C_{i}(\mathbf{r},t)}{\partial t}=D_{i}\nabla^{2}C_{i}(\mathbf{r},t)=0$ into (16a), we obtain that the following relation holds for the equilibrium concentration [31, Chapter 14]

[TABLE]

For instance, assuming that type-$A$ and type-$B$ molecules represent acid and base and the fluid is water, we have $\kappa_{f}=1.4\times 10^{11}\,\,1/(\text{Ms})$ and $\kappa_{b}=1.4\times 10^{-3}\,\,\text{M}/\text{s}$ at $25^{\circ}$C which leads to $\kappa_{w}=10^{-14}$ M2 where M, called molar, is a unit of concentration and represents the number of moles per every liter of solution in a three-dimensional environment555The concentration $C$ in molecule$\cdot$m-3 is related to the molar concentration $C^{\mathrm{mol}}$ in mol$\cdot$L-1 via $C=10^{3}n_{\mathrm{av}}C^{\mathrm{mol}}$ where $n_{\mathrm{av}}=6.02\times 10^{23}$ is the Avogadro constant [31]. [31, Chapter 14]. Therefore, when the concentrations of acid and base are identical, i.e., $C_{A}^{\mathrm{eq}}=C_{B}^{\mathrm{eq}}=10^{-7}$ M, the pH, defined by $\text{pH}=-\log_{10}(C_{A}^{\mathrm{eq}})$, is $7$. As can be seen from (18), regardless of whether or not molecules are released by the transmitter into the channel, some concentrations of type-$A$ and type-$B$ molecules are present in the environment which can be interpreted as environmental noise molecules. However, unlike non-reactive MC systems where the concentrations of signaling and noise molecules in the environment are typically assumed to be independent, for reactive MC systems, the concentrations of noise and signaling molecules are not independent and are jointly described by (16a).

As stated earlier, the general reaction-diffusion equations in (16a) have not yet been solved in closed form. The main difficulty in doing so arises from the coupling of the two equations and the non-linear term $\kappa_{f}C_{A}(\mathbf{r},t)C_{B}(\mathbf{r},t)$. Note that even if we assume that one of the variables, e.g. $C_{B}(\mathbf{r},t)$, is fixed, it is still challenging to solve (16aa) in terms of $C_{A}(\mathbf{r},t)$. Therefore, in the following, we derive a numerical method for solving (16a) in a computationally efficient manner. This is achieved by fully exploiting the characteristics of the considered system model.

IV-B Derivation of the CR

Let us assume that time is divided into a series of small intervals of length $\Delta t$. The main idea behind the proposed approach for computing the CR is to find the concentrations at the end of each time interval given the concentrations at the beginning of the time interval while exploiting the condition $\Delta t\to 0$. In particular, from the reaction-diffusion equations in (16a), we have

[TABLE]

where $\Delta C_{i}^{\mathrm{df}}(\mathbf{r},t+\Delta t)$ and $\Delta C_{i}^{\mathrm{rc}}(\mathbf{r},t+\Delta t)$ are concentration changes from time $t$ to time $t+\Delta t$ due to diffusion and (forward and backward) reaction, respectively. Note that, in general, the diffusion and reaction processes are coupled which makes the derivation of a closed-form expression for $\Delta C_{i}^{\mathrm{x}}(\mathbf{r},t+\Delta t),\,\,\mathrm{x}\in\{\mathrm{df},\mathrm{rc}\},$ intractable. Nevertheless, the following lemma specifies $\Delta C_{i}^{\mathrm{x}}(\mathbf{r},t+\Delta t)$ for $\Delta t\to 0$. First, let us formally introduce the following two sets of assumptions which are needed for a rigorous statement of Lemma 1: A1) We assume $C_{i}(\mathbf{r},t)\neq 0$ and $\nabla^{2}C_{i}(\mathbf{r},t)\neq 0$ hold for $\forall\mathbf{r},t$. A2) We assume $\mathcal{T}_{A}\cap\mathcal{T}_{B}$ is empty and that $\mathcal{T}_{i}\subset\mathcal{T}$ where $\mathcal{T}=\{t|t=m\Delta t-\varepsilon,\,\,m\in\mathbb{N}\}$, $\mathbb{N}$ is the set of natural numbers, and $\varepsilon$ is an arbitrary small positive real number satisfying $\varepsilon\ll\Delta t$. Note that assumption A1 is mild as it generally holds after the first release of molecules by the transmitter for any point $\mathbf{r}$ in the MC environment. Similarly, assumption A2 is not limiting since $\Delta t$ can be chosen to be sufficiently small to accommodate any desired release time pattern of molecules with high accuracy while satisfying A2.

Lemma 1

Under assumptions A1 and A2, the following asymptotic result holds for $\Delta C_{i}^{\mathrm{x}}(\mathbf{r},t+\Delta t)$:

[TABLE]

where $C_{i}^{\mathrm{df}}(\mathbf{r},t)$ and $C_{i}^{\mathrm{rc}}(\mathbf{r},t)$ are the type-$i$ molecule concentration assuming that in interval $[t,t+\Delta t]$, only diffusion and only reactions occur, respectively, and the other phenomenon is absent.

Proof:

The proof is given in Appendix C. ∎

The main result in the above lemma is that assuming $\Delta t\to 0$, we can decouple the impact of diffusion and reaction and compute the concentration change due to each phenomenon separately. More specifically, the above lemma indicates that the effect of the coupling of the involved phenomena on the concentration change has a negligible impact compared to the effect of the individual phenomena as $\Delta t\to 0$. The decoupling of reaction and diffusion equations in (20) can be seen as a form of general splitting methods used in numerical mathematics [32], [33], [34]. Based on this result, the concentration update rule for $C_{i}(\mathbf{r},t+\Delta t)$ for the considered MC system is given in the following proposition.

Proposition 2

Under assumptions A1 and A2, $C_{i}(\mathbf{r},t+\Delta t)$ is obtained as

[TABLE]

where $\bar{G}_{i}(\mathbf{r})=\sum_{t_{i}\in\mathcal{T}_{i}}N_{i}^{\mathrm{tx}}\delta_{\mathrm{d}}(\mathbf{r})\delta_{\mathrm{k}}(t+\Delta t-\varepsilon-t_{i})$ where $\delta_{\mathrm{k}}(\cdot)$ is the Kronecker delta function. Moreover, for the MC system under consideration, $C_{i}^{\mathrm{df}}(\mathbf{r},t+\Delta t)$ and $C_{i}^{\mathrm{rc}}(\mathbf{r},t+\Delta t)$ are given by (22a) on the top of the next page where $c_{1}(\mathbf{r},t)=C_{A}(\mathbf{r},t)-C_{B}(\mathbf{r},t)$, $c_{2}(\mathbf{r},t)=\sqrt{\kappa_{f}^{2}c_{1}^{2}(\mathbf{r},t)+4\kappa_{f}\kappa_{b}}$, $c_{3}(\mathbf{r},t)=C_{A}(\mathbf{r},t)+C_{B}(\mathbf{r},t)$, and $c_{4}(\mathbf{r},t)=\frac{c_{2}(\mathbf{r},t)-\kappa_{f}c_{3}(\mathbf{r},t)}{c_{2}(\mathbf{r},t)+\kappa_{f}c_{3}(\mathbf{r},t)}$.

Proof:

The proof is given in Appendix D. ∎

Algorithm 1 summarizes the simulation steps for CR computation using Proposition 2 where $T^{\max}$ is the maximum simulation time.

Remark 4

We note that as $\Delta t\to 0$, one can also use the following first order approximations for $C_{i}^{\mathrm{df}}(\mathbf{r},t+\Delta t)$ and $C_{i}^{\mathrm{rc}}(\mathbf{r},t+\Delta t)$ respectively:

[TABLE]

However, for these approximations to be valid, $\Delta t$ has to be chosen much smaller than when the exact expressions in Proposition 2 are adopted. Moreover, small values of $\Delta t$ lead to stability issues for computation of $\nabla^{2}C_{i}(\mathbf{r},t)$, see Section IV-D for a detailed discussion.

IV-C Simplifications for Special Cases

In the following, we further simplify the diffusion term $C_{i}^{\mathrm{df}}(\mathbf{r},t+\Delta t)$ and reaction term $C_{i}^{\mathrm{rc}}(\mathbf{r},t+\Delta t)$ for the considered system model.

Diffusion: The computational complexity of Algorithm 1 is largely determined by (22aa) since for each update, an $n$-dimensional integral has to be evaluated for each point $\mathbf{r}\in\mathbb{R}^{n}$ where $\mathbb{R}$ is the set of real numbers. Nevertheless, for the commonly adopted assumption of a point-source transmitter with impulsive release [11, 12, 16], the computation of $C_{i}^{\mathrm{df}}(\mathbf{r},t)$ can be significantly simplified using the following corollary.

Corollary 2

For the MC system under consideration, the concentrations of the molecules are symmetric with respect to origin and hence are only functions of variable $r\triangleq\|\mathbf{r}\|$. In this case, $C_{i}^{\mathrm{df}}(\mathbf{r},t+\Delta t)$ can be simplified as

[TABLE]

*where $W_{i}(r,\tilde{r})$ is given by *

[TABLE]

Here, $I_{0}(\cdot)$ is the zeroth order modified Bessel function of the first kind, $\sinh(\cdot)$ is the sine hyperbolic function, and $\cosh(\cdot)$ is the cosine hyperbolic function.

Proof:

The proof is given in Appendix E. ∎

Note that the $n$-dimensional integral in (22aa) is simplified to a one-dimensional integral in (27) which has to be evaluated for a one-dimensional space variable $r\in\mathbb{R}^{+}$, where $\mathbb{R}^{+}$ denotes the set of non-negative real numbers. In addition, the term $W_{i}(r,\tilde{r})$ in the integral does not depend on the molecule concentrations. Hence, we can evaluate it offline and use it for online concentration updates. In other words, the integral in (27) can be approximated by summation and multiplication operations.

Reaction: In certain cases, only the forward reaction (degradation) or the backward reaction (production) may be present. To this end, one must evaluate the limits of (22ab) and (22ac) for either $\kappa_{f}\to 0$ or $\kappa_{b}\to 0$. In particular, (22ab) and (22ac) have indeterminate forms at $\kappa_{f}=0$. Moreover, when $C_{A}(\mathbf{r},t)=C_{B}(\mathbf{r},t)$ holds, (22ab) and (22ac) have indeterminate forms at $\kappa_{b}=0$, too. The following corollary provides the results for these cases.

Corollary 3

For the considered MC system, the molecule concentrations when only the forward reaction is present, denoted by $C^{\mathrm{fr}}_{i}(\mathbf{r},t+\Delta t)$, and when only the backward reaction is present, denoted by $C^{\mathrm{br}}_{i}(\mathbf{r},t+\Delta t)$, are given by

[TABLE]

where $c_{5}(\mathbf{r},t)=\frac{C_{B}(\mathbf{r},t)}{C_{A}(\mathbf{r},t)}$.

Proof:

Due to the indeterminate form of $C^{\mathrm{rc}}_{i}(\mathbf{r},t+\Delta t)$ at $\kappa_{b}=0$, we use L’Hopital’s rule [35] to compute $C^{\mathrm{fr}}_{i}(\mathbf{r},t+\Delta t)=\underset{\kappa_{b}\to 0}{\lim}\,\,C^{\mathrm{rc}}_{i}(\mathbf{r},t+\Delta t)$ where $C^{\mathrm{rc}}_{i}(\mathbf{r},t+\Delta t)$ is given by (22ab) and (22ac). This leads to $C^{\mathrm{fr}}_{i}(\mathbf{r},t+\Delta t)$ given in (29a) and (30a) for the case when $c_{1}(\mathbf{r},t)\neq 0$. Note that these expressions for $C^{\mathrm{fr}}_{i}(\mathbf{r},t+\Delta t)$ are indeterminate at $c_{1}(\mathbf{r},t)=0$. In a similar manner, we apply L’Hopital’s rule to $C^{\mathrm{fr}}_{i}(\mathbf{r},t+\Delta t)$ given in (29a) and (30a) for $c_{1}(\mathbf{r},t)\neq 0$ in order to find the corresponding expressions when $c_{1}(\mathbf{r},t)=0$. For the case, when $\kappa_{f}=0$, it is easier to solve the original PDEs in (16a) which simplifies to $\frac{\partial C_{i}^{\mathrm{br}}(\mathbf{r},t)}{\partial t}=\kappa_{b}$ and has the immediate solution given in (31). This completes the proof. ∎

IV-D Discussion on Complexity

In the following, we compare the computational complexity of Algorithm 1 with that of the conventional numerical methods such as FDM. Most numerical methods in the literature rely on discretization of space and time to solve the reaction-diffusion equations in (16a) [14, 12]. In particular, for FDM, time and space are discretized into small intervals to approximate the differential operators in (16a). Let $N_{t}$ and $N_{r}$ denote the number of time steps and the number of space points per dimension used for discretization of time and space, respectively. The advantage of numerical methods such as FDM is their universality as they can be applied to general PDEs. However, for the approximation of the differential operators to be accurate, the adopted step sizes should be very small, i.e., $N_{t}$ and $N_{r}$ should be very large [14]. However, increasing $N_{t}$ and $N_{r}$ too much results in numerical instability for computation of the differential operators, e.g., $\nabla^{2}C_{i}(\mathbf{r},t)$. In fact, it is widely documented in the numerical mathematics literature that the overall accuracy and stability of FDM methods heavily depend on proper selection of the step sizes $N_{t}$ and $N_{r}$ [14, 12]. Moreover, the appropriate values of the step sizes are application-specific and there is no universal guideline that can be applied to all PDEs and system parameters. On the contrary, the proposed approach in Algorithm 1 is a hybrid method where time is discretized; however, the problem is solved analytically in (22a) and (27) with respect to the space variables. Since no differential operator is used in these update rules, Algorithm 1 does not face the stability issues of FDM methods. Moreover, compared to FDM methods, larger step sizes can be adopted in the proposed algorithm for the same overall accuracy. This is due to the fact that unlike FDM methods, which approximate the differential operators in (16a), we solve them exactly for the diffusion-only and reaction-only scenarios and errors occur only due to ignoring the coupling of reaction and diffusion in (16a). This makes the proposed algorithm much faster than pure numerical methods such as FDM and avoids the stability issues that arise in the FDM methods.

Remark 5

In this paper, we do not quantify the relative complexities of FDM methods and the proposed algorithm since this comparison depends on the specific adopted implementation of the FDM method and the values of the system parameters. However, in this remark, we discuss how the complexity of these methods scale in terms of $N_{t}$ and $N_{r}$. In particular, since the molecule concentrations have to be updated in each time step and for each space point, the computational complexity of FDM methods has order $\mathcal{O}(N_{t}N_{r}^{n})$ where $\mathcal{O}(\cdot)$ is the standard big-O notation. Moreover, if (16a) is rewritten in polar and spherical coordinates when $n=2$ and $3$, respectively, the computational complexity can be reduced to $\mathcal{O}(N_{t}N_{r})$ by exploiting symmetry with respect to $\|\mathbf{r}\|$. The complexity order of Algorithm 1 is given by $\mathcal{O}(N_{t}N_{r}^{2n})$. However, by exploiting the symmetry and using the update rule for $C_{i}^{\mathrm{df}}(\mathbf{r},t)$ given in (27), the complexity of Algorithm 1 is reduced to $\mathcal{O}(N_{t}N_{r}^{2}),\,\,n=1,2,3$. Furthermore, if one employs the first order approximations in (25a) as concentration update rules in Algorithm 1, complexity is further reduced to $\mathcal{O}(N_{t}N_{r})$. However, the implementation of the concentration update rules in (25a) suffers the same stability issues as the FDM methods. In fact, the main advantage of the proposed algorithm using the updates rules in (22a) and (27) is that it does not face the stability issues of FDM methods. Moreover, the required $N_{t}$ and $N_{r}$ are much smaller than those needed for FDM methods. Although for given $N_{t}$ and $N_{r}$, the complexity order of the proposed algorithm is higher than that of FDM methods. As a remedy to this issue, one may compute the integration in (22aa) over $\tilde{\mathbf{r}}\in\mathcal{R}(\mathbf{r})$, instead of over the entire space where \mathcal{R}(\mathbf{r})=\big{\{}\tilde{\mathbf{r}}\big{|}\|\tilde{\mathbf{r}}-\mathbf{r}\|\leq R\big{\}} and $R$ is a parameter referred to as integration radius. Thereby, the complexity order of the proposed algorithm reduces to $\mathcal{O}(N_{t}N_{r}^{n})$ and $\mathcal{O}(N_{t}N_{r})$ for concentrations that are asymmetric and symmetric with respect to the origin, respectively, i.e., the same complexity order as for FDM methods.

V Simulation Results

In this section, we first verify the accuracy of the proposed CR computation algorithm and the Poisson statistics assumed in (8) using the particle-based simulator described in Appendix F. Next, the effect of some important system parameters on the CR is studied, and finally, the performance of the proposed modulation schemes and detectors is evaluated in terms of the bit error rate (BER). For clarity of presentation, we focus on a three-dimensional environment and unless stated otherwise, the default values for the system parameters given in Table I are used. The time step size is $\Delta t=1\,\mu$s for Algorithm 1 and the particle-based simulator, respectively. Moreover, the simulation results shown in Figs. 3 and 4 are obtained by running the particle-based simulator $10^{4}$ times.

V-A Verification of the Proposed CR Computation Algorithm and the Poisson Statistics

In Fig. 3, we plot the expected number of molecules observed at the receiver, i.e., CR, where $\mathcal{T}_{A}=\{200\}$ and $\mathcal{T}_{B}=\{300\}$ $\mu$s and for initial condition $C_{A}^{\mathrm{init}}(\mathbf{r})=C_{B}^{\mathrm{init}}(\mathbf{r})=0$. We show numerical results obtained with Algorithm 1 and simulation results generated with the particle-based simulator introduced in Appendix F. From Fig. 3, we observe a perfect agreement between numerical and simulation results which validates the accuracy of the proposed numerical algorithm. In Fig. 3, we also plot the equilibrium concentration which can be computed based on Remark 3 as $\bar{y}_{i}^{\mathrm{eq}}=C_{i}^{\mathrm{eq}}V^{\mathrm{rx}}=0.5236$ where $C_{A}^{\mathrm{eq}}=C_{B}^{\mathrm{eq}}=\sqrt{\frac{\kappa_{b}}{\kappa_{f}}}$ and $V^{\mathrm{rx}}$ is the receiver volume. Note that we assume no molecules are present in the environment at $t=0$ and the transmitter also releases no molecule before $t=200\mu$s. Nevertheless, $\bar{y}_{i}(t)$ increases in interval $t\in[0,200\mu\text{s}]$ due to the backward reaction (production) in (7) and finally converges to $\bar{y}_{i}(t)=\bar{y}_{i}^{\mathrm{eq}}$ since the backward reaction and the forward reaction (degradation) in (7) reach an equilibrium, cf. Remark 3.

Next, we evaluate the accuracy of the Poisson statistical model for $y_{i}(t)$ introduced in (8), where the mean is obtained using Algorithm 1, by comparing it with the histogram of $y_{i}(t)$, which is obtained via particle-based simulation, i.e., the “exact distribution”. In particular, in Fig. 4, we show the PMF of $y_{i}(t)$ for a few different time instances, namely $t=200,300$ $\mu$s for type-$A$ molecules and $t=200,400$ $\mu$s for type-$B$ molecules. The same parameters as for Fig. 3 are adopted for Fig. 4. Note that the PMFs of $y_{A}(t)$ and $y_{B}(t)$ are identical at $t=200$ $\mu$s, i.e., before the release of molecules by the transmitter, due to the identical characteristics of the type-$A$ and type-$B$ molecules. Moreover, we emphasize that for the considered time instances, the mean of $y_{i}(t)$ assumes a small value $\bar{y}_{A}(t=200\mu\text{s})=\bar{y}_{B}(t=200\mu\text{s})\approx 0.5$, a moderate value $\bar{y}_{B}(t=400\mu\text{s})\approx 4$, and a large value $\bar{y}_{A}(t=300\mu\text{s})\approx 12$, and is given by the CR plotted in Fig. 3. As can be seen from Fig. 4, the Poisson statistics with the mean obtained by Algorithm 1 can accurately model the histograms obtained from particle-based simulation for the considered small, moderate, and large mean values. This validates the Poisson model for MC systems with reactive signaling and is in line with the results reported in the chemistry and physics literature [28, 29].

V-B Effect of System Parameters on CR

In this subsection, we study the impact of the diffusion coefficient, forward reaction rate, and backward reaction rate on the CR.

In Fig. 5, we plot the CR for type-$A$ molecules, $\bar{y}_{A}(t)$, versus time for $\mathcal{T}_{A}=\{200\}$ and $\mathcal{T}_{B}=\{300\}$ $\mu$s, and initial condition $C_{A}^{\mathrm{init}}(\mathbf{r})=C_{B}^{\mathrm{init}}(\mathbf{r})=0$. We consider a fixed $D_{A}=10^{-10}$ m${}^{2}\cdot$s*-1* and different $D_{B}=[0.1,0.5,1,2,10]\times 10^{-10}$ m${}^{2}\cdot$s*-1*. As can be seen from this figure, increasing the diffusion coefficient of the type-$B$ molecules shortens the CR. This is due to the fact that the cleaning type-$B$ molecules reach and cancel out the previously released type-$A$ molecules faster when $D_{B}$ is larger, i.e., the ISI reduction is more effective. Moreover, as expected, the equilibrium level of CR $\bar{y}_{A}^{\mathrm{eq}}$ is not affected by the value of $D_{B}$; however, a larger $D_{B}$ leads to a faster convergence to this equilibrium level.

Next, we study the effect of forward reaction rate constant $\kappa_{f}$ on the CR. We adopt $\kappa_{b}=0$ and $\kappa_{f}\in\{0,0.01,0.1,1,10,\infty\}\times 10^{-17}$ molecule${}^{-1}\cdot$m${}^{3}\cdot$s*-1* where the extreme values of $\kappa_{f}=0,\infty$ correspond to no reaction and instant reaction, respectively. Moreover, we assume $\mathcal{T}_{A}=\{0\}$ and $\mathcal{T}_{B}=\{100\}$ $\mu$s and initial concentrations $C_{A}^{\mathrm{init}}(\mathbf{r})=C_{B}^{\mathrm{init}}(\mathbf{r})=0$. In Fig. 6, we show the corresponding CR of the type-$A$ molecules, $\bar{y}_{A}(t)$, versus time. As can be seen from this figure, for the non-reactive scenario, $\kappa_{f}=0$, the CR is very wide which causes significant ISI for the considered nominal value of symbol duration, i.e., $T^{\mathrm{symb}}=200$ $\mu$s. On the other hand, as $\kappa_{f}$ increases, the CR is shortened and correspondingly ISI is reduced. For example, for the considered default system parameters, i.e., $\kappa_{f}=10^{-17}$, $T^{\mathrm{symb}}=200$ $\mu$s, and sampling time $T^{\mathrm{samp}}=100$ $\mu$s after the beginning of the symbol interval, ISI is reduced by a factor of $14$ compared to the non-reactive case. Furthermore, the reason why the CR $\bar{y}_{A}(t)$ does not suddenly drop to zero after the release of the type-$B$ molecules when $\kappa_{f}\to\infty$ is that the type-$B$ molecules need a finite time to diffuse and reach the type-$A$ molecules in the channel. In this case, the CR can be further shortened only by increasing $D_{B}$, cf. Fig. 5. Note that since we assume $\kappa_{b}=0$, the equilibrium level $\bar{y}_{A}^{\mathrm{eq}}$ is zero.

Finally, we study the effect of backward reaction rate constant $\kappa_{b}\in\{0,1,10,100\}\times 10^{25}$ molecule$\cdot$m${}^{-3}\\ \cdot$s*-1* on the CR when $\kappa_{f}=10^{-17}$ molecule${}^{-1}\cdot$m${}^{3}\cdot$s*-1*. Moreover, we assume $\mathcal{T}_{A}=\{200\}$ and $\mathcal{T}_{B}=\{300\}$ $\mu$s and the initial concentrations are zero. From Fig. 7, we observe that the main effect of increasing $\kappa_{b}$ is that the equilibrium level increases as well. This implies that more noise molecules are generated which can impair the performance of the MC system, cf. Remark 3. Interestingly, the channel dispersion does not significantly change with the values of $\kappa_{b}$. In other words, $\kappa_{b}$ impacts the level of noise molecules in the MC channel and does not considerably affect ISI.

In summary, Figs. 5-7 indicate that MC systems generally perform better by choosing molecules with large diffusion coefficients (for more effective ISI reduction as well as faster transmission), large forward reaction rates (for more effective ISI reduction), and small backward reaction rates (to decrease the number of noise molecules).

V-C Performance Evaluation

In this section, we evaluate the performance of the modulation and detection schemes proposed in this paper. We consider both 1TM and 2TM receivers. As benchmarks for non-reactive MC systems, we consider the following three modulation schemes: i) conventional OOK modulation with 1TM receiver [2], ii) MoSK modulation with 2TM receiver [22], and iii) the proposed OOK modulation with 2TM receiver where $\tau_{1}=0$666When the signaling molecules are not reactive, the proposed OOK modulation scheme reduces to sending the same data over two orthogonal channels using conventional OOK modulation which was previously studied in [36]. For a fair comparison, unlike [36] which employs $y(t)=y_{A}(t)-y_{B}(t)$ as the received signal, we perform joint ML detection based on $y_{A}(t)$ and $y_{B}(t)$. Note that based on the data processing inequality, optimal joint processing of $y_{A}(t)$ and $y_{B}(t)$ in general yields a better performance than the processing of any function of $y_{A}(t)$ and $y_{B}(t)$ including $y(t)$. Moreover, in the non-reactive case, choosing $\tau_{1}=0$ yields the best performance since both $y_{A}(t)$ and $y_{B}(t)$ attain their respective peak values at the same time.. For reactive MC systems, we consider the three modulation schemes introduced in Section II-A, namely conventional MoSK, the proposed OOK, and the proposed OSK. Note that the optimal and suboptimal detectors derived in Section III are applicable to all of the above binary modulation schemes in non-reactive and reactive MC systems.

In Fig. 8, we show CR $\bar{y}_{A}(t)$ and $\bar{y}_{B}(t)$ versus time for data sequence $\mathbf{s}=[0,1,0,1,0]$. We assume $\kappa_{f}=10^{-17}$ molecule${}^{-1}\cdot$m${}^{3}\cdot$s*-1* and $\kappa_{b}=10^{26}$ molecule$\cdot$m${}^{-3}\cdot$s*-1* which leads to a equilibrium level of $\bar{y}_{i}^{\mathrm{eq}}=1.65$ for the reactive MC system. Since this equilibrium level can be interpreted as the expected number of noise molecules within the receiver volume, cf. Remark 3, we assume that the non-reactive system is affected by the same fixed level of noise molecules. As can be seen from Fig. 8, for non-reactive signaling, $\bar{y}_{i}(t)$ increases over time which is due to the significant ISI from previous symbol intervals. On the contrary, for reactive signaling, the peak of $\bar{y}_{i}(t)$ remains almost constant. This is due to the fact that signaling molecules participate in a degradation reaction and cancel each other out. Moreover, employing reactive signaling molecules provides a huge potential for generating pulse shapes with desirable properties, i.e., causing minimum ISI. For the reactive systems, we observe from Fig. 8 that the CRs of the proposed OOK and OSK modulation schemes are narrower than those of conventional MoSK modulation which results in lower ISI. We emphasize that although, in this paper, we assumed an unbounded simulation environment, we expect that degradation via reactive molecules becomes even more important in a bounded environment where the accumulation of molecules can significantly contaminate the channel.

Now, we study the performance of the considered modulation and detection schemes in terms of BER. To this end, we simulate $10^{5}$ realizations of blocks of $K=10$ symbols and use $L=5$ for the proposed ML detector. We emphasize that, for simulation of $y_{A}$ and $y_{B}$, the memory of all previous symbol intervals is considered, whereas for the ML detector, only a memory length of size $L=5$ is assumed for simplicity.

We first focus on a reactive MC system. Fig. 9 shows the BER vs. the number of released molecules, $N^{\mathrm{tx}}_{A}=N^{\mathrm{tx}}_{B}\triangleq N^{\mathrm{tx}}$, for the genie-aided (GA) lower bound from Proposition 1 and Corollary 1, the ML-based detector with estimated ISI, and the suboptimal ISI-neglecting detector in (13) and (14). Fig. 9 a) and b) show the results for 1TM and 2TM receivers, respectively. The following observations can be made in Fig. 9:

• •

Let us first focus on the ML-based detector with estimated ISI. The ML-based detector performs very close to the corresponding genie-aided lower bound for all considered modulation schemes and both receiver types. For MoSK and OSK modulation, since information is encoded in both molecule types, they are expected to perform well for 2MT receivers. Fig. 9 supports this expectation. Interestingly, the proposed OSK modulation performs well even if a simple 1TM receiver is used whereas the performance of conventional MoSK significantly deteriorates. On the other hand, since information is encoded in type-$A$ molecules in OOK modulation777Although the proposed OOK modulation for reactive MC systems employs two types of molecules, the information is primarily encoded in the concentration of type-$A$ molecules. In fact, type-$B$ molecules are mainly used for ISI reduction, see Fig. 8., it is expected to perform well for both the 1TM receiver (that observes type-$A$ molecules) and the 2TM receiver. This is also consistent with the results shown in Fig. 9. Overall, we observe from Fig. 9 that for all considered modulation schemes, the 2TM receiver outperforms the 1TM receiver by exploiting the diversity gain that the observation of two molecule types offers. This additional gain of the 2TM receiver over the 1TM receiver is significant for MoSK and OSK modulation but is smaller for OOK modulation.

• •

Now, we focus on the suboptimal ISI-neglecting detector. For the conventional MoSK modulation, the performance of the suboptimal ISI-neglecting detector is not good for both 1TM and 2TM receivers. In fact, for MoSK modulation, the BER curves of the ISI-neglecting detector for both 1TM and 2TM receivers have an error floor at large $N^{\mathrm{tx}}$. The reason for this behavior is that this suboptimal scheme was proposed for the case when ISI is negligible whereas ISI is not sufficiently suppressed with conventional MoSK modulation. For the proposed OOK modulation, the ISI-neglecting detector performs well with respect to the corresponding genie-aided bound only if a 1TM receiver is used. For the 2TM receiver and OOK modulation, the ISI-neglecting detector in (13) does not perform well since this detector was designed for a modulation that encodes information in the concentration of two types of molecules which is not the case for OOK modulation. Finally, for the proposed OSK modulation, we observe that the ISI-neglecting detector has good performance relative to the genie-aided bound for both the 1TM and 2TM receivers. In fact, OSK modulation is well suited for 2TM receivers since at the sampling time, $\bar{y}_{A}^{(1)}(\mathbf{s})\gg\bar{y}_{B}^{(1)}(\mathbf{s})$ and $\bar{y}_{A}^{(0)}(\mathbf{s})\ll\bar{y}_{B}^{(0)}(\mathbf{s})$ hold. OSK modulation also performs well for 1TM receivers since at the sampling time, $\bar{y}_{A}^{(1)}(\mathbf{s})\gg\bar{y}_{A}^{(0)}(\mathbf{s})$ holds, see Fig. 8. These conditions hold due to effective ISI reduction for OSK modulation.

• •

In summary, among the considered modulation and detection schemes, OSK modulation with the suboptimal ISI-neglecting detector provides a favorable tradeoff between performance and complexity for both receiver types. Hereby, if a simpler 1TM receiver is adopted (compared to the 2TM receiver), the CR has to be acquired to determine the detection threshold $\gamma$. On the other hand, if a more complex 2TM receiver is adopted, the ISI-neglecting detector does not require any knowledge of the CR for its operation.

In Fig. 10, we compare the performance of reactive and non-reactive MC systems. In particular, Fig. 10 shows the BER vs. the number of released molecules, $N^{\mathrm{tx}}_{A}=N^{\mathrm{tx}}_{B}\triangleq N^{\mathrm{tx}}$, for the ML-based detector. We summarize our observations from Fig. 10 in the following:

• •

From Fig. 10, we see that the MC system with reactive signaling molecules has superior performance compared to the MC system with non-reactive signaling molecules for all considered modulation schemes and moderate-to-large $N^{\mathrm{tx}}$. This is due to the fact that for the adopted symbol duration, ISI is significantly reduced with reactive signaling; by contrast, for non-reactive signaling, the ISI is severe, particularly for large $N^{\mathrm{tx}}$, and thus limits performance. In fact, at large $N^{\mathrm{tx}}$, the slope of the BER curve for the proposed OOK and OSK modulation schemes is larger than that for the benchmark schemes which is due to their effective use of chemical reactions to reduce ISI.

• •

For OOK modulation with both reactive and non-reactive molecules, although the 2TM receiver outperforms the 1TM receiver, the slopes of the BER curves are similar for both receiver types. This behavior can be justified as follows. For non-reactive MC systems, the BER curve saturates for large $N^{\mathrm{tx}}$ due to severe ISI which ultimately limits the slope of the BER curves for both 1TM and 2TM receivers (and all considered modulation schemes). For reactive MC systems, type-$A$ molecules are the primary signaling molecules for the proposed OOK modulation and type-$B$ molecules are mainly used for cleaning the channel for the next release of type-$A$ molecules. Therefore, the observation of type-$B$ molecules at the receiver contains little information for the symbol sent in the current interval which explains the marginal gain of the 2TM receiver over the 1TM receiver in this case.

• •

For reactive MC systems, the slope of the BER curve for OSK modulation with 2TM receivers is larger than that of OOK modulation with 1TM and 2TM receivers. This is due to the effective use of two types of molecules in OSK modulation whose observation at the receiver provides an increased diversity gain (the slope of BER curve).

• •

Interestingly, for MoSK modulation in reactive MC systems, the slope of the BER curve is higher for small-to-moderate $N^{\mathrm{tx}}$ and decreases for large $N^{\mathrm{tx}}$. The reason for this behavior is attributed to high BER of those sequences that contain multiple consecutive ones or zeros. These sequences introduce a severe ISI which cannot be cope with in reactive MC systems using MoSK modulation.

The analytical investigation of the diversity gains of these schemes is an interesting topic for future research.

Recall that for the simulation results presented in Figs. 9 and 10, we assume blocks of length $K=10$ and generate $y_{i}(\mathbf{s}),\,\,i\in\{A,B\}$, without any CR truncation, i.e., the molecules observed at the tenth symbol interval may originate from the release in the first symbol interval. However, for detection, we take into account a channel memory length of $L=5$ for simplicity. In Figs. 11 and 12, we study the impact of parameters $L$ and $K$. In particular, in Fig. 11, we plot the BER vs. $L$ for $K=10$ and $N^{\mathrm{tx}}_{A}=N^{\mathrm{tx}}_{B}=5\times 10^{3}$. We observe from this figure that for the proposed OOK and OSK modulations, the ISI is negligible even for $L>3$ due to efficient ISI reduction enabled by chemical reactions. For the benchmark schemes, the impact of ISI is still observable for $L$ between $6$ and $8$ symbols. Overall, Fig. 11 implies that, for the proposed modulation schemes, the actual channel memory is smaller than the considered block length of $K=10$. Therefore, we expect that considering larger block lengths does not significantly affect the performance. This is validated in Fig. 12 where the BER is plotted vs. block length $K$ for $L=5$ and $N^{\mathrm{tx}}_{A}=N^{\mathrm{tx}}_{B}=5\times 10^{3}$. We observe from this figure that for the proposed OOK and OSK modulations in reactive MC systems, the BER remains almost constant as $K$ increases, which is due to the efficient ISI reduction inherent to these modulation schemes. For the benchmark schemes, the BER slightly increases as the block length increases. Therefore, considering $K>10$ does not change the BER of the proposed OOK and OSK modulation schemes whereas it may slightly deteriorate the performance of the benchmark schemes.

VI Conclusions and Future Work

We studied an MC system that employs two types of molecules for signaling. These molecules may participate in a reversible bimolecular reaction in the channel. We also considered two receiver types, namely 1TM and 2TM receivers. For this system, we proposed two binary modulation schemes, i.e., a modified OOK and OSK, and derived an optimal genie-aided ML detector and two suboptimal detectors, an ML-based detector and an ISI-neglecting detector. Knowledge of the CR was needed for all detectors except the ISI-neglecting detector for the 2TM receiver as well as for performance evaluation. However, deriving the CR of a reactive MC system is challenging since the underlying PDEs that describe the reaction-diffusion mechanism are coupled and non-linear. To address this issue, we developed a numerical algorithm for efficient computation of the CR. The accuracy of this algorithm has been validated via particle-based simulations. Moreover, simulation results confirmed the superior performance of reactive systems over non-reactive systems due to efficient ISI reduction enabled by chemical reactions. In particular, the BER performance of the proposed OOK and OSK modulations improves as the number of released molecules increases whereas the BERs of conventional OOK and MoSK modulations with non-reactive signaling molecules saturate due to severe ISI. In addition, for all considered modulation schemes and receiver types, the ML-based detectors yielded a similar performance as the corresponding genie-aided detectors. Furthermore, for OSK modulation, we observed that the ISI-neglecting detector provided a favorable tradeoff between performance and complexity and hence is a desirable candidate for practical MC systems with limited computational capabilities. The results of this paper can be further extended to non-binary signaling and other MC environments, e.g., considering other geometries and incorporating flow.

Appendix A Proof of Proposition 1

The log likelihood ratio (LLR) for problem (9) can be written as

[TABLE]

Due to the monotonicity of the logarithm, the ML problem in (9) can be rewritten as $\text{LLR}\overset{s=1}{\underset{s=0}{\gtreqless}}0$. Defining $\alpha(\mathbf{s})$ and $\beta(\mathbf{s})$ as in Proposition 1, we obtain (10) which concludes the proof.

Appendix B Proof of Corollary 1

The proof is similar to that given in Appendix A for 2TM receivers. In particular, the LLR for detection problem (11) can be written as

[TABLE]

The ML problem in (8) can be rewritten as $\text{LLR}\overset{s=1}{\underset{s=0}{\gtreqless}}0$ which in turn can be reformulated into $y_{A}\overset{s=1}{\underset{s=0}{\gtreqless}}\gamma(\mathbf{s})$ where $\gamma(\mathbf{s})$ is defined in Corollary 1. This completes the proof.

Appendix C Proof of Lemma 1

The concentration changes due to diffusion and reaction given in (19) can be written in form of

[TABLE]

where $f_{i}^{\mathrm{df}}(\mathbf{r},\tilde{t})=D_{i}\nabla^{2}C_{i}(\mathbf{r},\tilde{t})$ and $f_{i}^{\mathrm{rc}}(\mathbf{r},\tilde{t})=\kappa_{b}-\kappa_{f}C_{A}(\mathbf{r},\tilde{t})C_{B}(\mathbf{r},\tilde{t})$. Let us define $\bar{f}_{i}^{\mathrm{df}}(\mathbf{r},\tilde{t})=D_{i}\nabla^{2}C_{i}^{\mathrm{df}}(\mathbf{r},\tilde{t})$ where $C_{i}^{\mathrm{df}}(\mathbf{r},\tilde{t})$ denotes the concentration of type-$i$ molecules assuming during interval $\tilde{t}\in[t,t+\Delta t]$ diffusion is present and reaction is absent. Similarly, let us define $\bar{f}_{i}^{\mathrm{rc}}(\mathbf{r},\tilde{t})=\kappa_{b}-\kappa_{f}C_{A}^{\mathrm{rc}}(\mathbf{r},\tilde{t})C_{B}^{\mathrm{rc}}(\mathbf{r},\tilde{t})$ where $C_{i}^{\mathrm{rc}}(\mathbf{r},\tilde{t})$ denotes the concentration of type-$i$ molecules assuming during interval $\tilde{t}\in[t,t+\Delta t]$ reaction is present and diffusion is absent. For mathematical rigorousness, we distinguish the following two cases:

Case 1: In this case, we assume that the transmitter releases no molecules within interval $[t,t+\Delta t]$, i.e., $(\mathcal{T}_{A}\cup\mathcal{T}_{B})\cap[t,t+\Delta t]$ is empty. Thereby, by defining $\epsilon_{i}^{\mathrm{x}}(\mathbf{r},\tilde{t})\triangleq f_{i}^{\mathrm{x}}(\mathbf{r},\tilde{t})-\bar{f}_{i}^{\mathrm{x}}(\mathbf{r},\tilde{t})$, we have

[TABLE]

Next, we show that as $\Delta t\to 0$, $B_{i}^{\mathrm{x}}(\Delta t)$ becomes negligible compared to $A_{i}^{\mathrm{x}}(\Delta t)$. Under assumption A1, i.e., $C_{i}(\mathbf{r},t)\neq 0$ and $\nabla^{2}C_{i}(\mathbf{r},t)\neq 0$, we obtain $\bar{f}_{i}^{\mathrm{x}}(\mathbf{r},t)\neq 0$ which leads to $\underset{\Delta t\to 0}{\lim}A_{i}^{\mathrm{x}}(\Delta t)\sim o(\Delta t)$. Recall that $\epsilon_{i}^{\mathrm{x}}(\mathbf{r},\tilde{t})$ characterizes the error in decoupling the reaction and diffusion phenomena at time $\tilde{t}>t$. Therefore, for $B_{i}^{\mathrm{x}}(\Delta t)$, we have $\epsilon_{i}^{\mathrm{x}}(\mathbf{r},t)=0$ which leads to $\underset{\Delta t\to 0}{\lim}B_{i}^{\mathrm{x}}(\Delta t)\sim o(\Delta t^{p})$ for some $p>1$. In summary, we have $\underset{\Delta t\to 0}{\lim}\frac{B_{i}^{\mathrm{x}}(\Delta t)}{A_{i}^{\mathrm{x}}(\Delta t)}\sim o(\Delta t^{p-1})\to 0$.

Case 2: In this case, we assume the transmitter releases molecules within interval $[t,t+\Delta t]$, i.e., $(\mathcal{T}_{A}\cup\mathcal{T}_{B})\cap[t,t+\Delta t]$ is non-empty. For this case, the conditions on set $\mathcal{T}_{i}$ introduced in assumption A2 are needed such that the argument in Case 1 is also applicable to Case 2. In particular, under condition $\mathcal{T}_{i}\subset\mathcal{T}$, without ambiguity, each release by the transmitter is taken into account only in one of the time slots. Moreover, due to condition $\varepsilon\ll\Delta t$, the impact of diffusion and reaction on molecule concentration in interval $[t+\Delta t-\varepsilon,t+\Delta t]$ is negligible compared to that in $[t,t+\Delta t-\varepsilon)$. Furthermore, assuming that $\mathcal{T}_{A}\cap\mathcal{T}_{B}$ is empty is a reasonable assumption since simultaneous release of both types of molecules leads to their degradation via the forward reaction. These conditions enable us to directly apply the proof for $\underset{\Delta t\to 0}{\lim}\frac{B_{i}^{\mathrm{x}}(\Delta t)}{A_{i}^{\mathrm{x}}(\Delta t)}\sim o(\Delta t^{p-1})\to 0$ provided for Case 1 to Case 2.

In summary, under assumptions A1 and A2, we can conclude that as $\Delta t\to 0$, the relative effect arising from the coupling of diffusion and reaction becomes negligible, and the overall concentration change is simply a superposition of the concentration changes due to the individual phenomena. This concludes the proof.

Appendix D Proof of Proposition 2

Eq. (2) is obtained by applying (20) from Lemma 1 in (19). Under assumption A1, we directly obtain $\bar{G}_{i}(\mathbf{r})\triangleq\int_{\tilde{t}=t}^{t+\Delta t}G_{i}(\mathbf{r},t+\Delta t)\mathrm{d}\tilde{t}=\sum_{t_{i}\in\mathcal{T}_{i}}N_{i}^{\mathrm{tx}}\delta_{\mathrm{d}}(\mathbf{r})\delta_{\mathrm{k}}(t+\Delta t-\varepsilon-t_{i})$. In the following, we derive expressions for $C_{i}^{\mathrm{df}}(\mathbf{r},t+\Delta t)$ and $C_{i}^{\mathrm{rc}}(\mathbf{r},t+\Delta t)$ as a function of initial conditions at time $t$, respectively. Without loss of generality and in order to simplify the notations, we assume the initial time as zero and obtain the concentrations for $t>0$, denoted by $C_{i}^{\mathrm{df}}(\mathbf{r},t)$ and $C_{i}^{\mathrm{rc}}(\mathbf{r},t)$.

Diffusion: The following diffusion-only PDE has to be solved for $t>0$ and arbitrary initial conditions $C_{i}(\mathbf{r},0)$ at initial time $t=0$

[TABLE]

For one-dimensional diffusion, this problem is solved in [13, Chapter 1.7] by transforming it into the Laplace domain with respect to space variable $\mathbf{r}$. The generalization of this result to the $n\in\{2,3\}$ dimensional environment given in (22aa) is straightforward. Therefore, we skip the derivation due to space constraints and refer the readers to [13, Chapter 1.7].

Reversible Bimolecular Reaction: The PDE for the reaction-only case is given by

[TABLE]

Subtracting (37ab) from (37aa), we obtain

[TABLE]

which has the solution $C_{A}^{\mathrm{rc}}(\mathbf{r},t)-C_{B}^{\mathrm{rc}}(\mathbf{r},t)=c_{1}(\mathbf{r})$ where $c_{1}(\mathbf{r})$ is a constant with respect to $t$ and is set to $c_{1}(\mathbf{r})=C_{A}^{\mathrm{rc}}(\mathbf{r},0)-C_{B}^{\mathrm{rc}}(\mathbf{r},0)$ to satisfy the initial conditions at $t=0$. Substituting $C_{B}^{\mathrm{rc}}(\mathbf{r},t)=C_{A}^{\mathrm{rc}}(\mathbf{r},t)-c_{1}(\mathbf{r})$ into (37aa) leads to

[TABLE]

This equation is in form of a Bernoulli equation (or also in form of a Riccati equation) which can be solved by rewriting (40) as

[TABLE]

Integrating both sides of the above equation and using [37]

[TABLE]

where $\Delta=\sqrt{b^{2}-4ac}$, we obtain

[TABLE]

where $c_{2}(\mathbf{r})=\sqrt{\kappa_{f}^{2}c_{1}^{2}(\mathbf{r})+4\kappa_{f}\kappa_{b}}$ and $\tilde{c}_{4}(\mathbf{r})$ is a constant with respect to time. Rewriting the above equation leads to

[TABLE]

Using the initial condition $C_{A}^{\mathrm{rc}}(\mathbf{r},0)$, we obtain $\tilde{c}_{4}(\mathbf{r})=\frac{1}{c_{2}(\mathbf{r})}\mathrm{ln}\left(\frac{c_{2}(\mathbf{r})+\kappa_{f}c_{3}(\mathbf{r})}{c_{2}(\mathbf{r})-\kappa_{f}c_{3}(\mathbf{r})}\right)$ with $c_{3}(\mathbf{r})=C_{A}^{\mathrm{rc}}(\mathbf{r},0)+C_{B}^{\mathrm{rc}}(\mathbf{r},0)$. This leads to $C_{A}^{\mathrm{rc}}(\mathbf{r},t)$ given in (22ab) after defining $c_{4}(\mathbf{r})\triangleq\exp(c_{2}(\mathbf{r})\tilde{c}_{4}(\mathbf{r}))=\frac{c_{2}(\mathbf{r})-\kappa_{f}c_{3}(\mathbf{r})}{c_{2}(\mathbf{r})+\kappa_{f}c_{3}(\mathbf{r})}$ for notational simplicity. Using $C_{B}^{\mathrm{rc}}(\mathbf{r},t)=C_{A}^{\mathrm{rc}}(\mathbf{r},t)-c_{1}(\mathbf{r})$, we straightforwardly obtain $C_{B}^{\mathrm{rc}}(\mathbf{r},t)$ given in (22ac).

The above results are concisely given in Proposition 2 which concludes the proof.

Appendix E Proof of Corollary 2

Because of the geometrical symmetry of the problem, the concentration for the MC system under consideration is only a function of $r\triangleq\|\mathbf{r}\|$. In the following, we present the proof for one-, two-, and three-dimensional environments, respectively.

One-Dimensional Environment: In this case, we obtain (27) by substituting $n=1$ in (22aa) as

[TABLE]

where for equality $(a)$, we used $\tilde{r}=|\tilde{r}_{1}|$. Defining $W_{i}(r,\tilde{r})$ as in (28) leads to (27) for $n=1$.

Two-Dimensional Environment: For $n=2$, the simplification of $C_{i}^{\mathrm{df}}(\mathbf{r},t+\Delta t)$ requires transforming an integral from Cartesian coordinates to polar coordinates, i.e.,

[TABLE]

where $f(r_{1},r_{2})$ is an arbitrary function and $r\geq 0$ and $\theta\in[0,2\pi]$ are the variables of the polar coordinates. Moreover, without loss of generality, we consider $\mathbf{r}=(r,0)$ in order to simplify (22aa) as follows

[TABLE]

where in equality $(a)$, we used the identity $\int_{x=0}^{2\pi}\exp(c\cos(x))\mathrm{d}x=2\pi I_{0}(c)$ [37]. Defining $W_{i}(r,\tilde{r})$ as in (28) leads to (27) for $n=2$.

Three-Dimensional Environment: In a similar manner, the simplification of $C_{i}^{\mathrm{df}}(\mathbf{r},t+\Delta t)$ when $n=3$ requires transforming an integral from Cartesian coordinates to spherical coordinates, i.e.,

[TABLE]

where $f(r_{1},r_{2},r_{3})$ is an arbitrary function and $r\geq 0$, $\phi\in[0,\pi]$, and $\theta\in[0,2\pi]$ are the variables of the spherical coordinates. Here, without loss of generality, we consider $\mathbf{r}=(0,0,r)$ in order to simplify (22aa) as follows

[TABLE]

where for equality $(a)$, we used the identity $\int_{x=0}^{\pi}\sin(x)\exp(c\cos(x))\mathrm{d}x=\frac{2\sinh(c)}{c}$. Defining $W_{i}(r,\tilde{r})$ as in (28) leads to (27) for $n=3$ and completes the proof.

Appendix F Stochastic Simulation

Stochastic simulation is used to verify the accuracy of Algorithm 1 for CR computation and the Poisson statistics in (8). We refer the readers to [15, 38, 19, 17] for an overview of general stochastic simulation of reaction-diffusion systems. In the following, we explain the particular particle-based simulator that is used for performance verification in this paper.

For particle-based simulation, it is assumed that the spatial orientations and internal energy levels of the molecules fluctuate on time scales that are faster than the diffusion and reaction processes that are of interest [15]. Therefore, it suffices to track only the positions of the individual molecules during the simulation time. Let $\mathbf{r}_{i}(t)=(r_{1,i}^{(m)}(t),r_{2,i}^{(m)}(t),r_{3,i}^{(m)}(t))$ denote the coordinate of a specific type-$i$ molecule, indexed by $m$, at time $t$. Moreover, let $M_{i}$ denote a variable which tracks the total number of type-$i$ molecules generated in the simulation environment. In the following, we explain how we update the position of the molecules for the considered release, diffusion, and reaction mechanisms. For convenience of presentation, we focus on a three-dimensional environment.

F-A Transmitter Release

For instantaneous release from a point-source transmitter, we simply place $N^{\mathrm{tx}}_{i}$ type-$i$ molecules at the origin, i.e., $\mathbf{r}_{i}^{(m)}(t)=(0,0,0),\,\,m=M_{i}+1,\dots,M_{i}+N^{\mathrm{tx}}_{i}$, at any release time instant $t\in\mathcal{T}_{i}$. Variable $M_{i}$ is updated by $M_{i}+N^{\mathrm{tx}}_{i}$ after each release of type-$A$ molecules by the transmitter.

F-B Diffusion

According to Brownian dynamics, the position of each molecule at time $t+\Delta t$ is updated as [15]

[TABLE]

where $\Delta r_{1,i}^{(m)},\Delta r_{2,i}^{(m)},\Delta r_{3,i}^{(m)}\sim\mathcal{N}(0,1)$ and $\mathcal{N}(\mu,\sigma^{2})$ denotes a normal RV with mean $\mu$ and variance $\sigma^{2}$.

F-C Forward Reaction

The forward reaction in (7), i.e., $A+B\overset{\kappa_{f}}{\rightarrow}\varnothing$, is a second order bimolecular reaction. One approach for stochastic simulation of these reactions is based on the fact that a reaction occurs within interval $[t,t+\Delta t]$ when the reactant molecules are within a certain binding radius $\rho_{b}$ [15]. Therefore, in our simulation, we can simply remove one type-$A$ and one type-$B$ molecule if their distance is less than $\rho_{b}$. In other words, we remove each pair of molecules whose position variables $\mathbf{r}_{A}^{(m)}(t)$ and $\mathbf{r}_{B}^{(m^{\prime})}(t)$ satisfy $\|\mathbf{r}_{A}^{(m)}(t)-\mathbf{r}_{B}^{(m^{\prime})}(t)\|<\rho_{b}$. Unfortunately, the exact value of $\rho_{b}$ cannot be found analytically in general and depends on reaction rate constant $\kappa_{f}$ and the choice of $\Delta t$. Nevertheless, for the two special cases of $\Delta t\to 0$ and $\Delta t\to\infty$, the following simple relations are available [15, Eqs. (19) and (20)]

[TABLE]

where $\rho_{\mathrm{rms}}=\sqrt{2(D_{A}+D_{B})\Delta t}$ is the mutual root mean square step length of type-$A$ and type-$B$ molecules. Using the simplified formula for $\rho_{b}$ in (51), we obtain equivalent conditions for $\rho_{\mathrm{rms}}\ll\rho_{b}$ and $\rho_{\mathrm{rms}}\gg\rho_{b}$ as

[TABLE]

where $T^{\mathrm{crt}}=\frac{\kappa_{f}^{2}}{32\pi^{2}(D_{A}+D_{B})^{3}}$. In order to be able to use the closed-form relation in (51), we have to choose either large $\Delta t\gg\frac{9}{4}T^{\mathrm{crt}}$ or small $\Delta t\ll T^{\mathrm{crt}}$. This choice depends on the system under consideration and the total length of the simulation interval.

F-D Backward Reaction

The backward reaction in (7), i.e., $\varnothing\overset{\kappa_{b}}{\rightarrow}A+B$, is in form of a zeroth order reaction. Suppose that the simulation environment is a cube of volume $V=L^{3}$. Moreover, let RV $N(t)$ denote the number of times that the backward reaction occurs in a time interval $[t,t+\Delta t]$. Then, $N(t)$ follows a Poisson distribution [15]

[TABLE]

Here, we have to be careful about the positions of the type-$A$ and type-$B$ molecules that are generated via the backward reaction. In particular, if $\rho_{\mathrm{rms}}\ll\rho_{b}$ holds and we put the molecules generated by the backward reaction at the same location, then these molecules will directly participate in the forward reaction before they can diffuse away regardless of the value of $\rho_{b}$ [15, 19]. In order to avoid the automatic degradation of type-$A$ and type-$B$ molecules, the type-$A$ and type-$B$ product molecules are initially separated by a fixed distance, referred to as the unbinding radius $\rho_{u}$, which is larger than $\rho_{b}$ . Let $\mathbf{l}$ denote the center of the cube that defines the simulation environment. Then, the zeroth order reaction can be simulated by inserting each of the $N(t)$ molecules of type-$i$ at uniformly random positions $\mathbf{r}_{i}^{(m)}(t),\,\,m=M_{i}+1,\dots,M_{i}+N(t)$, obtained from [15]

[TABLE]

where $\Delta r_{1,A}^{(m)},\Delta r_{2,A}^{(m)},\Delta r_{3,A}^{(m)}\sim\mathcal{U}(-0.5,0.5)$ and $\mathcal{U}(a,b)$ is an RV uniformly distributed in the interval $[a,b]$. Moreover, we have $\Delta r_{1,B}^{(m)}=\sin(\theta^{(m)})\cos(\phi^{(m)}),\Delta r_{2,B}^{(m)}=\sin(\theta^{(m)})\sin(\phi^{(m)})$, and $\Delta r_{3,B}^{(m)}=\cos(\theta^{(m)})$ where $\theta^{(m)}\sim\mathcal{U}(0,\pi)$ and $\phi^{(m)}\sim\mathcal{U}(0,2\pi)$. On the other hand, if $\rho_{\mathrm{rms}}\gg\rho_{b}$ holds and diffusion is simulated after the backward reaction in the adopted simulator, it is very likely that the molecules diffuse out of the binding radius after one diffusion step. In this case, the value of the unbinding radius is not important and without loss of generality, we can choose $\rho_{u}=0$ which leads to $\mathbf{r}_{A}^{(m)}(t)=\mathbf{r}_{B}^{(m)}(t)$. Variable $M_{i}$ is updated by $M_{i}+N(t)$ at the end of this step.

Algorithm 2 summarizes the main steps required for particle-based simulation of the considered MC system.