Bayram Malcolm Egan

Cevdet Akdeniz

Bao Quoc Tang

Stochastic Reaction and Diffusion Systems in Molecular Communications: Recent Results and Open Problems

come   L'archive ouverte pluridisciplinaire

Introduction

A basic feature of a wide range of biochemical systems is the presence of chemical signals, which can be viewed as facilitating information exchange, or communication. A prototypical example is quorum sensing [START_REF] Miller | Quorum sensing in bacteria[END_REF], where bacteria cells coordinate their behavior based on local estimates of the quantity of proteins emitted by each cell into a common fluid environment. Other examples include inter-cell communication via calcium signaling [START_REF] Barros | Engineering calcum signaling of astrocytes for neural-molecular computing logic gates[END_REF] or in the development of the fruit fly Drosophila melanogaster via protein gradients [START_REF] Affolter | The Decapentaplegic morphogen gradient: from pattern formation to growth regulation[END_REF].

In 2005, the notion of molecular communications for synthetic or artificial devices was introduced in [4]. In contrast to chemical signals in natural biological systems, artificial molecular communication was developed in order to support communication between nanoscale devices, which may act as sensors or actuators. The potential applications of this idea in areas from healthcare to the environment has motivated a large body of literature over the intervening 15 years, including the monograph [START_REF] Nakano | Molecular Communication[END_REF].

It was quickly realized that a key challenge in molecular communications is uncertainty or randomness in the motion of the molecules comprising the chemical signal. Often, the motion of information-carrying molecules is driven by diffusion or convection within a fluid, and these molecules are larger than the other molecules comprising the fluid; that is, they are colloidal. It is well known from statistical physics that the motion of colloidal molecules depends on many factors including the composition of the fluid, the possibility of chemical reactions, the presence of external forces, molecular degradation, the presence of obstacles, and the structure of the fluid boundary [START_REF] Coffey | The Langevin Equation[END_REF].

One of the simplest stochastic models for molecular motion is the Wiener process, with non-interacting information carrying molecules, which lies in the general family of Smoluchowski equations that in turn lie in the family of Langevin diffusions [START_REF] Coffey | The Langevin Equation[END_REF][START_REF] Pavliotis | Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations[END_REF]. One advantage of this approach is that, under certain conditions, one can obtain explicit solutions, which are desirable for the design of modulation and detection schemes. Indeed, using the Wiener process model, there has been extensive work developing novel modulation and detection schemes, as well as parameter estimation and channel coding (e.g., [START_REF] Kuscu | Transmitter and receiver architectures for molecular communications: A survey on physical design with modulation, coding, and detection techniques[END_REF][START_REF] Kuran | A survey on modulation techniques in molecular communication via diffusion[END_REF][START_REF] Söldner | A survey of biological building blocks for synthetic molecular communication systems[END_REF][START_REF] Bi | A survey of molecular communication in cell biology: Establishing a new hierarchy for interdisciplinary applications[END_REF]). Techniques from information theory have also provided insights into the fundamental performance limits of these schemes (e.g., [START_REF] Srinivas | Molecular communication in fluid media: the additive inverse Gaussian noise channel[END_REF][START_REF] Akyildiz | An information theoretic framework to analyze communication systems based on statistical mechanics[END_REF][START_REF] Rose | Capacity bounds on point-to-point communication using molecules[END_REF]). Various boundary conditions have also been investigated in order to model reflection or absorption, which can capture the effect of certain kinds of chemical reactions (e.g., [START_REF] Deng | Modeling and simulation of molecular communication systems with a reversible adsorption receiver[END_REF]).

Nevertheless, the Wiener process does not capture all environments relevant for molecular communications, for example, the presence of external forces or general forms of chemical reactions between molecules in the system. For such cases, more general stochastic reaction-diffusion systems provide a powerful framework to develop implementations of fundamental signal processing algorithms and model biochemical systems.

Naturally, the additional complexity of these models makes it almost impossible to have any closed form solutions. Instead, recent mathematical advances in the characterization of reaction-diffusion systems (e.g., [START_REF] Fellner | The entropy method for reaction-diffusion systems without detailed balance: first order chemical reaction networks[END_REF]) play a key role in obtaining accurate and tractable approximations.

Main Contributions

This paper reviews recent contributions addressing these issues by introducing sophisticated stochastic models for the motion of information-carrying molecules and also exploiting these models for computation.

We start by reviewing two key families of stochastic models: Langevin diffusion and the reaction-diffusion master equation (RDME). The former model, in the form of the Smoluchowski equation, can be used to describe motion of a single molecule or multiple non-interacting molecules. In particular, for the Smoluchowski equation, when the external force is induced by a bounded scalar potential, a closed form expression is available for the equilibrium solution, which makes it straightforward to obtain the asymptotic statistics of the system. We also briefly describe the Euler-Maruyama method [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] to obtain realizations of sample paths (i.e., the realization of the trajectory for a single molecule) for the Smoluchowski equation. However, in general situations, where chemical reactions are either dominant or non-negligible, the Langevin diffusion model is not an accurate description for the dynamics of the system.

A first alternative is the chemical master equation (CME), a stochastic description of chemical reactions, whose (large time) behavior has been developed using tools from chemical reaction network theory [START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactors-I. the deficiency zero and deficiency one theorems[END_REF]. More precisely, when the chemical reaction network satisfies condition known as complex balance, then a stationary distribution can be determined explicitly using the complex balanced equilibrium. To capture both diffusion and chemical reactions, a second alternative called the reaction-diffusion master equation should be considered. The complexity of this model makes it more challenging to apply it to molecular communication. Nevertheless, there is a remarkable-but not, in general, completely rigorous-connection between this model and deterministic reaction-diffusion systems in the form of limit theorems [START_REF] Kurtz | Limit theorems for sequences of jump Markov processes[END_REF][START_REF] Arnold | Deterministic limit of the stochastic model of chemical reactions with diffusion[END_REF]. The qualitative and quantitative behavior of such deterministic reaction-diffusion systems have recently witnessed significant advances in the mathematics literature [START_REF] Fellner | The entropy method for reaction-diffusion systems without detailed balance: first order chemical reaction networks[END_REF]. This opens up the opportunity to utilize these recent results in molecular communication, and, at the same time, requires new rigorous, analytical studies of the relationship between deterministic reaction-diffusion systems and their stochastic counterpart. After laying out the underlying theoretical models, which are summarized in Fig. 1 we then overview their application in three contexts relevant for signal processing in molecular communication:

(i) Modulation and Detection: In order to design reliable modulation and detection algorithms, it is necessary to characterize the statistics for the quantity of molecules observed by a receiver. In this case, stochastic models provide the basis for analysis of these statistics.

(ii) Statistical Inference: For low complexity synthetic devices, molecular circuits-where computation is performed via chemical reactions-are a highly desirable means of implementing algorithms for statistical inference (e.g., the detection schemes in (i) or parameter estimation). Stochastic models, particularly stochastic chemical reaction systems, form a useful basis for designing these molecular circuits.

(iii) Modeling Biochemical Environments: In realistic molecular communication systems, informationcarrying molecules may interact with other biochemical processes in the environment. These processes may themselves be governed by stochastic models, and therefore it is a key issue to understand how the dynamics of the information-carrying molecules interact with the dynamics of such external processes.

We then discuss selected open problems when the channel, inference or the environment is modeled via stochastic reaction-diffusion systems. We highlight both engineering challenges and mathematical problems arising from system design.

Related Surveys

In recent years there have been a number of surveys on various aspects of molecular communication, ranging from channel models to applications in biology, medicine, and nanonetworking. An incomplete list of these surveys is summarized in Table 1.

Topic References

General Survey [START_REF] Farsad | A comprehensive survey of recent advancements in molecular communication[END_REF] Applications in Medicine [START_REF] Chude-Okonkwo | Molecular communication and nanonetwork for targeted drug delivery: a survey[END_REF][START_REF] Chahibi | Molecular communication for drug delivery systems: a survey[END_REF][START_REF] Felicetti | Applications of molecular communications to medicine: a survey[END_REF] Applications in Biology [START_REF] Bi | A survey of molecular communication in cell biology: Establishing a new hierarchy for interdisciplinary applications[END_REF][START_REF] Söldner | A survey of biological building blocks for synthetic molecular communication systems[END_REF][START_REF] Mcbride | The effect of loads in molecular communications[END_REF][START_REF] Kim | Redox is a global biodevice information processing modality[END_REF] Applications in Nanonetworking [START_REF] Nakano | Methods and applications of mobile molecular communication[END_REF][START_REF] Darchini | Molecular communication via microtubules and physical contact in nanonetworks: a survey[END_REF][START_REF] Yang | A comprehensive survey on hybrid communication in context of molecular communication and terahertz communication for body-centric nanonetworks[END_REF] Channel Models [START_REF] Jamali | Diffusive molecular communications with reactive molecules: Channel modeling and signal design[END_REF] Information Theoretic Limits [START_REF] Rose | Capacity bounds on point-to-point communication using molecules[END_REF][START_REF] Akyildiz | An information theoretic framework to analyze communication systems based on statistical mechanics[END_REF] Modulation and Equalization [START_REF] Kuscu | Transmitter and receiver architectures for molecular communications: A survey on physical design with modulation, coding, and detection techniques[END_REF][START_REF] Kuran | A survey on modulation techniques in molecular communication via diffusion[END_REF][START_REF] Gursoy | Towards high data-rate diffusive molecular communications: Performance enhancement strategies[END_REF] Synaptic Communication [START_REF] Veletić | Synaptic communication engineering for future cognitive brain-machine interfaces[END_REF] The vast majority of the existing surveys have focused on applications. A key novel aspect of this survey is its focus on methodology. In particular, in the first part of this survey, we overview a general family of stochastic models incorporating diffusion, external forces, and chemical reactions. This also contrasts with the surveys in [START_REF] Jamali | Diffusive molecular communications with reactive molecules: Channel modeling and signal design[END_REF][START_REF] Akyildiz | An information theoretic framework to analyze communication systems based on statistical mechanics[END_REF], which either emphasize deterministic channel models or a restricted family of Langevin diffusions. The other distinguishing aspect of this survey is in the second part, which overviews applications of stochastic models of chemical reactions and diffusion to inference via biological circuits and lab-on-a-chip devices, as well as the impact of external biochemical processes.

Organization

The paper is organized as follows. In Sec. 2, we overview stochastic models based on Langevin diffusion, while Sec. 3 is devoted to chemical master equations and reaction-diffusion master equations. In Sec. 4, we overview modulation and detection schemes applicable in general Smoluchowski and RDME models. In Sec. 5, we highlight methods for implementing statistical inference algorithms via stochastic models, which forms a basis for the construction of biological circuits. In Sec. 6, we overview recent work modeling the interaction between molecular communication systems and external biochemical processes. Sec. 7 outlines open problems and Sec. 8 concludes.

Notation

The notation in this paper is summarized in Table 2.

Langevin Diffusion

A basic source of randomness in the motion of colloidal molecules in a fluid is diffusion. Before turning to more complex models relevant for molecular communications, we first review the dynamics of a single molecule in a domain Ω ⊆ R 3 . The basic theory of diffusion for a single colloidal molecule was first developed by Einstein [START_REF] Einstein | On the movement of small particles suspended in stationary liquids required by the molecular-kinetic theory of heat[END_REF] in order to characterize Brownian motion and developed further by Smoluchowski [START_REF] Smoluchowski | Zur kinetischen theorie der Brownschen molekularbewegung und der suspensionen[END_REF] and Langevin [START_REF] Langevin | On the theory of Brownian motion[END_REF].

Dynamics of a Single Molecule

Let d ∈ {1, 2, 3} be the dimension of the system. The general form of the motion of a single molecule can be described by the differential equation, obtained via Newton's second law,

m d 2 x dt 2 = -γ(x) dx dt + F(x) + σ(x)ξ(t), (1)
where m is the mass of the molecule, x(t) ∈ R d is the particle's location at time t, γ(x) dx dt ∈ R d corresponds to viscous drag, F(x) ∈ R d is an external force (also known as convection), and σ(x)ξ(t) ∈ R d is a random perturbation due to collisions of the colloidal particle with the molecules comprising the fluid with σ(x) ∈ R d×d . Due to the randomness, (1) forms a stochastic differential equation, with solution depending on the statistical structure of the random perturbation.

The term F(x) corresponds to external forces. For example, in [START_REF] Lo | Modeling molecular communications in tubes with Poiseuille flow and Robin boundary condition[END_REF], flow is induced by a pressure gradient, known as Poiseuille flow. A common scenario, arising when electrostatic or gravitational fields act on the molecule, is that F(x) = ∇V (x), where V is known as a scalar potential. A widely used example in molecular communications [START_REF] Akyildiz | An information theoretic framework to analyze communication systems based on statistical mechanics[END_REF] is

V (x) = v 0 + a • x, v 0 ∈ R, a ∈ R 3
, where the force acting on the particle in each dimension is given by a i , i = 1, 2, 3, commonly known as drift. Another common choice is the quadratic potential, which in one dimension is given by V (x) = ax 2 , a ∈ R considered in [START_REF] Egan | Equilibrium signaling in spatially inhomogeneous diffusion and external forces[END_REF], arising when electrostatic steering is present [START_REF] Wade | Electrostatic steering and ionic tethering in enzyme-ligand binding: insights from simulations[END_REF]. A further example occurs in the presence of acoustic tweezers [START_REF] Takatori | Acoustic trapping of active matter[END_REF].

The random perturbation ξ(t), is assumed to satisfy the following basic properties [START_REF] Coffey | The Langevin Equation[END_REF]:

(i) ξ(t) is independent of the location x;

(ii) E[ξ(t)] = 0; (iii) E[ξ(t)ξ(t ) T ] = Iδ(t -t ),
where I is the identity matrix and δ(•) is the Dirac delta function. On the other hand, σ(x) ∈ R d×d is a deterministic matrix, which can introduce location dependent variations in correlation and variance of the random fluctuations. Often, it is assumed that σ(x) = σ(x)I, which means that the spatial components of the random perturbation are uncorrelated, potentially with location-dependent variance, often known as spatially inhomogeneous diffusion. In the vast majority of work on molecular communications and classical statistical mechanics, the term σ(x) = σ is independent of the location; however, this assumption can fail to hold in inhomogeneous fluid media [START_REF] Hughes | An introduction to inhomogeneous liquids, density functional theory, and the wetting transition[END_REF].

Making sense of the differential equation in (1) requires a non-standard calculus, often in the form of the Itô calculus [START_REF] Pavliotis | Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations[END_REF]. Roughly speaking, the starting point is to rewrite (1) as a stochastic integral equation. A rigorous definition of the random perturbation ξ(t) can be obtained via the Wiener process W(t) as [6]

W(t) = t 0 ξ(t )dt . (2)
In particular, the Wiener process is a continuous-time stochastic process, which can be viewed as a random function, satisfying

(i) W(0) = 0; (ii) W(t 2 ) -W(t 1 ) is independent of W(t 4 ) -W(t 3 ), with [t 1 , t 2 ] ∩ [t 3 , t 4 ] = ∅; (iii) W(t 2 ) -W(t 1
) is multivariate Gaussian with covariance matrix (t 2 -t 1 )I;

(iv) W(t) has continuous sample paths (i.e., realizations of (W(t), 0 ≤ T ) are continuous) almost surely.

For further properties of the Wiener process and its role in Langevin diffusion, we refer the reader to [START_REF] Coffey | The Langevin Equation[END_REF][START_REF] Pavliotis | Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations[END_REF].

An important special case of (1) arises when drag forces dominate inertial forces, which yields

γ(x) dx dt = F(x) + σ(x)ξ(t), (3)
where σ(x) = σ(x)I, known as the Smoluchowski equation. This situation often arises in molecular communications due to the fact that inertial forces are due to the mass of the particles, which are typically small.

∂ ∂t p(x, t|x 0 , t 0 ) = ∇ 2 (D(x)p(x, t|x 0 , t 0 )) -∇ • F(x) γ(x) p(x, t|x 0 , t 0 ) , (4)
where p(x, t|x 0 , t 0 ) is the probability density function for the location of the molecule at time t, given an initial location x 0 at time t 0 , and D(x) = σ(x) 2 2γ(x) 2 is the diffusion coefficient. Note that both the Smoluchowski equation in (3) and the Fokker-Planck equation in (4) are equivalent descriptions of the system dynamics.

This formulation allows for the diffusion coefficient to be spatially inhomogeneous; i.e., the diffusion coefficient depends on the location x ∈ R d .

In the special case where there are no external forces and the diffusion coefficient is constant, the Fokker-Planck equation for Brownian motion (commonly called Fick's second law) is obtained as [START_REF] Risken | The Fokker-Planck Equation: Methods of Solution and Applications[END_REF] ∂ ∂t p(x, t|x 0 , t 0 ) = D∇ 2 p(x, t|x 0 , t 0 ),

where D = σ 2 2γ 2 is the spatially homogeneous diffusion coefficient, independent of the location x, corresponding to constant coefficients σ, γ.

Solutions to the Fokker-Planck equations in (4) and ( 5) depend on boundary conditions. In the case of an infinite domain, conditions on the decay p(x, t|x 0 , t 0 ) → 0 as x → ∞ are imposed; for more details, see [START_REF] Schulten | Lectures in Theoretical Biophysics[END_REF]. In the case of finite domains (corresponding to finite volumes), which are common in the context of molecular communications, standard boundary conditions are reflective, reaction, Robin, or periodic.

Each of the boundary conditions represents a different physical behavior of molecules at the boundary.

We illustrate in Figure 2

) := n(x) • ∇(D(x)p(x, t|x 0 , t 0 )) - F(x) γ(x) p(x, t|x 0 , t 0 ) = 0, x ∈ ∂Ω, (6)
where n(x) is the unit outward normal vector on ∂Ω. The term j(x, t|x 0 , t 0 ) is known as the flux, and therefore ( 6) is interpreted as there is zero flux at the boundary. Another boundary condition is known as the homogeneous Dirichlet boundary condition, where

p(x, t|x 0 , t 0 ) = 0, x ∈ ∂Ω, (7)
which corresponds to the scenario where every particle is consumed at the boundary ∂Ω. It is also possible to combine the boundary conditions to yield the radiation boundary condition (or Robin boundary condition)

n(x) • j(x, t|x 0 , t 0 ) = kp(x, t|x 0 , t 0 ), x ∈ ∂Ω, (8)
where k ∈ [0, ∞). The fourth boundary condition that is often imposed is known as the periodic boundary condition, in case Ω is a bounded interval in R or a high dimensional torus T d . In one dimensional system defined on Ω = [0, b], it reads as

p(0, t|x 0 , t 0 ) = p(b, t|x 0 , t 0 ) ∀t > 0. (9)
In general, there do not exist closed-form expressions for p(x, t|x 0 , t 0 ) solving (4) under general boundary conditions. In special cases, such as for ( 5) subject to various boundary conditions, results can be found in [START_REF] Risken | The Fokker-Planck Equation: Methods of Solution and Applications[END_REF][START_REF] Ahmadzadeh | Comprehensive reactive receiver modeling for diffusive molecular communication systems: Reversible binding, molecule degradation, and finite number of receptors[END_REF][START_REF] Deng | Modeling and simulation of molecular communication systems with a reversible adsorption receiver[END_REF][START_REF] Deng | Analyzing large-scale multiuser molecular communication via 3-d stochastic geometry[END_REF].

Multiple Non-Interacting Molecules

Recall the Smoluchowski model in [START_REF] Affolter | The Decapentaplegic morphogen gradient: from pattern formation to growth regulation[END_REF]. A basic question is what are the statistics for how many molecules are present in a given compact subdomain Ω of Ω ⊂ R d at a given time? As will be discussed further in Section 4, this question is critical for receiver design in molecular communication. Suppose that all molecules are non-interacting; that is, their positions evolve independently. Physically, this means that the colloidal molecules do not collide with each other, or exert forces (e.g., electrostatic) on each other. Let X i (t) ∈ {0, 1}, i = 1, . . . , W indicate whether or not molecule i lies in the subdomain Ω at time t.

The total number of molecules in Ω at time t is then given by

X tot (t) = W i=1 X i (t), ( 10
)
where W is the total number of molecules in Ω. Since the X i (t) are independent by virtue of the absence of interactions, for W sufficiently large, the distribution of X tot (t) may be approximated via the central limit theorem [START_REF] Billingsley | Probability and Measure[END_REF]. In particular, a useful approximation for the distribution of X tot (t) is given by

Xtot (t) ∼ N (W p Ω , W p Ω ), (11)
where N (µ, σ 2 ) corresponds to the Gaussian distribution, with mean µ and variance σ 2 . The term p Ω denotes the probability a given colloidal molecule lies in Ω at time t. In general computing p Ω requires the solution of the Smoluchowski equation in [START_REF] Affolter | The Decapentaplegic morphogen gradient: from pattern formation to growth regulation[END_REF] or the Fokker-Planck equation in (4), which do not have a closed-form expressions for general functional forms of the external forces and boundary conditions.

For t sufficiently large, the Smoluchowski equation can admit simple behavior. This arises when an equilibrium solution exists, which means that p Ω converges to a quantity independent of the time t. Although the existence of an equilibrium solution must be postulated, as in the work of Einstein and Langevin, a sufficient condition in one-dimensional systems is the fluctuation-dissipation condition [START_REF] Schulten | Lectures in Theoretical Biophysics[END_REF][START_REF] Grassia | Dissipation, fluctuations, and conservation laws[END_REF], which is given by

∂ ∂x D(x) = F (x)(γ(x) -1 -D(x)β), x ∈ Ω, (12)
where β = 1 kBT , where k B is Boltzmann's constant and T is the temperature of the system. For two and three-dimensional systems, see [START_REF] Schulten | Lectures in Theoretical Biophysics[END_REF]. In the special case that the diffusion coefficient D(x) = D, x ∈ Ω, then the fluctuation-dissipation condition corresponds to the Einstein relation; namely,

D = 1 γβ . ( 13
)
Under the fluctuation-dissipation condition, the Fokker-Planck equation in (4) can be written as [START_REF] Schulten | Lectures in Theoretical Biophysics[END_REF]

∂ ∂t p(x, t|x 0 , t 0 ) = ∂ ∂x D(x) ∂ ∂x p(x, t|x 0 , t 0 ) - ∂ ∂x βF (x)D(x)p(x, t|x 0 , t 0 ). (14)
For a detailed derivation, we refer the reader to [42, Eq. 4.17].

If the external force is determined by a scalar potential function U , the force satisfies F = -∂U ∂x and U is bounded in [b l , b r ] and infinite elsewhere, then the equilibrium solution to the Fokker-Planck equation in (4) is given by [41, Sec. 5.2]

p ∞ (x) = Ze -βU (x) , x ∈ R (15)
where

Z = br b l e -βU (x) dx -1
is the normalizing constant. For d ∈ {2, 3} and a potential U (x), x ∈ R d , a similar result holds; namely,

p ∞ (x) = Ze -βU (x) , x ∈ R d , ( 16
)
where Z is the appropriate normalizing constant.

Note that scalar potentials are ubiquitous; for example, arising from electrostatic and gravitational fields.

We also highlight that the equilibrium solution p ∞ in ( 15) is remarkable in the sense that it only depends on β and the potential U . That is, given these parameters, it is straightforward to obtain the asymptotic statistics of the system. This is true even if the diffusion coefficient is spatially inhomogeneous. We note that for higher-dimensional systems, a similar equilibrium solution to (15) also holds [START_REF] Risken | The Fokker-Planck Equation: Methods of Solution and Applications[END_REF].

Simulation

In order to obtain realizations of the sample paths for the Smoluchowski equation in (3), a standard approach is known as the Euler-Maruyama method [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. Let a(t, x) ≡

F(x) γ(x) , b(t, x) ≡ σ(x) γ(x)
, and ∆W n be a realization drawn from N (0, ∆tI). The Euler-Maruyama method is detailed in Algorithm 1, which is analogous to the Euler method for deterministic partial differential equations. Boundary conditions can be treated via the methods in [START_REF] Lépingle | Euler scheme for reflected stochastic differential equations[END_REF].

Algorithm 1 The Euler-Maruyama Method 1: INPUT: a(t, x), b(t, x), number of iterations N .

2: Initialize: x(0)

3: For n = 1 : N 4: t n = t n-1 + ∆t 5: x n+1 = x n + a(t n , x n )∆t + b(t n , x n )∆W n

Stochastic Reaction-Diffusion Models

In the previous section, we overviewed Langevin models for diffusion and convection with a single molecule. In this section, we turn to the scenario where multiple molecules are present. In the case that no chemical reactions are present, this can be achieved via the Langevin models in combination with a central limit theorem or law of rare events. On the other hand, when chemical reactions are present, interactions between molecules must be considered. To this end, we first consider the case of well-mixed volumes where the chemical reactions can be exactly represented by the chemical reaction master equation. The relation of this stochastic description and deterministic reaction systems arising from chemical reaction network theory is also discussed. For systems with high heterogeneity, the occurrence of both reactions and diffusion are well captured by reaction-diffusion master equations (RDME). Note that the current setting of the RDME does not capture advection.

The Chemical Master Equation and Chemical Reaction Networks

Before introducing the RDME, we first overview the case where only reactions are present in a wellstirred container. Suppose that there are r 0 chemical reactions involving s 0 chemical species S 1 , . . . , S s0 .

The network of chemical reactions is given by

s0 i=1 ν ik S i → s0 i=1 ν ik S i , k = 1, . . . , r 0 , (17)
where ν ik , ν ik ∈ Z ≥0 are stoichiometric coefficients. For each k ∈ {1, . . . , r 0 }, the k-th reaction in [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] means that a complex consisting of ν ik molecules of S i , i = 1, . . . , s 0 , will produce, via reaction, a complex consisting of ν ik molecules of S i , i = 1, . . . , s 0 . We denote by ν k = [ν 1k , . . . , ν s0k ] T and ν k = [ν 1k , . . . , ν s0k ] T the (column) vectors of stoichiometric coefficients. As an illustration, suppose that the chemical species S 1 , S 2 , S 3 are present and reaction k is of the form

S 1 + 2S 2 → S 3 ,
then ν 1k = 1, ν 2k = 2, ν 3k = 0 and ν 1k = 0, ν 2k = 0, ν 3k = 1. With this notation, the reaction systems [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] can be compactly expressed as

ν k → ν k , k = 1, . . . , r 0 . (18)
At time t, the number of molecules for each species is described by a component of the random vector

X(t) = (X 1 (t), . . . , X s0 (t)) ∈ Z s0 ≥0 .
To account for the fact that reactions require individual molecules to be in close proximity, the rate that the k-th reaction occurs when the quantity of molecules for each species is given by x is denoted by λ k (x), called the propensity. If the k-th reaction occurs at time t, the new state is then

X(t) = X(t-) + ν k -ν k , ( 19
)
where t-is an infinitesimal time before t and X(t-) is the state of the system just before the reaction occurs.

The number of times the k-th reaction occurs by time t is given by the counting process [START_REF] Anderson | Stochastic Analysis of Biochemical Systems[END_REF] R

k (t) = Y k t 0 λ k (X(s))ds , (20)
where the Y k , k = 1, . . . , r 0 are independent unit rate Poisson processes. As such, the state of the system at time t is given by

X(t) = X(0) + r0 k=1 R k (t)(ν k -ν k ), (21)
which describes the dynamics of a stochastic chemical reaction system. The associated Kolmogorov forward equation is called the chemical master equation and is given by

dP (x, t) dt = r0 k=1 (λ k (x -ν k )P (x -ν k , t) -λ k (x)P (x, t)) , (22)
where P (x, t) is the probability of the state being x ∈ Z s0 ≥0 at time t, and λ k (x) is the propensity function for k-th reaction, which can be roughly interpreted as the probability per unit time that reaction k will occur given the system is in state x [START_REF] Gillespie | The chemical Langevin equation[END_REF].

As for the Smoluchowski equation, the evolution of the probability distribution function for the chemical master equation does not provide solutions in closed-form. It is therefore desirable to characterize equilibrium solutions, when they exist. It turns out that the conditions an equilibrium distribution exists is closely linked to the properties of an underlying deterministic chemical reaction system.

To illustrate the deterministic chemical reaction network framework (distinct from the above stochastic models) with a common example, we begin with the ubiquitous reversible 1 enzyme-activated chemical reaction network. This network consists of four chemical species: the enzyme E; the substrate S; the intermediate 1 To distinguish with irreversible enzyme reactions where there is no reverse reaction from the product back to the complex;

i.e. k b = 0 in [START_REF] Chahibi | Molecular communication for drug delivery systems: a survey[END_REF] complex ES; and the product P. The set of chemical species in this example is then S E = {E, S, ES, P}. In this system, there are four chemical reactions

E + S k1 k2 ES kcat k b E + P, (23)
where k 1 , k 2 , k cat , k b are positive reaction rate constants indicating the rate of the corresponding chemical reactions.

Denote by R E the set of all reactions in [START_REF] Chahibi | Molecular communication for drug delivery systems: a survey[END_REF], and [START_REF] Anderson | Stochastic Analysis of Biochemical Systems[END_REF], the concentrations of each species in the enzyme-activated system are governed by the following deterministic system of ordinary differential equations

K E = {k 1 , k 2 , k cat , k b }
d[E](t) dt = -k 1 [E](t)[S](t) -k b [E](t)[P](t) + k 2 [ES](t) + k cat [ES](t) d[S](t) dt = -k 1 [E](t)[S](t) + k 2 [ES](t) d[ES](t) dt = k 1 [E](t)[S](t) + k b [E](t)[P](t) -k 2 [ES](t) -k cat [ES](t) d[P](t) dt = -k b [E](t)[P](t) + k cat [ES](t), (24)
with non-negative initial conditions [E](0) = E 0 , [S](0) = S 0 , [ES](0) = ES 0 , and [P](0) = P 0 .

A general chemical reaction system, of which the example above is a special case, which is a chemical reaction network along with the corresponding reaction rate coefficients is defined as follows.

Definition 3.1. A chemical reaction system is the tuple (S, R, K) consisting of a set of chemical species

S = {S 1 , . . . , S s0 }, a set of reactions R = {ν → ν , ν, ν ∈ Z s0
≥0 }, and the rate function K : R → R + .

For general chemical reaction systems, the dynamics can be described as follows. Let z(t) ∈ R s0 be the vector consisting of concentrations of each chemical species at time t. Under mass-action kinetics, the chemical reaction system is then governed by

ż(t) = ν→ν ∈R k ν→ν z(t) ν (ν -ν), (25)
where z(t

) ν = z 1 (t) ν1 z 2 (t) ν2 • • • z s0 (t)
νs 0 and k ν→ν is the reaction rate constant associated to the reaction ν → ν ∈ R, where the notation ν → ν is understood as in [START_REF] Feinberg | Chemical reaction network structure and the stability of complex isothermal reactors-I. the deficiency zero and deficiency one theorems[END_REF].

A key property of chemical reaction networks is whether or not they admit a deterministic steady-state concentration; that is, for each species, the concentration converges to positive value as time t → ∞. It is important to note that these steady-states are not guaranteed to be unique nor even to exist. In fact, a longstanding open problem in the theory of chemical reaction networks is conditions under which the deterministic steady-state exists and the positive equilibrium is the only attracting point, known as the global attractor conjecture [START_REF] Yu | Mathematical analysis of chemical reaction systems[END_REF]. One such condition requires the notion of complex balanced chemical reaction systems.

Definition 3.2. A deterministic chemical reaction system (S, R, K) is said to be complex balanced if and only if there exists a point α ∈ R s0 >0 such that

ν→ν ∈R k ν→ν α ν (ν -ν) = ν →ν∈R k ν →ν α ν (ν -ν ). ( 26
)
Such a point α is called a positive complex balanced equilibrium.

One can show that the enzyme reaction system in ( 24) is complex balanced, since at equilibrium the reaction rate of each forward reaction is equal to the reaction rate of the corresponding backward reaction.

Returning to stochastic chemical reaction systems [START_REF] Arnold | Deterministic limit of the stochastic model of chemical reactions with diffusion[END_REF], an analogous question is whether a stationary distribution exists; that is, whether there existing a limiting distribution π satisfying

π(m) = lim t→∞ Pr (X(t) = x|X(0) = x 0 ) . ( 27
)
For a wide range of stochastic chemical reaction systems, the stationary distribution can be characterized; e.g., the following result in [51, Theorem 4.3].

Theorem 3.1. Let {S, R, K} be a weakly reversible2 deterministic chemical reaction system with rate constants K. Suppose that the deterministically modeled system is complex balanced with equilibrium concentration c ∈ R s0 >0 . Then, the stochastic system has a product form stationary distribution

π(x) = M c x x! , x ∈ Z s0 ≥0 , (28)
where M > 0 is a normalizing constant, and

x! = x 1 !x 2 ! • • • x s0 !.
For the case of two species reversible reaction S 1 S 2 , an application of Theorem 3.1 shows that the distribution of the number of molecules from a given species is binomially distributed [START_REF] Anderson | Product-form stationary distributions for deficiency zero chemical reaction networks[END_REF].

The time scale at which a stochastic chemical reaction system approaches an equilibrium state is highly dependent on environmental conditions (e.g., temperature) which affect reaction propensities and the initial quantity of each chemical species. One can obtain a reasonable estimate of the equilibriation time from the underlying deterministic chemical reaction network. For example, for certain enzyme reactions arising in bacterial chemotaxis, the equilibriation time is on the order of 0.2 s [START_REF] Silversmith | Kinetic characterization of catalysis by the chemotaxis phosphatase CheZ[END_REF].

In the context of molecular communications, the chemical master equation has been used in [START_REF] Pierobon | Noise analysis in ligand-binding reception for molecular communication in nanonetworks[END_REF] in order to model ligand-receptor binding. As will be detailed in Sec. 5, stochastic chemical reaction systems play a useful role in computation for signal processing.

Reaction-Diffusion Master Equation

Having defined stochastic chemical reaction systems, we now introduce, in addition, the presence of diffusion. While stochastic reaction-diffusion systems can be defined for high order chemical reactions, for 1

• • • i -1 i i + 1 • • • N S 1 a 1 i a 2 i S 2 M i (t) = [M 1 i (t), M 2 i (t)]
M i (t) = [M 1 i (t), M 2 i (t)]. M(t) = [M 1 (t), . . . , M i (t), . . . , M N (t)
] is governed by a reaction-diffusion master equation.

the purposes of exposition, we focus here on first-order reactions of the form

S 1 → S 2 S 2 → S 1 , (29)
which have been recently considered in the context of molecular communications [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF][START_REF] Akdeniz | Multi-level equilibrium signaling for molecular communication[END_REF][START_REF] Chou | Extended master equation models for molecular communication networks[END_REF].

In order to capture the effect of small quantities of each chemical species in the system, we consider a stochastic model for the kinetics. To formally describe the scenario, we introduce the following notation.

First, we divide the domain into N small voxels, each with a volume of V Vox as shown in Fig. 3. One main idea is that each voxel, well-mixing is assumed. Let M l i (t), l = 1, 2, i = 1, . . . , N denote the random variable for the number of molecules of species S 1 or S 2 in voxel i at time t. Denote M i (t) = [M 1 i (t), M 2 i (t)] as the state vector in voxel i and the matrix consisting of all state vectors as M(t) = [M 1 (t), . . . , M N (t)]. The probability that M(t) has value m at time t conditioned on there being a quantity of molecules m 0 in each voxel at time t = 0 is then denoted by P (m, t) = Pr(M(t) = m|M(0) = m 0 ), [START_REF] Jamali | Diffusive molecular communications with reactive molecules: Channel modeling and signal design[END_REF] where M(0) is the initial quantity of molecules of each species in each voxel.

Remark 3.1. If there is only a single voxel i = 1, then the system corresponds to a stochastic chemical reaction system with unimolecular reactions where M l 1 (t) corresponds to X l (t) in [START_REF] Farsad | A comprehensive survey of recent advancements in molecular communication[END_REF].

Since the reactions in ( 29) are unimolecular, it follows that in each reaction the number of molecules of the two species involved can only increase or decrease by one. Let 1 l i be the state where the number of molecules in all voxels is zero, except for species l in voxel i. That is, M(t) + 1 l i means that the number of molecules of species l in voxel i is increased by one.

A popular model for stochastic kinetics of molecules is the reaction-diffusion master equation (RDME) [START_REF] Isaacson | The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target[END_REF], also utilized in the context of molecular communications in [START_REF] Chou | Extended master equation models for molecular communication networks[END_REF]. In this model, the diffusive jump ratecorresponding to the rate of the underlying Markov process and not the transition probability-is denoted by κ l ij for each individual molecule of the l-th species moving from voxel j into voxel i, with κ l ii = 0, i = 1, . . . , N .

In particular, the probability per unit time that a molecule of S l diffuses from voxel j to voxel i at time t is given by κ l ij M l j (t).

In general, κ l ij depends on i, j and l; that is, the probability of a molecule diffusing between two voxels is not spatially homogeneous (diffusive jump rates vary from voxel to voxel). The diffusive jump rate out of the boundary voxels is zero, which can be interpreted as a reflective boundary condition 3 . We remark that spatial inhomogeneity of diffusion does not enforce any boundary constraints (other than on the reflective boundary). Nevertheless, if the diffusive jump rates near the boundary of the receiver are small, the flux of probability for molecules of S 1 into the receiver will be low. This behavior can model the impact of a membrane surrounding a receiver in a molecular communication system, which is more difficult for molecules to pass through than in the rest of the fluid.

In the case of mass-action kinetics and first-order reactions, the probability per unit time that a molecule of S l in voxel i reacts at time t is given by a l i M l i (t) with rate constants a l i . In general, the reaction rate is dependent on the voxel index; i.e., we allow for voxels i to have different reaction propensities. For small scale systems typical in molecular communications, the homogeneity in key parameters such as temperature can be assumed, which implies few variations in the reaction propensities within individual voxels. The net change of each chemical species due to the reaction with substrate S l is expressed via the vector

ν l = [ν l,1 , ν l,2 ] ∈ Z 2 ≥0 ,
where ν l,1 is the net gain of molecules of S 1 and ν l,2 is the net gain of molecules of S 2 . The term ν l 1 i indicates that M(t) changes by ν l in the i-th voxel.

In the RDME model, the probability distribution P(m, t) evolves according to the system of differential κ l ij (m l j + 1)P (m + 1 l j -1 l i , t) -κ l ji m l i P (m, t)

+ N i=1 2 l=1 a l i (m l i + 1)P (m -ν l 1 i , t) -a l i m l i P (m, t) , (31)
where the first triple sum corresponds to the net change in probability per unit time due to diffusion and the second double sum corresponds to the net increase in probability per unit time due to chemical reactions.

We highlight that implicitly in [START_REF] Gursoy | Towards high data-rate diffusive molecular communications: Performance enhancement strategies[END_REF], it is assumed that there is a finite, reflective boundary, no degradation of molecules, constant temperature, and no external forces (e.g. due to a linear potential inducing drift).

For more details, we refer the reader to [START_REF] Isaacson | The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target[END_REF][START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF]. A natural question is whether these assumptions can be relaxed. In [START_REF] Noel | Algorithm for mesoscopic advection-diffusion[END_REF], external forces were considered, with verification by numerical simulation. We postpone discussion of relaxing further the assumptions to Sec. 7.

The system of ordinary differential equations in [START_REF] Gursoy | Towards high data-rate diffusive molecular communications: Performance enhancement strategies[END_REF] corresponds to the Kolmogorov forward equation for a continuous-time Markov chain, that is, the evolution of the system state is Markovian. In our setting, due to the reversible reactions and diffusive jump rates detailed above, the Markov chain corresponding to the RDME is irreducible and positive recurrent. Therefore, the stationary distribution for (31) exists.

In the context of molecular communications, the use of the RDME was initiated by [START_REF] Chou | Extended master equation models for molecular communication networks[END_REF] and developed in 3 This is expressed mathematically as homogeneous Neumann boundary conditions.

several followup works [START_REF] Chou | Impact of receiver reaction mechanisms on the performance of molecular communication networks[END_REF][START_REF] Chou | A Markovian approach to the optimal demodulation of diffusion-based molecular communication networks[END_REF][START_REF] Chou | Maximum a-posteriori decoding for diffusion-based molecular communication using analog filters[END_REF][START_REF] Chou | Designing molecular circuits for approximate maximum a posteriori demodulation of concentration modulated signals[END_REF][START_REF] Awan | Generalized solution for the demodulation of reaction shift keying signals in molecular communication networks[END_REF][START_REF] Awan | Improving the capacity of molecular communication using enzymatic reaction cycles[END_REF]. The RDME has also been considered recently in the context of equilibrium signaling [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF][START_REF] Akdeniz | Multi-level equilibrium signaling for molecular communication[END_REF], which will be discussed further in Sec. 4.

For a sufficiently large number of molecules, the stationary distribution can in some cases be approximated via an underlying deterministic system. This approximation theory builds on results originally due to Kurtz for the chemical master equation [START_REF] Kurtz | Limit theorems for sequences of jump Markov processes[END_REF] and related density dependent jump Markov processes. For instance, in the work [START_REF] Arnold | Deterministic limit of the stochastic model of chemical reactions with diffusion[END_REF], limit theorems approximating the solution of the RDME to that of an underlying deterministic reaction-diffusion system have been developed for the case of a single chemical species that may be produced or annihilated. For the RDME in [START_REF] Gursoy | Towards high data-rate diffusive molecular communications: Performance enhancement strategies[END_REF], it may then be expected that for a sufficiently large number of molecules in the system, a good approximation for the concentrations of each species is given by the corresponding deterministic reaction diffusion system.

More precisely, denote by Ω ⊂ R d , d ≥ 1, a bounded domain comprising the communication channel, u l (x, t), l = 1, 2 is the concentration of S l at position x ∈ Ω and time t > 0. By applying the law of mass action, the deterministic reaction-diffusion system modeling the reversible reaction ( 29) reads as

             ∂ t u l -∇ • (D l (x)∇u l ) = 2 i=1 a l,i u i , x ∈ Ω, t > 0, l ∈ {1, 2}, D l (x)∇u l • n = 0, x ∈ ∂Ω, t > 0, l ∈ {1, 2}, u l (x, 0) = u l0 (x), x ∈ Ω, l ∈ {1, 2}, (32)
where the self-reaction rate a l,l is defined by a l,l = -a 3-l,l , for an appropriate choice of reaction rates a l,i and diffusion coefficients D l (x). The vector-valued function n(x) is the outer unit normal defined for x ∈ ∂Ω.

The homogeneous Neumann boundary condition D l (x)∇u l • n = 0 corresponds to a reflective boundary condition, and the initial concentations are assumed non-negative.

We remark that the unimolecular reaction [START_REF] Yang | A comprehensive survey on hybrid communication in context of molecular communication and terahertz communication for body-centric nanonetworks[END_REF] was chosen for the purposes of exposition. In general, one can consider a network of unimolecular reactions with arbitrary number of chemicals. The dynamical behaviour of the corresponding reaction-diffusion system has been completely characterized in [START_REF] Fellner | The entropy method for reaction-diffusion systems without detailed balance: first order chemical reaction networks[END_REF].

For the model studied in [START_REF] Arnold | Deterministic limit of the stochastic model of chemical reactions with diffusion[END_REF], it is possible to rigorously establish central limit theorems [START_REF] Kotelenez | Law of large numbers and central limit theorem for linear chemical reactions with diffusion[END_REF][START_REF] Blount | Comparison of stochastic and deterministic models of a linear chemical reaction with diffusion[END_REF], which facilitates the use of Gaussian approximations for the number of molecules in a given subdomain of the system. For the RDME in [START_REF] Gursoy | Towards high data-rate diffusive molecular communications: Performance enhancement strategies[END_REF], a Gaussian approximation has been justified empirically in [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF].

For the purpose of molecular communications, in [START_REF] Noel | Improving receiver performance of diffusive molecular communication with enzymes[END_REF] a deterministic approximation has been used

to approximate the expected number of molecules in a receiver in the presence of information molecule degrading enzymes in the fluid medium. In [START_REF] Farsad | A molecular communication system using acids, bases and hydrogen ions[END_REF], the deterministic system has been used to model the expected quantity of hydrogen ions when there are reactions between hydrogen ions and acids or bases.

Both of these deterministic models involve nonlinear reaction terms. A similar idea was introduced in [START_REF] Jamali | Diffusive molecular communications with reactive molecules: Channel modeling and signal design[END_REF][START_REF] Farahnak-Ghazani | On medium chemical reaction in diffusion-based molecular communication: A two-way relaying example[END_REF] exploiting a Poisson model with mean determined by an underlying deterministic reaction-diffusion system.

Simulation of the RDME in ( 31) can be performed via the Gillespie algorithm [START_REF] Gillespie | Stochastic simulation of chemical kinetics[END_REF], which has been widely used within molecular communications. Further variations are discussed in [START_REF] Isaacson | The reaction-diffusion master equation as an asymptotic approximation of diffusion to a small target[END_REF][START_REF] Elf | Mesoscopic reaction-diffusion in intracellular signaling[END_REF] and in [START_REF] Noel | Algorithm for mesoscopic advection-diffusion[END_REF] which incorporates the presence of drift into the RDME. --→ S 1 . Once the information molecules S 2 are released from the transmitter to transmit a message, the receiver (RX) observes these molecules until the system reaches the equilibrium state.

Modulation and Detection for Channels Governed by Stochastic Models

A basic problem in molecular communications is how to encode messages to be transmitted in the state of a set of molecules, known as modulation. On the receiver side, the challenge is then to detect which message was sent given the observations by the receiving device. The problem of modulation and detection has been a central question since the inception of molecular communications and there is now a wide body of work, recently surveyed in [START_REF] Kuran | A survey on modulation techniques in molecular communication via diffusion[END_REF]. In this section, we provide a brief review of the main approaches for modulation and detection, when molecular motion is governed by the RDME or Smoluchowski equation.

The general problem of maximum a posteriori detection can be written as

m = arg max m∈M Pr(N Rx |m), ( 33
)
where M is the finite set of possible messages and N Rx is the number of molecules observed at the receiver.

More concretely, suppose that the communication system is time slotted, with time slot duration T s . Consider the n-th time slot. Due to the previous n -1 transmissions, there are N Tx,l (nT s ), l = 1, 2 molecules of each species S l in the transmitter. At a time nT s + δ shortly after the beginning of the time slot, the transmitter produces a quantity of S 1 depending on the bit to be transmitted. In particular, for binary signaling with

M = {0, 1}, N Tx,1 (nT s + δ) =    N Tx,1 (nT s ) + ∆ s n = 1, N Tx,1 (nT s ) s n = 0, (34)
for δ > 0 a sufficiently small period of time, and ∆ > 0 corresponding to the number of molecules emitted with message m = 1 (not to be confused with the state m for the RDME). Zero molecules are emitted when m = 0.

Detection in General Channels

The key challenge in this setting is to reliably detect the transmitted binary symbol; i.e., whether m = 0 or m = 1. A standard method, known as concentration shift keying (CSK), is often defined as the emission of different quantities of molecules for each message; however, in practice, it is ubiquitous to restrict the receiver to sample at a time T s such that the quantity of molecules at the receiver is maximized in some sense (e.g., the expected number of molecules at the receiver is maximized). For the simplest diffusion model ( 5), the statistics for the number of molecules in the receiver, N Rx , can be computed in closed-form for an infinite domain. However, a characterization of the statistics is more challenging when boundary conditions or chemical reactions are present, although some results have been obtained in this setting (see e.g., [START_REF] Deng | Analyzing large-scale multiuser molecular communication via 3-d stochastic geometry[END_REF]).

Extensions to non-binary signaling schemes have been studied in, for example, [START_REF] Akdeniz | Multi-level equilibrium signaling for molecular communication[END_REF].

In the absence of chemical reactions, a number of approaches have been proposed; e.g., type based, timing-based modulation, and spatial-based. Type-based signaling is one of the most widespread, where multiple types of molecules can be used as information-carrying molecules. These benefits can be obtained by representing each symbol with different types of molecules as proposed in [START_REF] Kuran | Modulation techniques for communication via diffusion in nanonetworks[END_REF] and known as Molecule Shift Keying (MoSK). Alternatively, in Molecular Concentration Shift Keying (MCSK), it has been proposed to use different types of molecules for two consecutive transmissions to reduce interference [START_REF] Arjmandi | Diffusion-based nanonetworking: A new modulation technique and performance analysis[END_REF]. In [START_REF] Kabir | D-MoSK modulation in molecular communications[END_REF], different types of molecules are used as orthogonal channels and it is proposed to use multiple CSK channels simultaneously which is called as D-MoSK in the literature. Timing based approaches are based on representing the intended message by release time as done for binary signaling in [START_REF] Garralda | Diffusion-based physical channel identification in molecular nanonetworks[END_REF] and multi level signaling in [START_REF] Akdeniz | Position-based modulation in molecular communications[END_REF]. Spatial-based approaches can also be developed by utilizing MIMO concepts [START_REF] Koo | Detection algorithms for molecular MIMO[END_REF] and index modulation [START_REF] Gursoy | Index modulation for molecular communication via diffusion systems[END_REF]. In addition to these, hybrid approaches have also been proposed such as [START_REF] Gursoy | Pulse position-based spatial modulation for molecular communications[END_REF] and [START_REF] Gursoy | Concentration and position-based hybrid modulation scheme for molecular communications[END_REF].

In addition to these works and methods presented here, there are several other approaches available in the literature and we refer the interested reader to the surveys [START_REF] Kuscu | Transmitter and receiver architectures for molecular communications: A survey on physical design with modulation, coding, and detection techniques[END_REF] and [START_REF] Kuran | A survey on modulation techniques in molecular communication via diffusion[END_REF] for more details.

Detection in RDME Channels

On the other hand, few approaches have been developed for molecular motion governed by the RDME.

A recently proposed approach [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF][START_REF] Akdeniz | Multi-level equilibrium signaling for molecular communication[END_REF][START_REF] Egan | Equilibrium signaling in spatially inhomogeneous diffusion and external forces[END_REF], instead exploits the equilibrium behavior of the RDME in the presence of a reflective boundary-leading to a scheme known as equilibrium signaling, which is also applicable to a general family of Langevin diffusions [START_REF] Egan | Equilibrium signaling in spatially inhomogeneous diffusion and external forces[END_REF]. In equilibrium signaling, T s is chosen to be sufficiently large such that the statistics for the number of molecules in the receiver can be approximated in closed-form via [START_REF] Deng | Modeling and simulation of molecular communication systems with a reversible adsorption receiver[END_REF] or a Gaussian approximation for the RDME (see [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF]). An advantage of this approach is that the stationary distribution is independent of most system parameters (e.g., diffusion coefficients, presence of obstacles, and changes to the location of the transmitter or receiver). With the stationary distribution, it is then feasible to derive near-optimal detection rules [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF][START_REF] Egan | Equilibrium signaling in spatially inhomogeneous diffusion and external forces[END_REF].

To illustrate the idea behind equilibrium signaling, consider Fig. 4. Here, the transmitter seeks to send a message to the receiver by releasing information molecules of species S 2 through the communication channel be obtained analytically. Therefore, extensive simulations are typically required if conventional detection schemes are utilized. Furthermore, the detection methods will have different performance when the location of the transmitter or receiver, the fluid boundary, or obstacles varies. For example, if CSK is used for the channels given in Fig. 5, if the channel parameters (e.g., diffusion coefficients) are kept constant, we may expect different performance in the channels in Fig. 5(a) and 5(b) due to the differences in geometry. This is true even if the channels in This phenomenon is illustrated in Fig 6 . In this figure, three different communication schemes are evaluated: equilibrium signaling with two species S 1 , S 2 ; equilibrium signaling with one species S 1 ; and CSK at the peak time of the observed signal (CSK opt ). Note that for equilibrium signaling with one or two species, all parameters are kept the same including the communication time slot which is equal to the approximate equilibrium time and can be approximated as in [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF]. As expected and can also be observed in Fig 6, as the volume of the system and the receiver are kept constant, the performance of equilibrium signaling and equilibrium signaling with one species does not change for different scenarios.

On the other hand, using an optimum time slot for CSK communication can be obtained by selecting the time which maximizes the expected number of molecules as suggested in [START_REF] Llatser | Detection techniques for diffusion-based molecular communication[END_REF]. Therefore for all different scenarios-assuming that all channel parameters are known-it is possible to obtain an optimal detector for CSK. As can be seen in The main alternative to equilibrium signaling in the context of the RDME is known as reaction-shift keying (RSK) [START_REF] Awan | Generalized solution for the demodulation of reaction shift keying signals in molecular communication networks[END_REF], where a continuous-time molecular signal is transmitted and detected via chemical reactions. In common with equilibrium signaling (with multiple chemical species) is the use of chemical reactions in the transmitter and receiver; however, like CSK and MoSK, RSK does not exploit an equilibrium state and requires computation of the expected quantity of molecules in the receiver at a given time. As this requires a full solution of the equation of motion, good estimates of transmitter-receiver distances and diffusion coefficients are implicitly required. This again contrasts with equilibrium signaling-which also has a dramatically lower complexity-where these estimates are not required. Nevertheless, RSK does not require the system to reach an approximate equilibrium state and has advantages for systems with short symbol periods.

Inference via Stochastic Reaction-Diffusion Systems

It is widely known that chemical reaction networks, particularly genetic circuits, form a powerful means of performing computation. This has recently been developed in the context of molecular communications via microfluidic circuits [START_REF] Bi | A survey of molecular communication in cell biology: Establishing a new hierarchy for interdisciplinary applications[END_REF] and for encoding and decoding of channel codes [START_REF] Marcone | Parity-check coding based on genetic circuits for engineered molecular communication between biological cells[END_REF]. It is less commonly observed, but equally true, that stochastic reaction-diffusion can be viewed as stochastic molecular circuits for statistical inference. In this section, we review the use of stochastic reaction-diffusion systems to solve statistical inference problems.

Receiver Design

In molecular communications, a key question is how the detection rules discussed in Section 4 can be implemented. Due to the small scale and limited resources available to receivers, particularly for nanoscale devices, computation via biochemical processes is highly desirable. The main work in this direction is due to Chou and his collaborators [START_REF] Chou | Impact of receiver reaction mechanisms on the performance of molecular communication networks[END_REF][START_REF] Awan | Generalized solution for the demodulation of reaction shift keying signals in molecular communication networks[END_REF][START_REF] Chou | A Markovian approach to the optimal demodulation of diffusion-based molecular communication networks[END_REF][START_REF] Chou | Maximum a-posteriori decoding for diffusion-based molecular communication using analog filters[END_REF][START_REF] Awan | Improving the capacity of molecular communication using enzymatic reaction cycles[END_REF][START_REF] Riaz | Using spatial partitioning to reduce the bit error rate of diffusion-based molecular communications[END_REF], which exploits stochastic reaction-diffusion systems within the receiver.

A key example of this approach is in [START_REF] Chou | Designing molecular circuits for approximate maximum a posteriori demodulation of concentration modulated signals[END_REF], which considers a receiver consisting of receptors that can exist in two states: an inactive state X; and an active state X * . An inactive receptor can be activated by a signaling molecule, where the system is governed by the reactions

X + X g+ → S + X * X * g- → X, (38)
where g + and g -are the propensities (discussed in detail in Sec. 3) and S is a signaling molecule.

Detection is then based on the log posterior probability that symbol k has been sent given history X * (t), denoted by L k (t) = log P (k|X * (t)). Chou then showed that L k (t) obeys

dL k (t) dt = dx * (t) dt + log(E[n R (t)|k, X * (t)]) -g + (M -x * (t))E[n R (t)|k, X * (t)], (39)
where

+ log(σ k (t)) -g + (M -x * (t))E[n R (t)|k, X * (t)], (40)
with Z k (t) is initialized to the logarithm of the prior probability that symbol k is transmitted. Detection is then performed via k = argmax k Z k (t). While ( 40) can be computed numerically, a key contribution in [START_REF] Chou | Designing molecular circuits for approximate maximum a posteriori demodulation of concentration modulated signals[END_REF] shows that it can also be approximated via a molecular circuit based on the gene promotor circuit DCS2, which is found in the yeast Saccharomyces cerevisiae.

Parameter Estimation

One example of an inference problem relevant for molecular communications is to estimate the quantity of molecules of different types via non-specific receptors on the surface of the transmitter. In this setting, the number of each type of molecule produced by the receptors depends linearly on the number of molecules that bind with the receptors. In particular,

r = Ol, (41)
where r ∈ Z nR ≥0 is the number of each type of receptor molecule produced, which can be directly observed within the transmitting device. The vector l ∈ Z nL ≥0 is the number of each type of species in the environment that bind to receptors on the transmitter. The matrix O is an n R × n L observation matrix determining the production of receptor molecules given the number of molecules from the environment that bind with each receptor.

In order to estimate the quantity l, a Bayesian approach can be adopted, to yield samples from Pr(l = u|r).

Motivated by Theorem 3.1, consider the product-Poisson prior for the number of molecules in the environment that bind with the receptors, given by

Pr(l = u) = nL i=1 e -qi qui i u i ! , (42)
where l ∈ Z nL ≥0 , r ∈ Z nR ≥0 , q = [q 1 , . . . , qnL ] T are the positive rate parameters. In [START_REF] Virinchi | A stochastic molecular scheme for an artificial cell to infer its environment from partial observations[END_REF], it has been shown that the distribution of the posterior Pr(l|r) with the product-Poisson prior can be sampled from the quantities of molecules of species X 1 , . . . , X n at a sufficiently large time t.

This approach has been adopted within molecular communications for signal processing in [START_REF] Egan | Biological circuits for detection in MoSK-based molecular communication[END_REF][START_REF] Egan | Estimation and optimization for molecular communications with a coexistence constraint[END_REF].

Indeed, stochastic chemical reaction networks provide a powerful framework for computation. In particular, chemical reaction networks have been shown to implement expectation-maximization [START_REF] Virinchi | A reaction network scheme which implements the EM algorithm[END_REF] and Boltzmann machines (binary stochastic neural networks) [START_REF] Poole | Chemical Boltzmann machines[END_REF][START_REF] Virinchi | A reaction network scheme which implements the EM algorithm[END_REF]. This work focuses on the problem where it is desirable to maximize the likelihood Pr(x|θ), where x is data and θ is the parameter to be estimated. A key assumption is that x is drawn from an exponential family and only observations s = Ax are available; that is, the data x is a latent variable and the observations s are obtained via a linear transformation A. We note that a similar problem arises in the context of compressed sensing, which has been widely applied in various communication systems (see e.g., [START_REF] Leinonen | Compressed sensing with applications in wireless networks[END_REF]).

Stochastic Biological Environments

In molecular communication systems embedded in a biological environment, there may be interactions between information-carrying molecules and external biological organisms or competition for genetic resources in the process of encoding or decoding information. This problem, addressed in [START_REF] Egan | Strategies for coexistence in molecular communication[END_REF][START_REF] Akdeniz | A reactive signaling approach to ensure coexistence between molecular communication and external biochemical systems[END_REF][START_REF] Mcbride | The effect of loads in molecular communications[END_REF][START_REF] Egan | Coexistence in molecular communications[END_REF][START_REF] Gyorgy | Limitations and trade-offs in gene expression due to competition for shared cellular resources[END_REF][START_REF] Loscri | On the interaction between a nanoparticulate system and the human body in body area networks[END_REF] and often called coexistence, must account for the stochasticity due to small numbers of information carrying molecules or, for example, limited ribosomes during DNA translation [START_REF] Gyorgy | Limitations and trade-offs in gene expression due to competition for shared cellular resources[END_REF]. As such, stochastic reactiondiffusion systems can play an important role.

The problem of coexistence is how to ensure that information-carrying molecules do not significantly perturb an external biological system, where the external biological system is modeled via a stochastic chemical reaction system, illustrated in Fig. 7 for two types of signaling molecules. In [START_REF] Mcbride | The effect of loads in molecular communications[END_REF][START_REF] Gyorgy | Limitations and trade-offs in gene expression due to competition for shared cellular resources[END_REF], components of a biological system are modeled via deterministic chemical reaction systems. There are two basic reasons why the interactions between a molecular communication system and an external biological system may be better modeled via stochastic models:

(i) the motion of the information-carrying molecules is governed via Langevin or RDME models;

(ii) when a biological system has relatively few molecules of a given species, the internal dynamics are better captured via stochastic, typically chemical reaction network, models.

A general framework to establish information theoretic limits of molecular communication such that perturbations of the biological system are limited was developed in [START_REF] Egan | Coexistence in molecular communications[END_REF]. At its heart, this approach placed constraints on the statistics for the number of information-carrying molecules that bind with an external biological system. As such, techniques from covert communications [START_REF] Bloch | Covert communication over noisy channels: A resolvability perspective[END_REF] can be adapted to derive scaling laws for the number of messages that can be reliably communicated. An important conclusion is for a code blocklength n (corresponding to n symbols), the number of messages scales with √ n when reliability and coexistence constraints are imposed. This is in contrast with the case where the biological system is not present, where via Shannon's noisy channel coding theorem [START_REF] Shannon | A mathematical theory of communication[END_REF], the number of messages that can be reliably communicated scales with n. This suggests that coding in the presence of external biochemical systems may need to be adapted; namely, utilizing longer codes or transmission of less information.

From the perspective of practical modulation schemes, in [START_REF] Akdeniz | A reactive signaling approach to ensure coexistence between molecular communication and external biochemical systems[END_REF] an approach based on reactive signaling [START_REF] Jamali | Diffusive molecular communications with reactive molecules: Channel modeling and signal design[END_REF] was proposed. In particular, the number and time of release of two signaling molecules that can react with each other was optimized subject to a constraint on the expected number of molecules that pass through an external biological system.

Open Problems

As detailed in the previous sections, there has been significant progress in developing molecular communication systems in the presence of and exploiting stochastic models. Nevertheless, each scheme has some limits on its applicability. In this section, we highlight several open problems that we believe will enhance the range of applications that can be supported via molecular communications.

Detection

As noted in Sec. 4, a key challenge is to derive reliable detection rules for various modulation schemes in the presence of channels governed by stochastic models. While the simplest form of Langevin diffusionnamely standard Brownian motion-has been widely studied, there is very little work considering inhomogeneous diffusion and external forces.

The main work in this direction is [START_REF] Egan | Equilibrium signaling in spatially inhomogeneous diffusion and external forces[END_REF], based on equilibrium signaling; however, a key limitation is the requirement that information-carrying molecules do not degrade. Indeed, if molecules degrade, the equilibrium distribution for the number of molecules is concentrated around zero; that is, no informationcarrying molecules remain in the system at the sampling time. A potential means of extending the work in [START_REF] Egan | Equilibrium signaling in spatially inhomogeneous diffusion and external forces[END_REF] to the case of molecule degradation is to instead study quasi-equilibrium distributions. While quasiequilibrium distributions have been considered in various chemical reaction systems [START_REF] Anderson | Stochastic analysis of biochemical reaction networks with absolute concentration robustness[END_REF], they have not yet been studied in the context of molecular communications. A similar issue also arises for the equilibrium signaling scheme in [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF], in the context of the RDME model.

Another challenge related to equilibrium signaling is the choice of the sampling time. In order to select the sampling time, it is necessary to determine when the system is "close" to equilibrium. In [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF], this choice was based on analysis of an approximate underlying deterministic system. In order to carry out this analysis for general RDME models, a characterization of the optimal rate of convergence to equilibrium should be investigated. While bounds on the equilibrium convergence rate are known for fairly general deterministic reaction-diffusion systems [START_REF] Fellner | The entropy method for reaction-diffusion systems without detailed balance: first order chemical reaction networks[END_REF], the optimal rate has not yet been obtained and remains a key open problem.

An important limitation of the vast majority of work on molecular communications is the assumption that information-carrying molecules do not interact. Such an interaction can arise due to electrostatic forces; e.g., if the molecules are ionized. In order to derive reliable detection rules in this scenario it appears necessary to draw techniques from the field of interacting particle systems [START_REF] Liggett | Interacting Particle Systems[END_REF].

Finally, when the fluid media is, for example, turbulent, Langevin diffusion is not an appropriate model.

Instead, anomalous diffusion models are more appropriate. While there is some work in this direction [START_REF] Mai | Event detection in molecular communication networks with anomalous diffusion[END_REF][START_REF] Trinh | Molecular communication with anomalous diffusion in stochastic nanonetworks[END_REF][START_REF] Trinh | Molecular communication in H-diffusion[END_REF][START_REF] Cao | Anomalous diffusion in molecular communication[END_REF][START_REF] Chouhan | Molecular communication in fractional diffusive channel[END_REF][START_REF] Huang | Molecular communication via subdiffusion with a spherical absorbing receiver[END_REF], the intractability of the solutions to the resulting generalized Fokker-Planck equations introduces new difficulties for the design of detection rules.

Stochastic Molecular Circuits

At present, there has been very limited use of stochastic models to implement molecular circuits for the purpose of signal processing. As transmitting and receiving devices typically require even limited capabilities to perform coding and decoding, as well as parameter estimation, there are several opportunities for new applications of stochastic models in these settings. On the one hand, there are extensions of the work in [START_REF] Chou | Designing molecular circuits for approximate maximum a posteriori demodulation of concentration modulated signals[END_REF] exploiting RDME models for detection. On the other hand, exploiting the capabilities of stochastic chemical reaction systems, as in the work in [START_REF] Virinchi | A stochastic molecular scheme for an artificial cell to infer its environment from partial observations[END_REF][START_REF] Virinchi | A reaction network scheme which implements the EM algorithm[END_REF], appears to be a fruitful direction for future work.

A specific open problem is the development of molecular circuits for detection and channel estimation, which can account for non-specific receptors. Indeed, nanoscale networks are expected to exploit devices that are enclosed by a membrane and require a mechanism to transfer information-carrying molecules into a receiving. Such receptors can bind to information-carrying molecules, and potentially also other chemical species in the environment. The techniques developed in [START_REF] Virinchi | A stochastic molecular scheme for an artificial cell to infer its environment from partial observations[END_REF][START_REF] Virinchi | A reaction network scheme which implements the EM algorithm[END_REF] are ideally suited to coping with this problem and an initial study has been carried out in [START_REF] Egan | Biological circuits for detection in MoSK-based molecular communication[END_REF].

External Biochemical Systems

Despite the possibility of information-carrying molecules interacting with external biochemical systems, there has been very limited work investigating how these interactions can be mitigated and the impact on the performance of the communication system. This is particularly true when the internal stochastic dynamics of the biochemical systems are considered.

A starting point to begin characterizing the performance reduction can be found in [START_REF] Egan | Strategies for coexistence in molecular communication[END_REF], which exploited a fundamental connection between the coexistence problem and covert communications to derive information theoretic limits of communication. However, this work assumed that the biological system does not have memory; that is, the biological system is "reset" after each symbol is transmitted. As realistic biological systems are likely to have memory of their previous state (due to the Markovian nature of their dynamics), addressing this issue seems important in order to understand the true limits of molecular communication systems in the presence of external biochemical systems. A further extension would be to consider the "finite blocklength regime" as considered in [START_REF] Polyanskiy | Channel coding rate in the finite blocklength regime[END_REF][START_REF] Lin | Weak flip codes and their optimality on the binary erasure channel[END_REF].

Another open issue is the development of practical signaling strategies to minimize the impact of molecular communication systems on their environment. One approach is to build on recent work studying the impact of multiple receivers on reliability [START_REF] Lu | The effect of two receivers on broadcast molecular communication systems[END_REF]. Another angle is to view the molecular communication link as a perturbation of the reaction-diffusion system governing the external biological system. Here, variations on perturbed global attractor theorems in chemical reaction network theory may play a useful role [START_REF] Gopalkrishnan | A scheme for molecular computation of maximum likelihood estimators for loglinear models[END_REF].

Distributed Spectroscopy

Many proposed engineering applications of molecular communications are in highly uncertain environments; e.g., in vitro applications in the human body. Nevertheless, in some settings it is possible to design the communication channel in the sense that molecular motion can be controlled. A key example of such a setting is in microfluidic systems, where the stochastic model governing motion can be adapted to the engineering application. A basic question is then how to optimize the channel?

For deterministic models of microfluidic channels, there have been several works investigating how channel parameters affect the dynamics [START_REF] Leo | Communications and switching in microfluidic systems: pure hydrodynamic control for networking labs-on-a-chip[END_REF][START_REF] Hamidović | Information encoding in droplet-based microfluidic systems: First practical study[END_REF][START_REF] Farsad | On-chip molecular communication: analysis and design[END_REF][START_REF] Bicen | Shannon meets Fick on the microfluidic channel: diffusion limit to sum broadcast capacity for molecular communication[END_REF][START_REF] Krishnaswamy | Time-elapse communication: bacterial communication on a microfluidic chip[END_REF][START_REF] Galluccio | Capacity of a binary droplet-based microfluidic channel with memory and anticipation for flow-induced molecular communications[END_REF][START_REF] Bicen | System-theoretic analysis and least-squares design of microfluidic channels for flow-induced molecular communication[END_REF]. For channels governed by stochastic models, how to tailor the channel to an engineering application has not been widely investigated. One exception is the recent work in [START_REF] Akdeniz | A molecular communication scheme to estimate the state of biochemical processes on a lab-on-a-chip[END_REF][START_REF] Akdeniz | Molecular communication for equilibrium state estimation in biochemical processes on a lab-on-a-chip[END_REF], where channel design for the purpose of spectroscopy was investigated.

Indeed, a key challenge in molecular biology is to understand how environmental conditions affect the dynamics of biochemical processes. To answer this questions in practice, it is necessary to measure the concentration of chemical species produced by the biochemical process under observation-known as an assay.

However, it is not always straightforward to perform such measurements in situ or in vivo as the spectroscopy measurement devices require special environmental conditions; for example, in force spectroscopy utilizing atomic force microscopy, well-defined, flat and non-reactive substrates are required [START_REF] Ferrari | Interaction of cement model systems with superplasticizers investigated by atomic force microscopy, zeta potential, and adsorption measures[END_REF]. In this case, it may desirable to extract samples from the process and perform spectroscopy elsewhere.

One approach to performing state estimation of a biochemical process in a distributed fashion has recently been proposed in [START_REF] Akdeniz | A molecular communication scheme to estimate the state of biochemical processes on a lab-on-a-chip[END_REF]. This approach is illustrated in Fig. 8, where a biochemical process under observation contained in a chamber, connected to spectroscopy chamber-where the quantity of a given chemical species can be directly observed-via a microfluidic channel. A key open problem is how to design the joint sensing, communication and estimation system tailored to particular spectroscopy methods. Such design relies heavily on the behavior of a RDME analogous. In particular, instead of a molecular communication detection problem (as considered in Sec. 7.1), it is necessary to estimate an unknown quantity of a given chemical species. An even more challenging variant is to design such systems for in vivo applications, perhaps building on the Raman spectroscopy-based methods in [START_REF] Guo | Cooperative Raman spectroscopy for real-time in vivo nano-biosensing[END_REF].

Conclusions

Stochastic models for reaction and diffusion play a critical role in developing molecular communication systems. This is due to the need to characterize the statistics for the number of molecules that are observed by a receiver, to implement signal processing via molecular circuits, and to mitigate the impact on external biochemical processes. In this paper, we have surveyed general stochastic models and their applications in recent work. By outlining a number of open problems, we hope to encourage further developments to address outstanding issues of both an engineering and a mathematical nature.

Figure 1 :

 1 Figure 1: An overview of the stochastic models of reaction and diffusion considered in this paper. Models lower down in the hierarchy are special cases of models higher up.

Figure 2 :

 2 Figure 2: Solutions at time t = 2 for diffusion equation ∂tp(x, t) = D∂xxp(x, t), x ∈ Ω = (0, 1), D = 0.01 subject to different boundary conditions. For the homogeneous Dirichlet boundary condition, there are no molecules at the boundary for all time.For the homogeneous Neumann boundary condition, the flux of molecules at the boundary, in this case it is simply the slope of the solution state, is zero. For the periodic boundary condition, the number of molecules at the left and right boundary are always identical. Finally, for Robin boundary condition, the flux at the boundary is proportional to probability density of the molecule, i.e. the slope is steeper if the probability density is larger.

Figure 3 :

 3 Figure 3: RDME model. The environment is divided into N voxels whose volumes are equal and V Vox . The molecules in ith voxel at time t are represented asM i (t) = [M 1 i (t), M 2 i (t)]. M(t) = [M 1 (t), . . . , M i (t), . . . , M N (t)] is governed by a reaction-diffusion master equation.

Figure 4 : 1 -

 41 Figure 4: Communication Channel Model for equilibrium signaling. Transmitter (Tx), uses ∆ new S 1 molecules to produce information molecules S 2 via S 1 a 1 --→ S 2 and S 2 a 2

  Fig 5(b) and 5(c) have the same fluid boundary and obstacles, since the distance between the transmitter and the receiver are different. On the other hand, in equilibrium signaling, the channel response only depends on the volume of the channel and the receiver. As such, the performance in the three different channels in Fig. 5 is the same.

Fig 6 ,

 6 these curves differ for each scenario and equilibrium signaling can still outperform standard implementations of CSK.

Figure 5 :

 5 Figure 5: Channel with different geometries. Blue boxes indicate obstacles.

Figure 6 :

 6 Figure 6: Performance comparison of CSK and equilibrium signaling for the three scenarios presented in Fig. 5. Note that for conventional CSK, the performances differ in different channels as expected while in equilibrium signaling since the volumes are kept constant, same BER values are observed.

Figure 7 :

 7 Figure 7: Molecular communication in the presence of a biological system [92].

Figure 8 :

 8 Figure 8: Distributed spectroscopy model in[START_REF] Akdeniz | A molecular communication scheme to estimate the state of biochemical processes on a lab-on-a-chip[END_REF][START_REF] Akdeniz | Molecular communication for equilibrium state estimation in biochemical processes on a lab-on-a-chip[END_REF].

Table 1 :

 1 Summaries of Surveys on Molecular Communication

Table 2 :

 2 Notation.

	Variable	Definition			Variable	Definition
	N	Number of voxels in the system.	F(x)	external force
	γ(x)	Drag parameter in Langevin dif-	m	mass of a molecule
		fusion.
	σ(x)	Covariance matrix of the random	x(t)	location of a particle at time t
		perturbation in Langevin diffu-
		sion.
	S 1 , S 2	Chemical species.		E[•]	Expectation operator.
	a l i , l = 1, 2.	Reaction rate constants.	∂ t	Partial derivative with respect to
						t.
	ν l,k	Quantity of each species l pro-	∂ xx	Second partial derivative with re-
		duced or removed in reaction k.		spect to x.
	N Rx j (t)	Number of S j molecules in the	∇ 2 f (x)
		receiver at time t.		= ∇ • ∇f (x)	Laplacian of a function f : R d →
						R.
	T s	Communication time interval.	V Vox	volume of a voxel
	∆	Number	of	transmitted	S m n	Binary sequence with length n in
		molecules for each bit 1 trans-		which m bits are 1
		mission.
	D

M i (t) State vector of voxel i in time t. ∇f Gradient of a function f : R d → R. κ l ij Diffusive jump rate. ∇ • f Divergence of a vector field f : R d → R d . l (x) Spatially dependent diffusion coefficient for substance S l d Spatial dimension of the system

  such behavior for the simple case of the diffusion equation in one-dimension

	∂ t p(x, t) = D∂ xx p(x, t) for x ∈ Ω := (0, 1), D > 0, subject to either (a): Homogeneous Dirichlet boundary
	condition p(0, t) = p(1, t) = 0, (b): Homogeneous Neumann boundary condition p (0, t) = p (1, t) = 0,
	(c): Periodic boundary condition p(0, t) = p(1, t), or (d): Robin boundary condition p (1, t) = kp(1, t) and
	-p (0, t) = kp(0, t) for some k > 0.
	For the general Fokker-Planck equation (4) in a bounded domain Ω ⊂ R d with boundary ∂Ω, the boundary
	conditions in the case of reflective boundaries (or the Neumann boundary condition) are given by
	n(x) • j(x, t|x 0 , t 0

  the set of reaction rate constants. The tuple (S E , R E , K E ) is called a chemical reaction network. In order to characterize this network, it is also necessary to consider suitable chemical kinetics (e.g., mass-action or Michaelis-Menten).

	Let [E](t), [S](t), [ES](t), [P](t) denote the concentration of each chemical species at time t. Under the standard
	assumption of mass-action kinetics

  x * (t) denotes the number of molecules X * at time t and n R (t) is the number of signaling molecules in the receiver. Computing the term E[n R (t)|k, X * (t)] is computationally challenging and therefore in[START_REF] Chou | Designing molecular circuits for approximate maximum a posteriori demodulation of concentration modulated signals[END_REF], itis approximated via σ k (t) = E[n R (t)|k]. This leads to the demodulation filter

	dZ k (t) dt	=	dx * (t) dt

A weakly reversible chemical reaction system means that for any sequence of reactions from one chemical complex to another, there exists a reverse sequence of reactions to return to the original complex.

by producing molecules via the first order chemical reactions S 1 a 1 -→ S 2 and S 2 a 2 -→ S 1 . In equilibrium signaling, the receiver ideally observes the quantity of molecules of S 1 , S 2 in its domain at the sampling time T s → ∞. In practice, as shown in [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF], the required time t * to reach the equilibrium state of the system with a tolerable error, can be much much smaller and such time can be derived in terms of channel parameters. Once, T s is chosen as t * , the transmitter communicates by releasing and the receiver observes the quantity of molecules in its domain.

One can view equilibrium signaling as lying in the family of CSK schemes, although-as noted earlier-it differs from how such schemes have been designed in practice. We also remark that equilibrium signaling can be used to increase the data rate for multi-level signaling [START_REF] Akdeniz | Multi-level equilibrium signaling for molecular communication[END_REF]. The scheme is also highly robust to uncertainties in diffusion coefficients, the shape of the boundary, and locations of the transmitter and receiver [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF] To illustrate the robustness of equilibrium signaling to uncertainties in the channel, suppose that a sequence of bits, s 1 , . . . , s n , over a period of n sampling intervals is sent. We can represent this sequence as S m n to show that this sequence has m transmissions (each transmission corresponds to a bit 1). Further, let N Rx,1 (nT s |S m n ) and N Rx,2 (nT s |S m n ) denote the number of molecules of species S 1 and S 2 , respectively, observed by the receiver at the end of the n-th symbol period (i.e., at time (n + 1)T s ), given the transmitted sequence S m n . Under these circumstances with sufficiently large ∆ as shown in [START_REF] Akdeniz | Equilibrium signaling: Molecular communication robust to geometry uncertainties[END_REF], the distribution of molecules in the receiver can be well approximated by

where µ r,1 , µ r,2 > 0 are known constants, dependent on the volume of the enclosing container and not the specific geometry, and N (µ, σ 2 ) denotes the Gaussian law with mean µ and variance σ 2 . In particular,

where V tot = N V Vox is the total volume of the system. Moreover,

Therefore, unlike other signaling schemes, in equilibrium signaling, there is no need to know exact geometry of the channel or the distance between the transmitter and the receiver. The only required information to model the signal is the volume of the channel and the volume of the receiver.

To show the performance of equilibrium signaling, we focus on the scenario of a bounded fluid domain and the presence of obstacles as in the three different channels presented in Fig. 5. In order to derive detection rules in these channels, it is necessary to evaluate the channel impulse responses, which cannot