Diffusive Molecular Communications with Reactive Molecules: Channel Modeling and Signal Design
Vahid Jamali Nariman Farsad Robert Schober Andrea Goldsmith





Diffusive Molecular Communications with Reactive Molecules: Channel Modeling and Signal Design
    
Vahid Jamali, Student Member IEEE, Nariman Farsad, Member IEEE, 
 Robert Schober, Fellow IEEE,  and Andrea Goldsmith, Fellow IEEE 


This paper has been presented in part at IEEE ICC 2018 <cit.>.
This work was supported in part by the German Research Foundation under Project SCHO 831/7-1, in part by the Friedrich-Alexander
	University Erlangen-Nürnberg under the Emerging Fields Initiative, in part by the STAEDTLER Foundation, in part by the NSF Center for Science of Information under grant NSF-CCF-0939370, and in part by German Academic Exchange Service (DAAD).
V. Jamali and R. Schober are with the Institute for Digital Communications,
Friedrich-Alexander University (FAU), Erlangen D-91058, Germany (e-mail:
vahid.jamali@fau.de; robert.schober@fau.de).
N. Farsad and A. Goldsmith are with the Electrical Engineering Department,
Stanford University, Stanford, CA 94305 USA (e-mail: nfarsad@stanford.edu;
andrea@ee.stanford.edu). 

    Received 18 December 2018 / Accepted 26 February 2019
===============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




 
MCmolecular communication
ISIinter-symbol interference
MLmaximum likelihood
CRchannel response
PDEpartial differential equation
FDMfinite difference method
BERbit error rate
RVrandom variable
PMFprobability mass function
LLRlog likelihood ratio



This paper focuses on 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 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 ML detector for general binary modulation is derived under the assumption that the detector has perfect knowledge of the ISI-causing sequence. The performance of this genie-aided ML detector yields an upper bound on the performance of any practical detector. 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 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. Simulation results  obtained using this CR computation algorithm confirm the effectiveness of the proposed modulation and detection schemes. Moreover, these results 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. 


 
Diffusive molecular communications, reactive signaling, detector design, channel modeling, particle-based simulation, and inter-symbol interference.
 




§ INTRODUCTION

Recent advances in biology, nanotechnology, and medicine have given rise to the need for communication in nano/micrometer scale dimensions <cit.>. 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 <cit.>. 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 <cit.> and non-biological <cit.> 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 <cit.>. It has been shown in <cit.> 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 <cit.>, 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 <cit.>. These advantages of the MC system in <cit.> 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 PDE 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  <cit.>, see <cit.> for a recent overview of various channel models that have been developed so far for MC systems. For instance,  in <cit.>, 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   <cit.>.  This approach was applied to MC in <cit.> where the authors employed a 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  <cit.>. However, the computational complexity of these numerical and simulation methods is very high.  In <cit.>, 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 <cit.>. 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-dimensional[One- 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 <cit.> and propagation in a low-depth petri dish or thin tissues can be approximated as two-dimensional <cit.>.]) 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 <cit.>, and avoids their inherent numerical stability issues <cit.>.  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 <cit.>,  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 <cit.> in several directions. First, only 2TM receivers were considered in <cit.>, 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 <cit.>. Finally, unlike <cit.>, 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 <ref>, the considered system is presented.  The proposed optimal and suboptimal detectors are derived in Section <ref> and the proposed algorithm for computing the CR is introduced in Section <ref>. Simulation results are reported in Section <ref>, and conclusions are drawn in Section <ref>.



§ SYSTEM MODEL


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







 §.§ Transmitter


We assume a point-source transmitter located at the origin of the n-dimensional Cartesian coordinate system, i.e., 𝐫=(r_1,…,r_n)=(0,…,0),  n∈{1,2,3}.  The transmitter employs two types of molecules for signaling, namely type-A and type-B molecules. In particular, the transmitter releases i type-i molecules into the channel at time instances t∈𝒯_i, i∈{A,B}. By properly defining 𝒯_i, different modulation schemes can be realized. Let s[k]∈{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) <cit.>, a proposed variant of on-off keying (OOK), and a new modulation scheme named order shift keying (OSK), see Fig. <ref>. 

Conventional MoSK:
For binary one, s[k]=1, the transmitter releases A 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 B type-B molecules at the beginning of the symbol interval and no type-A molecules. For this modulation scheme, we obtain 

lll 
𝒯_A={t|t=(k-1),  ∀k  and   s[k]=1} 

𝒯_B={t|t=(k-1),  ∀k  and   s[k]=0}, 

where  denotes the length of a symbol interval. 

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

lll 
𝒯_A={t|t=(k-1),  ∀k  and  s[k]=1} 

𝒯_B={t|t=(k-1)+τ_1,  ∀k  and  s[k]=1}. 

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 <ref>, this modulation scheme  considerably reduces ISI compared to conventional OOK modulation. A similar idea was introduced in <cit.> where the release patterns of type-A and type-B molecules were studied for efficient pulse shaping. However, the analysis in <cit.> 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 A type-A molecules at the beginning of the symbol interval and B type-B molecules  τ_1 seconds after the start of the symbol interval. In a similar manner, for binary zero, s[k]=0, the transmitter releases B type-B molecules at the beginning of the symbol interval and A type-A molecules  τ_0 seconds after the start of the symbol interval. In particular, we choose τ_1 (τ_0) as the peak of the CR assuming instantaneous release of only A (B) type-A (type-B) molecules at t=0. For this modulation scheme, we obtain 

lll 
𝒯_A={t|t=(k-1)+(1-s[k])τ_0,  ∀k} 

𝒯_B={t|t=(k-1)+ s[k]τ_1,  ∀k}. 

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 c̅_i(t),  i∈{A,B}, denote the expected concentration of type-i molecules at the receiver  when only N_i^tx type-i molecules are released at time t=0. Then, we choose τ_1=argmax_t c̅_A(t) and τ_0=argmax_t c̅_B(t).	



 §.§ Channel


We assume an unbounded n-dimensional environment for n∈{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·s^-1, respectively, and may participate in the following bimolecular reaction

lll 
A+B κ_bκ_f⇌ ∅,


where κ_f and κ_b denote the forward reaction rate constant in molecule^-1·m^n·s^-1,  n∈{1,2,3}, and the backward reaction rate constant in molecule·m^-n·s^-1,  n∈{1,2,3}, respectively. Moreover, symbol ∅ denotes chemical species which are of no interest for communication. Note that (<ref>) includes the reactions considered in <cit.>. Moreover, if  type-B molecules  represent enzymes and only type-A molecules are used for signaling, (<ref>) includes the degradation reaction in <cit.> when the enzyme concentration is constant everywhere and is much larger than the concentration of the type-A molecules such that the reaction in (<ref>) does not change the enzyme concentration. 



 §.§ 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 <ref>,  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 At and Bt 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 it,  i∈{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 memory[Theoretically, 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 𝐬∈{0,1}^L-1 denote the vector of the L-1 previously transmitted symbols. We also refer to 𝐬 as the ISI-causing sequence[For notational simplicity, we drop the symbol index k in the remainder of the paper.]. Therefore, given 𝐬 and s, the number of  type-i molecules counted at the receiver at sample time t_s is modelled as <cit.>

lll 
y_i ∼𝒫(si1+(1-s)i0),  i∈{A,B},


where 𝒫(λ) denotes a Poisson RV with mean λ. Moreover, is is it_s under the condition that the symbol in the current symbol interval is s∈{0,1} and the ISI-causing sequence is 𝐬∈{0,1}^L-1. We note that the Poisson distribution assumed for y_i in (<ref>) is an approximation which has been shown to be accurate for reaction-diffusion processes in the chemistry and physics literature <cit.>. In Section <ref>, we will validate the Poisson model   in (<ref>)  using the particle-based simulator developed in arXivAppendix F.<cit.>.

Note that due to the reaction process, the CR of the considered MC system is is a non-linear function of the transmitted data symbols s and 𝐬. Therefore, we cannot simply compute the CR for one shot transmission and use convolution to capture the effect of  ISI <cit.>. 
In particular, to fully characterize the average behavior of the system, one has to compute the CR for both symbol hypotheses s∈{0,1} and all 2^L-1 possible ISI-causing sequences. In Section <ref>, we derive optimal and suboptimal detectors assuming is,  ∀ s,𝐬, is known. Then, in Section <ref>, we derive an efficient numerical algorithm for computation of CR is,  ∀ s,𝐬.



§ 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.



 §.§ 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

lll 
   = s∈{0,1}argmax  (y_A,y_B|s,𝐬) 

   (a)= s∈{0,1}argmax   f_𝒫(y_A|s,𝐬)f_𝒫(y_B|s,𝐬),


where (·) denotes probability and f_𝒫(x)=λ^x e^-λ/x! is the PMF of a Poisson RV with mean λ. Equality (a) follows from the fact that conditioned on is,  i∈{A,B}, and (s,𝐬), RVs y_A and y_B are independent. 
The optimal detector for 2TM receivers is given in the following proposition.


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

lll 
= 
1,      if   y_A ≥α(𝐬) y_B + β(𝐬) 

0,      otherwise,



where α(𝐬)=1/χ(𝐬)log(B0/B1), 
β(𝐬)=1/χ(𝐬)(A1+B1-A0-B0), and 
χ(𝐬) = log(A1/A0).


The proof is given in Appendix <ref>.


Proposition <ref> 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., α(𝐬) and β(𝐬), depend on the CR is and the ISI-causing sequence 𝐬.

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

lll 
   = s∈{0,1}argmax  (y_A|s,𝐬) = s∈{0,1}argmax   f_𝒫(y_A|s,𝐬).


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


The genie-aided ML detector as a solution of (<ref>) is given by

lll 
= 
1,      if   y_A ≥γ(𝐬) 

0,      otherwise,



where 
γ(𝐬)=(A1-A0)/log(A1/A0).


The proof is given in Appendix <ref>.


Corollary <ref> shows that the optimal detection rule for 1TM receivers is a simple threshold detector.  The detection threshold, γ(𝐬), depends on the CR As and the ISI-causing sequence 𝐬.


We note that the detectors in Proposition <ref> and Corollary <ref> are not limited to the specific modulation schemes introduced in Section <ref> 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  is, of course.


The optimal detectors in Proposition <ref> and Corollary <ref> provide performance upper bounds since they require  knowledge of the ISI-causing sequence 𝐬 which is not available in practice. Therefore, in the following, we propose some practical suboptimal detectors.



 §.§ Suboptimal Detectors

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



  §.§.§ ML-Based Detectors
 Here, we assume that the receiver employs the detectors in Proposition <ref> or Corollary <ref> 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 <ref> 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.



  §.§.§ 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 <ref> and Corollary <ref> assuming an ISI-free channel and propose to use them when the ISI is sufficiently small. In particular, for an ISI-free channel, parameters α(𝐬),  β(𝐬), and γ(𝐬)  do not depend on previous symbols and are denoted by α, β, and γ, respectively.  This simplifies the detector in Proposition <ref> for 2TM receivers to

lll 
ŝ= 
1,      if   y_A ≥y_B 

0,      otherwise,
 


and the detector in Corollary <ref> for 1TM receivers to

lll 
ŝ= 
1,      if   y_A ≥γ

0,      otherwise,
 


where in (<ref>), we further assume symmetry with respect to the two molecule types, i.e., D_A=D_B and N_A=N_B hold[Note that if D_A=D_B and A=B hold, we obtain  A1=B0≜ a and A0=B1≜ b. This leads to α=1/log(a/b)log(a/b)=1 and β=1/log(a/b)(a+b-b-a)=0 in Proposition <ref>.]. 


Note that the detector in (<ref>) does not need knowledge of the CR which makes it suitable for MC systems with limited computational capabilities. However, for the detector in (<ref>), the computation of γ requires knowledge of the CR, see e.g. <cit.> for the optimal value of γ for an ISI-free MC channel.  On the other hand, the advantage of the detector in (<ref>) over that in (<ref>) is that the corresponding receiver is simpler because it needs to  be able to count only one type of molecule.




§ 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.




 §.§ Problem Statement

Let A𝐫t and  B𝐫t denote the concentrations of type-A and type-B molecules in point 𝐫 and at time t.  Considering a passive receiver, it is obtained as

lll 
it = 𝐫∈𝒱^rx∫ i𝐫t d𝐫,  i∈{A,B},


where 𝒱^rx is the set of points within the receiver space.  Concentrations A𝐫t and  B𝐫t  can be found using the following reaction-diffusion equations 
<cit.>

lll 
∂A𝐫t/∂t = 
   D_A ∇^2 A𝐫t 

     - κ_f A𝐫tB𝐫t + κ_b+ A𝐫t  

∂B𝐫t/∂t = 
   D_B ∇^2 B𝐫t 

    - κ_f A𝐫tB𝐫t + κ_b + B𝐫t,       

where ∇^2=∂^2/∂ r_1^2+⋯+∂^2/∂ r_n^2,  n∈{1,2,3}, is the Laplace operator. Moreover,  i𝐫t= ∑_t_i∈𝒯_iiδ_d(𝐫,t - t_i),  i∈{A,B}, represents the concentration of type-i molecules that are released by the transmitter into the channel where δ_d(𝐫,t)=δ_d(r_1)⋯δ_d(r_n)δ_d(t),  n∈{1,2,3}, and δ_d(·) is the Dirac delta function. In addition, without loss of generality, we assume  initial condition ϵ→ 0lim  i𝐫-ϵ=i𝐫 where i𝐫 is the spatial concentration of type-i molecules before transmission starts at t=0. 



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., i𝐫≜ C_i^eq. Therefore, substituting  ∂i𝐫t/∂ t = 
D_i ∇^2 i𝐫t =0 into (<ref>), we obtain that the following relation holds for the equilibrium concentration <cit.>

lll 
C_A^eqC_B^eq = κ_b/κ_f≜κ_w.


For instance, assuming that type-A and type-B molecules represent acid and base and the fluid is water, we have κ_f=1.4× 10^11  1/(Ms) and κ_b=1.4× 10^-3  M/s at 25^∘C which leads to κ_w=10^-14 M^2 where M, called molar, is a unit of concentration and represents the number of moles per every liter of solution in a three-dimensional environment[The concentration C in molecule·m^-3 is related to the molar concentration C^mol in mol·L^-1 via C= 10^3 n_avC^mol where n_av=6.02× 10^23 is the Avogadro constant <cit.>.] <cit.>.  Therefore, when the concentrations of acid and base are identical, i.e., C_A^eq=C_B^eq=10^-7 M, the pH, defined by pH = -log_10(C_A^eq), is 7.  As can be seen from (<ref>), 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 (<ref>).


As stated earlier, the general reaction-diffusion equations in (<ref>) 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 κ_f A𝐫tB𝐫t. Note that even if we assume that one of the variables, e.g. C_B(𝐫,t), is fixed, it is still challenging to solve (<ref>a) in terms of C_A(𝐫,t).  Therefore, in the following, we derive a numerical method for solving (<ref>) in a computationally efficient manner. This is achieved by fully exploiting the characteristics of the considered system model. 



 §.§ Derivation of the CR


Let us assume that time is divided into a series of small intervals of length Δ 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 Δ t→ 0. In particular, from the reaction-diffusion equations in (<ref>), we have

lll 
i𝐫t+Δt =     
i𝐫t + ∫_t̃=t^t+Δt G_i(𝐫,t+Δt)dt̃ 

    +ΔC_i^df(𝐫,t+Δt)∫_t̃=t^t+Δt  D_i ∇^2 C_i(𝐫,t̃) dt̃  

   +  ΔC_i^rc(𝐫,t+Δt)∫_t̃=t^t+Δt   (κ_b- κ_f C_A(𝐫,t̃)C_B(𝐫,t̃) ) dt̃,   
  

where  Δ C_i^df(𝐫,t+Δ t) and Δ C_i^rc(𝐫,t+Δ t) are concentration changes from time t to time t+Δ 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 Δ C_i^x(𝐫,t+Δ t),  x∈{df,rc}, intractable. Nevertheless, the following lemma specifies Δ C_i^x(𝐫,t+Δ t) for Δ t→ 0. First, let us formally introduce the following two sets of assumptions which are needed for a rigorous statement of Lemma <ref>: A1) We assume C_i(𝐫,t)≠ 0 and ∇^2 C_i(𝐫,t)≠ 0 hold for ∀𝐫,t. A2) We assume 𝒯_A∩𝒯_B is empty and that 𝒯_i⊂𝒯 where  𝒯={t|t=mΔ t - ε,  m∈ℕ}, ℕ is the set of natural numbers, and ε is an arbitrary small positive real number satisfying ε≪Δ t. Note that assumption A1 is mild as it generally holds after the first release of molecules by the transmitter for any point 𝐫 in the MC environment. Similarly, assumption A2 is not limiting since Δ t can be chosen to be sufficiently small  to accommodate any desired release time pattern of molecules with high accuracy while satisfying A2.


 Under assumptions A1 and A2, the following asymptotic result holds for  Δ C_i^x(𝐫,t+Δ t):

lll  
Δt→0lim  ΔC_i^x(𝐫,t+Δt) = C_i^x(𝐫,t+Δt) - C_i(𝐫,t), x∈{df,rc},  

where C_i^df(𝐫,t) and C_i^rc(𝐫,t) are the type-i molecule concentration assuming that in interval [t,t+Δ t], only diffusion and only reactions occur, respectively,  and the other phenomenon is absent. 


The proof is given in Appendix <ref>.


The main result in the above lemma is that assuming Δ t→ 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  Δ t→ 0. The decoupling of reaction and diffusion equations in (<ref>) can be seen as a form of general splitting methods used in numerical mathematics <cit.>, <cit.>, <cit.>. Based on this result, the concentration update rule for i𝐫t+Δ t for the considered MC system is given in the following proposition.



Under assumptions A1 and A2, i𝐫t+Δ t is obtained as

lll 
	i𝐫t+Δt 

	= G̅_i(𝐫)  + C_i^df(𝐫,t+Δt) + C_i^rc(𝐫,t+Δt)- i𝐫t,


where G̅_i(𝐫)= ∑_t_i∈𝒯_iiδ_d(𝐫)δ_k(t+Δ t-ε-t_i) where δ_k(·) is the Kronecker delta function. Moreover, for the MC system under consideration, C_i^df(𝐫,t+Δ t) and C_i^rc(𝐫,t+Δ t) are given by (<ref>) on the top of the next page  
where c_1(𝐫,t)=C_A(𝐫,t)-C_B(𝐫,t), c_2(𝐫,t)=√(κ_f^2c_1^2(𝐫,t)+4κ_fκ_b),  c_3(𝐫,t)=C_A(𝐫,t)+C_B(𝐫,t), and c_4(𝐫,t)=c_2(𝐫,t)-κ_fc_3(𝐫,t)/c_2(𝐫,t)+κ_fc_3(𝐫,t). 



The proof is given in Appendix <ref>.



Algorithm <ref> summarizes the simulation steps for CR computation using Proposition <ref> where T^max is the maximum simulation time.






We note that as Δ t→ 0, one can also use the following first order approximations for C_i^df(𝐫,t+Δ t) and C_i^rc(𝐫,t+Δ t) respectively:

rll
C_i^df(𝐫,t+Δt)     = C_i(𝐫,t) +  D_i ∇^2 C_i(𝐫,t) Δt, 

C^rc_i(𝐫,t+Δt)    = C_i(𝐫,t) + (κ_b- κ_f C_A(𝐫,t)C_B(𝐫,t))Δt. 

 However, for these approximations to be valid, Δ t has to be chosen much smaller than when the exact expressions in Proposition <ref> are adopted. Moreover, small values of Δ t lead to stability issues for computation of ∇^2 C_i(𝐫,t), see Section <ref> for a detailed discussion.  




 §.§ Simplifications for Special Cases


In the following, we further simplify the diffusion term C_i^df(𝐫,t+Δ t) and reaction term C_i^rc(𝐫,t+Δ t) for the considered system model.

Diffusion: The computational complexity of Algorithm <ref> is largely determined by (<ref>a) since for each update, an n-dimensional integral has to be evaluated for each point 𝐫∈ℝ^n where ℝ is the set of real numbers.  Nevertheless, for the commonly adopted assumption of a point-source transmitter with impulsive release <cit.>,  the computation of C_i^df(𝐫,t)  can be significantly simplified using the following corollary.
 

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≜𝐫. In this case,  C_i^df(𝐫,t+Δ t) can be simplified as

rll 
C_i^df(𝐫,t+Δt)   
   =  1/(4πD_i Δt)^n/2 
∫_r̃=0^∞ C_i(r̃,t) 
W_i(r,r̃) dr̃,  

where W_i(r,r̃) is given by

lll 
W_i(r,r̃) =    2 
exp(-r̃^2+r^2/4D_i Δt) 

   ×
cosh(rr̃/2D_iΔt),    if   n=1

πr̃ I_0(rr̃/2D_i Δt),    if   n=2

4πD_ir̃ Δt/rsinh(rr̃/2D_iΔt),    if   n=3
.


Here, I_0(·) is the zeroth order modified Bessel function of the first kind, sinh(·) is the sine hyperbolic function, and   cosh(·) is the cosine hyperbolic function.


The proof is given in Appendix <ref>.


Note that the n-dimensional integral in (<ref>a) is simplified to a one-dimensional integral in (<ref>) which has to be evaluated for a one-dimensional space variable r∈ℝ^+, where ℝ^+ denotes the set of non-negative real numbers. In addition, the term W_i(r,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 (<ref>) 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 (<ref>b) and (<ref>c) for either κ_f→ 0 or κ_b → 0. In particular, (<ref>b) and (<ref>c) have indeterminate forms at κ_f=0. Moreover, when C_A(𝐫,t)=C_B(𝐫,t) holds, (<ref>b) and (<ref>c) have indeterminate forms at κ_b=0, too. The following corollary provides the results for these cases.


For the considered MC system, the molecule concentrations when only the forward reaction is present, denoted by C^fr_i(𝐫,t+Δ t), and when only the backward reaction is present, denoted by C^br_i(𝐫,t+Δ t), are given by

rll
C^fr_A(𝐫,t+Δt)     =κ_b→0lim    C^rc_A(𝐫,t+Δt) 

   = 
c_1(𝐫,t)/1-c_5(𝐫,t)exp(-κ_f c_1(𝐫,t)Δt),     if   c_1(𝐫,t)≠0 
 
1/κ_f t +1/C_A(𝐫,t),    if  c_1(𝐫,t) = 0
  

C^fr_B(𝐫,t+Δt)     = κ_b→0lim    C^rc_B(𝐫,t+Δt) 

   = 
c_1(𝐫,t)c_5(𝐫,t)exp(-κ_f c_1(𝐫,t)Δt)/1-c_5(𝐫,t)exp(-κ_f c_1(𝐫,t)Δt),     if   c_1(𝐫,t)≠0 
 
1/κ_f t +1/C_B(𝐫,t),    if  c_1(𝐫,t) = 0
      

C_i^br(𝐫,t+Δt)      = κ_f→0lim    C^rc_i(𝐫,t+Δt) 

   = C_i(𝐫,t) + κ_b Δt,  i∈{A,B}, 


where c_5(𝐫,t)=C_B(𝐫,t)/C_A(𝐫,t). 


Due to the indeterminate form of C^rc_i(𝐫,t+Δ t) at κ_b=0, we use L'Hopital's rule <cit.> to compute   C^fr_i(𝐫,t+Δ t)=κ_b→ 0lim   C^rc_i(𝐫,t+Δ t) where C^rc_i(𝐫,t+Δ t) is given by (<ref>b) and (<ref>c). This leads to  C^fr_i(𝐫,t+Δ t) given in (<ref>) and (<ref>) for the case when c_1(𝐫,t)≠ 0. Note that these expressions for C^fr_i(𝐫,t+Δ t) are indeterminate at c_1(𝐫,t)= 0. In a similar manner, we apply L'Hopital's rule to C^fr_i(𝐫,t+Δ t) given in (<ref>) and (<ref>) for c_1(𝐫,t)≠ 0 in order to find the corresponding expressions when c_1(𝐫,t)= 0. For the case, when κ_f=0, it is easier to solve the original PDEs in (<ref>) which simplifies to ∂ C_i^br(𝐫,t)/∂ t=κ_b and has the immediate solution given in (<ref>). This completes the proof.




 §.§ Discussion on Complexity


In the following, we compare the computational complexity of Algorithm <ref> 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 (<ref>) <cit.>. In particular, for FDM, time  and space are discretized into small intervals to approximate the differential operators in (<ref>). 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 PDE. 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 <cit.>. However, increasing N_t and N_r too much results in numerical instability for computation of the differential operators, e.g., ∇^2 C_i(𝐫,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 <cit.>. Moreover, the appropriate values of the step sizes are application-specific and there is no universal guideline that can be applied to all PDE and system parameters.  On the contrary, the proposed approach in Algorithm <ref> is a hybrid method where time is discretized; however, the problem is solved analytically in (<ref>) and (<ref>) with respect to the space variables. Since no differential operator is used in these update rules, Algorithm <ref> 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 (<ref>), 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 (<ref>).  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.


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 𝒪(N_tN_r^n) where 𝒪(·) is the standard big-O notation. Moreover, if (<ref>) is rewritten in polar and spherical coordinates when n=2 and 3, respectively, the computational complexity can be reduced to 𝒪(N_tN_r) by exploiting symmetry with respect to 𝐫. The complexity order of Algorithm <ref> is given by 𝒪(N_tN_r^2n). However, by exploiting the symmetry and using the update rule for C_i^df(𝐫,t) given in (<ref>), the complexity of Algorithm <ref>  is reduced to 𝒪(N_tN_r^2),  n=1,2,3. Furthermore, if one employs the first order approximations in (<ref>) as concentration update rules in Algorithm <ref>,  complexity is further reduced to 𝒪(N_tN_r). However, the implementation of the concentration update rules in (<ref>) suffers the same stability issues as the FDM methods. In fact, the main advantage of the proposed algorithm using the updates rules in (<ref>) and (<ref>) 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 (<ref>a) over 𝐫̃∈ℛ(𝐫), instead of over the entire space where ℛ(𝐫)={𝐫̃|𝐫̃-𝐫≤ R} and R is a parameter referred to as integration radius. Thereby, the complexity order of the proposed algorithm reduces to 𝒪(N_tN_r^n) and 𝒪(N_tN_r) for concentrations that are asymmetric and symmetric with respect to the origin, respectively, i.e., the same complexity order as for FDM methods. 

	








§ SIMULATION RESULTS


In this section, we first verify the accuracy of the proposed CR computation algorithm and the Poisson statistics assumed in (<ref>) using the particle-based simulator described in arXivAppendix F.<cit.>. 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 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  Δ t=1 μs for Algorithm <ref> and the particle-based simulator, respectively. Moreover, the simulation results shown in Figs. <ref> and <ref> are obtained by  running the particle-based simulator  10^4 times.





 §.§ Verification of the Proposed CR Computation Algorithm and the Poisson Statistics



In Fig. <ref>, we plot the expected number of molecules observed at the receiver, i.e., CR, where 𝒯_A={200} and 𝒯_B={300} μs and for initial condition C_A^init(𝐫)=C_B^init(𝐫)=0. We show numerical results obtained with Algorithm <ref>  and simulation results generated with the particle-based simulator introduced in Appendix <ref>. From Fig. <ref>, we observe a perfect agreement between numerical and simulation results which validates the accuracy of the proposed numerical algorithm. In Fig. <ref>, we also plot the equilibrium concentration which can be computed based on Remark <ref> as y̅_i^eq=C_i^eqV^rx=0.5236 where C_A^eq=C_B^eq=√(κ_b/κ_f) and V^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μs. Nevertheless, y̅_i(t) increases in interval t∈[0,200μs] due to the backward reaction (production) in (<ref>) and finally converges to y̅_i(t)=y̅_i^eq since the backward reaction and the forward reaction (degradation) in (<ref>) reach an equilibrium, cf. Remark <ref>. 









Next, we evaluate the accuracy of the Poisson statistical model for y_i(t) introduced in (<ref>), where the mean is obtained using Algorithm <ref>, 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. <ref>, we show the PMF of y_i(t) for a few different time instances, namely t=200,300 μs for type-A molecules and t=200,400 μs for type-B molecules. The same parameters as for Fig. <ref> are adopted for Fig. <ref>.  Note that the PMFs of y_A(t) and y_B(t) are identical at t=200 μ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 y̅_A(t=200μs)=y̅_B(t=200μs)≈ 0.5, a moderate value y̅_B(t=400μs)≈ 4, and a large value y̅_A(t=300μs)≈ 12, and is given by the CR plotted in Fig. <ref>. As can be seen from Fig. <ref>, the Poisson statistics with the mean obtained by Algorithm <ref> 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 <cit.>. 



 §.§ 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. <ref>, we plot the CR for type-A molecules, y̅_A(t),  versus time for 𝒯_A={200} and 𝒯_B={300} μs, and initial condition C_A^init(𝐫)=C_B^init(𝐫)=0. We consider a fixed D_A=10^-10 m^2·s^-1 and different D_B=[0.1, 0.5, 1, 2, 10]× 10^-10 m^2·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 y̅_A^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 κ_f on the CR. We adopt κ_b=0 and κ_f∈{0,0.01,0.1,1,10,∞}× 10^-17 molecule^-1·m^3·s^-1 where the extreme values of κ_f=0,∞ correspond to no reaction and instant reaction, respectively.  Moreover, we assume 𝒯_A={0} and 𝒯_B={100} μs and initial concentrations C_A^init(𝐫)=C_B^init(𝐫)=0. In Fig. <ref>, we show the corresponding CR of the type-A molecules, y̅_A(t),  versus time. As can be seen from this figure, for the non-reactive scenario, κ_f=0, the CR is very wide which causes significant ISI for the considered nominal value of symbol duration, i.e., T^symb=200 μs. On the other hand, as κ_f increases, the CR is shortened and correspondingly ISI is reduced. For example, for the considered default system parameters, i.e., κ_f=10^-17, T^symb=200 μs, and sampling time T^samp=100 μ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 y̅_A(t) does not suddenly drop to zero after the release of the type-B molecules when κ_f→∞ 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. <ref>. Note  that since we assume κ_b=0, the equilibrium level y̅_A^eq is zero.



Finally, we study the effect of backward reaction rate constant κ_b∈{0,1,10,100}× 10^25 molecule·m^-3
·s^-1 on the CR when κ_f=10^-17 molecule^-1·m^3·s^-1. Moreover, we assume 𝒯_A={200} and 𝒯_B={300} μs and the initial concentrations are zero. From Fig. <ref>, we observe that the main effect of increasing κ_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 <ref>. Interestingly, the channel dispersion does not significantly change  with the values of κ_b. In other words, κ_b impacts the level of noise molecules in the MC channel and does not considerably affect ISI.  


In summary, Figs. <ref>-<ref> 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).

















 §.§ 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 <cit.>, ii) MoSK modulation with 2TM receiver <cit.>, and iii) the proposed OOK modulation with 2TM receiver where τ_1=0[When 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 <cit.>. For a fair comparison, unlike <cit.> 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 τ_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 <ref>, namely conventional MoSK, the proposed OOK, and the proposed OSK.  Note that the optimal and suboptimal detectors derived in Section <ref> are applicable to all of the above binary modulation schemes in non-reactive and reactive MC systems. 


In Fig. <ref>, we show CR y̅_A(t) and y̅_B(t) versus time for data sequence 𝐬=[0,1,0,1,0]. We assume κ_f=10^-17 molecule^-1·m^3·s^-1 and κ_b=10^26 molecule·m^-3·s^-1 which leads to a equilibrium level of y̅_i^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 <ref>, we assume that the non-reactive system is affected by the same fixed level of noise molecules. As can be seen from Fig. <ref>,  for non-reactive signaling, 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 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. <ref> 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. <ref> shows the BER vs. the number of released molecules, N^tx_A = N^tx_B≜ N^tx, for the genie-aided (GA) lower bound from Proposition <ref> and Corollary <ref>, the ML-based detector with estimated ISI, and the suboptimal ISI-neglecting detector in (<ref>) and (<ref>). Fig. <ref> a) and b) show the results for 1TM and 2TM receivers, respectively. The following observations can be made in Fig. <ref>:



  * 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. <ref> 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 modulation[Although 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. <ref>.], 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. <ref>. Overall,  we observe from Fig. <ref> 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^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 (<ref>) 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,  A1≫B1 and A0≪B0 hold. OSK modulation also performs well for 1TM receivers since at the sampling time, A1≫A0 holds, see Fig. <ref>. 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 γ. 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. <ref>, we compare the performance of reactive and non-reactive MC systems. In particular, Fig. <ref> shows the BER vs. the number of released molecules, N^tx_A = N^tx_B≜ N^tx, for the ML-based detector. We summarize our observations from Fig. <ref> in the following:



  * From Fig. <ref>, 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^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^tx, and thus limits performance. 
In fact, at large N^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^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^tx and decreases for large N^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. <ref> and <ref>, we assume blocks of length K=10 and generate y_i(𝐬),  i∈{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. <ref> and <ref>, we study the impact of parameters L and K. In particular, in Fig. <ref>, we plot the BER vs. L for K=10 and N^tx_A = N^tx_B = 5× 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. <ref> 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. <ref> where the BER is plotted vs. block length K for L=5 and N^tx_A = N^tx_B = 5× 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.	











§ 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.





§ PROOF OF PROPOSITION <REF>
 

The LLR for  problem (<ref>)  can be written as

rll 
LLR=   log(f_𝒫(y_A|s=1,𝐬)f_𝒫(y_B|s=1,𝐬)/f_𝒫(y_A|s=0,𝐬)f_𝒫(y_B|s=0,𝐬)) 

=    log((y̅_A^(1)(𝐬))^y_A e^-y̅_A^(1)(𝐬) (y̅_B^(1)(𝐬))^y_B e^-y̅_B^(1)(𝐬)/(y̅_A^(0)(𝐬))^y_A e^-y̅_A^(0)(𝐬) (y̅_B^(0)(𝐬))^y_B e^-y̅_B^(0)(𝐬)) 

=    y_A log(y̅_A^(1)(𝐬)/y̅_A^(0)(𝐬))
-y_B log(y̅_B^(0)(𝐬)/y̅_B^(1)(𝐬)) 
    -y̅_A^(1)(𝐬)-y̅_B^(1)(𝐬)+y̅_A^(0)(𝐬)+y̅_B^(0)(𝐬).


Due to the monotonicity of the logarithm, the ML problem in (<ref>) can be rewritten as LLRs=1s=0⋛ 0.  Defining α(𝐬) and β(𝐬)  as in Proposition <ref>, we obtain (<ref>) which concludes the proof.



§ PROOF OF COROLLARY <REF>
 

	The proof is similar to that given in Appendix <ref> for 2TM receivers. In particular, the LLR for detection problem (<ref>) can be written as
	
	lll 
		LLR   =log(f_𝒫(y_A|s=1,𝐬)/f_𝒫(y_A|s=0,𝐬)) = log((y̅_A^(1)(𝐬))^y_A e^-y̅_A^(1)(𝐬) /(y̅_A^(0)(𝐬))^y_A e^-y̅_A^(0)(𝐬)) 

		   = y_A log(y̅_A^(1)(𝐬)/y̅_A^(0)(𝐬)) -y̅_A^(1)(𝐬)+y̅_A^(0)(𝐬).
	
	
	The ML problem in (8) can be rewritten as LLRs=1s=0⋛ 0 which in turn can be reformulated into y_A s=1s=0⋛γ(𝐬) where  γ(𝐬) is defined in Corollary 1. This completes the proof.



§ PROOF OF LEMMA <REF>



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

lll 
ΔC_i^x(𝐫,t+Δt) =  ∫_t̃=t^t+Δt f_i^x(𝐫,t̃) dt̃,   x∈{df,rc},


where f_i^df(𝐫,t̃)=D_i ∇^2 C_i(𝐫,t̃) and f_i^rc(𝐫,t̃)=κ_b- κ_f C_A(𝐫,t̃)C_B(𝐫,t̃). Let us define f̅_i^df(𝐫,t̃)=D_i ∇^2 C_i^df(𝐫,t̃) where C_i^df(𝐫,t̃)  denotes the concentration of type-i molecules assuming during interval t̃∈[t,t+Δ t] diffusion is present and reaction is absent. Similarly, let us define f̅_i^rc(𝐫,t̃)=κ_b- κ_f C_A^rc(𝐫,t̃)C_B^rc(𝐫,t̃) where C_i^rc(𝐫,t̃)  denotes the concentration of type-i molecules assuming during interval t̃∈[t,t+Δ 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+Δ t], i.e., (𝒯_A∪𝒯_B)∩[t,t+Δ t] is empty. Thereby, by defining ϵ_i^x(𝐫,t̃)≜ f_i^x(𝐫,t̃)-f̅_i^x(𝐫,t̃), we have 

lll 
ΔC_i^x(𝐫,t+Δt) 
= ≜A_i^x(Δt)∫_t̃=t^t+Δt f̅_i^x(𝐫,t̃)dt̃ 
+≜B_i^x(Δt)∫_t̃=t^t+Δt ϵ_i^x(𝐫,t̃) dt̃.


Next, we show that as Δ t → 0, B_i^x(Δ t) becomes negligible compared to A_i^x(Δ t). Under assumption A1, i.e., C_i(𝐫,t)≠ 0 and ∇^2 C_i(𝐫,t)≠ 0, we obtain f̅_i^x(𝐫,t)≠ 0 which leads to Δ t → 0lim A_i^x(Δ t)∼ o(Δ t). Recall that ϵ_i^x(𝐫,t̃) characterizes the error in decoupling the reaction and diffusion phenomena at time t̃>t. Therefore, for B_i^x(Δ t), we have ϵ_i^x(𝐫,t)=0 which leads to Δ t → 0lim B_i^x(Δ t)∼ o(Δ t^p) for some p>1. In summary, we have Δ t → 0limB_i^x(Δ t)/A_i^x(Δ t)∼ o(Δ t^p-1)→ 0.

Case 2: In this case, we assume the transmitter releases molecules within interval [t,t+Δ t], i.e., (𝒯_A∪𝒯_B)∩[t,t+Δ t] is non-empty. For this case, the conditions on set 𝒯_i introduced in assumption A2 are needed such that the argument in Case 1 is also applicable to Case 2.  In particular, under condition 𝒯_i⊂𝒯, without ambiguity, each release by the transmitter is taken into account only in one of the time slots.  Moreover, due to condition ε≪Δ t, the impact of diffusion and reaction on molecule concentration in interval [t+Δ t - ε,t+Δ t] is negligible compared to that in [t,t+Δ t - ε). Furthermore, assuming that 𝒯_A∩𝒯_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 Δ t → 0limB_i^x(Δ t)/A_i^x(Δ t)∼ o(Δ t^p-1)→ 0 provided for Case 1 to Case 2.  

In summary, under assumptions A1 and A2, we can conclude that as Δ t → 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. 



§ PROOF OF PROPOSITION <REF>


Eq. (<ref>) is obtained by applying (<ref>) from Lemma <ref> in (<ref>). Under assumption A1, we directly obtain  G̅_i(𝐫)≜∫_t̃=t^t+Δ t G_i(𝐫,t+Δ t)dt̃= ∑_t_i∈𝒯_iiδ_d(𝐫)δ_k(t+Δ t-ε-t_i). In the following, we derive expressions for C_i^df(𝐫,t+Δ t) and C_i^rc(𝐫,t+Δ 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^df(𝐫,t) and C_i^rc(𝐫,t). 

Diffusion: 
The following diffusion-only PDE  has to be solved for t>0 and arbitrary initial conditions  C_i(𝐫,0) at initial time t=0

lll 
∂C_i^df(𝐫,t)/∂t    = D_i ∇^2 C_i^df(𝐫,t).


For one-dimensional diffusion, this problem is solved in <cit.> by transforming it into the Laplace domain with respect to space variable 𝐫. The generalization of this result to the n∈{2,3} dimensional environment given in (<ref>a) is straightforward. Therefore, we skip the derivation due to space constraints and refer the readers to <cit.>.


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

lll 
∂C_A^rc(𝐫,t)/∂t    = -κ_f C_A^rc(𝐫,t)C_B^rc(𝐫,t) +κ_b 

∂C_B^rc(𝐫,t)/∂t    = -κ_f C_A^rc(𝐫,t)C_B^rc(𝐫,t)  +κ_b. 

Subtracting (<ref>b) from (<ref>a), we obtain

lll 
∂( C_A^rc(𝐫,t)-C_B^rc(𝐫,t) ) /∂t = 0,


which has the solution C_A^rc(𝐫,t)-C_B^rc(𝐫,t)=c_1(𝐫) where c_1(𝐫) is a constant with respect to t and is set to c_1(𝐫)=C_A^rc(𝐫,0)-C_B^rc(𝐫,0) to satisfy the initial conditions at t=0. Substituting C_B^rc(𝐫,t)=C_A^rc(𝐫,t)-c_1(𝐫) into (<ref>a) leads to

lll 
∂C_A^rc(𝐫,t)/∂t = - (κ_f(C_A^rc(𝐫,t))^2-κ_f c_1(𝐫)C_A^rc(𝐫,t)-κ_b).


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

lll 
∂C_A^rc(𝐫,t)/κ_f(C_A^rc(𝐫,t))^2-κ_f c_1(𝐫)C_A^rc(𝐫,t)-κ_b
= -∂t.


Integrating both sides of the above equation and using <cit.>

lll 
∫dx/ax^2+bx+c=1/Δlog(Δ-2ax-b/Δ+2ax+b),


where Δ=√(b^2-4ac), we obtain

lll 
- t + c̃_4(𝐫)=
1/c_2(𝐫)log(c_2(𝐫)-2κ_f C_A^rc(𝐫,t)+κ_fc_1(𝐫)/c_2(𝐫)+2κ_f C_A^rc(𝐫,t)-κ_fc_1(𝐫))


where c_2(𝐫)=√(κ_f^2c_1^2(𝐫)+4κ_fκ_b) and c̃_4(𝐫) is a constant with respect to time. Rewriting the above equation leads to

lll 
C_A^rc(𝐫,t) = 

c_2(𝐫)+κ_fc_1(𝐫)-(c_2(𝐫)-κ_fc_1(𝐫))
exp(-c_2(𝐫)(t-c̃_4(𝐫)))/2κ_f[1+exp(-c_2(𝐫)(t-c̃_4(𝐫)))].

Using the initial condition C_A^rc(𝐫,0), we obtain c̃_4(𝐫)=1/c_2(𝐫)ln(c_2(𝐫)+κ_fc_3(𝐫)/c_2(𝐫)-κ_fc_3(𝐫)) with c_3(𝐫)= C_A^rc(𝐫,0)+C_B^rc(𝐫,0). This leads to C_A^rc(𝐫,t) given in (<ref>b) after defining c_4(𝐫)≜exp(c_2(𝐫)c̃_4(𝐫))=c_2(𝐫)-κ_fc_3(𝐫)/c_2(𝐫)+κ_fc_3(𝐫) for notational simplicity. Using C_B^rc(𝐫,t)=C_A^rc(𝐫,t)-c_1(𝐫), we straightforwardly obtain C_B^rc(𝐫,t) given in (<ref>c).  

The above results are concisely given in Proposition <ref> which concludes the proof.








§ PROOF OF COROLLARY <REF>
 

Because of the geometrical symmetry of the problem, the concentration for the MC system under consideration  is only a function of  r≜𝐫. In the following, we present the proof for one-, two-, and three-dimensional environments, respectively. 

One-Dimensional Environment: In this case, we obtain (<ref>) by substituting n=1 in (<ref>a) as

lll 
C_i^df(𝐫,t+Δt)  

 =  1/(4πD_i Δt)^1/2 
∫_r̃_1
C_i(r̃_1,t) exp(-|r̃_1|^2+r^2-2rr̃_1/4D_i Δt)dr̃_1, 

 = 1/(4πD_i Δt)^1/2 
∫_r̃_1=-∞^0
C_i(r̃_1,t) exp(-|r̃_1|^2+r^2-2rr̃_1/4D_i Δt)dr̃_1, 

  + 1/(4πD_i Δt)^1/2  ∫_r̃_1=0^∞
C_i(r̃_1,t) exp(-|r̃_1|^2+r^2-2rr̃_1/4D_i Δt)dr̃_1, 

(a)= 1/(4πD_i Δt)^1/2 
∫_r̃=0^∞
C_i(r̃,t) 2exp(-r̃^2+r^2/4D_i Δt)
cosh(r̃r/2D_i Δt)dr̃, 

where for equality (a), we used r̃=|r̃_1|. Defining W_i(r,r̃) as in (<ref>) leads to  (<ref>) for n=1. 


Two-Dimensional Environment: For n=2, the simplification of C_i^df(𝐫,t+Δ t) requires transforming an integral from  Cartesian coordinates to polar coordinates, i.e., 

lll 
∫_r_1∫_r_2 f(r_1,r_2) dr_1dr_2 
= ∫_r∫_θ f(rcos(θ),rsin(θ))  r  dr dθ,       


where f(r_1,r_2) is an arbitrary function and r≥ 0 and θ∈[0,2π] are the variables of the polar coordinates. Moreover, without loss of generality, we consider 𝐫=(r,0) in order to simplify (<ref>a) as follows

lll 
C_i^df(𝐫,t+Δt) 

 =  1/4πD_i Δt 
∬_𝐫̃
C_i(𝐫̃,t) exp(-𝐫̃^2+r^2-2rr̃_1/4D_i Δt)d𝐫̃, 

 = 1/4πD_i Δt 
∫_r̃=0^∞ C_i(r̃,t) 
r̃exp(-r̃^2+r^2/4D_i Δt) 

         ×∫_θ̃=0^2πexp(rr̃cos(θ̃)/2D_i Δt) 
dθ̃ dr̃,       

 (a)= 1/4πD_i Δt
∫_r̃=0^∞ C_i(r̃,t) 
2πr̃
exp(-r̃^2+r^2/4D_i Δt)
I_0(rr̃/2D_i Δt)  dr̃,      


where in equality (a),  we used the identity ∫_x=0^2πexp(ccos(x))dx = 2π I_0(c) <cit.>. Defining W_i(r,r̃) as in (<ref>) leads to  (<ref>) for n=2. 

Three-Dimensional Environment: In a similar manner, the simplification of C_i^df(𝐫,t+Δ t) when  n=3 requires transforming an integral from  Cartesian coordinates to spherical coordinates, i.e., 

lll 
∫_r_1∫_r_2∫_r_3 f(r_1,r_2,r_3) dr_1dr_2dr_3 

= ∫_r∫_ϕ∫_θ f(rsin(ϕ)cos(θ),rsin(ϕ)sin(θ),rcos(ϕ)) 

        ×r^2 sin(ϕ) dr
dϕdθ,   

where f(r_1,r_2,r_3) is an arbitrary function and r≥ 0, ϕ∈[0,π], and θ∈[0,2π] are the variables of the spherical coordinates. Here, without loss of generality, we consider 𝐫=(0,0,r) in order to simplify (<ref>a) as follows

lll 
 C_i^df(𝐫,t+Δt)  

=  1/(4πD_i Δt)^3/2 
∭_𝐫̃
C_i(𝐫̃,t) exp(-𝐫̃^2+r^2-2rr̃_3/4D_i Δt)d𝐫̃, 

 = 1/(4πD_i Δt)^3/2 
∫_r̃=0^∞ C_i(r̃,t) 
r̃^2exp(-r̃^2+r^2/4D_i Δt) 

  ×∫_ϕ̃=0^π ∫_θ̃=0^2π
sin(ϕ̃)exp(rr̃cos(ϕ̃)/2D_i Δt) 
dθ̃ dϕ̃dr̃,       

 (a)= 1/(4πD_i Δt)^3/2 
∫_r̃=0^∞ C_i(r̃,t) 
8πD_i r̃Δt/r 

  ×exp(-r̃^2+r^2/4D_i Δt)
sinh(rr̃/2D_iΔt) dr̃, 


where for equality (a), we used the identity ∫_x=0^πsin(x)exp(ccos(x))dx = 2sinh(c)/c. Defining W_i(r,r̃) as in (<ref>) leads to  (<ref>) for n=3 and completes the proof.

arXiv
	
		
	

§ STOCHASTIC SIMULATION

	Stochastic simulation is used to verify the accuracy of Algorithm <ref> for CR computation and the Poisson statistics in (<ref>). We refer the readers to <cit.> 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 <cit.>. Therefore, it suffices  to track only the positions of the individual molecules  during the simulation time. Let 𝐫_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.
	
	

 §.§ Transmitter Release

	For instantaneous release from a point-source transmitter, we simply place N^tx_i type-i molecules at the origin, i.e., 𝐫_i^(m)(t) = (0,0,0),  m=M_i+1,…,M_i+N^tx_i, at any release time instant t∈𝒯_i. Variable M_i is updated by M_i+N^tx_i after each release of type-A molecules by the transmitter.
	
	

 §.§ Diffusion
 
	According to Brownian dynamics, the position of each molecule at time t+Δ t is updated as <cit.>
	
	lll 
		𝐫_i^(m)(t+Δt) = 𝐫_i^(m)(t) + √(2D_iΔt) (Δr_1,i^(m),Δr_2,i^(m),Δr_3,i^(m)),
	
	
	where Δ r_1,i^(m),Δ r_2,i^(m),Δ r_3,i^(m)∼𝒩(0,1) and 𝒩(μ,σ^2) denotes a normal RV with mean μ and variance σ^2.
	
	
	

 §.§ Forward Reaction

	
	The forward reaction in (<ref>), i.e., A+Bκ_f→∅, 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+Δ t] when the reactant molecules are within a certain binding radius ρ_b <cit.>.  Therefore, in our simulation, we can simply remove one type-A and one type-B molecule if their distance is less than ρ_b. In other words, we remove each pair of molecules whose position variables 𝐫_A^(m)(t) and 𝐫_B^(m')(t) satisfy 𝐫_A^(m)(t)-𝐫_B^(m')(t)<ρ_b. Unfortunately, the exact value of ρ_b cannot be found analytically in general  and depends on reaction rate constant κ_f and the choice of Δ t. Nevertheless, for the two special cases of Δ t→ 0 and Δ t→∞, the following simple relations are available <cit.>
	
	lll 
		κ_f = 
			4πρ_b(D_A+D_B),     if  ρ_rms≪ρ_b 

			4πρ_b^3/(3Δt),     if  ρ_rms ≫ρ_b
		
	
	
	where ρ_rms=√(2(D_A+D_B)Δ t) is the mutual root mean square step length of type-A and type-B molecules. Using the simplified formula for ρ_b in (<ref>), we obtain equivalent conditions for ρ_rms≪ρ_b and ρ_rms≫ρ_b as
	
	lll 
		
			ρ_rms≪ρ_b →Δt ≪T^crt

			ρ_rms ≫ρ_b →Δt ≫9/4 T^crt,
		
	
	
	where T^crt=κ_f^2/32π^2(D_A+D_B)^3. In order to be able to use the closed-form relation in (<ref>), we have to choose either large Δ t ≫9/4 T^crt or small Δ t ≪ T^crt. This choice depends on the system under consideration and the total length of the simulation interval.
	
	
	

 §.§ Backward Reaction

	The backward reaction in  (<ref>), i.e., ∅κ_b→ 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+Δ t]. Then, N(t) follows a Poisson distribution <cit.>
	
	lll 
		N(t) = 𝒫(Vκ_bΔt).
	
	
	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 ρ_rms≪ρ_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 ρ_b <cit.>.  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 ρ_u, which is larger than ρ_b .
	Let  𝐥 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 𝐫_i^(m)(t),  m=M_i+1,…,M_i+N(t),  obtained from <cit.>
	
	lll 
		𝐫_A^(m)(t)    = 𝐥 +  L(Δr_1,A^(m),Δr_2,A^(m),Δr_3,A^(m)) 

		𝐫_B^(m)(t)    = 𝐫_A^(m)(t) + ρ_u(Δr_1,B^(m),Δr_2,B^(m),Δr_3,B^(m) ),
	
	where Δ r_1,A^(m),Δ r_2,A^(m),Δ r_3,A^(m)∼𝒰(-0.5,0.5) and 𝒰(a,b) is an RV uniformly distributed in the interval [a,b]. Moreover, we have Δ r_1,B^(m)=sin(θ^(m))cos(ϕ^(m)),Δ r_2,B^(m)=sin(θ^(m))sin(ϕ^(m)), and Δ r_3,B^(m)=cos(θ^(m)) where θ^(m)∼𝒰(0,π) and ϕ^(m)∼𝒰(0,2π). On the other hand, if ρ_rms≫ρ_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 ρ_u=0 which leads to 𝐫_A^(m)(t) = 𝐫_B^(m)(t). Variable M_i is updated by M_i+N(t) at the end of this step. 
	
	Algorithm <ref> summarizes the main steps required for particle-based simulation of the considered MC system.
	
	
	
	
	




IEEEtran


arXiv
	


[
    < g r a p h i c s >
]Vahid Jamali (S'12) received the B.S. and M.S. degrees (honors) in electrical engineering from the
	K. N. Toosi University of Technology, Tehran,
	Iran, in 2010 and 2012, respectively, and
	the Ph.D. degree (with distinctions) from
	the Friedrich-Alexander-University (FAU) of
	Erlangen-Nürnberg, Erlangen, Germany,
	in 2019. In 2017, he was a Visiting Research Scholar with Stanford
	University, CA, USA. He is currently a Postdoctoral Fellow
	with the Institute for Digital Communication, FAU. His 
	research interests include wireless and molecular communications,
	Bayesian inference and learning, and multiuser information theory.
	
	Dr. Jamali has served as a member of the Technical Program
	Committee for several IEEE conferences. He received several
	awards, including the Exemplary Reviewer Certificates from IEEE Communications Letters
	in 2014 and the IEEE Transactions on Communications in 2017 and 2018, the Best Paper Award from
	the IEEE International Conference on Communications in 2016,
	the Doctoral Scholarship from the German Academic Exchange
	Service (DAAD) in 2017, the Goldener Igel Publication Award from
	the Telecommunications Laboratory (LNT), FAU, in 2018, and
	he was the Winner of the Best 3 Minutes Ph.D. Thesis (3MT)
	Presentation from the IEEE Wireless Communications and Networking Conference in 2018. He is currently an Associate Editor of IEEE Communications Letters.



[
    < g r a p h i c s >
]Nariman Farsad (S’07-M'15) is currently a Postdoctoral Fellow with the Department of Electrical Engineering at Stanford University, where he is a recipient of Natural Sciences and Engineering Research Council of Canada (NSERC) Postdoctoral Fellowship. His research interests are on emerging communication technologies such as molecular communication, and on improving the performance of communication systems through machine learning and deep learning. 
	
Dr. Farsad has won the second prize in 2014 IEEE ComSoc Student Competition: Communications Technology Changing the World, the best demo award at INFOCOM’2015, and was recognized as a finalist for the 2014 Bell Labs Prize. He has been an Area Associate Editor for IEEE Journal of Selected Areas of Communication–Special Issue on Emerging Technologies in Communications, and a Technical Reviewer for a number of journals including IEEE Transactions on Signal Processing, and IEEE Transactions on Information Theory. He was also a member of the Technical Program Committees for the ICC’2015, ICC’2018, BICT’2015, GLOBCOM’2015, GLOBCOM’2016, and GLOBECOM’2017. 


[
    < g r a p h i c s >
]Robert Schober (S'98, M'01, SM'08, F'10) was born in Neuendettelsau, Germany, in 1971. He received the Diplom (Univ.) and the Ph.D. degrees in electrical engineering from the Friedrich-Alexander-University of Erlangen-Nuremberg (FAU), Germany, in 1997 and 2000, respectively. From May 2001 to April 2002 he was a Postdoctoral Fellow at the University of Toronto, Canada, sponsored by the German Academic Exchange Service (DAAD). From 2002 to 2011, he was a Professor and Canada Research Chair at the University of British Columbia (UBC), Vancouver, Canada. Since January 2012 he is an Alexander von Humboldt Professor and the Chair for Digital Communication at FAU. His research interests fall into the broad areas of Communication Theory, Wireless Communications, and Statistical Signal Processing.
	
	Dr. Schober received several awards for his work including the 2002 Heinz Maier–Leibnitz Award of the German Science Foundation (DFG), the 2004 Innovations Award of the Vodafone Foundation for Research in Mobile Communications, the 2006 UBC Killam Research Prize, the 2007 Wilhelm Friedrich Bessel Research Award of the Alexander von Humboldt Foundation, the 2008 Charles McDowell Award for Excellence in Research from UBC, a 2011 Alexander von Humboldt Professorship, a 2012 NSERC E.W.R. Stacie Fellowship, and a 2017 Wireless Communications Recognition Award. Furthermore, he has been listed as a Highly Cited Researcher by Clarivate Analytics. Robert is a Fellow of the Canadian Academy of Engineering and a Fellow of the Engineering Institute of Canada. From 2012 to 2015, he served as Editor-in-Chief of the IEEE Transactions on Communications. Currently, he is the Chair of the Steering Committee of the IEEE Transactions on Molecular, Biological and Multiscale Communication and serves on the Editorial Board of the Proceedings of the IEEE. Furthermore, he is a Member at Large of the Board of Governors and a Distinguished Lecturer of the IEEE Communications Society.


[
    < g r a p h i c s >
]Andrea Goldsmith (S'90-M'93-SM'99-F'05) received
	the B.S., M.S., and Ph.D. degrees in electrical
	engineering from the University of California at
	Berkeley. She is the Stephen Harris Professor in the School
	of Engineering and a Professor of Electrical Engineering
	at Stanford University. She was previously
	on the faculty of Electrical Engineering at Caltech.
	Her research interests are in information theory
	and communication theory, and their application to
	wireless communications and related fields. She cofounded
	and served as Chief Scientist of Plume WiFi, and previously cofounded
	and served as CTO of Quantenna Communications, Inc. She has also
	held industry positions at Maxim Technologies, Memorylink Corporation, and
	AT&T Bell Laboratories. She is author of the book Wireless Communications
	and co-author of the books MIMO Wireless Communications and Principles of
	Cognitive Radio, all published by Cambridge University Press, as well as an
	inventor on 28 patents.
	
	Dr. Goldsmith has received several awards for her work, including the
	IEEE ComSoc Edwin H. Armstrong Achievement Award, as well as Technical
	Achievement Awards in Communications Theory and in Wireless Communications,
	the National Academy of Engineering Gilbreth Lecture Award,
	the IEEE ComSoc and Information Theory Society Joint Paper Award, the
	IEEE ComSoc Best Tutorial Paper Award, the Alfred P. Sloan Fellowship,
	the WICE Technical Achievement Award, and the Silicon Valley/San Jose
	Business Journals Women of Influence Award. She has served on the Steering
	Committee for the IEEE Transactions on Wireless Communications
	and as editor for the IEEE Transactions on Information Theory,
	the Journal on Foundations and Trends in Communications and Information
	Theory and in Networks, the IEEE Transactions on Communications,
	and the IEEE Wireless Communications Magazine. She participates actively
	in committees and conference organization for the IEEE Information Theory
	and Communications Societies and has served on the Board of Governors for
	both societies. She has also been a Distinguished Lecturer for both societies,
	served as President of the IEEE Information Theory Society in 2009, founded
	and chaired the student committee of the IEEE Information Theory society,
	and chaired the Emerging Technology Committee of the IEEE Communications
	Society. At Stanford she received the inaugural University Postdoc Mentoring
	Award, served as Chair of Stanford’s Faculty Senate and for multiple terms as a
	Senator, and currently serves on its Budget Group, Committee on Research, and
	Task Force on Women and Leadership. She is a Fellow of Stanford University.



