Channel Modeling for Diffusive Molecular Communication - A Tutorial Review
Vahid Jamali, Arman Ahmadzadeh, Wayan Wicke, Adam Noel, and Robert, Schober
TL;DR
This tutorial review comprehensively discusses mathematical channel models for diffusive molecular communication, covering deterministic, statistical, time-varying, simulation-driven, and experimental models, highlighting challenges and future research directions.
Contribution
It provides a unified framework for various MC channel models, integrating biological, chemical, and physical phenomena, and discusses advanced scenarios like mobile nanomachines.
Findings
Deterministic models predict expected molecule counts at receivers.
Statistical models characterize observed molecule variability.
Models for mobile MC systems enable dynamic application analysis.
Abstract
Molecular communication (MC) is a new communication engineering paradigm where molecules are employed as information carriers. MC systems are expected to enable new revolutionary applications such as sensing of target substances in biotechnology, smart drug delivery in medicine, and monitoring of oil pipelines or chemical reactors in industrial settings. As for any other kind of communication, simple yet sufficiently accurate channel models are needed for the design, analysis, and efficient operation of MC systems. In this paper, we provide a tutorial review on mathematical channel modeling for diffusive MC systems. The considered end-to-end MC channel models incorporate the effects of the release mechanism, the MC environment, and the reception mechanism on the observed information molecules. Thereby, the various existing models for the different components of an MC system are…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 57
Figure 57
Figure 5| Ref | Transmitter | Physical Channel | Receiver | Dim | (32) | |||||||||||||
| Surface | PG | REM | Diff | Adv | Geo | D&P | Passive | Active | SP | Joi | Sep | |||||||
| w/o UCA | w/ UCA | Sto | Det | |||||||||||||||
| (37), (38) | [33] | point | ✓ | ✓ | ✓ | 3D | ||||||||||||
| (41) | [73] | point | ✓ | 1st-order | 3D | ✓ | ||||||||||||
| [74, Eq. (13)] | [74] | point | ✓ | 1st-order | 3D | ✓ | ||||||||||||
| [64, Eq. (12)] | [64] | point | ✓ | (25) | 1st-order | 3D | ✓ | |||||||||||
| Num | [61] | point | ✓ | 2nd-order | 3D | ✓ | ||||||||||||
| [60, Eq. (29)] | [60] | point | ✓ | (25) | 2nd-order | 3D | ✓ | |||||||||||
| Num | [81] | trans | ✓ | ✓ | ✓ | 3D | ✓ | |||||||||||
| (44) | [82] | trans/refl | ✓ | ✓ | 1st-order | 1D / 3D | ✓ | |||||||||||
| Num | [83] | refl | ✓ | 1st-order | 3D | ✓ | ||||||||||||
| Num | [84] | refl | ✓ | ✓ | ✓ | 3D | ✓ | |||||||||||
| Num | [85, 86, 87] | trans | ✓ | ✓ | ✓ | ✓ | 3D | ✓ | ||||||||||
| (50), (53) | This article | point | ✓ | uni | ✓ | 3D | ||||||||||||
| (54) | [37] | point | ✓ | uni | ✓ | ✓ | 3D | |||||||||||
| [47, Eq. (15)] | [47] | point | ✓ | uni | (25) | ✓ | 3D | ✓ | ||||||||||
| (– ‣ • ‣ III-D), (59), (60), (62) | [54] | trans | point / | ✓ | lam | ✓ | 3D | |||||||||||
| [5] | trans | ✓ | ✓ | 2nd-order | 3D | ✓ | ||||||||||||
| [80, Eq. (8)] | [80] | point | ✓ | ✓ | 2nd-order | 3D | ✓ | |||||||||||
| (– ‣ • ‣ III-D), (– ‣ • ‣ III-D) | [58] | point | ✓ | (29a) | ✓ | 3D | ✓ | |||||||||||
| [79, Eq. (12)] | [79] | point | ✓ | uni | 1st-order | 1D | ✓ | |||||||||||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Channel Modeling for Diffusive Molecular Communication – A Tutorial Review
Vahid Jamali, Arman Ahmadzadeh, Wayan Wicke,
Adam Noel, and Robert Schober
Friedrich-Alexander University of Erlangen-Nürnberg, Germany
University of Warwick, UK ⋆ Co-first authors.
Abstract
Molecular communication (MC) is a new communication engineering paradigm where molecules are employed as information carriers. MC systems are expected to enable new revolutionary applications such as sensing of target substances in biotechnology, smart drug delivery in medicine, and monitoring of oil pipelines or chemical reactors in industrial settings. As for any other kind of communication, simple yet sufficiently accurate channel models are needed for the design, analysis, and efficient operation of MC systems. In this paper, we provide a tutorial review on mathematical channel modeling for diffusive MC systems. The considered end-to-end MC channel models incorporate the effects of the release mechanism, the MC environment, and the reception mechanism on the observed information molecules. Thereby, the various existing models for the different components of an MC system are presented under a common framework and the underlying biological, chemical, and physical phenomena are discussed. Deterministic models characterizing the expected number of molecules observed at the receiver and statistical models characterizing the actual number of observed molecules are developed. In addition, we provide channel models for time-varying MC systems with moving transmitters and receivers, which are relevant for advanced applications such as smart drug delivery with mobile nanomachines. For complex scenarios, where simple MC channel models cannot be obtained from first principles, we investigate simulation-driven and experimentally-driven channel models. Finally, we provide a detailed discussion of potential challenges, open research problems, and future directions in channel modeling for diffusive MC systems.
Index Terms:
Molecular communications, diffusion, flow, reaction, end-to-end CIR, statistical model, simulation-driven models, and experiment-driven models.
I Introduction
Wireless communication networks have permeated throughout modern society, but existing systems are constrained by where conventional radio frequency technologies can be deployed. There are emerging applications where wireless communication could be a vital component, but where conventional implementations would be unsafe or impractical. An alternative approach that has received increasing attention within the communications research community over the last decade is molecular communication (MC), where molecules are employed as the information carriers111We note that, in this paper, we use the terms “molecule” and “particle” interchangeably.. MC was first proposed for the design of synthetic communication networks in [1]. The topic has received steady growth since the seminal survey on nanonetworks in [2], which are networks of devices with nanoscale functional components. An attractive feature of MC is its ubiquitous deployment in natural biochemical and biophysical systems, which lends credibility to its potential for biological applications such as targeting substances, smart drug delivery, and designing lab-on-a-chip systems [3]. Furthermore, MC could be deployed in industrial settings, including the monitoring of chemical reactors and nanoscale manufacturing, or for larger activities such as monitoring the emission of pollutants or the transport of oil [4].
Motivated by natural MC systems, several different mechanisms have been considered for MC in the literature including free diffusion [5, 6], gap junctions [7], molecular motors [8], and bacterial motors [9]; see Fig. 1. In particular, diffusion is referred to as the random movement of small particles suspended in a fluid medium as a result of their collisions with other particles in the fluid. Diffusion is one of the dominant propagation mechanisms in nature including communication inside cells and between cells, e.g., in quorum sensing among bacteria and in the synaptic cleft between neurons. Gap junctions enable another form of communication between cells where the molecules pass through small channels connecting the cytosol of neighboring cells. Calcium signaling is an example of this form of MC that is used by adjacent cells to regulate a large number of cellular processes including fertilization, proliferation, and death of mammalian cells [7, 10]. Molecular motors enable a form of active transportation of large signaling molecules via a special rail-like infrastructure, e.g., actin or microtubule filaments [11]. The motor moves along the rail by using repeated cycles of coordinated binding and unbinding of its legs to the rail. This type of MC is used for intracellular communication among organelles inside a cell [8]. Finally, bacterial motors enable another kind of active transport where the bacteria can pick up large signaling molecules, e.g., deoxyribonucleic acid (DNA), and move in a specific direction, e.g., due to a food concentration gradient, using their tiny propellers (known as flagella) [9].
Diffusion-based MC, sometimes in combination with advection and chemical reaction networks (CRNs), has been the prevalent approach considered in the literature thus far; see [3, Table 4]. The main advantages of diffusion-based MC include that, unlike gap junction-based MC, special infrastructure is not needed, and unlike motor-based MC, external energy for propagation of the signaling molecules is not required. Moreover, the simplicity of diffusion makes it an attractive propagation scheme, especially for ad hoc networks of devices with limited computational resources. Hence, in this tutorial, we focus on diffusion-based MC, where we also consider environments with advection and CRNs.
I-A Scope
A fundamental aspect in the analysis of any communication system is characterizing and understanding the physical layer, i.e., the channel between the transmitter and the receiver. Generally, we rely on a channel model that is simple and yet sufficiently accurate for us to design, analyze, and operate a system that communicates over the channel. A complete view of a diffusion-based channel includes the release of molecules from a transmitter, their propagation in the fluid environment, and the reception mechanism at the receiver. While there is rich historical literature on the physics of diffusion and characterizing expected diffusion in environments of different shapes, cf. e.g., [12, 13], the communications research community has expanded these models to account for the behavior of the end-to-end system, for the inclusion of non-diffusive phenomena that play important roles in biophysical systems, and for the statistics of molecular behavior.
Recent surveys, in particular [14, 3], have provided excellent qualitative summaries of diffusive MC and included some of the most common channel models available at the time of writing. A more complete mathematical treatment of diffusion-based modeling of MC can be found in [15]. However, there have been significant advances in channel modeling in the years since the publication of [15], and also since the most recent major survey of models in [3]. In particular, non-diffusive effects that can be coupled with diffusion, such as advective flow and chemical reaction kinetics, have been integrated in many channel models to make them more practical and more accurate.
Due to the rapid growth in channel models, it has become difficult for an interested researcher to enter the MC field and become familiar with the state-of-the-art in diffusion-based channel modeling. It has also become more challenging for practitioners in this field to stay up to date. The aim of this tutorial review is to satisfy both audiences. We present a detailed and rigorous mathematical development of diffusive MC channel models. We seek to provide a useful comprehensive reference on channel models that is both approachable for an audience that is new to the field and also convenient for active practitioners to assess and select a model. To do so, we begin with a review of the underlying fundamental laws that govern diffusive MC channels and show how they are used in the literature to derive the channel impulse responses (CIRs) of different MC systems. In addition, we present different deterministic and statistical models developed for the observation signal at the receiver. We also discuss the complementary roles of simulations and physical experiments to both support analytical modeling and to provide data-driven models when simple analytical models cannot be readily obtained.
I-B Contributions
In this tutorial review, we make the following contributions:
By taking a mathematically rigorous approach, we first provide a tutorial on the underlying phenomena from biology, chemistry, and physics, and their effect on the components of MC systems. Specifically, we start with Fick’s laws of diffusion and build towards the general advection-reaction-diffusion equation. We discuss the common assumptions and special cases that enable the general equation to be solved for the CIR in closed form. 2. 2.
We review the major end-to-end channel models in the diffusive MC literature including the effects of release mechanisms, the physical channel, and reception mechanisms. In particular, we include the relevant classical models from the physical sciences literature, as well as a comprehensive presentation of the models that have been developed and the equations that have been derived within the communications engineering community over the last few years. 3. 3.
We present a unified definition for the observed signal at a receiver. The unified definition encompasses both timing and counting receivers and helps to better understand the basic assumptions that have been made to arrive at the well-known signal models used in the MC literature and how they relate to each other. Then, we focus on counting receivers and derive signal models relevant for different time scales. We further generalize these models to account for interfering noise molecules and inter-symbol interference (ISI). Finally, we study the correlation between the received signals observed at different time scales. 4. 4.
We discuss the integral role of simulations and experiments, in particular to gain insight from a data-driven model when closed-form solutions for the CIR are not readily available. We also describe how to implement simple stochastic simulations as well as how to derive an example data-driven model based on experimental data.
For clarity of presentation, the focus of the channel models presented in this work is on a single communication link between one transmitter and one receiver. Many of the envisioned applications of diffusive MC systems will depend on many links within a network of devices. While there have been a number of relevant contributions that consider the propagation of signals over multiple links, such as via relaying and cooperative detection (cf. e.g., [16, 17, 18, 19, 20]), these models can often be decomposed into a superposition of individual links. However, it is important to note that we cannot always directly apply single link analysis to multiple-link systems. In particular, we must exercise caution when there are multiple non-transparent devices (such as reactive receivers) in the system that molecules can collide or react with. The presence of one such device will impact the signal received at any receiver. Some special cases can still be treated in closed form, such as having two absorbing receivers placed on either side of a transmitter in [21], but otherwise we need to resort to data-driven models, such as the simulation of multiple absorbing receivers in [22].
I-C Organization
The rest of this tutorial review is organized as follows, and also summarized in Table I to show how the content of Sections II-V is connected. We review the fundamental physical principles that govern diffusion-based MC systems in Section II. In particular, we model diffusion, advection, and chemical reactions, which lead to a general advection-reaction-diffusion partial differential equation (PDE) to describe the spatio-temporal variation in molecule concentrations.
In Section III, we discuss the components of MC systems and their effect on the end-to-end CIR. Our definition of the end-to-end channel includes the physical and chemical properties of the transmitter and receiver, as well as the fluid medium in which they are located. A table to summarize the reviewed CIRs is also provided.
In Section IV, we present a unified definition for the diffusive signal observed at the receiver. We focus on counting receivers to derive deterministic and statistical signal models that are valid for different time scales. We also consider the impact of interfering noise, including the interference caused by repeated transmissions at the transmitter, and the correlation among received signals sampled at different time scales.
We discuss simulation- and experimental-driven models in Section V. We describe the different physical scales for simulating diffusion-based systems, including existing simulation platforms for each scale, and discuss how to implement simple stochastic simulations. We review a selection of experimental platforms and discuss features that cannot be readily captured via modeling or simulation. Thereby, we propose to employ physically-motivated parametric models and neural networks whose parameters are found using experimental data. One experimental system is presented as a case study for data-driven modeling.
We end this tutorial review with a discussion of future work and open challenges in Section VI before presenting our conclusions in Section VII.
II Fundamental Governing Physical Principles in MC Systems
In this section, we review the fundamental laws that govern the propagation of molecules. In particular, we mathematically model the impact of diffusion, advection, and reaction on the spatio-temporal distribution of molecules. This modeling is essential for the development of channel models. A solid understanding of these phenomena is needed to develop intuition for molecule propagation in diffusive MC systems. Furthermore, in Section III, we will use the mathematical tools introduced in this section for the derivation of the CIR for several different diffusive MC systems.
II-A Free Diffusion
Molecules in a fluid environment, such as a liquid or a gas, are affected by thermal vibrations and collisions with other molecules. The resulting movement of the molecules is a purely random without any preferred direction and is referred to as a random walk or Brownian motion. Let denote a vector specifying the position of the -th molecule in three-dimensional (3D) Cartesian coordinates at time . Thereby, the random walk is modeled by [23, Eqs. (1.3) and (1.21)]
[TABLE]
where is the time step size and in [m2s*-1*] is the diffusion coefficient of the -th molecule. Moreover, denotes a multivariate Gaussian random variable (RV) with mean vector and covariance matrix , represents a vector whose elements are all zeros, and is the identity matrix. The diffusion coefficient determines how fast the molecule moves. The larger the diffusion coefficient, the larger the average displacement of the molecule in a given time interval. The value of the diffusion coefficient depends on the environment as well as the shape and the size of the particle. For large spherical particles, the diffusion coefficient can be determined based on the so-called Einstein relation [24, Eq. (4.15)]
[TABLE]
where is the Boltzmann constant, is the temperature in kelvin, is the (dynamic) viscosity of the fluid ( for water at ), and is the radius of the particle. Note that larger particles have a smaller diffusion coefficient and are hence less affected by diffusion.
Remark 1
Besides the ideal free diffusion with constant diffusion coefficient discussed above, there are also other types of diffusion. For instance, in contrast to the typical free diffusion where the mean squared displacement (MSD) is linearly proportional to time, i.e., , in anomalous diffusion, the MSD follows a nonlinear relation, i.e., where . Sub-diffusion occurs when and can be used to model diffusion inside biological cells where the presence of the organelles does not allow ideal free diffusion to take place [25]. Super-diffusion occurs when and can be used to model active cellular transport processes [26]. Moreover, in (1), we assumed the diffusion coefficient to be constant. However, the diffusion coefficient may depend on the local concentration of the molecules [27]. For the constant diffusion coefficient assumption to hold, the temperature and viscosity of the environment are assumed to be uniform and constant and all solute molecules (dissolved molecules) are assumed to be locally dilute everywhere, i.e., the number of solute molecules is sufficiently small everywhere. These assumptions allow us to ignore potential collisions between solute molecules such that the diffusion coefficient does not vary with the local concentration [27, 28]. We refer the readers to [29] for the study of diffusion with non-constant diffusion coefficients. Another example of a complex diffusion process is the diffusion of protons in water. Here, the movement of the protons is a combination of ideal free diffusion and the so-called structural diffusion where protons hop from one water molecule to the next. Nevertheless, it has been shown in [30] that proton transport can be well approximated by free diffusion with an effective diffusion coefficient.
We let denote the concentration of the solute molecules, i.e., the average number of solute molecules per unit volume, at coordinate and time . The random movement of molecules due to diffusion, described by (1), leads to variation in across time and space that obeys Fick’s second law of diffusion222Fick’s first law of diffusion relates the diffusive flux, denoted by , to the concentration as .
[TABLE]
where is the Laplace operator, e.g., in Cartesian coordinates. The PDE in (3) can be solved for simple initial conditions (ICs) and simple boundary conditions (BCs). In the following, we consider a simple example, namely diffusion in an unbounded 3D environment with an impulsive point release, which has been the most widely studied case in the MC literature due to its simplicity [31, 32, 3, 33, 34, 35, 36, 37, 38, 39]. In the remainder of this paper, we denote the solutions of the considered PDEs by .
Example 1** (Diffusion in an Unbounded 3D Environment with Impulsive Point Release)**
Consider a 3D diffusion process with instantaneous release of solute molecules from at time . To obtain , we have to solve (3) with the following initial and boundary conditions
[TABLE]
where , and is the Dirac delta function. Solving (3) with and yields [23, Eq. (2.8)]
[TABLE]
In Fig. 2, the molecule concentration [molecules/m3] in (6) is plotted versus time [s] at distance with nm for an initial release of molecules with m2/s from the origin at time . From Fig. 2, we observe that first increases with time, which is due to the non-zero propagation time that the molecules need to reach , before it decreases since the molecules diffuse away. Moreover, as distance increases, the peak of the concentration decreases since the molecules are spread over a larger volume. Furthermore, the time when the concentration peak occurs, denoted by , increases with distance.
The assumption of an unbounded environment is accurate when the actual boundaries of the system are far away from the region of interest (i.e., from transmitter and receiver), such that the impact of the boundaries on the diffusing molecules can be neglected. In the following, we present an example where the effect of the boundaries cannot be neglected.
Example 2** (Diffusion in an Unbounded Straight Duct with Impulsive Release from Cross-Section)**
We assume a straight duct333A duct is a pipe, tube, or channel which carries a liquid or gas. channel with circular cross-section and for convenience, we employ cylindrical coordinates, i.e., with , where denotes the radius of the circular cross-section of the duct. We assume that, at the time of release, , the molecules are uniformly distributed across the cross-section at . Therefore, we have the following initial and boundary conditions
[TABLE]
where enforces the reflection of the molecules at the wall, i.e., a fully reflective wall is assumed. Solving (3) with , , and yields [40]
[TABLE]
As can be seen from (2), does not depend on variables and due to the symmetry of the initial condition and the environment with respect to and . This model can be used to characterize the propagation of molecules in blood vessels as is necessary for drug delivery applications of MC in the cardiovascular system [41, 42, 43, 44, 45].
II-B Advection
Besides diffusion, advection is another fundamental mechanism for solute particle transport in a fluid environment. In the following, we first specify how mass transport by advection affects a single solute particle. Subsequently, we distinguish between two types of advection, namely drift and fluid flow, and give the particle velocity vector for some example cases. Moreover, we present the advection equation which describes the change in molecule concentration due to advection. Finally, we introduce the advection-diffusion equation, which captures the joint impact of diffusion and advection, and characterize the relative importance of diffusion and advection.
In general, transport by advection can be described by a velocity vector which generally may depend on position and time . When considering the movement of the -th particle at position due to advection, its position at time can be modeled by
[TABLE]
where should be small enough such that the velocity vector is constant between and . Next, we discuss what may cause the velocity vector and what form it may take.
II-B1 Velocity Vector Field
Transport by advection can be mediated by different physical mechanisms which we categorize as force-induced drift and bulk flow [46, 47].
Force-Induced Drift: Advection can be caused by external forces acting on the particles but not on the fluid containing the particles. An external force can be modeled by force vector which describes the force on a particle at position at time . These external forces can be electrical, e.g., if the particles are ions, or magnetic, e.g., if the particles are magnetic nanoparticles, or gravitational, e.g., if the particles have sufficient mass, or a combination of forces [47, 48]. When the force is not too large, the velocity vector can be determined from the corresponding force by Stokes’ law via [49, Eq. (2.65)]
[TABLE]
where is a proportionality constant referred to as the friction coefficient. The friction coefficient can be related to the diffusion coefficient via . In other words, using (2), we obtain . Force may vary with time (e.g., for ions if the electric field changes over time) and space (e.g., for magnetic nanoparticles, the magnetic force generally decreases rapidly with increasing distance from the magnet) [47, 48].
Bulk Flow: If the particle movement is induced by the movement of the fluid, then the resulting transport by advection is referred to as flow. Flow can be encountered in many MC environments such as blood vessels and microfluidic channels [50]. In MC, we typically have dilute particle suspensions, where the flow velocity is independent from the particle concentration. Thereby, the velocity vector will depend on space if there are boundaries or obstacles in the environment, e.g., in a duct, the flow velocity is typically largest in the center and smallest at the boundary where the fluid is subject to friction. The flow may also depend on time, e.g., in a blood vessel the flow is generated by the periodic contractions of the heart.
Remark 2
Although both flow and external force cause the particles to drift, which can be modeled by (11), they may require quite different considerations. For instance, any object in the environment influences the velocity vector caused by bulk flow since the flow may not be able to penetrate the object and has to go around the object. On the other hand, the drift velocity vector caused by an external force is not necessarily influenced by objects in the environment.
Flow can be also categorized into two classes, namely turbulent and laminar flow. In particular, when the variations of the flow velocity, over space and/or time, are stochastic, e.g., due to rough surfaces and high flow velocities [51], we refer to the flow as turbulent. If the flow is not turbulent, it is referred to as laminar. For flow in a bounded environment of effective length and with an effective velocity of , the Reynolds number can be used as a criterion for predicting laminar or turbulent flow and is given by [51, Eq. (1.24)]
[TABLE]
where is the kinematic viscosity [] of the fluid444Kinematic viscosity is related to (dynamic) viscosity according to where [] is the fluid density.. For example, for flow in a straight pipe with circular cross-section of radius , the flow can be assumed to be laminar and turbulent for and , respectively, where [51]. For microfluidic settings, typically and hence laminar flow can be assumed [49]. For most blood vessels, holds and hence the blood flow is typically laminar [52, 53]. Only in large arteries such as the aorta (the largest artery in the human body), the Reynolds number can be in the range and thereby blood flow exhibits turbulent behavior [53].
Generally, for a given environment, the flow velocity vector as a function of space and time can be determined by solving the so-called Navier-Stokes equation with appropriate boundary conditions, see e.g. [49, Eq. (5.22)]. Let us review two special cases of , which have been widely studied in the MC literature [37, 3, 46, 54, 55] and are also considered in Section III.
Example 3** (Uniform and/or Constant Advection)**
For uniform advection, the velocity vector is constant across space but can be time-dependent, i.e., [55]. For advection by flow in an unbounded environment, uniform flow solves the Navier-Stokes equation and hence can be physically plausible. Moreover, for advection by drift, uniform drift is applicable when the corresponding force vector does not depend on space, see (12). As a special case, the velocity vector may be constant across both space and time, i.e., . Due to its simplicity, advection with constant velocity is the most widely-studied advection model in the MC literature [37, 3, 46].
Example 4** (Steady Flow in an Infinite Straight Duct with Circular Cross-Section)**
For this example, we concentrate on advection by fluid flow because force-induced drift is completely specified by (12). In particular, in this case, the flow velocity vector in cylindrical coordinates can be obtained as [51, Eq. (4.134)]
[TABLE]
where is the center velocity. The flow described in (14) is laminar and can be interpreted as follows. For a given , the flow velocity in (14) is constant but increases from the boundary where towards the center where , i.e., for each , we can think of a circular layer within the duct that slides along its neighboring layers with a constant velocity. The velocity vector in (14) is known as the Poiseuille flow profile and is a common model for the flow in blood capillaries [54].
While for other environments and boundary conditions the velocity vector can still in principle be obtained from the Navier-Stokes equation, it is often not possible to do so analytically. In these cases, the Navier-Stokes equation can be solved by numerical algorithms that are well-established in computational fluid dynamics [51].
II-B2 Advection Equation
Given , the change in concentration with respect to time due to advective transport is modeled by the following PDE, which is referred to as the advection equation or continuity equation [49, Eq. (4.14)]
[TABLE]
where denotes the gradient operator and denotes the inner product of two vectors and . In general, (15) cannot be readily solved for a given velocity vector and numerical methods have to be employed [56]. Nevertheless, for the velocity vectors in Examples 3 and 4, (15) can be solved as shown in the following.
Example 5
*Assuming initial condition at , the advection equation (15) has the following solution for *
[TABLE]
We note that while the solutions in (5) appear similar, they are actually fundamentally different. In particular, for constant uniform flow and uniform flow (space-independent flow profiles), the initial concentration is simply translated to a different position without changing its shape. However, for Poiseuille flow (space-dependent flow profile), the concentration generally spreads in space over time depending on the initial concentration.
II-B3 Advection-Diffusion Equation
In many application scenarios, such as drug delivery via the capillary networks [41, 42, 43, 44, 45], advection and diffusion are both present in the MC environment. Thereby, the combined effect of both advection and diffusion is characterized by the following PDE known as the advection-diffusion equation
[TABLE]
Similar to diffusion equation (3), (17) cannot be solved analytically for general velocity vectors and general boundary and initial conditions. In the following, we first provide the solution of (17) for constant uniform flow in an unbounded environment. Subsequently, we quantify the relative impact of diffusion over advection by introducing the notions of Péclet number and dispersion factor.
Example 6
Consider an unbounded 3D environment with instantaneous release of solute molecules at at time . Solving (17) with initial condition in (4), boundary condition in (5), and constant uniform velocity vector yields [37, Eq. (18)]
[TABLE]
In Fig. 3, we show molecule concentration [molecules/m3] in (18) versus time [s] at nm for initial release of molecules with m2/s from at , and flow velocity vector with m/s. From Fig. 3, we observe that as the flow velocity increases, the concentration peak increases and decreases. This is mainly due to the fact that the flow is in the same direction as the point where the concentration is measured, i.e., parallel flow is considered. Parallel flow can considerably enhance the coverage of a diffusion-based MC system, e.g., in blood vessels. Moreover, by increasing , the tail of over time is decreased, which is useful for ISI reduction in MC systems [57, 37].
Relative Importance of Advection over Diffusion for Molecule Transport: Advection and diffusion can both displace and transport molecules, albeit in different ways. An important question is under what conditions is one more effective than the other. The Péclet number, denoted by , can be used to answer this question. Let us assume a velocity vector with strength and transport over a distance which is referred to as the characteristic length. The Péclet number quantifies the ratio of time required for particles to be transported by diffusion over distance (which is proportional to ) with the time required for particles to be transported by advection over distance (given by ). This ratio is given by [49, Eq. (4.44)]
[TABLE]
Note that is a dimensionless number. If holds, diffusion dominates advection and the spreading of molecules is almost isotropic despite a weak biased transport in the direction of the flow. In this case, the solution of the diffusion equation (3) provides an accurate estimate of the molecule concentration. On the other hand, if holds, advection dominates diffusion and is the main cause for molecule transport. In this case, the advection equation (15) can be solved to obtain an accurate estimate of the molecule concentration. Finally, for , molecule transport is sensitive to both diffusion and advection and the advection-diffusion equation in (17) should be solved.
Relative Importance of Advection over Diffusion for Dispersion: Let us consider a straight duct with a circular cross-section, see Examples 2 and 4, where advection is the main transport mechanism along the duct. In other words, holds where denotes the Péclet number for transport along the -axis, is the effective flow velocity in the duct (see (14)), and is the desired transport length along the -axis. In this case, we are interested in studying the dispersion (spatial spreading) of individual particles across the cross-section over the time when transport along the -axis occurs. In particular, one may distinguish between the following two extreme regimes, namely the non-dispersive and dispersive regimes:
i) Non-dispersive regime: Here, particles do not considerably diffuse across the cross-section while being transported by advection. Therefore, each particle is simply transported along the -axis by advection with a velocity strength that depends on the radial position of the particle, , according to (14). We note that although the dispersion of individual particles is negligible in this regime, the shape of the concentration profile varies over time since the flow has a different effect at different radial positions, i.e., particles closer to the center of the duct travel faster.
ii) Dispersive regime: In the dispersive regime, particles fully diffuse across the cross-section while also being transported along the -axis by advection. In addition to the dispersion across the cross-section, there is also dispersion along the -axis, due to the combined impact of diffusion and advection with space-dependent flow profile (14).
In the following, we mathematically quantify the dispersive and non-dispersive regimes in terms of system parameters, i.e., , , , and . We choose the characteristic length as the distance over which the velocity vector changes (usually a fraction of ). Moreover, we define as the corresponding dimensionless normalized distance with respect to characteristic distance . Then, we can compare the characteristic time required for particles to be transported by advection over distance (given by ) with the time required for diffusion over distance (which is proportional to ). To compare these two time scales, we can define a dispersion factor as
[TABLE]
Here, signifies that there is not enough time for particles to diffuse across the cross-section while being transported by advection over distance , i.e., we are in the non-dispersive regime. On the other hand, for , diffusion causes considerable dispersion across the cross-section, which in turn causes significant dispersion along the -axis due to space-dependent flow velocity (14), i.e., we are in the dispersive regime. In other words, in terms of the Péclet number , we have non-dispersive and dispersive regimes if and hold, respectively.
Fig. 4 illustrates different dispersion regimes for a 3D straight duct. For clarity of presentation, we only show those particles where the -component of their position lies in interval . As can be seen from Fig. 4, for , the positions of the particles simply follow the velocity profile in (14) whereas for , particles are significantly dispersed in the environment.
II-C Chemical Reactions
Another important phenomenon affecting the propagation of signaling molecules in diffusive MC systems is chemical reactions. On the one hand, chemical reactions may occur naturally in MC environments and their impact must be taken into account for communication design. On the other hand, chemical reactions have been exploited in the MC literature to achieve certain objectives, such as ISI reduction [58, 59, 36, 37] and ligand-based reception modeling [60, 61]. Therefore, in the following, we first review general chemical reactions, the corresponding reaction equations, and examples of reactions widely considered in the MC literature. Subsequently, we study the joint impact of all three phenomena discussed in this section, namely diffusion, advection, and reaction, on the propagation of the molecules and solve the corresponding advection-reaction-diffusion equation for a simple example.
II-C1 Reaction Equation
Consider a general reaction of the form [62, Eq. (13)]
[TABLE]
where are reactant molecules, is the set of reactant molecules, are product molecules, is the set of product molecules, and are non-negative integers, and is the reaction rate constant. Let and denote the concentration of type- and type- molecules at coordinate and time , respectively. Reactions locally change the concentration of particles over time which is described by the following PDEs, known as reaction equations
[TABLE]
where denotes the reaction rate function, which depends on the reaction rate constant and the concentrations of the reactant molecules. The reaction rate function has the following general form, known as the rate law [63, Eq. (9.2)]
[TABLE]
where is the order of the reaction with respect to type- reactant molecules and typically takes an integer value (but in principle may also assume real values). The overall reaction order is defined as [63, 40]. Note that the units of reaction rate function and reaction rate constant are and \frac{1}{\text{s}}\big{(}\frac{\text{molecule}}{\text{m}^{3}}\big{)}^{1-\sum_{I\in\mathcal{I}}\varepsilon_{I}}, respectively.
In the following, we present three important classes of reactions, namely unimolecular degradation, bimolecular reactions, and enzymatic reactions, which can all play important roles in MC systems [37, 58, 64, 65, 66]. In particular, degradation is a natural characteristic of some types of molecules and its effect has to be accounted for in communication design, see Section III-D and [37, 64]. Bimolecular reactions can be used to analayze ligand-receptor binding [60, 61] and reactive signaling [59, 67]. In addition, enzymatic reactions have been studied in the MC literature for the purpose of ISI reduction [58, 66].
Example 7** (Unimolecular Degradation)**
This reaction is used to describe the degradation of a desired type of molecule, e.g., type , into a new type of molecule, denoted by , which is of no interest for the considered communications. In fact, unimolecular degradation is often used as a first-order approximation of more complex reactions such as bimolecular and enzymatic reactions, see Examples 8 and 9. Unimolecular degradation is modeled by [63, Ch. 9]
[TABLE]
where [\frac{1}{\text{s}}\big{(}\frac{\text{molecule}}{\text{m}^{3}}\big{)}^{1-\varepsilon_{A}}] is the reaction rate constant, is the reaction rate function, and is the reaction order. In the MC literature, first-order reactions are used to model degradation, i.e., [37, 64]. However, depending on the speed of reaction, higher and lower order reactions may be relevant, e.g., zero-order () or second-order (Type-I) () reactions [63, Ch. 9]. Assuming an initial condition at , (22a) has the following solution for
[TABLE]
where . Note that the speed of molecule concentration decay is hyperbolic for second-order degradations, which is faster than the exponential decay for first-order degradations, which in turn is faster than the linear decay for zero-order degradations. Nevertheless, for sufficiently large , for second-order degradations is larger than that for first-order degradations, whereas , holds for zero-order degradations.
Example 8** (Bimolecular Reactions)**
Some reactions may involve the interaction of two reactant chemical species, e.g., and , to produce product molecule(s), e.g., . For instance, in [60], the activation of ligand receptors via signaling molecules was modeled by a second-order bimolecular reaction. Moreover, in [59] and [67], acids and bases were used as reactive signaling molecules to reduce ISI. Acids and bases cancel each other out to produce salt and water. This process is modeled by a second-order bimolecular reaction. In particular, the second-order (Type-II) bimolecular reaction is given by [68]
[TABLE]
where is the forward reaction rate constant [], [] is the backward reaction rate constant, and is the reaction rate function. The PDEs corresponding to (27) are nonlinear and challenging to solve. However, after introducing some approximations, in Section III, we use (27) to derive the CIRs of MC systems affected by bimolecular reactions. Moreover, let us assume and that the concentration of type- molecules is sufficiently large such that the reaction in (27) does not considerably change over time, i.e., . In this case, the bimolecular reaction in (27) can be approximated by the first-order unimolecular reaction in (25) with [60].
Example 9** (Enzymatic Reactions)**
For typical scenarios, the speed of natural degradation might be too slow compared to the desired time scale of communication. In this case, enzymes can be used to accelerate the reaction process. Enzymes, denoted by , are specific proteins that bind to the desired molecule (also referred to as the substrate), and lower the activation energy needed for a reaction to occur. Enzymatic degradations are modeled by the following reactions [58, Eq. (1)]
[TABLE]
where is an intermediate chemical species and is the product molecule. Moreover, [], [], and [] denote the reaction rate constants of the forward, backward555The forward and backward reaction rate constants are also referred to as binding and unbinding reaction rate constants, respectively., and degradation reactions, respectively. As can be seen from (28), the enzyme molecules are not consumed in the reaction process. The following set of PDEs, known as Michaelis-Menten kinetics, describe the evolution of the concentrations of the participating molecules
[TABLE]
Solving the above system of coupled and nonlinear PDEs is challenging. Let us consider very fast degradation reactions, i.e., , slow backward reactions, i.e., , and that the concentration of enzyme molecules is much larger than the concentration of type- molecules. In this case, the formation of intermediate molecules does not last long and hence, we obtain . In [58], it was shown that under the aforementioned assumptions, the enzymatic reaction in (28) can be approximated by the first-order unimolecular reaction in (25) with reaction rate constant .
II-C2 Advection-Reaction-Diffusion Equation
Next, we consider the joint effects of diffusion, drift, and reactions. For simplicity, we focus on a single molecule type and drop the corresponding subscript. In this case, the general advection-reaction-diffusion equation is given by the following PDE [23, 69]
[TABLE]
where and hold if the considered molecule is the product and the reactant of the reaction, respectively. Solving (32) for general initial and boundary conditions is again difficult for most practical MC environments. Hence, in the following, we make some simplifying assumptions that enable us to solve (32) in closed form for one example scenario [58].
Example 10
*Let us assume the impulsive release of molecules at time by a point source located at into an unbounded 3D environment, i.e., initial condition in (4) and boundary condition in (5) hold. Moreover, we assume uniform flow and the first-order degradation reaction in (25), i.e., and . Based on these assumptions, (32) has the following closed-form solution [70, 71] *
[TABLE]
In Fig. 5, the molecule concentration [molecules/m3] is shown versus time [s] at nm for an initial release of molecules from and at , m2/s, flow velocity m/s, and 1/s. This figure shows that as the degradation rate constant increases, the concentration peak decreases, which is not desirable for an MC system, in general. However, the tail of the concentration for large fades away much faster for larger degradation rates, which was exploited for ISI reduction in [58].
III Component Modeling
In this section, we review the existing component models for the transmitter, receiver, and physical channel of diffusive MC systems. To this end, in Section III-A, we first define the end-to-end CIR of single-link diffusive MC systems, and discuss the relevant mechanisms of each component and their impact on the end-to-end CIR. We use the CIR to characterize the components of MC systems, since the impulse response fully characterizes the behaviour of linear systems, and linearity is commonly assumed in the MC literature. Subsequently, in Sections III-B, III-C, and III-D, we review the existing models that have been developed by taking into account the impact of the receiver, transmitter, and physical channel on the end-to-end CIR, respectively. Finally, in Section III-E, we provide a summary table of all reviewed end-to-end CIR models.
III-A Channel Impulse Response
In this subsection, we first briefly discuss the relevant mechanisms that characterize the functionalities of the transmitter and receiver, and the phenomena and impairments that occur in the physical channel of diffusive MC systems. Then, we provide a formal definition of what we refer to as the end-to-end channel of diffusive MC systems and we show how the CIR corresponding to the end-to-end channel can be obtained using the tools introduced in Section II.
Similar to traditional communication systems, the end-to-end chain of diffusive MC systems consists of three components, namely the transmitter, the physical channel, and the receiver; see Fig. 6. Each of these components has unique features and responsibilities, which are outlined below; see also Fig. 7.
- •
Transmitter: The transmitter is responsible for the encoding and modulation of information bits. In MC, the information is typically encoded in the number, type, or time of release of signaling molecules. Furthermore, the transmitter has to generate the signaling molecules, (e.g. by CRNs inside the transmitter), store the signaling molecules, (e.g. in vesicles), and control their release into the physical channel.
- •
Physical Channel: The physical channel is the environment in which the signaling molecules move and propagate once they leave the transmitter. In diffusive MC systems, the movement of signaling molecules, at its most basic level, is described by the diffusion process. However, during the course of diffusion, the random walk of signaling molecules may be affected by several other factors and noise sources such as advection, CRNs degrading the signaling molecules, environment geometry, and obstacles inside the physical channel, see Section II.
- •
Receiver: Signaling particles that reach the vicinity of the receiver can be observed and processed by the receiver to extract the information that is necessary for performing detection and decoding. The reception mechanism of the receiver may include the following functionalities, depending on its structure: i) external sensory units for detecting the presence of signaling molecules, membrane receptors of cells in nature, or sensing component(s) of macro-scale receivers such as the alcohol sensor in [57] and the magnetic coils of the susceptometer in [72]; ii) internal relaying and interface components to convey and convert the measurements of the sensory unit into quantitites suitable for detection and decoding of the information bits. For instance, in nature, this task is performed by the CRNs inside cells, which are referred to as downstream signaling pathways [10]. Downstream signaling pathways may be driven by activated receptors or directly by signaling molecules that passively enter the cells.
In the following, we formally define the end-to-end channel to study the reviewed CIR models in a unified manner.
Definition 1** (End-to-end Channel)**
We define the end-to-end channel as the effective channel that not only includes the physical channel but also the impact of the physical and chemical properties of the transmitter and receiver, including the effects of signaling molecule generation, release mechanisms, sensory units, and internal receiver components.
Note that our definition of the end-to-end channel does not include the coding, modulation, detection, and decoding operations that the transmitter and receiver may perform; see also Fig. 6. This definition of the end-to-end channel is analogous to that in traditional wireless communication systems, where the antennas, power amplifiers, and filters of the transmitter and receiver are also included in the model for the wireless end-to-end channel. The input to the end-to-end channel is the signal representing the modulated information symbol, which we also refer to as the stimulation signal. The stimulation signal can be an electrical (voltage or current), magnetic, mechanical, optical, chemical, or temperature signal. The output of the end-to-end channel is referred to as the observed signal and should be in a form that is suitable for the subsequent detection and decoding operations. Depending on the structure of the receiver, the observed signal can be either a number of output molecules or any secondary signal derived from the output molecules. In particular, output molecules may represent: i) signaling molecules that can passively enter the receiver; ii) absorbed molecules that hit the receiver surface; or iii) activated receptors. Furthermore, the secondary signal derived from output molecules may be an electrical signal, e.g., the output voltage or output current of the alcohol sensor in [57]. In the following, for the definition of the CIR of the end-to-end channel, we emphasize that we consider the number of the output molecules as the observed signal, as it is commonly assumed in the MC literature, although our definition can be easily extended to other forms of the observed signal.
Definition 2** (Channel Impulse Response)**
We define the CIR of the end-to-end channel, denoted by , as the probability of observation of one output molecule at time at the receiver when the transmitter is stimulated in an impulsive manner at time .
We note that defining the CIR as a probability has several advantages. In particular, it facilitates the definition of the received signal in Section IV. There, we propose a general received signal model that takes into account both the arrival time and the numbers of observed output molecules. As will be shown in Section IV, both of these quantities can be readily obtained from the probability of observation of one output molecule.
In our definition of the CIR, the quantitative meaning of the term observation depends on the type of receiver and is defined for each considered receiver model in detail in the next subsection, e.g., for passive receivers the observed signal is defined as the number of signaling molecules inside the receiver, while for reactive receivers it is defined as the number of activated receptor molecules. Furthermore, we assume that the transmitter stimulation is an impulsive input that either controls the opening and closing of the signaling molecule reservoir or drives the CRNs inside the transmitter responsible for the generation of the signaling molecules.
In this section, we assume that the parameters of the considered MC system are constant, i.e., the end-to-end CIR is time-invariant. In the following, we refer to the signaling molecules as molecules. The following phenomena may affect the propagation of the molecules, and as a result, :
Particle generation: Generation of the molecules is performed, e.g., by the CRNs inside the transmitter. 2. 2.
Release mechanism: The release mechanism can be chemical, electrical, or mechanical and controls the release of the molecules into the physical channel. 3. 3.
Diffusion: Diffusion refers to the propagation of molecules by Brownian motion. 4. 4.
Degradation and production: CRNs may degrade or produce molecules in the physical channel. 5. 5.
Advection: Advection may affect the transportation of the molecules in the physical channel. 6. 6.
Geometry: Potentially, the geometry of the individual components of the end-to-end channel can influence the propagation of signaling molecules. 7. 7.
Receptor kinetics: Receptor kinetics affect the interaction of the molecules with the receptors of the sensory unit at the receiver. 8. 8.
Signaling pathways: The signaling pathways transducing the observed molecules into secondary signal affect the received signal.
In order to obtain for a specific MC system, one has to solve the advection-reaction-diffusion equation (32) or a simplified version thereof, depending on the MC system under consideration, with the appropriate initial and boundary conditions. The initial conditions of the system capture the initial states of the CRNs, the time of production of the molecules, and the location of the produced molecules. The boundary conditions capture the physical and chemical properties of the components of the end-to-end channel. As discussed in the previous section, the solution to this system of PDEs does not exist in closed-form for many environments. However, as we will see in the remainder of this section, in the MC literature, different approximations have been developed to arrive at approximate yet meaningful solutions for that can still capture the main effects and phenomena of the end-to-end channel. These approximate models focus on one of the components of the MC system and make simplifying assumptions about the other two. Accordingly, we will consider such receiver, transmitter, and channel centric models in the following three subsections.
III-B Receiver Models
In this section, we review some of the existing end-to-end CIR models that focus particularly on the properties of the receiver, while simplifying assumptions for the transmitter and MC environment are made. The reception mechanism of the receiver can be categorized into two classes: i) passive reception, where the receiver does not impede the movement of signaling molecules; and ii) active reception, where the receiver may affect the movement of signaling molecules either by their absorption on its surface, or by chemically reacting with them via receptors (and thereby forming ligand-receptor complexes) embedded in the receiver surface. For active reception, both mechanisms can be described by a form of chemical reaction. Moreover, the received signaling molecules may be converted via signaling pathways into secondary molecules, which can later be used for detection or decoding of the information. In nature, cells have diverse types of signaling pathways, each of which is responsible for relaying a particular type of measurement taken in the extracellular space to the organelles in the cytosol, which ultimately causes a response by the cell. For more information on the signaling pathways in natural cells, we refer the interested reader to [10].
For the CIR models considered in the following, we adopt rather simple models for the transmitter and the physical channel. Specifically, we assume that the transmitter is a point that releases one molecule instantaneously upon stimulation at time at location , where denotes the location of the center of the transmitter; see Section III-C for more details on the point transmitter model. In other words, a point transmitter implicitly implies that upon stimulation, the signaling molecule is immediately produced and enters the physical channel. We denote the location of the center of the receiver by , and the distance between the center of the transmitter and the center of the receiver by . Furthermore, for the physical channel, we consider an unbounded environment affected only by diffusion noise; see Section III-D for more complex MC environments.
Passive receiver: Passive receivers (also referred to as transparent receivers or perfect monitoring receivers) employ passive reception mechanisms and are commonly considered in the MC literature, see e.g. [31, 32, 3, 33, 34, 35, 36, 37, 38, 39]. In particular, signaling molecules in the vicinity of the receiver can enter and leave the receiver via free diffusion; see e.g. Fig. 8a). The passive receiver model is a good approximation for the diffusion of small uncharged molecules such as ethanol, urea, and oxygen. These molecules can enter and leave a cell by passive diffusion across the plasma membrane [10]. A passive receiver model is also valid for the experimental system in [72], where the susceptometer that serves as the receiver does not impede the movement of the magnetic nanoparticles passing through it. For passive receivers, the set of all points inside the volume of the receiver, , constitutes the sensing area, and the number of molecules in constitutes the observed signal. Let denotes the number of molecules that the transmitter releases. Since we are interested in computing CIR , i.e., the probability that a molecule released by the transmitter at is observed at the receiver at time , we set . Moreover, we use the notation which can be interpreted as the PDF of a molecule released by the transmitter at with respect to at time . In other words, is the probability that the molecule is observed in volume at time . Since we focus on linear systems, solving with and solving for are related as . For the considered MC system with a point transmitter and unbounded environment, the CIR of a passive receiver can be obtained by first finding from (3) with the following initial and boundary conditions
[TABLE]
Given the solution of (3), , can be written as
[TABLE]
The solution of the integral in (36) can be readily obtained when the receiver is sufficiently far away from the transmitter, i.e., is very large relative to the largest dimension of the receiver. In this case, a common approach, which is referred to as the uniform concentration assumption (UCA), is to approximate everywhere inside the volume of the receiver by its value at the center of the receiver, i.e., . This leads to the following simple expression for , [31, 32, 3, 33, 34, 35, 36, 37, 38, 39]
[TABLE]
where is a constant denoting the volume of the receiver. We note that (37) is valid independent of the geometry of the receiver. Specifically, the UCA is one of the most useful approximation methods in the MC literature, since it directly relates the solution of (3), (17), and (32) to the CIR of the corresponding system. Thus, many results in the rich literature on solving PDEs, see [13], can be used to obtain the CIR in MC systems with passive receivers under the UCA.
The problem of solving (36) may become cumbersome when the receiver is close to the transmitter. In this case, the solution of the integral depends on the geometry of the receiver and the UCA cannot be hold. It has been shown in [33, Eq. (27)] that for a spherical passive receiver with radius , is given by
[TABLE]
where denotes the error function. Eq. (37) provides an accurate approximation for (38) if [33].
Remark 3
We refer the interested reader to [33] for an analytical expression for for a passive receiver with rectangular geometry.
Fully-absorbing Receiver: For fully-absorbing receivers [73, 74, 64, 75, 22, 76, 77] (also referred to as perfect sinks), unlike the passive receiver model, the physical and chemical properties of the receiver geometry are taken into account. In particular, the signaling molecules that reach the receiver via diffusion are absorbed as soon as they hit the receiver surface, see Fig. 8b). The sensing area of a fully-absorbing receiver is defined as all points on the surface of the receiver, , and the observed signal is the number of absorbed molecules during an infinitesimally small time . Here, a useful quantity that facilitates the derivation of is the rate of absorption of the molecule, which we denote by . Given , we have . Now, to derive , we first have to solve (3) with (34), (35), and the following boundary condition that models the absorption of the molecule on the surface of the receiver
[TABLE]
where in a spherical coordinate system, , for a spherical receiver with radius located at the origin of the coordinate system, i.e., , we have . Given , i.e., the solution of (3) with , , and , is given by [78, Eq. (3.106)]
[TABLE]
In [73], for a spherical absorbing receiver is introduced to the MC community and is calculated as [73, Eq. (22)]
[TABLE]
Another quantity of interest is the probability that a given molecule is absorbed by time , , which can be obtained as
[TABLE]
where is the complementary error function.
Remark 4
Alternatively, when the receiver counts the number of absorbed molecules during observation window , can be defined as
[TABLE]
Remark 5
For a fully-absorbing receiver, it is implicitly assumed that the whole surface of the receiver is fully-absorbing. The extension of this model to the case where the receiver surface is partially covered by fully absorbing receptor patches is considered in [74]. Moreover, the extension of the fully-absorbing receiver to take the impact of degradation and production noise into account, is considered in [64].
Remark 6
We note that one of earliest CIR models taking the absorption of particles in a 1D diffusion channel with uniform drift into account is proposed in [79]. There, a closed-form expression is derived for the probability of the time of absorption of the signaling molecules.
Reactive Receiver: Large or polar signaling molecules cannot passively diffuse through the membrane of cells and are detected by external receptors embedded in the cell membrane. In particular, the diffusive signaling molecules that reach the cell may participate in a reversible bimolecular second-order reaction with receptor protein molecules on the cell surface and form a ligand-receptor complex molecules; see e.g. Fig. 8c). The ligand-receptor interaction can be modelled as shown in (27) with binding reaction rate constant in [] and unbinding reaction rate constant in []. For such reactive receivers, the sensing area is that part of the receiver surface that is covered by receptors, denoted by , and the number of activated receptors constitute the received signal. We refer the interested reader to [60] for a closed-form CIR expression for reactive receivers.
Remark 7
In [61], a reactive receiver with an infinite number of receptor molecules covering the whole surface of the recevier, (i.e., a homogenous receiver surface, which is a special case of [60]), was considered and the corresponding CIR was numerically evaluated. Furthermore, in the MC literature, first steps to take the impact of ligand-receptor interaction on the CIR into account are made in [5] and [80]. There, for the evaluation of , the diffusion equation and the reaction equation are solved separately, unlike [60, 61] where a coupled diffusion-reaction equation is considered.
Remark 8
The fully-absorbing receiver is a special case of the reactive receiver when the whole surface of the receiver is covered with infinitely many molecules, , and . In this case, reaction equation (27) becomes a pseudo first-order reaction of the form , with binding reaction rate constant , where now corresponds to the number of absorbed molecules. However, implies that any collision of a signaling molecule with the receiver surface leads to the formation of a molecule, i.e., the reaction is deterministic. We refer the interested reader to [60] where it is shown how the CIR of the reactive receiver, under the above assumptions, is simplified to the CIR of the fully-absorbing receiver.
Remark 9
A receiver model that, unlike the CIR models reviewed in this section so far, also accounts for the impact of the signaling pathways, is proposed in [81]. In that model, two simple approximate signaling pathways, modeled via first-order and second-order CRNs, are considered. The CIR model in [81] is derived based on a mesoscopic modeling approach; see Section V for more details on mesoscopic modeling.
III-C Transmitter Models
In this section, we review some of the existing end-to-end CIR models developed in the MC literature that mainly focus on the properties of the transmitter. The main features of the transmitter that can potentially affect the end-to-end CIR include: i) the geometry of the transmitter, i.e., the volume, boundaries, and shape of the transmitter [2, 14]; ii) the particle generation via chemical reactions, which can take different forms ranging from a simple zero-order production reaction to more complex CRNs that take several aspects of molecule generation into account including, e.g., energy consumption via hydrolization of adenosine triphosphate (ATP) molecules [10]; and iii) the release mechanism controlling the release of the molecules into the physical channel. In particular, after production, the molecules can leave the transmitter either passively, for instance by passive diffusion through channels or gates embedded in the hull of the transmitter, or actively, for example via pumps integrated in the hull of the transmitter. In nature, passive and active transportation occur in cells via ion channels and transporters, respectively, see [10]. In the following, we study transmitter models that partially take the effects of the geometry, release mechanisms, and particle generation into account.
Point Transmitter: The point transmitter is the most widely used transmitter model in the MC literature mainly due to its simplicity, see [3]. However, this model takes none of the above mentioned features into account. In particular, the point transmitter, as the name suggests, is modelled as a zero-dimensional point, i.e., the impact of the geometry of a physical transmitter is not included in the model; see Fig. 9a). Furthermore, it is commonly assumed that the molecules are produced instantaneously and enter the physical channel immediately. These assumptions imply that the effects of the particle generation and the release mechanism on are neglected.
Volume Transmitter: Unlike point transmitters, where all molecules are generated at the same location, volume transmitter models take the transmitter geometry into account by assuming that the molecules are initially distributed over the transmitter volume666We note that, here, the term “volume” is generic and may refer to a volume or a surface in a 3D space, a surface or a line in a 2D space, and a line in a 1D space. [82]; see Fig. 9b). This leads to more realistic models since, in reality, signaling molecules are physical quantities that occupy space. However, volume transmitter models assume that the molecules are generated instantaneously, and that the surface of the transmitter is transparent and does not impede the diffusion of the molecules. With these two assumptions, volume transmitters neglect the effect of the particle generation and the impact of the release mechanisms. Let us, for the moment, denote the CIR models obtained for a point transmitter model, e.g., (38), (37), (41), by . Then, employing the principle of superposition and assuming a uniform particle distribution over the volume of the transmitter, , the CIR of the corresponding volume transmitter can be written as [82, Eq. (12)]
[TABLE]
where denotes the volume of the transmitter.
Remark 10
In [82], (44) is solved numerically for a 3D spherical transmitter and both passive and fully-absorbing receivers. Furthermore, in [82], closed-form expressions are given for corresponding one-dimensional scenarios.
One useful approximation of (44) can be obtained when the transmitter is sufficiently far away from the receiver. Then, the distance of any point inside the transmitter to the receiver can be approximated by and (44) simplifies to
[TABLE]
We refer the interested reader to [82], where the accuracy of the above approximation has been investigated for several environments.
Remark 11
The analytical transmitter models in [82] assume that the molecules are generated throughout . [82] and [83] simulated a volume transmitter model where the molecules are generated on the surface of a reflective spherical transmitter. In [83], a parametric model is proposed for the CIR of an MC system employing the considered transmitter and a fully-absorbing receiver. A machine learning approach is used to obtain the parameters of the parametric model.
Ion-Channel Based Transmitter: Ion-channel based (IC) transmitters are considered in [84] to model the effect of the release of the signaling molecules into the physical channel. IC transmitters take both the transmitter geometry and the release mechanism into account. In particular, IC transmitters are modelled as spherical objects with ion-channels embedded in their membrane; see Fig. 9c). The opening and closing of the ion-channels is controlled via a so-called gating parameter such as a voltage or a ligand. When the gating parameter is applied, e.g., in the form of a voltage across the transmitter membrane, the ion-channels open and the molecules can leave the transmitter via passive diffusion. The impact of the particle generation is neglected in [84]. In particular, it is assumed that the molecules diffuse with different diffusion coefficients inside and outside the transmitter. In [84] an expression is derived for the average rate of signaling molecules entering the physical channel upon transmitter stimulation. Furthermore, an approximate solution for the CIR of the corresponding end-to-end channel is provided under the conditions that the entire surface of the transmitter is covered by a large number of open ion-channels and that the signaling molecules diffuse with the same diffusion coefficient inside and outside the transmitter. There, a passive receiver under the UCA and an unbounded environment is considered. Then, the CIR is approximated as [84, Eq. (42)]
[TABLE]
where denotes the hyperbolic sine function. In fact, (46) is actually the CIR of a volume transmitter, since the assumption of having many open ion-channels is equivalent to assuming that the entire surface of the transmitter is a transparent membrane.
Remark 12
None of the transmitter models reviewed so far consider the impact of the particle generation via CRNs inside the transmitter. This is mainly due to the fact that to take the particle generation into account, a coupled reaction-diffusion equation has to be solved, which is a challenging task. Neverdeless, the effect of particle generation is studied in [85, 86, 87]. There, a common methodology for solving the corresponding reaction-diffusion equation is to adopt mesoscopic models and numerically solve the problem.
III-D Physical Channel Models
In this section, we review some of the existing end-to-end CIR models that emphasize the phenomena or impairments of the physical channel. In diffusive MC systems, the signaling molecules that enter the physical channel may be affected by several factors and noise sources besides diffusion, including: i) advection that can be constructive or destructive depending on the direction and strength of the velocity vector; ii) the geometry of the physical channel, e.g., bounded or unbounded environments, constraining the dispersion of the particles; and iii) degradation and production of molecules. For the CIR models reviewed in this section, in order to be able to focus on how is affected by the phenomena in the physical channel, we adopt the point or volume transmitter model and the passive receiver model.
Bounded Diffusion Channels: The CIR models reviewed previously were obtained under the assumption of an unbounded physical channel. Now, we focus on CIR models that assume a more elaborate physical channel geometry. To determine for bounded physical channels, generally, one has to solve diffusion equation (3) with appropriate boundary conditions reflecting the physical and chemical properties of the geometry of the channel. Unfortunately, for many practical geometries, simple and insightful solutions of (3) do not exist. Thus, approximations are needed to model practical geometries. In the following, we focus on a class of bounded physical channels that are referred to as duct channels. In particular, we consider duct channels with rectangular and circular cross sections; see Fig. 10. These two duct channels are of particular importance since they approximate the geometry of microfluidic channels and blood vessels, respectively.
- •
Rectangular Duct Channel: For a rectangular duct channel with dimensions , fully reflective walls, a point transmitter at , and a receiver at , the CIR can be obtained by solving (3) for with and the following boundary conditions
[TABLE]
where and capture the reflection of the molecule on the duct walls. Since, for the considered geometry, the diffusion of the molecule in one Cartesian coordinate does not influence its diffusion in the other coordinates, we can write . Now, solving (3) for , , and with , , and , respectively, and considering a passive receiver under the UCA, can be obtained as follows [13, Eq. (14.4.4)]
[TABLE]
- •
Circular Duct Channel: For a circular duct channel with dimensions in cylindrical coordinates, fully reflective walls, a point transmitter at , and a receiver at , the CIR can be derived by solving (3) with (34) and the following boundary conditions
[TABLE]
Here, models the reflection of the molecule at the wall of the duct with radius . Employing the same technique as for rectangular duct channels, using [13, Eq. (14.13.7)] and considering a passive receiver under the UCA, can be obtained as follows
[TABLE]
where the summation in is over the positive roots of . Here, denotes the -th order Bessel function of the first kind and denotes its derivative.
The CIR expressions (50) and (53) indicate that even for simple bounded environments, the solution to may be quite complicated and difficult to interpret. To gain more insight, in Fig. 11, we compare the CIR of an unbounded physical channel to that of a circular duct channel for system parameters , , receiver radius , and . Fig. 11 shows that when duct radius is small, the CIR of the duct channel is much larger than the CIR of the unbounded channel, i.e., for a given time it is more likely to observe the signaling molecule. This is because when is small, the signaling molecule is reflected more frequently on the duct walls which increases its chance of being observed at the receiver compared to the unbounded case where the molecules can diffuse away. However, for large , the CIR of the duct channel approaches the CIR of the unbounded environment, i.e., the CIR of the unbounded channel provides a good approximation for the CIR of a large bounded circular duct channel.
Remark 13
We note that the necessary condition for the validity of the UCA developed for passive receivers in unbounded channels, i.e., , is not applicable for bounded physical channels. However, we expect that as , the accuracy of (50) and (53) improves.
Advection Channels: Next, we consider physical channels in which the signaling molecules experience advection in addition to their random walk. In particular, for the CIR models reviewed in this section, for mathematical tractability, we consider advection processes with a time-invariant velocity field, i.e., we assume .
- •
Uniform constant advection: In this case, the magnitude and the direction of the velocity field are identical at any point in space, i.e., , where is the set of all points in the 3D Cartesian system, see Example 3 in Section II. Vector can be effectively decomposed into two components, a parallel component and an orthogonal component with respect to . Let us assume, without loss of generality, a point transmitter at and a passive receiver located at , such that . Then, , and we can write .
- –
Unbounded Channel with UCA: For an unbounded channel and a passive receiver under the UCA, can be obtained by solving advection-diffusion equation (17). Using the method of moving reference frame, i.e., assuming that the reference of the coordinate system is moving with , it can be readily verified that can be obtained from (37) as [37, Eq. (18)]
[TABLE]
Eq. (54) can be also directly obtained from (18) after setting , multiplying with , and using the , and mentioned above.
- –
Unbounded channel without UCA: For the case when the UCA does not hold, , and , can be solved numerically. However, it is shown in [37] that for the special case of , can be obtained from (38) after substituting with .
- –
Bounded channel with UCA: In this case, i.e., when we have bounded channels such as duct channels, and for the general case where , , we cannot apply the technique of moving reference frame in the dimensions of the coordinate system where the physical channel is bounded. Thus, has to be directly evaluated via the corresponding advection-diffusion equation. However, for the special case of , after substituting with , the corresponding CIR s of the rectangular and circular duct channels can be obtained from (50) and (53), respectively.
- –
Bounded channel without UCA: In this case, the general form of depends on the geometries of the bounded physical channel and the passive receiver. However, for a rectangular duct channel, a rectangular passive receiver, and , , an analytical expression for is derived in [47]. We note that in [47], it is assumed that and are a fluid velocity field and a drift velocity caused by a magnetic field, respectively. However, the derived expression for is valid independent of the origin of and . Furthermore, in [47], the case of partially absorbing duct channel walls is also considered.
- •
Laminar flow: In this case, we only focus on bounded channels, and in particular on circular duct channels, since laminar flow arises in bounded environments. Thus, we consider given in (14). For the CIR models reviewed here, we distinguish between point and volume transmitter models with axial position , and consider the passive receiver model with the following dimensions in cylindrical coordinates ; see Fig. 12. In particular, we distinguish between two cases, namely the dispersion regime () and the flow dominant regime (), see (20).
- –
Dispersion regime with UCA: In the dispersion regime, holds in (20). As a result, the time required for transportation of the molecule in the direction via flow, , is much larger than , which is the characteristic time for diffusion of the molecule over distance . This fact has two immediate consequences: i) by the time that the released molecule reaches the receiver, it experiences the average flow velocity, i.e., , due to its fast diffusion across the cross section; ii) there is no difference between point and uniform release and only depends on . Thus, the corresponding advection-diffusion equation in three dimensional space can be effectively approximated by its one dimensional equation with effective velocity and effective diffusion coefficient as follows
[TABLE]
where is the Aris-Taylor effective diffusion coefficient and can be obtained as [50, Eq. (4.6.35)]
[TABLE]
Solving (55) with the UCA approximation, , and the following initial condition for uniform release across the cross section
[TABLE]
leads to [54, Eq. (11)]
[TABLE]
- –
Dispersion regime without UCA: In this case, can be obtained by taking the integral of the solution of (55) over the volume of the receiver, which leads to [54, Eq. (13)]
[TABLE]
where denotes the Gaussian Q-function.
Remark 14
We note that the accuracy of both (– ‣ • ‣ III-D) and (59) depends on the value of in (20). For example, by increasing and , decreases and the accuracy of the dispersion regime improves; see [54].
- –
Flow dominant regime with volume transmitter: In this case, i.e., , the impact of diffusion is negligible. Thus, the signaling molecules do not have sufficient time to disperse across the cross section of the duct channel before they arrive at the receiver. As a result, a particle released at radial position is observed approximately at the same radial position at the receiver. Thus, we have to distinguish between the volume and point transmitter models. For the case of uniform release, can be approximated as [54, Eq. (16)]
[TABLE]
where
[TABLE]
In (61), and are the times when the parabolic velocity field first hits and leaves the receiver volume, respectively.
- –
Flow dominant regime with point transmitter: For the case of a point transmitter, when the molecule is released away from the receiver but within the and coordinates defining the dimensions of the receiver, i.e., at and , we observe the molecule with certainty if
[TABLE]
where if .
Remark 15
In the MC literature, first steps towards the extension of some of the CIR models reveiwed for the advection channel to more complex networks of interconnected bounded duct channels are provided in [88, 42, 45, 89]. For example, in [88], the Aris-Taylor effective diffusion coefficient approximation is employed to calculate the end-to-end CIR of multiple interconnected blood vessels for drug delivery applications. Furthermore, in [45], the uniform advection model is adopted to model blood vessel networks for abnormality detection applications in biological systems.
Degradation Channels: In degradation channels, the arrival of the signaling molecules is affected by their possible degradation and production. In this case, can be obtained by solving the diffusion-reaction equation (32) (with ) given appropriate initial and boundary conditions. However, the solution to (32) depends very much on the structure of the corresponding CRN described by reaction rate function . Often, the reaction terms in are highly nonlinear and coupled, which makes the problem of solving (32) very challenging. Here, in order to arrive at mathematically tractable and insightful results describing the general behaviour of degradation channels, we focus on two particular forms of degradation and production noise, namely first-order degradation and enzymatic degradation; see Examples 7 and 9, and [58, 18, 60, 64].
- •
First-order degradation: Let us for the moment denote the CIR expressions developed in the previous sections by . It can be shown that , for a physical channel with first-order degradation reaction of the form (25) and reaction rate constant can be readily obtained from when the following assumptions hold. A1) The signaling molecule is affected by the degradation reaction uniformly and equally throughout the entire end-to-end channel, and A2) it is not involved in any other CRN from stimulation time until observation time . In this case, we can write
[TABLE]
In (63), the term captures the surviving probability of the signaling molecule, which is a monotonically decreasing function of time, i.e., as increases, it becomes more likely that the signaling molecule is degraded in the channel. As a result, for degradation channels, at any instant, is smaller than the corresponding CIR without degradation.
Assumption A2) holds for all CIR models presented so far except the CIR model for the reactive receiver in Section III-B. This is because for reactive receivers, the signaling molecule is involved in a ligand-receptor kinetic reaction (27), in addition to the degradation reaction, and may experience several binding/unbinding events before reception time . We refer the interested reader to [60], where a closed-form expression is derived for the CIR of an MC system with a reactive receiver and first-order degradation channel.
- •
Enzymatic degradation: The impact of enzymatic degradation reactions in the channel is studied for passive receiver and point transmitter models in [58]. Enzymatic reactions include a second-order reaction, and as a result, in order to obtain , we have to solve (32) with , , and driven by (29a). However, this system of nonlinear equations does not facilitate a closed-form solution for . As a result, in [58], several approximate solutions were proposed for the CIR of an MC system with point transmitter, passive receiver, and unbounded environment as follows:
- –
If the concentration of signaling molecules, , and the concentration of the intermediate species , , remain constant over time and space, then is described only by (28a) and has the following solution
[TABLE]
- –
If and , then the total concentration of the enzyme molecules remains constant and . Then, a lower bound for is obtained via (– ‣ • ‣ III-D) after setting .
- –
Another useful approximation is obtained by assuming that is constant [58]. This is a valid assumption when and . Then, as explained in Section II, the enzymatic reaction in (28) can be approximated by the first-order unimolecular reaction in (25) and the corresponding can be written as
[TABLE]
where denotes the total concentration of the enzyme molecules, including both bound and unbound enzyme. We refer the interested reader to [58] for verification of the accuracy of the proposed approximate expressions for .
Remark 16
Signaling molecules of different types may also degrade each other. For instance, the MC testbed in [59] uses acids and bases as signaling molecules that can participate in a bimolecular reaction and cancel each other out, cf. (27). Unfortunately, the underlying PDEs that describe the bimolecular reaction are coupled and nonlinear and closed-form expressions for the CIR are not available. In [67], a numerical method was developed which decouples reaction and diffusion in each time slot and computes the channel response in an iterative manner.
III-E Summary of End-to-End CIR Models
To conclude this section, we provide a summary of the reviewed CIR models in Table II. For conciseness of presentation, we use the following abbreviations. For the transmitter, “” indicates whether a point transmitter is assumed or a volume transmitter releasing the molecules from its volume () or surface (). In the later case, we also specify whether the surface is transparent (“trans”) or reflective (“refl”). Furthermore, “PG” and “REM” specify whether particle generation and the release mechanism are taken into account, respectively. For the physical channel, “Diff” and “Adv” denote whether diffusion and advection processes are taken into account, respectively. In the case of advection, “uni” and “lam” specify uniform and laminar advections in the physical channel. “Geo” specifies whether a bounded (“”) or unbounded (“”) geometry is considered. “D&P” indicates degradation and production reactions of signaling molecules in the physical channel. For the receiver, “” indicates whether the volume of the receiver, , constitutes the sensing area of the receiver or a surface, , or a partial surface, . Furthermore, “passive” and “active” refer to the reception mechanism of the receiver. In the later case, “Det” and “Sto” specify whether the corresponding reaction in active reception is considered deterministically or stochastically, respectively. “SP” indicates whether signaling pathways are considered. Moreover, “Dim” denotes the dimension of the considered end-to-end channel. “Num” indicates that the CIR is obtained numerically. Whenever possible, we also refer to the equation number of the corresponding CIR . Finally, whenever the reaction-diffusion equation was involved for obtaining , we emphasize whether the reaction and diffusion processes are solved separately (“Sep”) or jointly (“Joi”). Thus, the reaction process can potentially describe an active reception mechanism, D&P reaction, or both.
IV Received Signal Modeling
In this section, we provide mathematical models for the signals used for estimation of the system parameters and detection of the transmitted data by MC receivers. To this end, we first present a unified signal representation for MC systems. Next, we introduce three time scales for the signal observed at the receiver, and subsequently, we provide signal models for each of these time scales. In addition, we generalize these models to account for the interfering noise molecules in the environment. Subsequently, time-slotted communication is considered and a corresponding signal model is developed which accounts for the impact of ISI. Finally, the correlation of the signals received at different time instants is discussed for the considered time scales.
The models that we present in this section are general in the sense that they apply to all MC systems discussed in Section III. More specifically, these models only depend on the CIR within the considered observation window or at the considered sampling times. We note that for most MC environments, derivation of the CIR in closed form, as was done for specific cases in Section III, is challenging. In Section V, we present numerical and simulation methods to obtain the CIR of more complex MC systems. In addition, in practical MC systems, the transmitter may send known pilot symbols that enable the receiver to estimate the CIR from its observations (e.g., see [38] and [90] for channel parameter estimators for MC systems). The models developed in this section are applicable for analytically derived, simulated, and estimated CIRs.
IV-A Unified Signal Definition
In the MC literature, different physical quantities have been modeled as the received signal. Important examples include i) the number of molecules observed at a given time within the volume of a transparent receiver [31, 32, 33, 34, 35, 36, 37, 38, 39], ii) the number of molecules bound at a given time to the receptors of a reactive receiver [60, 61, 74], iii) the accumulated number of molecules observed by a fully-absorbing receiver within a given observation time window [91, 73, 92, 93, 94], and iv) the arrival times of the molecules at a fully-absorbing receiver [95, 96, 97, 98, 99, 100]. In the following, we first provide a unified definition of the received signal of general MC receivers including the aforementioned special cases. Since the presented general signal model is difficult to analyze, subsequently, we introduce the concept of counting receivers, which are widely considered in the literature and allow for simple mathematical modeling. The main purpose for introducing a general representation of the received signal is to highlight the basic assumptions that have been made to arrive at specific signal models used in the literature and to unveil the connections between different signal models.
IV-A1 Generalized Receivers
Since different molecules of the same type are indistinguishable for the receiver, the most detailed information that the receiver could access at a given time is the arrival (and departure, if relevant) times of the molecules at (or from) the receiver up to that time. We use this fact to introduce a unified representation of the received signal of general MC receivers. For mathematical rigor, let us first formally distinguish between two types of receivers, namely recurrent and non-recurrent receivers.
Definition 3
If a given molecule can be observed by the receiver at most once, then the receiver is referred to as non-recurrent; otherwise, it is referred to as recurrent.
Transparent and reactive receivers with unbinding rate are recurrent since a given molecule can be observed multiple times by the receiver. On the other hand, fully-absorbing receivers and reactive receivers with are non-recurrent since after a molecule has been observed at the receiver, it cannot be observed again. For non-recurrent receivers, the time instants at which the molecules are observed constitute the most general signal representation. Let us define
[TABLE]
as the vector containing the arrival times , of all molecules observed by time in an ascending order. We note that both the number of molecules observed by time , i.e., , and their arrival times , are RVs. On the other hand, for recurrent receivers, in addition to , we need to keep track of the molecules that have been un-observed, i.e., have left the receiver. To this end, let us define
[TABLE]
as the vector containing the departure times , of all molecules that have left the receiver by time in an ascending order. We note that by the above formulation, non-recurrent receivers can be seen as a special case of recurrent receivers where . In summary, and constitute a complete and unified representation of the received signal of MC receivers. As will be shown in the following, different notions of received signal used in the MC literature can be seen as special cases of and .
IV-A2 Timing-based Receivers
In the MC literature, timing channels have been used as a model for non-recurrent receivers [95, 98, 99, 96, 97]. Let denote the vector containing the release times of the molecules by the transmitter and let be the vector containing the respective arrival times of the molecules at the receiver. Thus, is related to according to [95, 96, 97, 98, 99, 100]
[TABLE]
where is a vector containing the random delays between the release of the molecules by the transmitter and their observation at the receiver. Moreover, it is typically assumed that the release, propagation, and reception of molecules are independent from each other, which we refer to as the independent molecule behavior assumption [98, 99]. Based on this assumption, the elements in are independent and identically distributed and assume only non-negative real values. For an unbounded 1D environment, the random observation delay follows a Levy distribution if no flow is present [97] and the inverse Gaussian distribution if flow in the direction of the receiver is present [95].
We note that, in practice, is not available at the receiver since i) different molecules of the same type are indistinguishable by the receiver and ii) out of the total number of released molecules, only molecules are observed by time . In fact, is the actual observation signal available to the receiver. To arrive at a model for , we introduce the following definitions and assumptions. Let us assume that molecules are released by the transmitter within interval and their release times are collected in . Since the molecules observed at the receiver are indistinguishable, we do not know which molecules out of the total released molecules have been observed. In general, there are at most possibilities for selecting (observed) molecules from the (released) molecules. Therefore, we define , as a vector which contains the -th possible order index of the observed molecules. Moreover, let and denote the probability density function (PDF) and cumulative density function (CDF) of RV at , respectively. We note that due to causality, , has to hold where RV denotes the random delay of a given molecule. Following a similar framework as developed in [98, 99], the PDF of the observation vector conditioned on the molecule release time vector , denoted by , is obtained as
[TABLE]
where equality holds when all molecules are released at time zero. The above formulation provides a general framework for modeling the arrival times of non-recurrent receivers. Unfortunately, (IV-A2) cannot be easily simplified and its generalization to recurrent receivers or the cases when interfering noise molecules or ISI are present is cumbersome. In fact, the results reported in [95, 96, 97, 98, 99, 100] are valid for non-recurrent receivers when ISI and interfering noise molecules do not exist. In addition, perfect synchronization is a key underlying assumption for most timing channels considered in the literature [95, 96, 97, 98, 99, 100] and hence the performance of timing receivers is very sensitive to synchronization errors. Therefore, in the remainder of this section, we consider special receivers, namely molecule counting receivers, whose signal is a function of and only. Molecule counting receivers are widely adopted in the literature and the corresponding received signal lends itself to more tractable models and analysis.
IV-A3 Counting Receivers
These receivers consider the number of observed molecules as the received signal. In general, the receiver may count the number of observed molecules multiple times, which is referred to as a multi-sample detector [35, 32, 31, 46, 101, 83]. Let denote the received signal at sample time , where is the sample interval. To formally characterize , we distinguish two types of counting receivers, namely arrival-counting and observation-counting receivers.
Definition 4
If a receiver counts the number of molecules that have arrived within the observation window at its reception site, i.e., , then it is referred to as an accumulative-molecule-counting (AMC) receiver, whereas if the receiver counts the number of molecules that are observed at a given time at its reception site, i.e., , then it is referred to as an instantaneous-molecule-counting (IMC) receiver.
In general, there are four types of receivers based on the recurrent/non-recurrent and AMC/IMC classifications. In the following, we present the different counting receivers used in the MC literature as special cases of these four categories:
Non-Recurrent Accumulative-Molecule-Counting (nR-AMC) Receivers: The signal in this case is where . For instance, for fully-absorbing receivers, denotes the number of molecules that have arrived within interval [91, 73, 92, 93, 94].
Recurrent Accumulative-Molecule-Counting (R-AMC) Receivers: The signal in this case is where . Although the mathematical form looks identical to that for nR-AMC receivers, the modeling for R-AMC receivers is much more cumbersome since one molecule might be counted multiple times within the observation window . Furthermore, we note that the expected number of observed molecules for R-AMC receivers is larger than that for nR-AMC receivers because some molecules may be counted multiple times.
Recurrent Instantaneous-Molecule-Counting (R-IMC) Receivers: The signal in this case is where . For instance, for transparent receivers, denotes the number of molecules within the receiver volume at time [31, 32, 33, 34, 35, 36, 37, 38, 39], and for reactive receivers, is the number of molecules bound to the receiver’s receptors at time [60, 61, 74].
Non-Recurrent Instantaneous-Molecule-Counting (nR-IMC) Receivers: The signal in this case is where and . We note that since the received molecules do not leave the receiver, is a non-decreasing function of time.
In the remainder of this section, we focus on the modeling of for R-IMC receivers as a function of CIR , i.e., the probability of a molecule being observed at time seconds after its release by the transmitter; see Section III. This model is also valid for nR-AMC and nR-IMC receivers if is substituted by the probability of observing a molecule within intervals and , respectively, after its release by the transmitter at time , cf. (43). Modeling of for R-AMC receivers is cumbersome due to the possibility of counting a given molecule multiple times within the observation window. This type of signal is relevant, e.g., for ligand-based receivers when a ligand molecule can activate the receptors on the receiver surface multiple times. However, this problem has not yet been studied in the MC literature and is a potential topic for future research. Finally, in the following, we assume that the sampling interval is sufficiently large such that consecutive samples are statistically independent. Therefore, we drop index in Sections IV-B, IV-C, and IV-D for simplicity. How large should be chosen to guarantee sample independence will be discussed in Section IV-F.
IV-B Three Time-Scale Signal Representation
Let us define as the number of molecules observed at the receiver seconds after its release is stimulated by the transmitter at time . Then, can be modeled as
[TABLE]
where denotes the mean of the signal for a fixed set of channel parameters, denotes the random fluctuations around the mean (e.g., caused by diffusion), and denotes expectation. We note that the channel parameters may also change over time; however, this is over a scale that is much slower than the signal variations. In other words, the mean of the signal, , varies over time due to diffusion, advection, and reactions in the channel, but it also varies over the larger time scale due to variations of system parameters such as temperature, viscosity, and the distance between a mobile transmitter and receiver [102, 103, 65]. In summary, we have variations on three time scales in : Time Scale 1) variations of around its mean , i.e., noisy fluctuations , Time Scale 2) variations of the signal mean over observation time , which are slower than the variations of , and Time Scale 3) variations of over the release time , which are slower than the variations of with respect to . For instance, for typical MC systems at microscale, e.g., cell-to-cell communication, the noisy fluctuations vary on the order of a few s, the variations of the signal mean over time are on the order of tens or hundreds of s, and the change in the parameters, e.g., due to the mobility of the nodes, can be on time scales much larger than ms [60]. Fig. 13 illustrates versus for three values of . The aforementioned three time scales are illustrated in this figure: 1) the actual received signals, , are denoted by colored solid lines; 2) the black dashed lines denote the signal mean ; and 3) the variations of the signal due to changes in the system parameters over time scale are represented by different colors.
Remark 17
The three time-scale signal representation for MC is analogous to a similar signal representation in wireless communications. In particular, in a wideband wireless communication system impaired by additive white Gaussian noise (AWGN) and fading, the AWGN is analogous to the random fluctuations of the signal in MC, the CIR of the wireless communication channel is analogous to the signal mean in MC, and the variations of the CIR over time (due to the movement of the transmitter and/or receiver) are analogous to the time-variant signal mean in MC [104].
IV-C Signal Models
In the following, we first derive the expected number of molecules observed at the receiver, which we refer to as the deterministic model of the received signal. Subsequently, we derive statistical models of the received signal that capture the random fluctuations of the observed molecules. Finally, we study time-variant channels and derive stochastic models to capture the effect of the time variance.
IV-C1 Deterministic Models
In Section III, we derived the CIR which can be interpreted as the probability of a molecule released at time being observed at the receiver at time . Let us define as the probability of a molecule released by the transmitter at time being observed at the receiver at time . In the following, we first assume a time-invariant MC channel which leads to . Then, in Section IV-C3, we analyze the impact of time variance of the channel. Following the independent molecule behavior assumption [28, 91], the expected number of molecules observed at the receiver at time due to the release of molecules by the transmitter at time is readily obtained as
[TABLE]
For a given set of system parameters, the expected behavior is non-random and we have a deterministic signal model. Thus, each of the CIR expressions derived in Section III constitutes a deterministic representation of the respective MC system.
Remark 18
The independent molecule behavior assumption has to hold for (71) to be valid. However, for some MC systems, this assumption does not hold. For instance, if the number of receptors on the surface of a reactive receiver is finite, becomes a nonlinear function of the released molecules and cannot be described by the simple linear expression in (71). This effect is known as receptor occupancy [60]. In these cases, has to be found for a given either numerically or via simulation, cf. Section V-A for a detailed discussion on simulation methods.
Remark 19
The deterministic model in (71) assumes an impulsive release of molecules at time by the transmitter. In general, the transmitter may release the molecules continuously over a finite time interval of length . Let denote the release function satisfying and . Then, the expected number of molecules observed at the receiver at time due to the release of molecules by the transmitter with release function is given by
[TABLE]
We note that (72) reduces to (71) for . In the remainder of this section, we focus on impulsive release, as this is typically assumed in the MC literature.
IV-C2 Statistical Models
In the following, we develop statistical models for the number of molecules observed at the receiver as a function of for time-invariant MC channels.
Binomial Model: Based on the independent molecule behavior assumption and since any given molecule released by the transmitter is either observed by the receiver or not, a binary state model applies and the number of observed molecules follows the Binomial distribution with trials and success probability , i.e.,
[TABLE]
where represents a Binomial distribution with parameters and denoting the number of trials and the success probability, respectively. Under the Binomial model, the probability mass function (PMF) of , denoted by , is given by
[TABLE]
for . Unfortunately, the Binomial distribution considerably complicates the analysis of MC systems. Therefore, in the following, we present two approximations of the Binomial model with better mathematical tractability.
Gaussian Model: If the expected number of molecules observed at the receiver, i.e., , is sufficiently large, then we can apply the central limit theorem (CLT) and approximate by a Gaussian RV with mean and variance identical to that of the Binomial RV. This leads to
[TABLE]
Under the Gaussian model, the PDF of , denoted by , is given by
[TABLE]
The Gaussian distribution is much more amenable to analysis than the Binomial distribution. However, the basic assumption behind the applicability of the Gaussian distribution, namely large , may not hold in MC systems. In fact, although the number of released molecules can be quite large, the expected number of observed molecules can be very small. Moreover, Gaussian RVs are continuous and can assume non-integer and negative values, which contradicts the true nature of RV as discrete and non-negative.
Poisson Model: For the case when the number of trials is large and the mean of the Binomial RV is small, the Binomial distribution can be well approximated by a Poisson distribution with the same mean , i.e.,
[TABLE]
where represents the Poisson distribution with parameter denoting the mean of the RV. Under the Poisson model, the PMF of , denoted by , is given by
[TABLE]
In fact, assuming is fixed, the proof simply follows from [40]
[TABLE]
where for equality we used and [105].
Comparison: In order to quantify the accuracy of the Gaussian and Poisson approximations, we define the root mean square error (RMSE) between the the approximated Gaussian and Poisson CDFs, denoted by , and the Binomial CDF, denoted by , as [91, 93, 106]
[TABLE]
In Fig. 14, the RMSE between the approximate Gaussian and Poisson CDFs and the Binomial CDF versus is shown for . We observe from this figure that by increasing , the accuracy of the Poisson model deteriorates whereas the accuracy of the Gaussian model improves, which is consistent with the respective assumptions that led to their derivation. Moreover, as increases, the Gaussian model becomes more accurate whereas this is not true for the Poisson model if is very small. In fact, for small and , both the Binomial and Poisson distributions approach the binary distribution, i.e., either zero or one molecule is observed and the probability of observing more than one molecule becomes negligible. Since for typical MC systems, the value of is expected to be much smaller than , the Poisson model is generally a more accurate model. Nevertheless, the fact that the accuracy of the Gaussian model increases with increasing makes it a suitable model for macroscale applications when is potentially very large. Moreover, the Gaussian model is attractive for asymptotically high signal-to-noise ratio (SNR) analysis. These observations are consistent with the results reported in [91].
IV-C3 Time-Variant Models
Until now, we have assumed time-invariant MC channels where the channel parameters are fixed. Hence, and consequently were only functions of . Now, we consider time-variant MC channels where and are in general functions of both and . More specifically, we study the impact of system parameter variations on the mean received signal . In principle, each of the system parameters such as , , the physical and chemical properties of the boundaries of the end-to-end channel, the reaction rates of the involved CRNs, and and can potentially vary over time, which in turn leads to a variation of . For instance, the diffusion coefficient appears in the expressions for for all diffusive MC systems, and consequently in . As we can see from (2), changes in the fluid environment channel parameters, e.g., the viscosity or temperature, will result in a change in . In fact, the impact of variations in on for a point transmitter and passive receiver in 1D is investigated in [103]. The authors in [65] consider a point transmitter with impulsive release, a passive receiver with the UCA, and an unbounded 3D environment with uniform flow and uniformly distributed enzymes, cf. (36) and (10). There, the impact of Gaussian variations in the diffusion coefficient, flow velocity, and enzyme concentration is modeled by a parametric model where the parameters of the model are obtained via curve fitting. The impact of the mobility of a point transmitter and a passive receiver on the CIR is studied in [102] and a stochastic model for is derived. Similarly, a stochastic model for mobile MC systems with a point transmitter and fully-absorbing receiver is derived in [94]. We note that mobile transceivers are relevant for many envisioned applications of synthetic MC systems such as targeted drug delivery and health monitoring [42, 43, 44, 45]. Therefore, in the following, we focus on diffusive mobile MC systems and review some of the results reported in [102].
We assume a point transmitter, a passive receiver with the UCA approximation, and an unbounded diffusive channel without advection. Furthermore, we model the mobility of transmitter and receiver via 3D diffusion, since diffusion is a common cause of mobility and can also be used to model more elaborate movements such as cell migrations and bacteria chemotaxis [107, 108]. In particular, we denote the diffusion coefficients of transmitter and receiver by and , respectively, and their corresponding locations at time by and , respectively. Then, it can be shown that follows a Gaussian distribution [102, Eq. (2)]
[TABLE]
where is an effective diffusion coefficient capturing the relative motion of transmitter and receiver, see [109, Eq. (10)]. Then, given (37), the CIR of the end-to-end channel can be rewritten as
[TABLE]
where and is the effective diffusion coefficient capturing the relative motion of the signaling molecules and the receiver, see [109, Eq. (8)]. The movement of the receiver affects both (81) and (82) via and , respectively, as long as its movement with respect to the transmitter and the signaling molecules is accounted for. For any given , is a stochastic process with RVs . In the following, we analyze the mean, the variance, and an approximate expression for the PDF of .
Mean: Let denote the distance between the transmitter and receiver at , i.e., . Given (81) and (82), the mean of the time-variant channel, denoted by , can be evaluated as [102, Eq. (14)]
[TABLE]
where in , the expectation is taken with respect to the RV . As we expected, is a function of , because of the mobility of transmitter and receiver. As a result, is a non-stationary stochastic process. Moreover, due to the assumption of an unbounded environment, on average the transmitter and receiver diffuse away from each other. Therefore, when at least one of the transceivers is mobile, i.e., , we obtain as .
Variance: The variance of , denoted by , is given by
[TABLE]
where the second order moment is obtained as [102, Eq. (21)]
[TABLE]
We note that as , which is due to the fact that as . On the other hand, it can be shown that the normalized variance as . In other words, the normalized variance increases as increases. This in turn implies that due to the random walk, the uncertainty that we have about increases as increases.
Approximate PDF: In the following, we present the approximated PDF of the considered time-variant channel with mobile transceivers and refer the interested reader to [102] for the exact expressions of the CDF and PDF. In particular, it is shown in [102] that when holds for any , then the PDF of the CIR can be accurately approximated via a log-normal distribution [102, Eq. (29)]
[TABLE]
with
[TABLE]
where and denote the mean and the variance of the log-normal distribution. Given (86), the PDF of , denoted by , can be written as
[TABLE]
The above stochastic model can be used for the design and performance analysis of time-variant MC systems. For instance, (87) was used in [109] to compute the expected error probability of a mobile MC system when the knowledge of CIR used for detection becomes gradually outdated due to the mobility of the transceivers. Moreover, in [65], a stochastic channel model was used to develop non-coherent detectors. In contrast to the stochastic model in (87) for mobile MC systems, it was shown in [65] that the Gamma distribution is a good fit for (Gaussian) variations in the diffusion coefficient, flow velocity, and enzyme concentration in a non-mobile MC system.
IV-D Interfering Noise Molecules
In the previous section, we have considered statistical models for the number of molecules observed at the receiver due to the release of signaling molecules by the transmitter. However, MC systems may be impaired by noise molecules that are not released by the transmitter but originate from interfering natural or synthetic sources. In the following, we introduce statistical models to account for the number of noise molecules that are observed at the receiver. Since information and noise molecules are indistinguishable, the receiver treats the total numbers of observed signaling molecules, denoted by , and interfering noise molecules, denoted by , as the received signal , i.e.,
[TABLE]
To derive a statistical model for , we focus on a passive receiver. Similar arguments apply for other receiver types. We make the following assumptions. A1) Let denote the expected number of noise molecules observed within the receiver volume at a given sample time . We assume that the value of is constant over observation time . Nevertheless, may change over larger time scale due to variations in the system parameters such as the temperature, cf. Section IV-B and Section IV-C3. A2) It is further assumed that the observation of one noise molecule at the receiver is independent from the observations of other noise molecules. A3) Finally, we assume that the expected number of noise molecules observed within a given volume in space is proportional to the size of that volume.
Based on assumptions A1-A3, the statistics of the observed noise molecules is Poisson following the law of rare events (LRE) [110]. In particular, suppose the receiver volume is divided into subvolumes where . Thus, can be interpreted as the probability that one noise molecule is observed in one of these subvolumes at the receiver. The probability that two noise molecules are simultaneously observed in one subvolume becomes negligible for large . Therefore, the number of noise molecules observed over the entire volume of the receiver follows a Binomial distribution with trials and success probability . Consequently, since is a free variable, one can assume such that the Binomial distribution approaches the Poisson distribution , cf. (78). In summary, under assumptions A1-A3, we obtain .
Remark 20
The choice of the Poisson distribution for the number of environmental noise molecules observed at the receiver, , can be further justified from an information-theoretic perspective [111]. Let us define RV where denotes the coordinates of the -th noise molecule observed at the receiver and we drop argument of in for notational simplicity. In particular, the maximum entropy distribution for corresponds to a Poisson distribution for the number of observed noise molecules . Therefore, the most random noise under assumptions A1-A3 is Poisson noise, i.e., a worst-case scenario. To see this, let denote the distribution of RV . Using the chain rule, we have where is the conditional distribution of given and denotes the distribution of . For maximum entropy, should be a uniform distribution across the receiver volume. Substituting this result in , we obtain that has to be the Poisson distribution to maximize the entropy of [111, Appendix 8].
We note that the Poisson distribution approaches a Gaussian distribution for . Therefore, for very noisy environments, the approximation becomes valid.
Remark 21
More detailed models for the number of noise molecules observed at the receiver can be found in [71] where it was assumed that the noise molecules originate from external sources that continuously release molecules or transmitters of other communication links which use the same type of molecule for signaling. There, the expected number of noise molecules observed at the receiver at time after the noise sources start releasing molecules, denoted by , was derived as a function of the system parameters such as the distances between the noise sources and the receiver. It was shown that asymptotically as , converges to a constant value, i.e., , which is consistent with assumption A1 made earlier in this section.
IV-E Continuous Transmission
The statistical models developed so far are appropriate for one-shot transmission. Nevertheless, in most communication systems, the transmitter may send multiple symbols consecutively to the receiver. To develop a model valid for continuous transmission, we consider a time-slotted communication system where one symbol is transmitted in each time slot, also referred to as a symbol interval, of length . We focus on concentration shift keying (CSK) modulation where the transmitter releases molecules at the beginning of the -th symbol interval to convey information symbol [3]. We assume synchronous transmission and that the receiver counts the number of observed molecules multiple times in each symbol interval with sampling interval [101]. Because of the memory of the MC channel, ISI occurs. To take this into account, we assume that the MC channel has a memory of symbol intervals, i.e., the ISI in symbol interval originates from the symbols transmitted in the previous symbol intervals. We further take into account that communication may be impaired by noise molecules that originate from interfering natural or synthetic sources. Finally, we assume that the MC channel parameters remain unchanged for the considered observation window, and hence, we drop argument in , , and for notational simplicity. In the following, we first provide the signal model for a general case and subsequently simplify it for extreme SNR regimes to obtain further insight.
IV-E1 General Case
Let denote the total number of molecules observed at the receiver for sample in symbol interval , i.e., where and is given in (70). Then, following the discussion in Section IV-C2, can be accurately modeled as a Poisson RV, i.e.,
[TABLE]
where . Moreover, we used the superposition property of Poisson RV s, i.e., if and are two independent Poisson RV s with means and , respectively, then is also a Poisson RV with mean [112]. Alternatively, defining , one can obtain the following more familiar additive signal model
[TABLE]
where denotes the diffusion noise and denotes the interfering noise molecules. Here, we use the notation when where , i.e., is a Poisson RV whose mean has been subtracted.
When the expected numbers of information and interfering noise molecules are large, one may use the Gaussian model for the number of observed molecules, i.e., where . One can also write in the form of (90) where for the Gaussian model, we have and . We note that unlike for the AWGN channel in conventional wireless communication, the Gaussian diffusion noise in MC is signal dependent.
Remark 22
In (90), we distinguish between two types of additive noise, namely , which originates from signaling molecules, and , which originates from external interfering noise molecules. We note that the randomness of and can be attributed to the random Brownian motion of the signaling and noise molecules, respectively. In addition to the aforementioned noises, other types of noises may be present. For instance, in a reactive receiver, the noisy measurements of the activated receptors, caused by the randomness of diffusion and ligand-receptor interactions, may be relayed by signaling pathways to the interior of the receiver (e.g. a cell), which may add extra noise [5, 106]. We refer to this noise as counting noise to contrast it with the diffusion noise.
IV-E2 Simplifications for Extreme SNR Regimes
In the following, we further simplify the model in (90) for two asymptotic SNR regimes, namely the diffusion-noise-limited and interference-limited regimes. To do so, we first formally define SNR as [113]
[TABLE]
where denotes the expected number of signaling molecules received at the sampling time. In the following, we focus on the ISI-free channel, i.e., , and a single-sample detector. Therefore, we drop indices and for notational simplicity.
Remark 23
One approach to obtain an approximately ISI-free channel is to choose a sufficiently large symbol interval such that the CIR practically fully decays to zero within one symbol interval. In such a case, the transmission rate may be severely reduced which may lead to an inefficient system design. Fortunately, it has been shown in the literature that reactions can be beneficial for ISI mitigation [58, 59, 67]. In particular, enzymes [58] and reactive signaling molecules, such as acid and base molecules [59, 67], may be used to speed up the decay of the CIR as a function of time, which would increase the accuracy of the assumption of an ISI-free channel, see Fig. 5. For instance, in [67], a reactive signaling MC system was assumed where the transmitter employs different molecules that react with each other, e.g., acids and bases. Then, after the release of the signaling molecules (e.g., an acid), the transmitter may release so-called cleaning molecules (e.g., a base). It is shown that the resulting CIR is considerably shortened, which makes the ISI-free channel an accurate model, see Fig. 15. Moreover, the peak of the received signal remains unchanged since the cleaning molecules are released after the peak is observed at the receiver.
Diffusion-Noise-Limited SNR Regime: In this case, we assume holds. Thus, the model in (90) simplifies to
[TABLE]
where and hold for the Poisson and Gaussian models, respectively. The SNR in this case is obtained as .
Interference-Limited SNR Regime: In this case, we assume holds. Thus, the model in (90) simplifies to
[TABLE]
where and hold for the Poisson and Gaussian models, respectively. The SNR in this case is obtained as . We note that this special case yields a signal-independent (Gaussian) model as it is widely adopted in conventional wireless communications.
Finally, we note that it may be necessary to use a combination of both of the above special cases for analysis of MC systems. For instance, for a simple on-off keying (OOK) modulation, i.e., , the interference noise molecules are dominant for bit whereas the diffusion noise is dominant for bit . This is schematically illustrated in Fig. 16 where it can be observed that the noise power for symbol (diffusion-noise-limited regime) is larger than that for symbol (interference-limited regime). A similar observation has also been reported for photon-counting receivers in optical wireless communications where the shot noise at the receiver has two components, one generated by the laser transmitter (analogous to diffusion noise) and one generated by the ambient background light (analogous to interfering noise molecules) [114].
IV-F Time Correlation
In the following, we discuss the signal correlation with respect to observation time scale and release time scale .
IV-F1 Sample Correlation
In (89) and (90), we assume that the number of molecules counted at different time instants within one symbol interval or in different symbol intervals are independent from each other. However, this assumption holds only if the sampling interval is chosen large enough such that the independence of consecutive samples is guaranteed. In [46], the mutual information between two samples and was numerically computed and the minimum spacing needed to ensure independence between consecutive samples was found such that the corresponding mutual information is below some threshold. Since the mutual information between two samples is difficult to derive in closed form, one may consider the Pearson correlation coefficient [115] among two consecutive samples instead, i.e.,
[TABLE]
where equality follows from the fact that under both Poisson and Gaussian statistics, the variance of is . The cross-correlation term depends on the specific adopted receiver type. Note that by definition, holds. Typically, the sample times and should be separated such that falls below a certain threshold, denoted by , i.e., should be large enough such that holds. In Fig. 17, we show the absolute correlation versus where denotes the peak of the expected received signal. As can be seen from this figure, the correlation decreases as increases. Moreover, as an example, we choose the value of the threshold as . One can observe from Fig. 17 that as the diffusion coefficient of the molecules increases, the minimum sample spacing needed to ensure decreases.
IV-F2 Mean Correlation
Recall that if the system parameters change, the mean signal varies over time scale . In a similar manner as for sample correlation, one can define a correlation factor between the mean signals at time and as follows
[TABLE]
For the case when transmitter and receiver mobility are the cause of the variations in , cf. Section IV-C3, and are given by (IV-C3) and (84), respectively. Moreover, the cross-correlation, denoted by , for two arbitrary times and is derived in [102, Eq. (19)]
[TABLE]
where for compactness , , , , and are respectively defined as
[TABLE]
In order to quantify the time variations of the end-to-end MC channel, we define the coherence time, , as the minimum time for which falls below a certain threshold value , i.e., [102]
[TABLE]
The coherence time of the channel is a metric which determines the time over which the channel does not change substantially. As such, the particular choice of depends on the application of interest. Future applications of synthetic MC systems that are more robust to CIR variations can assume smaller values of , whereas applications that are more sensitive to CIR variations may require larger values of . For example, typical values of reported in the conventional wireless communications literature span the range from to [116, 117, 118]. Smaller values of are often employed for resource allocation problems, while larger values of are used for channel estimation. In Fig. 18, we show the absolute correlation versus for different scenarios of transmitter and receiver mobility, i.e., . As can be seen from Fig. 18, the channel mean decorrelates as increases. Moreover, assuming a fixed threshold , the coherence time decreases as the diffusion coefficients of the transmitter and receiver increase.
V Simulation- and Experiment-Driven Models
The analytical results presented thus far in this tutorial have focused on tractable solutions based on the underlying physical principles of advection, reaction, and diffusion. In order to arrive at these results, we often had to make assumptions that simplify the physical transmitter, receiver, and channel. However, this approach has limitations. Assumptions are generally constrained by specific channel parameters or the conditions for which they accurately apply. For example, we can assume that an environment’s outer boundary is unbounded if it is sufficiently larger than the signaling range; see Fig. 11 and [119]. Similarly, we can model the locally-varying concentration due to a molecule source as uniform if we are observing from a distance that is sufficiently far from the source; see the UCA in Section III and [33].
Sometimes we are able to relax assumptions and still maintain analytical tractability, cf. Section III. When this occurs, we can define a reliable rule of thumb that dictates explicit conditions under which the assumption can be satisfied to some degree of accuracy. For instance in [33], it was shown that the simplified CIR with UCA was within 2 % of the ideal CIR for most of the time of interest if the radius of a spherical receiver was no more than 15 % of the distance from the molecule source to the center of the receiver. However, we generally do not have the option to relax assumptions for analytical tractability while maintaining sufficient accuracy. Furthermore, we might encounter a channel with complex or novel phenomena where we do not yet know what suitable assumptions might be.
In the absence of reliable analytical results, we must rely on data-driven approaches to model a channel. Such approaches can also be used to help verify analytical results. This section reviews simulation and experimental approaches for generating data. Simulations can provide an efficient means for channel modeling, even in the presence of complex and coupled physical phenomena. Reliable experimental data may be preferred, but can be time-consuming and expensive to obtain.
V-A Simulation-Driven Models
Simulations of reaction-diffusion systems can be performed over a range of physical scales. As such, there are a range of simulation classes available, which we summarize in Fig. 19 and also discussed in [119]. We refer to these classes as continuum simulations, mesoscopic simulations, microscopic simulations, and molecular dynamics simulations. Generally, each class is suitable for a particular scale. Not surprisingly, there is an inherent trade-off between the physical resolution of a simulation and the computational resources (whether measured in time or memory) that are required to simulate it. The continuum approach is most suitable for macroscale systems. Both the microscopic and mesoscopic approaches can be appropriate for microscale systems. The molecular dynamics approach is most suitable for systems at the nanoscale and smaller. While the microscopic approach has been the most common simulation method within the MC research community, here we discuss all four approaches, their relevance, and also the potential to combine them in a single simulation. For the microscopic and mesoscopic approaches, we also describe how to implement a simple simulation.
V-A1 Continuum Simulations
When the physical scale of a simulation, including the number of molecules, is sufficiently large, then the evolution of the system can be directly described using the corresponding spatio-temporal PDEs, see e.g. (32). We refer to these as continuum simulations. Specifically, finite element analysis is used to spatially partition the system into a grid (see Fig. 19a)), and the system is simulated over a sequence of time steps. The molecule concentrations at each node in the grid are updated in every time step according to the differential equations that describe the phenomena. The updated concentrations are always non-negative real values. Popular commercial solvers that follow this approach include COMSOL Multiphysics [120] and ANSYS [121].
Unless the differential equations are stochastic or explicit noise sources are introduced, then continuum simulation of a system is deterministic. Generally, the accuracy depends on the resolution of the grid and the size of the time step; more accurate simulations can be performed by increasing the grid resolution and decreasing the size of the time step. However, as the nodes in the grid become increasingly close, the number of molecules associated with each node decreases. When the molecule concentrations get sufficiently small, it becomes more appropriate to consider integer numbers of molecules instead of continuous-valued concentrations. Thus, we next discuss mesoscopic simulations.
V-A2 Mesoscopic Simulations
Like continuum simulations, mesoscopic simulations also partition the system into a grid. The resulting containers are commonly referred to as subvolumes or voxels. However, instead of tracking continuous molecule concentrations, mesoscopic modeling counts discrete numbers of molecules in each subvolume; see Fig. 19b). Instead of deterministically solving the system’s set of PDEs, mesoscopic simulations proceed by stochastically generating event times, where each event is the occurrence of a chemical reaction or a molecule’s transition between two subvolumes. The key assumptions of this scheme are that i) molecules within a given subvolume are uniformly distributed, and ii) the solvent molecules in a subvolume can be treated as a homogeneous continuum that is in thermal equilibrium. When these assumptions are satisfied, the mesoscopic approach can simulate chemical reactions exactly (in a statistical sense, as proven in [122]). Furthermore, if the subvolume sizes are appropriately chosen, advection-reaction-diffusion systems can also be simulated exactly; see [123, 124]. An important constraint is that the subvolume size (e.g., cube length) must be much smaller than both and , where is the dimension of the subvolume, is the characteristic time of the fastest reaction in the system, is the diffusion coefficient of the largest corresponding reactant, and is the magnitude of the flow velocity. If this is not satisfied, then we cannot safely assume that the subvolumes are well-stirred (i.e., that the molecules are uniformly distributed).
Simple Implementation: A basic implementation of a mesoscopic simulation with equal-sized subvolumes (of length ) is as follows. Let be the number of molecules of the -th type that are in the -th subvolume. Events are associated with propensities. The propensity of a transition of a molecule of the -th type to diffuse from the -th subvolume to the -th subvolume, where these two subvolumes are adjacent and share a face, is [125, Eq. (1.6)]
[TABLE]
where is the diffusion coefficient of the -th molecule type. The propensity of the -th chemical reaction in the -th subvolume is [126, Eq. (6)]
[TABLE]
for zeroth-, first-, and second-order reactions, respectively, where the order corresponds to the number of reactants. is the corresponding reaction rate constant, is the subvolume volume, and and are the corresponding numbers of reactant molecules. For the entire system, the total propensity is then
[TABLE]
and we can simulate the time of the next event in the system by generating exponential random variable
[TABLE]
where is a random number uniformly distributed between 0 and 1. We can determine which of the possible events occurred by tossing a weighted die, where the likelihood of each event is proportional to its associated propensity. Once the event is determined, we update the molecule counts, update the corresponding propensities, and repeat the process to find the next event.
Remark 24
We note that there are mathematically equivalent but more computationally efficient implementations, particularly when updating propensities. These include Gibson and Bruck’s Next Reaction Method; see [127]. Furthermore, different accuracy-efficiency trade-offs can be introduced to provide more flexible scalability. For example, tau-leaping can be used to execute multiple events in a constant time step, where “tau” refers to the time step size; see [128]. Tau-leaping enables a transition between continuum and mesoscopic simulations; if the number of events during one “leap” is sufficiently large, then it can be treated as a deterministic value. As long as the propensities do not significantly change between time steps, then tau-leaping’s computational efficiency gains can be made with minimal losses in accuracy.
V-A3 Microscopic Simulations
In some sense, microscopic simulations are the dual of the mesoscopic approach. Whereas the (non-leaping) mesoscopic approach is continuous over time and discrete over space, the common microscopic approach implementation is discrete over time and continuous over space; see [69]. Instead of relying on well-stirred subvolumes, microscopic simulations track every molecule individually (i.e., particle-based simulation); see Fig. 19c). Nevertheless, they still assume that the solvent is a continuum of molecules, which means that the diffusion of the molecules of interest is still governed by a diffusion coefficient.
Simple Implementation: A basic implementation of a microscopic simulation with flow and first-order reactions in the propagation environment is as follows. In each time step , every molecule is tested for every possible first-order reaction. If there is only one potential reaction, and the associated reaction rate is , then the corresponding reaction probability is [69, Eq. (14)]
[TABLE]
If a coin flip with this probability is successful, then the molecule is converted to the corresponding reaction product. After all of the possible reactions have been tested, the remaining molecules are diffused along every available dimension by adding a displacement of towards each dimension of the Cartesian coordinate system, cf. (1). The realizations are independent for every molecule and along every dimension. Furthermore, if the environment has a bulk flow with a component along a particular dimension, then every molecule should have an additional displacement of along that dimension, cf. (11). Diffusion should be unimpeded, unless there are boundaries in the environment. For example, if a molecule crosses a solid reflective surface, then the coordinate that is normal to the surface is reverted to its value before diffusion. If a molecule crosses an absorbing surface, then it should be consumed by the absorbing reaction.
Due to their simplicity and their suitability for simulations over a range of nanometers to micrometers, microscopic simulations have been common for cellular systems and also specifically for MC systems. Mature tools from the physical chemistry community include Smoldyn (see [69, 129]). Microscopic tools that have been developed specifically for the MC community include BiNS2 [130], N3Sim [131], MUCIN [92], and AcCoRD [119].
V-A4 Molecular Dynamics Simulations
At a more precise scale, solvent molecules and their interactions with solute molecules and with each other can be modeled in detail; see Fig. 19d). These are molecular dynamics simulations, and they might account for intermolecular forces (including those imposed by charge potentials) and collision dynamics. One such example is the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS); see [132]. Due to the very large density of molecules to be considered, molecular dynamics simulations are best suited for very small systems, e.g., on a nanoscale. For instance, molecular dynamic simulations can be used to study how the conformation of receptor proteins change after binding to a specific molecule. Thus, they have generally not been applied to study MC systems.
Remark 25** (Hybrid Simulations)**
The aforementioned discussion of simulation classes has emphasized their suitability for simulations over different physical scales. However, a particular system might have multiple scales of interest. In order to avoid constraining the entire simulation by the most granular approach needed, hybrid simulation tools have sought to integrate different classes within a single simulation. One approach has been to combine microscopic and mesoscopic models, using hybrid interfaces such as that proposed in [125] and later implemented in Smoldyn (see [133]) and AcCoRD (see [119]). Other examples include the integration of the molecular dynamics solver LAMMPS with a continuum model (see [134]), and the integration of the continuum solver Virtual Cell with the microscopic approach in Smoldyn (see [135]).
Example 11** (Example Simulation)**
We complete our discussion of simulations with a brief demonstration. We consider an extension of the bounded rectangular-duct channel discussed in Section III-D that has no readily available analytical channel response. Nevertheless, we can simulate the system. The environment is a microfluidic system where two chambers are connected via a long pipe, as shown and described in Fig. 20. We place molecules uniformly within one of the chambers (i.e., a cube). These molecules can diffuse out through the rectangular pipe and into the other chamber and no flow is considered. The second chamber has a perfectly-absorbing surface and we count the number of molecules that are absorbed. We assume that the receiver counts the number of molecules absorbed by time . Therefore, the receiver can be classified as nR-AMC, i.e., non-recurrent and accumulative-molecule-counting, with received signal , cf. Section IV-A. A realization of this system is simulated using a microscopic approach in the AcCoRD simulator with a simulation time step of ms. The number of absorbed molecules is plotted in Fig. 21 for different pipe lengths. We see that all released molecules are absorbed within about s for the shortest pipe length (i.e., m). As the distance between the two chambers increases, fewer molecules get absorbed within the same time. We note that one can obtain the CIR of this system, i.e., defined as the probability of a molecule being absorbed at the receiver in interval after its release by the transmitter at , by simulating the system for many realizations and averaging the result, i.e., .
V-B Experimentally-driven Models
In the previous subsection, we have seen how elaborate simulations can be used for scenarios where it is difficult or even impossible to derive an analytical model based on physical principles. However, for practical systems, we may also face situations where even simulation of certain phenomena is challenging. In fact, even complex simulation methods typically cannot account for all characteristics of a real experimental environment. In the following, we first highlight some of the unique characteristics of two existing experimental platforms for MC [57, 136] that cannot be easily modeled or simulated. Subsequently, to cope with the aforementioned challenges, we present a general data-driven modeling approach which is then applied to an example experimental system.
V-B1 Challenges in Modeling Existing Experimental Systems
Several experimental systems exist for demonstrating MC. These testbeds include both non-biological systems [57, 59, 137, 138, 139, 140, 141, 72] and biological systems [142, 143, 144, 145, 136]. To show the need for experimentally-driven models, we review the challenges of channel modeling for two of these testbeds.
Example 12** (Non-biological Testbed [57])**
The first experimental MC system was presented in [57] and is based on spraying and detecting alcohol in open space. In [146], it was shown that a simple model based on diffusion and the flow generated by a fan cannot accurately explain the measurements obtained from the testbed in [57] due to system nonlinearities whose exact cause is not known. For example, the spray that is used for releasing the chemicals may not produce consistently-sized droplets in the spray stream across different trials, the flow may show turbulent behavior that is difficult to model, and the receiver sensor is prone to long response and recovery times.
Example 13** (Biological Testbed [136])**
The biological MC testbed reported in [136] converts an electrically controlled optical signal into a chemical signal. In particular, for this testbed, E. coli bacteria were genetically modified to incorporate light-driven proton pumps in their cell membranes. Upon a light stimulus, the modified bacteria then pump protons out into the environment which increases the proton concentration outside the bacteria. The resulting proton concentration was measured by a pH sensor playing the role of the receiver. Although complex models were developed in the biology literature for describing the proton release rate of proton pumps as a function of a given induced optical intensity [147, 148, 149], they typically do not account for all of the dynamics inherent to living cells. In fact, the growth, dying, and varying living conditions of the bacteria due to constant exposure to light may impact the channel model of the MC system in [136] and cannot be easily captured analytically or via simulation.
We note that similar inherent randomness and nonlinearities as discussed for the two examples above also exist for other experimental testbeds [59, 137, 138, 139, 140, 141, 72, 142, 143, 144, 145] and are challenging to model analytically or even simulate since their exact cause is unknown.
V-B2 Data-Driven Model
To address the aforementioned shortcomings of analytical and simulation models, we propose to employ data-driven models to account for the unpredictable randomness and nonlinearities of real MC systems. The basic idea behind these models is to select an appropriate parametric model and choose its corresponding parameters to fit the measurement data. In the following, we describe two different approaches for selecting a suitable parametric model.
Physically-Motivated Parametric Models: Here, the model is chosen based on physics’ first principles. For instance, in [146], a mathematical model is developed for the testbed in [57] which is based on the solution to the advection-diffusion equation with uniform flow, cf. (18). Nevertheless, the parameters of the original analytical model were modified to fit the model to the experimental data. As another example, we consider the system in [72]. The model uses magnetic nanoparticles in duct flow that are detected upon moving through a coil enclosing the duct. Here, the parametric model is based on laminar flow, cf. (14), and depends on the initial distribution of the particles released across the cross-section of the duct. The adopted parametric model was then shown to accurately model the complex advection-diffusion process in the duct after fitting its parameters to the measurement data. In general, after choosing the parametric model, standard curve fitting toolboxes can be employed to find the model parameters. One common approach is to use the parameter set that minimizes the mean square error between the model and the measurement data [146, 72, 136].
Blind Models based on Neural Networks: In the absence of an appropriate physically-motivated model, an alternative option is to employ blind models based on neural networks to jointly learn the model and its parameters [150, 151]. One suitable network architecture for this purpose is the generative adversarial network (GAN) which is able to generate a model that creates artificial data very similar to the measurement data [152]. The advantage of such blind parametric models is that they can be universally applied to general MC systems, whereas physically-motivated models have to be carefully chosen according to the MC system under consideration. On the other hand, the parameters of a physically-motivated model have physical meaning, which is not the case for the parameters of a trained neural network. The other challenge of channel modeling based on neural networks is that they typically require much more experimental data than parametric models to construct the model. This is not surprising since without domain knowledge, the number of parameters to be learned for a neural network is much larger than that for a parametric model.
V-B3 Example of an Experimentally-driven Model
In order to further familiarize the reader with the main steps of developing an experimentally-driven channel model and to highlight some peculiarities that may arise, we present the modeling methodology for the biological testbed in [136] in some detail, cf. Example 13 and Fig. 22.
Simple Physically-Motivated Parametric Model: In order to arrive at an analytical model, the following assumptions are made in [136]. It is assumed that the bacteria (i.e., the transmitter) are uniformly distributed in their container and that all bacteria are subject to the same light stimulus at the same intensity because the bacteria suspension is continuously stirred. It is further assumed that the bacteria begin and stop pumping protons (i.e., signaling molecules) instantly when the light is activated and deactivated, respectively. Furthermore, it is assumed that the bacteria take up protons in a passive manner, i.e., protons are consumed by the bacteria which lowers the measured proton concentration. Finally, it is assumed that the pH measuring device (i.e., the receiver) is passive and that its presence does not change the proton concentration or the behavior of the bacteria. These assumptions do not strictly hold but are reasonable in consideration of the size of the setup, the pumping speed of a proton pump, and the characteristics of the bacteria [136]. We note that counting the individual molecules observed at the receiver, , might be a reasonable assumption for nanomachines; however, for experimental testbeds such as [136], computing is not feasible. Hence, in [136], the proton concentration obtained from the measured pH was considered as the received signal and was modeled as
[TABLE]
where is the expected proton concentration and is a random additive noise. It was shown that follows a Gaussian distribution and this was justified using the CLT since consists of different types of noises including diffusion (counting) noise, pH sensor circuitry noise, and the noise inherent to the biological machinery of the bacteria. We assume that the light stimulates the bacteria over a time interval (i.e., a rectangular input pulse) since an impulsive stimulus (i.e., a delta input) does not effectively stimulate the bacteria to release a sufficient number of protons into the MC channel, see Remark 19. Under the aforementioned assumptions, the expected proton concentration , depending on whether the light is on () or off (), can be obtained as
[TABLE]
where is the initial concentration at starting time , is the saturation concentration, and is a time constant. The parameters of this model are , , , and , which are found using nonlinear least square error minimization to fit the measurement data. For example, in Fig. 23a), we consider constant illumination for minutes followed by darkness. The corresponding measurement signal and the fitted model using (106) are shown in Fig. 23b) using blue and black lines, respectively. As expected, the concentration increases upon illumination and decreases quickly in darkness. Nevertheless, the model in (106) fails to accurately follow the measurement data. In fact, there exists an additional persistent decreasing bias in the measurement signal which is not anticipated by the saturation model in (106).
Enhanced Parametric Model: The assumptions made to arrive at (106) do not account for the dynamics inherent to living cells. As such, cells can be growing in number or dying, or their fidelity can change. Motivated by the observations from the measurement data in Fig. 23b), it was suggested in [136] to enhance the model (106) with a simple additive linear offset as follows
[TABLE]
where is a parameter controlling the slope of the bias. The extended model is then given by where . From Fig. 23b), we can observe that the enhanced model, shown in red, fits the measurement data well. This example shows that further modification of a model that was obtained solely based on physical principles may be needed to arrive at an appropriate parametric model for an experimental system.
VI Challenges and Directions for Future Work
MC is still in its early stages of development and our understanding of MC channels is still quite limited. In the following, we review some potential challenges and open research problems which have to be addressed for successful deployment of MC systems.
Particle Generation and Signaling Pathways: Although the impact of CRNs in the physical channel (e.g. degradation reactions) on the CIR of the MC systems has been studied, see Section III-D and [58, 59, 36, 37, 60, 61, 67], less attention has been dedicated to the analysis of the influence of the CRNs at the transmitter and receiver. Such CRNs include particle generation reaction networks at the transmitter and signaling pathways at the receiver. In the MC literature, there are some preliminary works that have studied the impact of particle generation reaction networks and also simplified signaling pathways; e.g., see [85, 81]. However, the corresponding models are derived for a mesoscopic modeling approach. Analytical CIR models that take the impact of these CRNs into account are crucial for system design and hence constitute an interesting research challenge.
Turbulent Flow: The majority of the CIR models in the MC literature for advection channels have been developed based on the assumption of a uniform or laminar flow velocity field, cf. Section III-D. However, for several MC environments, such as large arteries (e.g. the aorta) and macroscale environments (e.g. oil pipe lines), flow may exhibit turbulent behavior [53]. In particular, turbulence can occur when the MC channel is non-homogeneous, e.g., due to the presence of obstacles in the physical channel. Therefore, studying and analyzing advection channels with a turbulent velocity field is an important open research problem. First results towards analyzing turbulent flow for MC systems were reported in [137].
Sample Correlation: Multiple-sample detectors are used in the MC literature to improve the detection performance [35, 32, 31, 58, 37, 46, 90, 101, 83]. It is typically assumed that different samples are statistically independent from each other. However, this assumption holds only when the sampling interval is chosen large enough such that the independence of consecutive samples is ensured. In Section IV-F1, we have numerically evaluated the correlation among consecutive samples, and in [46], the mutual information between consecutive samples is numerically evaluated. We note that sample correlation significantly depends on the type of receiver, e.g., a recurrent, non-recurrent, AMC, or IMC receiver, and its physical and chemical properties, e.g., the size, the number of receptors, and the reaction rate constants of the binding and unbinding reactions for a reactive receiver. Therefore, a careful study of sample correlation and the minimum sampling interval needed to ensure sample independence for different receiver types is essential for the applicability and performance analysis of the multiple-sample detectors proposed in the literature.
Complex Networks: In this tutorial, we focused on a single one-way communication link from a transmitter to a receiver. This is the simplest communication architecture and hence the basis for more complex network typologies. We note that although multi-node networks can often be decomposed into a superposition of individual links, there are certain scenarios where such a decomposition is invalid. For instance, if multiple reactive receivers are in the environment, the presence of each receiver will impact the signal received at any other receiver, see e.g. [21] and [22]. Moreover, in Section III, we considered MC environments with simple boundary and initial conditions in order to derive analytical channel models. However, some important MC environments, such as the cardiovascular system, are quite complex and cannot be fully modeled based on physics’ first principals. One approach is to develop simulation environments for such complex networks, see [130] and Section V-A. Nevertheless, for system design, it is desirable to have simple yet sufficiently accurate analytical models for complex multi-node networks. Developing such analytical models constitutes an important future research topic, see [153, 42, 143, 75] for some related works.
Microscale and Macroscale Models: MC systems have numerous potential applications which range from targeted drug delivery and health monitoring for microscale systems to communication in oil pipelines or chemical reactors and environmental monitoring for macroscale systems. Nevertheless, most of the current literature has targeted microscale applications and the available models are typically developed for microscale MC environments. However, macroscale and microscale MC systems may require quite different considerations. For instance, the number of molecules needed for communication at macroscale is typically much larger than that needed for communication at microscale. Moreover, while at microscale, molecules can be counted at the receivers (e.g. via ligand-receptors), at macroscale, receivers usually measure a quantity that is a function of the molecule concentration (e.g. a pH sensor was used in [59, 72] and mass spectroscopy was used in [154, 140]). In summary, the development of channel models for macroscale MC systems is an important and interesting topic for future research. We refer the interested reader to [4, 155, 154, 140, 72, 146] for preliminary results on channel modeling for macroscale MC systems.
Generally-Accepted and Experimentally-Verified Models: Over the past years, several non-biological experimental testbeds [57, 59, 137, 138, 139, 140, 141, 72] and biological experimental testbeds [142, 143, 144, 145, 136] have been developed to demonstrate MC. Most of these experimental testbeds were developed as proofs-of-concept for human-designed MC and have not led to the development of general mathematical MC channel models. However, for the advancement of MC research, it will be crucial to specify generally-accepted test channels with corresponding experimentally-verified mathematical channel models. Then, researchers in the MC community can use these established models for the design and performance analysis of newly-developed communication schemes.
VII Conclusions
This paper provided a comprehensive tutorial review of the diffusive MC channel models available in the literature. To this end, we first presented the underlying fundamental laws that govern diffusion, advection, and chemical reactions in MC channels and constitute the essential mathematical tools from biology, chemistry, and physics required for the development of MC channel models. Subsequently, we reviewed the main end-to-end channel models reported in the diffusive MC literature and showed how they were developed from basic physical principles. The reviewed end-to-end channel models included the joint effects of release mechanisms, the physical channel, and reception mechanisms. Moreover, we provided a unified definition for the received signal that included the representation obtained by both timing and counting receivers as special cases. Furthermore, for counting receivers, we derived signal models relevant for different time scales. We generalized these models to account for interfering noise molecules and ISI and studied the correlation among the received signals observed at different time scales. In addition, simulation-driven and experimentally-driven channel models were investigated for complex scenarios where simple MC channel models cannot be obtained from basic physical principles. Finally, we provided a discussion of challenges, open research problems, and future directions for channel modeling of diffusive MC systems.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] S. Hiyama, Y. Moritani, T. Suda, R. Egashira, A. Enomoto, M. J. Moore, and T. Nakano, “Molecular Communication,” in Proc. NSTI Nanotech , May 2005, pp. 391–394.
- 2[2] I. Akyildiz, F. Brunetti, and C. Blazquez, “Nanonetworks: A New Communication Paradigm,” Comput. Net. , vol. 52, pp. 2260–2279, Apr. 2008.
- 3[3] N. Farsad, H. Yilmaz, A. Eckford, C. Chae, and W. Guo, “A Comprehensive Survey of Recent Advancements in Molecular Communication,” IEEE Commun. Surveys Tutorials , vol. 18, no. 3, pp. 1887–1919, third quarter 2016.
- 4[4] L. P. Giné and I. F. Akyildiz, “Molecular Communication Options for Long Range Nanonetworks,” Computer Netw. , vol. 53, no. 16, pp. 2753 – 2766, 2009.
- 5[5] M. Pierobon and I. Akyildiz, “Noise Analysis in Ligand-Binding Reception for Molecular Communication in Nanonetworks,” IEEE Trans. Sig. Process. , vol. 59, no. 9, pp. 4168–4182, Sep. 2011.
- 6[6] S. Kadloor, R. Adve, and A. Eckford, “Molecular Communication Using Brownian Motion With Drift,” IEEE Trans. Nano Biosci. , vol. 11, no. 2, pp. 89–99, Jun. 2012.
- 7[7] T. Nakano, T. Suda, T. Koujin, T. Haraguchi, and Y. Hiraoka, Molecular Communication Through Gap Junction Channels . Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, pp. 81–99.
- 8[8] M. Moore, A. Enomoto, T. Nakano, R. Egashira, T. Suda, A. Kayasuga, H. Kojima, H. Sakakibara, and K. Oiwa, “A Design of a Molecular Communication System for Nanomachines Using Molecular Motors,” in Proc. Pervasive Comput. Commun. Workshops , Mar. 2006, pp. 6 pp.–559.