Using spatial partitioning to reduce the bit error rate of diffusion-based molecular communications

TL;DR

This paper extends demodulation techniques for diffusion-based molecular communications to spatially partitioned receivers, demonstrating that such partitioning can effectively reduce the bit-error rate in these systems.

Contribution

It introduces a new ODE-based demodulation method for multi-voxel receivers and shows how spatial partitioning improves communication reliability.

Findings

Spatial partitioning reduces bit-error rate.

Extended ODE for multi-voxel receiver demodulation.

Analysis of different receiver volumes and shapes.

Abstract

This work builds on our earlier work on designing demodulators for diffusion-based molecular communications using a Markovian approach. The demodulation filters take the form of an ordinary differential equation (ODE) which computes the log-posteriori probability of observing a transmission symbol given the continuous history of receptor activities. A limitation of our earlier work is that the receiver is assumed to be a small cubic volume called a voxel. In this work, we extend the maximum a-posteriori demodulation to the case where the receiver may consist of multiple voxels and derive the ODE for log-posteriori probability calculation. This extension allows us to study receiver behaviour of different volumes and shapes. In particular, it also allows us to consider spatially partitioned receivers where the chemicals in the receiver are not allowed to mix. The key result of this paper is that spatial partitioning can be used to reduce bit-error rate in diffusion-based molecular communications.

Tables

Table 1. TABLE I : Notation and Chemical Symbols

Symbols	Notation and Value
$w$	Dimension of one side of a voxel
$D$	diffusion coefficient of signalling molecules
$d$	Inter-voxel diffusion constant for signalling molecules
$d_r$	Inter-voxel diffusion constant for \ceeX^*
$g_{+}$	Propensity function constant for the activation reaction (1a)
$g_{-}$	Propensity function constant for deactivation reaction (1)
$P$	Number of receiver voxels
\ceeS, \ceeX, \ceeX^*	Chemical symbols for signalling molecule, inactivate receptor and active receptor
$X_p\left(t\right)$ ($X_p^{\ast }\left(t\right)$)	Number of active (inactive) receptors in the $p$-th receiver voxel at time $t$
$N_{R,p}\left(t\right)$	Number of signalling molecules in the $p$-th receiver voxel at time $t$

Equations

$$\frac{d{L}_k(t)}{dt}=$$

$$\frac{d{Z}_k(t)}{dt}={\left[\frac{d{X}^{\ast }(t)}{dt}\right]}_{+}\log ({\sigma }_k(t))-g_{+}(M-X^{\ast }(t)){\sigma }_k(t)$$

$$\frac{d{L}_k(t)}{dt}=$$

$$\frac{Z_k(t)}{dt}=$$

$$Z_{k,p}(t)=$$

$$\frac{d{L}_k(t)}{dt}=$$

$$\frac{Z_k(t)}{dt}\approx$$

$$\frac{d{\bar{x}}_1^{\ast }(t)}{dt}=$$

$$\frac{d{\beta }_{k,p}(t)}{dt}=$$

$$\frac{d{\tilde{x}}_1^{\ast }(t)}{dt}=$$

$$\sqrt{{\hat{k}}_{+}{\bar{n}}_{R,1}(t){\bar{x}}_1(t)}{\gamma }_{\mathrm{ar},rv1}(t)+\sqrt{k_{-}{\bar{x}}_1^{\ast }(t)}{\gamma }_{dr,RV1}(t)+\sqrt{d_r{\bar{x}}_1^{\ast }(t)}{\gamma }_{d,X\ast ,1to2}(t)+\sqrt{d_r{\bar{x}}_2^{\ast }(t)}{\gamma }_{d,X\ast ,2to1}(t)$$

$${\left[\frac{d{X}_p^{\ast }(t)}{dt}\right]}_{+}=$$

$$\sqrt{{\hat{k}}_{+}{\bar{n}}_{R,1}(t){\bar{x}}_1(t)}{\gamma }_{\mathrm{ar},rv1}(t)+\sqrt{d_r{\bar{x}}_2^{\ast }(t)}{\gamma }_{d,X\ast ,2to1}(t).)$$

$$\mathbf{P}[X_1^{\ast }(t+\Delta t),X_2^{\ast }(t+\Delta t)|k,{\mathcal{X}}_1^{\ast }(t),{\mathcal{X}}_2^{\ast }(t)]=$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t)+1,X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t))Q_{1,a}+$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t)-1,X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t))Q_{1,d}+$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t),X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t)+1)Q_{2,a}+$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t),X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t)-1)Q_{2,d}+$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t),X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t))Q_0$$

$$Q_{1,a}=g_{+}(M_1-X_1^{\ast }(t))E[N_{R,1}(t)|k,{\mathcal{X}}_1^{\ast }(t),{\mathcal{X}}_2^{\ast }(t)]\Delta t$$

$$Q_{1,d}=g_{-}{X}_1^{\ast }(t)\Delta t$$

$$Q_{2,a}=g_{+}(M_1-X_1^{\ast }(t))E[N_{R,2}(t)|k,{\mathcal{X}}_1^{\ast }(t),{\mathcal{X}}_2^{\ast }(t)]\Delta t$$

$$Q_{2,d}=g_{-}{X}_2^{\ast }(t)\Delta t$$

$$Q_0=1-(Q_{1,a}+Q_{1,d}+Q_{2,a}+Q_{2,d})$$

$$\frac{d{L}_k(t)}{dt}=$$

$$\frac{d{L}_k(t)}{dt}=\underset{\Delta t\to 0}{\lim }\frac{\log (\mathbf{P}[X_P^{\ast }(t+\Delta t)|k,{\mathcal{X}}_P^{\ast }(t)])}{\Delta t}$$

$$\frac{d{L}_k(t)}{dt}=$$

$$g_{+}(M_p-X_P^{\ast }(t))\mathbf{E}[N_p(t)|k,{\mathcal{X}}_R^{\ast }(t)])$$

$$\mathbf{P}[X_1^{\ast }(t+\Delta t),X_2^{\ast }(t+\Delta t)|k,{\mathcal{X}}_1^{\ast }(t),{\mathcal{X}}_2^{\ast }(t)]=$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t)+1,X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t))Q_{1,a}+$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t)-1,X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t))Q_{1,d}+$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t),X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t)+1)Q_{2,a}+$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t),X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t)-1)Q_{2,d}+$$

$$\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t)-1,X_2^{\ast }(t+\Delta t)=X_2^{\ast }(t)+1)$$

$$Q_{1\to 2}+\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t)+1,X_2^{\ast }(t+\Delta t)=$$

$$X_2^{\ast }(t)-1)Q_{2\to 1}+\delta (X_1^{\ast }(t+\Delta t)=X_1^{\ast }(t),X_2^{\ast }(t+\Delta t)=$$

$$X_2^{\ast }(t))Q_0$$

$$Q_{1,a}=g_{+}E[X_1(t)N_{R1}(t)|k,{\mathcal{X}}_1^{\ast },{\mathcal{X}}_2^{\ast }(t)]\Delta t$$

$$Q_{1,d}=g_{-}{X}_1^{\ast }(t)\Delta t$$

$$Q_{2,a}=g_{+}E[X_2(t)N_{R2}(t)|k,{\mathcal{X}}_1^{\ast },{\mathcal{X}}_2^{\ast }(t)]\Delta t$$

$$Q_{2,d}=g_{-}{X}_2^{\ast }(t)\Delta t$$

$$Q_{1\to 2}=d_r{X}_1^{\ast }(t)\Delta t$$

$$Q_{2\to 1}=d_r{X}_2^{\ast }(t)\Delta t$$

Full text

I Introduction

Molecular communication is an emerging field which focuses on realizing communication between nano-scale devices [1, 2] and especially the internet of bio-nano things [3]. A key characteristic of molecular communication is the use of molecules as the information carrier. This paper considers diffusion-based molecular communications [4][5].

The receiver is an important component in any communications system. We can divide the techniques for improving receiver performance of diffusion-based molecular communications into two categories. The first category of work uses signal processing techniques. Some examples in this category are: the paper [6] designs a receiver based on minimum mean square error method; [7] uses multiple samples per symbol and maximum likelihood method to design a receiver; [8] designs a matched filter to maximise the signal-to-interference-plus-noise ratio; [9] studies the design of estimate-and-forward relay nodes; and, [10] uses data fusion of decisions from multiple receivers.

The second category of work uses physical or chemical properties to improve communications performance. Some examples are: [11] uses enzyme in the medium to reduce interference; [12] uses flow to improve the performance of the weighted sum detector at the receiver; [13] uses chemical reactions in the transmission medium to reduce interference.

This paper proposes to use the physical mechanism of spatial partitioning to improve communication performance. Our work is inspired by the fact that receptors on the cell membrane are organized into spatially separated clusters of receptors [14]. In this paper, we assume that the receiver uses receptors which can be activated by the signalling molecules. We propose to segregate these receptors into a number of clusters separated spatially. We consider two configurations: partitioned and mixed. In the partitioned configurations, spatial isolation is perfect and receptors cannot move between clusters; however, in mixed configurations, which can be considered as imperfect isolation, receptors can diffuse between the clusters. For both configurations, we derive the maximum-a-posteriori (MAP) demodulator by leveraging our earlier work on designing demodulators for diffusion-based molecular communications using a Markovian approach [15][16].

This paper also removes a limitation of our earlier work in [15][16] where the receiver is assumed to be a small cubic volume called a voxel. In this paper, we assume that the receiver consists of multiple voxels and the union of these voxels define the shape of the receiver. This method of modelling a 3-dimensional volume is similar to that in finite difference method. In this paper, we extend the method of demodulator design in [15][16] to the multi-voxel case.

The contributions of this paper are:

• •

We propose to use spatial partitioning to reduce the bit error rate (BER) of diffusion-based molecular communications.

• •

We derive the MAP demodulator for both the partitioned and mixed configurations. This derivation of the MAP demodulator applies to a receiver with multiple voxels.

• •

We derive approximations of the derived MAP demodulators. We show that the approximate demodulator for the partitioned case has a lower BER than that of the mixed configuration. We also show that our approximate demodulator for the mixed configuration can offer gradual degradation in BER from the partitioned configuration which means that it can deal with imperfect partitioning.

• •

We derive a method to analytically estimate the BER of the approximate demodulator.

The rest of the paper is organized as follows. Section II discusses the related work. In Section III we present the background and a summary of our previous work in [15][16]. Section IV presents the MAP demodulator for the partitioned and mixed configurations. Section V presents an analytical method to calculate the BER of the approximate demodulators. Section VI presents numerical studies on comparing the performance of the partitioned and mixed configurations. Finally, Section VII concludes the paper.

II Related work

A number of insightful surveys have been written on the topic of molecular communications, see [1, 2, 17, 18, 4].

There are three main components of a molecular communications system, namely transmitter, propagation medium and receiver. For molecular communications transmitters, many different modulation methods have been proposed. These include Molecule Shift Keying, Frequency Shift Keying, on-off keying and Concentration Shift Keying [19, 20, 21]. In this paper, our transmitter uses Reaction Shift Keying (RSK) which is proposed in our earlier work [16], [15], [22] [23]. In RSK, the transmission symbols are generated by a set of chemical reactions and are characterised by the time-varying concentration of the signalling molecules.

For the modelling of the propagation medium, various models have been proposed in literature. A common model is to assume that space is continuous, see e.g. [24, 25, 26, 11, 13]. In this paper, we assume that the medium is divided into voxels [27, 28, 29]. The voxel based setting allows us to model the molecular communication networks, which consist of both diffusion and chemical reactions, by using the reaction-diffusion master equation (RDME) [30].

In Section I, we discussed two categories of work on improving receiver performance. Another method to classify different receivers that have been used in molecular communication is: those that use chemical reactions at the receiver, e.g. [31, 32, 33, 34, 35], and those that do not, e.g. [36, 37, 38, 21]. In this work, we assume that the receiver uses chemical reactions at the receiver.

The demodulation problem considered in this paper can be considered as using data fusion to combine the measurements coming from multiple voxels in order to carry out demodulation. Data fusion has previously been considered in molecular communications in [10]. The paper [10] assumes that there are multiple receivers; each receiver makes a hard decision for ON-OFF keying and a fusion centre makes an overall decision based on the decisions from the receivers. This paper uses a different method for data fusion. For our proposed demodulator, each voxel computes an approximate log-posteriori probability and we add these probabilities up to obtain the log-posteriori probability of the receiver.

In this paper, we consider using spatial partitioning to reduce the BER in molecular communications. Spatial partitioning has been studied in biophysics, e.g. the authors of [14] study how spatial partitioning can reduce noise, and the paper [39] studies the impact of receptor clustering on the noise in cell signalling. In the area of chemical based computation, the authors of [40] use partitioning to improve the performance of a molecular computing system.

In our previous paper [41] we show that the partitioning can be used to reduce the signal variance in the receiver. However, no receiver mechanisms have been proposed. This paper derives the MAP demodulators for both the partitioned and mixed configurations.

III Background and summary of previous Work

This section provides a summary of our previous work in [15, 16] on using a Markovian approach to derive a MAP demodulator. This section is divided into two subsections. Section III-A presents the modelling framework while Section III-B presents the MAP demodulator.

III-A Modelling framework

In [15, 16], we consider a molecular communication system with a transmitter and a receiver inside a fluid propagation medium. Fig. 2 shows an overview of the system. We assume that the transmitter and the receiver communicate with one type of signalling molecules denoted by \ceeS. We will now discuss the modelling of the system components in further details.

III-A1 Propagation medium

Signalling molecules diffuse from transmitter to receiver through a propagation medium. We model the propagation medium as a rectangular prism. We divide the medium into voxels. This is the same as applying a finite difference spatial discretization to a volume. Fig. 2 shows a 2-dimensional projection of a medium consisting of 4-by-4-by-1 voxels. In order to facilitate the description later on, we index each voxel in this example by using integers 1,..,16 which are shown in the top-right corner of each voxel. We model the diffusion of the signalling molecules by using spatially discrete jumps between neighbouring voxels. For example, a signalling molecule in Voxel 1 in Fig. 2 can diffuse to any of its neighbouring voxels, which are Voxels 2 and 5. The double headed arrows in Fig. 2 show the direction of diffusion. We assume the diffusion of the signalling molecules are independent each other.

We assume the medium is homogeneous with a constant diffusion coefficient $D$ for the signalling molecules. By applying a finite difference discretization to the 3-dimensional diffusion equation [30], it can be shown that within an infinitesimal time $\Delta t$, the probability that a signalling molecule will jump from a voxel to a neighbouring voxel is $\frac{D}{w^{2}}\;\Delta t$ where $w$ denotes the length a voxel edge.

Lastly, our model can be used to model two types of boundary conditions: reflecting boundary condition where signalling molecules are not allowed to leave the medium; or, absorbing boundary condition where signalling molecules may leave the medium forever, e.g. the single headed arrow in Voxel 4 in Fig. 2 shows that signalling molecules may leave the medium.

III-A2 Transmitter

The transmitter is assumed to occupy only one voxel. The role of the transmitter is to produce the signalling molecules for the transmission symbols. We assume that the transmitter uses RSK. This means that each symbol corresponds to a time-varying concentration profile of signalling molecules produced by a set of chemical reactions, Fig. 2 depicts a transmitter which can produce $K$ different transmission symbols. Once the signalling molecules have been produced by the transmitter, they are free to diffuse in the medium.

III-A3 Receiver

In [15, 16], we assume that the receiver occupies one voxel, which is an assumption which we will change in Section IV of this paper. In this section, we follow the assumption in [15, 16].

We divide the operation of the receiver into two blocks, which we will refer to as the front-end and back-end blocks, see Fig. 2. The front-end block consists of a molecular circuit, which is another name for a set of chemical reactions. In [15], the front-end block is assumed to be a reversible ligand-receptor binding while in [16], the front-end block can be any molecular circuit. For example, the following molecular circuit is considered in [42]:

[TABLE]

where $g_{+}$ and $g_{-}$ are propensity function constants; and \ceeX and \ceeX^* are, respectively, the inactive and active forms of a species. Reaction (1a) is an activation reaction where the signalling molecule \ceeS turns inactive \ceeX into active \ceeX^*, while Reaction (1) is a deactivation reaction.

Note that Table I contains a list of commonly used constants and chemical symbols for easy referral.

One or more chemical species in the front-end molecular circuit are chosen as the output species. The molecular counts of the output species over time are the output signals of the front-end, which are fed into the back-end as the input signals, see Fig. 2. The back-end of the receiver is the demodulator whose aim is to infer the symbol that the transmitter has sent by using the input signals to the demodulator, i.e. the molecular counts of the output species of the front-end molecular circuit.

III-B MAP demodulator

This section summarises the steps of deriving the MAP demodulator in our Markovian framework [15, 16]. For illustration, we assume that the front-end molecular circuit of the receiver is given by Reactions (1). We further assume that the molecules \ceeX and \ceeX^* can only be found in the receiver voxels, and they are uniformly distributed in the voxel. We also assume that the total number of \ceeX and \ceeX^* molecules is a constant $M$. We designate \ceeX^* as the output species. Let $X(t)$ and $X^{*}(t)$ denote, respectively, the number of $X$ and $X^{*}$ molecules at time $t$. Note that both $X(t)$ and $X^{*}(t)$ are piece-wise constant because they are molecular counts.

We model the transmitter, medium and receiver front-end by using a RDME111Note that RDME assumes that space is discretised into voxels and time is continuous, therefore RDME is compatible with the discretisation of the medium discussed in Section III-A1. There is another rationale why we choose RDME. For the design of molecular communications systems, we need a stochastic model that can model systems with both diffusion and chemical reactions. There are three main classes of such models: Smoluchowski equation [43], RDME and the Langevin equation [44]. The Smoluchowski equation is based on particle dynamics. It is a fine grained model but hard to work with analytically. Both RDME and Langevin are easier to work with analytically but master equation has a finer scale and granularity compared to the Langevin equation [45]. Therefore we choose to use RDME which allows us to use the Markovian theory for analysis and is at the same time a finer grained model. Note that there is some recent work in combining the voxel-based approach (also known as the mesoscopic approach) and the Smoluchowski equation (also known as the microscopic approach) in simulating systems with both reactions and diffusion, see [46] and [47].. Note that RDME is a specific type of continuous-time Markov process (CTMP). This means the signal $X^{*}(t)$ is a realization of a CTMP. The RDME assumes that the species in each voxel is well-mixed. For RDME, the position of each molecule is only known to the spatial scale of a voxel; in other words, the exact position (say, in terms of $x$, $y$, $z$ co-ordinates) of each molecules is not known. This means that for RDME, the state of a molecular system is characterised by using the count of each possible type of molecules in each voxel.

Since \ceeX^* has been designated as the output species of the front-end molecular circuit, the signal $X^{*}(t)$ is available to the demodulator. For the derivation of the demodulation filter, we assume that at time $t$, the data available to the demodulation filter is $X^{*}(\tau)$ for all $\tau\in[0,t]$. We use ${\cal X}^{*}(t)$ to denote the continuous-time history of $X^{*}(t)$ up to time $t$. We use the Bayesian framework for demodulation. Let ${\mathbf{P}}[k|{\cal X^{*}}(t)]$ denote the posteriori probability that symbol $k$ has been sent given the history ${\cal X^{*}}(t)$. Instead of working with ${\mathbf{P}}[k|{\cal X^{*}}(t)]$, we will work with its logarithm. Let $L_{k}(t)=\log({\mathbf{P}}[k|{\cal X^{*}}(t)])$. Note that the log-posteriori probability diverges in continuous time, but we are able to compute a shifted version of it. For ease of reference, we will simply refer to the shifted version as the log-posteriori probability, We can use the method in [16] to show that we can compute $L_{k}(t)$ by using the following ordinary differential equation (ODE):

[TABLE]

where $[\xi]_{+}=\max(\xi,0)$ and ${\mathbf{E}}[N_{R}(t)|k,{\cal X^{*}}(t)]$ is the estimation of the mean number of signalling molecules in the receiver voxel. Recall that $X^{*}(t)$ is the count of \ceeX^* molecules at time $t$. It is therefore a piecewise constant signal and the derivative $\left[\frac{dX^{*}(t)}{dt}\right]_{+}$ is a series of Dirac deltas situated at the times at which the count of \ceeX^* changes.

The computation of ${\mathbf{E}}[N_{R}(t)|k,{\cal X^{*}}(t)]$ requires extensive computation because it requires an optimal Bayesian filtering problem to be solved222 In order to solve the optimal Bayesian filtering problem, we need to track the probability of each possible system state. For our problem, the system state is the count of each possible type of molecule in each voxel. The total number of possible states grows rapidly due to combinatorial explosion. This explanation also applies to subsequent problems considered in this paper.. In [15], we propose to replace the hard-to-compute ${\mathbf{E}}[N_{R}(t)|k,{\cal X^{*}}(t)]$ by ${\mathbf{E}}[N_{R}(t)|k]$, which is the mean number of signalling molecules if Symbol $k$ is transmitted. Let $\sigma_{k}(t)$ denote ${\mathbf{E}}[N_{R}(t)|k]$. With the proposed replacement, Eq. (1b) becomes:

[TABLE]

where $Z_{k}(t)$ is an approximation of $L_{k}(t)$. We can interpret Eq. (1n) as using $\sigma_{k}(t)$ as an internal model or prior knowledge for demodulation. The use of internal model is fairly common in signal processing and communications, e.g. in matched filtering.

Note that Eq. (1b) is the optimal solution for the demodulation problem. The replacement of ${\mathbf{E}}[N_{R}(t)|k,{\cal X^{*}}(t)]$ by ${\mathbf{E}}[N_{R}(t)|k]$ means that Eq. (1n) is a approximate solution. We show in [15] that $Z_{k}(t)$ from Eq. (1n) approximates $L_{k}(t)$ in Eq. (1b) well.

To make the decision at time $t$, the demodulator decide symbol $\hat{k}$ is transmitted if $\hat{k}={\arg\max}_{k=0,...,K-1}Z_{k}(t)$. Also, $Z_{k}(0)$ is initialized to the logarithm of the prior probability that the transmitter sends Symbol $k$. The structure of the demodulator is depicted in Fig. 4. The demodulator consists of $K$ filters computing $Z_{k}(t)$ and a maximum block to determine the transmission symbol.

IV Receiver with multiple voxels

This section derives the demodulator when the receiver consists of multiple voxels. We will consider both the partitioned and mixed configurations.

IV-A Multi-voxel receiver demodulation problem

We assume the receiver consists of $P$ voxels where $P>1$. Without loss of generality, we assume these $P$ voxels form a connected volume. We will use $p$ to index the voxels in the receiver, where $p=1,\ldots,P$. We will explain the demodulation problem using an example.

Example 1

In this example, we assume the molecular circuit (1) is present in all the $P$ receiver voxels. We will continue to use the RDME to model the system. The RDME framework requires us to distinguish between the \ceeS, \ceeX and \ceeX^ in different voxels. We will use \ceeS_p, \ceeX_p and \ceeX^_P to denote the \ceeS, \ceeX and \ceeX^ in the $p$-th receiver voxel where $p=1,...,P$. The $P$ molecular circuits in the receiver voxels are:*

[TABLE]

We designate \ceeX^_P, for $p=1,\ldots,P$, as the output species of the receiver. Let $X_{p}(t)$ and $X^{*}_{p}(t)$ denote respectively the number of \ceeX_p and \ceeX^_P molecules at time $t$.

Let also ${\cal X}^{*}_{P}(t)$ be the continuous history of the species \ceeX^_P up to and including time $t$. Our goal is to determine the posteriori probability ${\bf P}[k|{\cal X}^{*}_{1}(t),\ldots,{\cal X}^{*}_{P}(t)]$ that the $k$-th transmitter symbol has been sent given the histories ${\cal X}^{*}_{1}(t)$, $\ldots$, ${\cal X}^{*}_{P}(t)$. $\Box$*

Note that it is possible to generalise the above examples in a few different ways: we can use molecular circuits other than (1); we can use different molecular circuits in different receiver voxels; we can use different choices of output species. We note that the methodology that we have developed to derive the demodulator can deal with all these generalisations.

We make the assumption that the species in the front-end molecular circuits (e.g. \ceeX and \ceeX^ in (1)) can only be found in the receiver. Also these species are confined within the receiver and cannot get into the propagation environment. In the text below, we will to refer to them as the receiver species.*

IV-B Mixed and partitioned configurations

Inspired by the study **[14]** which shows that spatial partitioning of the receptors on cell membrane can be used to reduce noise, we will consider two different configurations for the receiver. For illustration, we assume the receiver consists of 4 voxels arranged in 2-by-2 configuration as shown in Fig. 4.

In the first configuration, which is illustrated in Fig. 4a, a receiver species is allowed to diffuse from a receiver voxel to another receiver voxel. We will refer to this configuration as mixed as the receiver species is allowed to mix among the receiver voxels.

In the second configuration, which is illustrated in Fig. 4b, a receiver species cannot move from a receiver voxel to another. In this case, we can assume the receiver voxel wall is a selective membrane which prevent certain species from moving between voxels. We will refer to this configuration as partitioned.

We note that the voxel framework is well suited to model the partitioned and mixed configurations. This is because, for the voxel framework, it is possible to choose the value of the diffusion coefficient of a species between any two voxels. If a species is not allowed to pass between the interface of two voxels, then its diffusion coefficient for that interface is zero.

We now discuss how the two configurations can be realised. We assume that the cross-sectional area of the signalling molecules is small while that of the receiver species is large. We assume that the boundaries between receiver voxels are made of membrane with pores on it. If we choose the size of the pores to be large enough for the signalling molecules to pass through unhindered but small enough to prevent the receiver species from passing through, then we obtain the partitioned configuration. If we increase the size of the pores to allow the receiver species to pass through them with a non-zero probability, then we get the mixed configuration. A possible physical implementation is based on structural DNA nanotechnology **[48, 49]** which can create 3-dimensional wireframes with pores of different shapes and dimension, e.g. **[50]** demonstrated a design with pores down to approximately 6nm by 6nm. In nature, many signalling molecules and receptors are protein. The cross-sectional area of large protein molecules can be 10nm or more, while that of small protein molecules is a few nanometres **[51]**.

We will now derive the MAP demodulator for the partitioned configuration.

IV-C MAP demodulator for the partitioned configuration

In this section, we will derive the MAP demodulator for the demodulation problem illustrated in Example 1. In particular, in this section, we consider the partitioned configuration where the species \ceeX and \ceeX^ are not allowed to move among the receiver voxels. This means the total number of \ceeX and \ceeX^* molecules in a receiver voxel remains constant. We will use $M_{p}$ to denote the total number of \ceeX and \ceeX^* molecules in the $p$-th receiver voxel.*

The aim of the demodulation problem is to determine the posteriori probability ${\bf P}[k|{\cal X}^{*}_{1}(t),\ldots,{\cal X}^{*}_{P}(t)]$ that the $k$-th transmitter symbol has been sent given the histories ${\cal X}^{*}_{1}(t)$, $\ldots$, ${\cal X}^{*}_{P}(t)$. A key step in the derivation is to compute the probability ${\bf P}[X^{*}_{1}(t+\Delta t),X^{*}_{2}(t+\Delta t),\ldots,X^{*}_{P}(t+\Delta t)|k,{\cal X}^{*}_{1}(t),\ldots,{\cal X}^{*}_{P}(t)]$ which predicts the counts of the output species based on their histories. This is an optimal Bayesian filtering problem. The structure of this filtering problem is identical to the one considered in **[16]**. In **[16]**, we consider a filtering problem where the state vector evolves according to a CTMP. In that problem, only some elements of the state vector can be observed and the aim is to predict the future values of the observable elements from their past histories. We can therefore apply the method in **[16]** to the demodulation problem considered in this paper. We remark that although the presentation in **[16]** states that the output species come from one voxel, it is the structure of the filtering problem that really matters; in other words, it is not important whether the output species come from one or many voxels.

In order to simplify the notation, we will use the shorthand ${\cal X}^{*}_{R}(t)$ to denote the histories ${\cal X}^{*}_{1}(t)$, $\ldots$, ${\cal X}^{*}_{P}(t)$ from all the receiver voxels; note that the subscript R in ${\cal X}^{*}_{R}(t)$ is short for receiver. Let $L_{k}(t)=\log({\bf P}[k|{\cal X}^{*}_{R}(t)]$ be the log-posteriori probability that the $k$-th symbol has been sent given the histories. In Appendix A, we show that $L_{k}(t)$ evolves according to:

[TABLE]

where ${\mathbf{E}}[N_{R,p}(t)|k,{\cal X}^{*}_{R}(t)]$ is the estimated number of signalling molecules in the $p$-th receiver voxel given $k$ and the histories ${\cal X}^{*}_{R}(t)$.

Although the term ${\mathbf{E}}[N_{R,p}(t)|k,{\cal X}^{*}_{R}(t)]$ can be computed by solving an optimal Bayesian filtering problem, its computational complexity is high. Another issue, which is found in the multi-voxel receiver case but not the single-voxel receiver case, is that the estimation of the mean number of signalling molecule in the $p$-th receiver voxel via ${\mathbf{E}}[N_{R,p}(t)|k,{\cal X}^{*}_{R}(t)]$ requires the history in the $p$-th receiver voxel as well as other receiver voxels. In other words, the computation of ${\mathbf{E}}[N_{R,p}(t)|k,{\cal X}^{*}_{R}(t)]$ requires communications between receiver voxels, which adds burden to the receiver. We follow our earlier work **[15]**[16]** and replace ${\mathbf{E}}[N_{R,p}(t)|k,{\cal X}^{*}_{R}(t)]$ with the prior knowledge ${\mathbf{E}}[N_{R,p}(t)|k]$, which will be denoted by $\alpha_{p}(t)$. We will use a numerical example in Section VI-B to demonstrate that this approximation is accurate. With this replacement, an approximate demodulator is:

[TABLE]

We can interpret $Z_{k,p}(t)$ in Eq. (1or) as the approximate log-posteriori probability computed using the measurements $X^{*}_{p}(t)$ from the $p$-the voxel. These approximate log-probabilities are then added up or fused in Eq. (1oq) to obtain the approximate log-posteriori probability of the receiver.

IV-D MAP demodulator for the mixed configuration

This section considers the demodulation for the mixed configuration. For the mixed configuration, we assume that the \ceeX and \ceeX^ molecules are allowed to diffuse between the receiver voxels but they are not allowed to leave the receiver. Let $D_{r}$ be the diffusion coefficient of \ceeX^* in the receiver voxels. Note that we do not specify the diffusion coefficient of \ceeX because it does not enter the ODE for computing log-posteriori probability.*

If \ceeX and \ceeX^ are not allowed to diffuse, i.e. for the partitioned case as discussed in Section IV-B, then $X_{p}(t)+X^{*}_{P}(t)$ is a constant for all receiver voxels $p=1,\ldots,P$. However, this condition no longer holds for the mixed configuration. This means that $M_{p}$ are not parameters of the mixed configuration.*

The aim of the demodulation problem is to determine the posteriori probability ${\bf P}[k|{\cal X}^{*}_{1}(t),\ldots,{\cal X}^{*}_{P}(t)]$ that the $k$-th transmitter symbol has been sent given the histories ${\cal X}^{*}_{1}(t)$, $\ldots$, ${\cal X}^{*}_{P}(t)$. We again use ${\cal X}^{*}_{R}(t)$ to denote the histories.

In Appendix B, we derive the optimal demodulation filter:

[TABLE]

where $A_{p}(t)$ is a sequence of Dirac deltas situated at the times where activation reactions (1a) occur in voxel $p$ according to the given histories and ${\mathbf{E}}[X_{p}(t)N_{R,p}(t)|k,{\cal X}^{*}_{R}(t))]$ is the estimated mean of the product of $X_{p}(t)$ and $N_{R,p}(t)$ by using the histories.

We now consider how the terms in (1os) can be approximated. We first use an example to point out that $A_{p}(t)$ is different from $\left[\frac{dX^{*}_{p}(t)}{dt}\right]_{+}$. We illustrate this by assuming that the receiver consists of two neighbouring receiver voxels, which we will refer to as RV1 and RV2. Fig. 5 shows two sample trajectories for $X^{*}_{1}(t)$ and $X^{*}_{2}(t)$, which are of the number of \ceeX^ molecules in RV1 and RV2. The figure also depicts the signals $A_{1}(t)$ and $\left[\frac{dX^{*}_{1}(t)}{dt}\right]_{+}$ for RV1. At time $t_{1}$, an activation reaction occurs in RV1 and $X^{*}_{1}(t)$ is increased by 1 and therefore a Dirac delta is situated at time $t_{1}$ in both $A_{p}(t)$ and $\left[\frac{dX^{*}_{P}(t)}{dt}\right]_{+}$. We remark that $X^{*}_{2}(t)$ remains unchanged at $t_{1}$, which is a fact that we will use later. At time $t_{2}$, an \ceeX^* molecule moves from RV2 to RV1. This means $X^{*}_{1}(t)$ is increased by 1 while $X^{*}_{2}(t)$ is decreased by 1. We note that at time $t_{2}$, there is a Dirac delta in $\left[\frac{dX^{*}_{1}(t)}{dt}\right]_{+}$ but not $A_{1}(t)$.*

The above example shows that Dirac deltas in $A_{1}(t)$ occur at the times at which $X^{*}_{1}(t)$ is increased by 1 while $X^{*}_{2}(t)$ remains the same. This means the computation of $A_{1}(t)$ requires both the trajectories of $X^{*}_{1}(t)$ and $X^{*}_{2}(t)$. In general, the computation of $A_{p}(t)$ requires $X^{*}_{p}(t)$ as well as the histories of \ceeX^ in the neighbouring voxels. This implies that the computation $A_{p}(t)$ requires communications between neighbouring voxels, which is not realistic. An alternative is to change the activation reaction (1a) so that every time when an activation reaction occurs, an extra molecule is produced to indicate its occurrence. However, there are a number of open issues with this approach: how to ensure this extra molecule is local both in time and space, and to make sure the continual supply of this extra molecule. In this paper, we will approximate $A_{p}(t)$ by $\left[\frac{dX^{*}_{P}(t)}{dt}\right]_{+}$. We are conscious that $\left[\frac{dX^{*}_{P}(t)}{dt}\right]_{+}$ has more Dirac deltas than $A_{p}(t)$ and these extra Dirac deltas are due to the movement of \ceeX^* into receiver voxel $p$ from its neighbouring receiver voxels. We will demonstrate in Section VI that if \ceeX^* diffuses slowly, then the error is small.*

In order to avoid solving an optimal Bayesian filtering problem, we will replace ${\mathbf{E}}[X_{p}(t)N_{R,p}(t)|k,{\cal X}^{*}_{R}(t))]$ in (1os) by ${\mathbf{E}}[X_{p}(t)N_{R,p}(t)|k]$, which is denoted by $\beta_{k,p}(t)$. This replacement also removes the need for the voxels to communicate with each other. We will show that this approximation is accurate in Section VI-E. With these approximations, Eq. (1os) becomes:

[TABLE]

This will be the approximation demodulation filter for the $k$-th symbol for the mixed configuration.

Remark 1

Although we have presented the results assuming that the molecular circuit is given by Reactions (1). Our methodology can be used for any molecular circuit. The key to deriving the demodulator is to compute the counterpart of the probability ${\bf P}[X^{*}_{1}(t+\Delta t),X^{*}_{2}(t+\Delta t),\ldots,X^{*}_{P}(t+\Delta t)|k,{\cal X}^{*}_{1}(t),\ldots,{\cal X}^{*}_{P}(t)]$ when different molecular circuits are used. In words, this probability can be stated as: $\mathbf{P}[$ counts of all the output species from all receiver voxels at time $t+\Delta t|k$, histories of the counts of all the output species from all receiver voxels up to time t$]$. We can use the method in [16] to determine this probability. $\Box$

Remark 2

It is possible for a receiver to have both mixed and partitioned configurations. The derivation of the filters for this case is straightforward. We use the summand of the RHS of Eq. (1op) for those voxels that are partitioned and that of Eq. (1os) for those voxels that are mixed. After that we sum up the log-posteriori probability contributed by each voxel.

We mentioned earlier that there is no loss of generality to consider all $P$ voxels being connected. If the receiver does not form one connected component, then we can compute the log-posteriori probability for each connected component and sum them up. $\Box$

Remark 3

We would like to discuss how the approximate demodulation filters may be implemented. For the partitioned case, we first point out that (1op) has the same form as (1n). We have recently derived a molecular circuit implementation of (1n) for the case where the transmission symbols are concentration-modulated signals in [42]. By leveraging these results, we have recently presented a molecular circuit implementation for the partitioned case in [52]. A work-in-progress is to implement the approximate demodulation filters for the mixed configuration.

We have so far assumed that the channel is stationary. An interesting open research problem is to consider a time-varying communication channel. This would require us to adapt the molecular circuit implementation of the demodulation filters to the channel. This appears to be a challenging problem but we want to point out that molecular circuit implementation of some relevant tools do exist, e.g. least-squares estimation [53] and feedback control for tracking [54].

$\Box$**

V Analytical approximation of BER

This section presents an analytical method to compute the BER of the sub-optimal demodulators for the partitioned and mixed configurations. In order to analytically compute the BER, we need to derive the joint probability distribution of $Z_{0}(t)$, $Z_{1}(t)$, … , $Z_{K-1}(t)$, where $Z_{k}(t)$ is the approximate log-posteriori probability given in (1or) or (1ot). This in turn requires us to compute the probability distribution of the molecular counts of the species in the system, whose exact solution is given by the solution of the RDME. However, for systems with nonlinear reaction rates, which is the case for our system, it is not possible to solve the RDME analytically because of the moment closure problem **[55]** which refers to the issue that the computation of the lower order moments requires the values of the higher order moments. An alternative method is the Linear Noise Approximation (LNA) **[30, 56]** which uses Gaussian distribution to approximate the distribution of the molecular counts when the system is in the steady state. The LNA is known to be a good approximation when the number of the molecules in the system is large. We will use LNA in this section.

In this section, we use a system which consists of three voxels, to illustrate the method. The voxels are arranged in a line and are numbered as 1, 2 and 3. The transmitter voxel is Voxel 3 and the two receiver voxels are in Voxels 1 and 2. Receiver voxel 1 (resp. 2) is in Voxel 1 (2). Note that LNA is based on concentration rather than molecular count. Earlier, we have used the notation $X_{p}(t)$, $X^{*}_{p}(t)$, $N_{R,p}(t)$ and $N_{i}(t)$ to denote molecular counts. Here we will use the corresponding lowercase notation — $x_{p}(t)$, $x^{*}_{p}(t)$, $n_{R,p}(t)$ and $n_{i}(t)$ — to denote the concentration. Let $V=w^{3}$ be the volume of a voxel, then we have $x_{p}(t)=\frac{X_{p}(t)}{V}$ etc. Also, when concentration is used, the reaction rates are expressed in terms of reaction rate constant rather than propensity function constants. For the activation reaction (1a), we denote the reaction rate constant by $\hat{k}_{+}$ and it can be shown that the reaction rate constant $\hat{k}_{+}$ and propensity function constant $k_{+}$ are related by $\hat{k}_{+}=Vk_{+}$. For the inactivation reaction (1) and diffusion between voxels, the rate constants are the same as the propensity function constants so we will not define new notation for them. We begin by explaining how the mean quantities can be computed. We will focus on the mixed configuration because it is the more complicated case and it provides us with insights on the property of the approximate demodulator.

V-A Mean molecular counts and reference signal

We first explain how the mean concentration of all species in the mixed configuration can be computed. For notation, we will add an overhead bar, e.g. $\bar{x}^{*}_{1}(t),\bar{X}^{*}_{1}(t)$ etc, to denote the mean quantity.

The derivation is based on the reaction rate law **[57]**. We will use the concentration of \ceeX^ in receiver voxel 1 as an example. The concentration of this species is affected by two reactions and two diffusion events: (i) activation reaction at a rate of $\hat{k}_{+}\bar{n}_{R,1}(t)\bar{x}_{1}(t)$; (ii) deactivation reaction at a rate of $k_{-}\bar{x}^{*}_{1}(t)$; (iii) diffusion of \ceeX^* from receiver voxel 1 to receiver voxel 2 at a rate of $d_{r}\bar{x}^{*}_{1}(t)$; (iv) diffusion of \ceeX^* from receiver voxel 2 to receiver voxel 1 at a rate of $d_{r}\bar{x}^{*}_{2}(t)$. Therefore the rate of change of $\bar{x}^{*}_{1}(t)$ is:*

[TABLE]

We can use the same method to write down the ODEs that govern the evolution of all species in the system. The full system of ODEs is in Appendix C. The solution of this system of ODEs will give us the mean concentration. The mean species count can then be computed by multiplying the mean concentration by the voxel volume. For the mixed configuration, we require $\beta_{k,p}(t)={\mathbf{E}}[X_{p}(t)N_{R,p}(t)|k]$ and we propose to approximate this by

[TABLE]

V-B Second order statistics

This section explains how the second order statistics of the molecular counts can be computed. We begin with using a simple example to explain the intuition behind the computation. Consider a chemical reaction that takes place at a constant mean rate of $r$ reactions per unit volume per unit time, then the number of reactions $\rho$ that occurs in a time interval $\tau$ in a volume $V$ is Poisson distributed with mean $Vr\tau$. By using Gaussian approximation, we can write $\rho$ as $Vr\tau+\sqrt{Vr\tau}\xi$ where $\xi$ is a standard Gaussian random variable and we have used the fact that the variance of a Poisson random variable is equal to its mean. One key idea of the LNA is to approximate the reaction rate using this method. However, additional approximations are required because the reaction rate are nonlinear and one also needs to approximate the square root term. Let us use $x^{*}_{1}(t)$ as an example. For LNA, the second order statistics of $x^{*}_{1}(t)$ is computed from $\sqrt{V}\widetilde{x}^{*}_{1}(t)$ where $\widetilde{x}^{*}_{1}(t)$ is a Gaussian random variable whose evolution is governed by:

[TABLE]

where $\gamma_{\rm ar,rv1}(t)$, $\gamma_{dr,RV1}$, $\gamma_{d,X*,1to2}(t)$, $\gamma_{d,X*,1to2}(t)$ and $\gamma_{d,X*,2to1}(t)$ are continuous-time standard Gaussian white noise associated with the reactions. Note that linearisation has been applied to the reaction rate $\hat{k}_{+}n_{R,1}(t)x^{*}_{1}(t)$ in the first line of (1ow). One can see that, if the mean concentrations are given, then (1ow) is a linear stochastic differential equation (SDE) in $\widetilde{x}^{*}_{1}(t)$, $\widetilde{n}_{R,1}(t)$ etc. We can use similar method to write down the SDEs for other species in the system. With this set of SDEs, we can compute the covariance matrix for vectored Gaussian distribution of the species concentrations by solving a Lyapunov differential equation **[30, 56]**.

V-C Mean and covariance of $Z_{k}(t)$

Having obtained the mean and covariance of the Gaussian distribution associated with concentrations of the species in the system, we can now compute the joint probability distribution of $Z_{0}(t)$, $Z_{1}(t)$, … , $Z_{K-1}(t)$. We can see from (1or) and (1ot) that $Z_{k}(t)$ is linear in the stochastic quantities, i.e. $\left[\frac{dX^{*}_{p}(t)}{dt}\right]_{+}$ and $X^{*}_{P}(t)$, therefore the $Z_{k}(t)$’s are jointly Gaussian distributed. The key problem here is to estimate the statistical property of $\left[\frac{dX^{*}_{p}(t)}{dt}\right]_{+}$. Since the increment of $X^{*}_{p}(t)$ can be due to an activation reaction or the diffusion of an \ceeX_ into the voxel, we model this term by using the SDE:*

[TABLE]

By combining this with the Lyapunov differential equation for calculating the covariance of the concentration, we can compute the joint probability distribution of $Z_{0}(t)$, $Z_{1}(t)$, … , $Z_{K-1}(t)$.

When the number of symbols $K$ is 2, the BER can be computed as follows. If Symbol 0 is transmitted, bit error occurs if $Z_{0}(t)<Z_{1}(t)$. Since $Z_{0}(t)$ and $Z_{1}(t)$ are jointly Gaussian distributed, the random variable $Y(t)=Z_{0}(t)-Z_{1}(t)$ is Gaussian distributed, therefore the BER equals to the probability of $Y(t)<0$. The computation of the BER for Symbol 1 is similar.

VI Numerical Examples

This section presents numerical results to understand the performance of the partitioned and mixed configurations. We first describe the methodology and then the results.

VI-A Methodology

We assume the propagation medium has the shape of a rectangular prism divided into voxels with edge length $w$. Different numerical studies will use different medium size and we will also study the impact of the parameter $w$, so we will specify their numerical values when the present the studies.

The diffusion coefficient $D$ of the propagation medium is 1 $\mu$m2s*-1**. Unless otherwise stated, we assume an absorbing boundary condition where signalling molecules may leave the surface of boundary voxel at a rate $\frac{d}{50}$.*

The transmitter is assumed to use $K=2$ symbols. Unless otherwise stated, each symbol is represented by an emission pattern which is generated by a chemical reaction of the form:

[TABLE]

where $r_{k}$ is the rate of production of signalling molecules \ceeS when symbol $k$ is transmitted. We assume that symbols 0 and 1 cause, respectively, 10 and 40 signalling molecules to be generated per second on average by the transmitter. The reaction rate constants for the reactions (1oa) is 0.005 $\mu$m3s*-1** and for reaction (1o) is 1s**-1**. Note that the diffusion and reaction parameters are the same as those we used in our previous work [15], and are realistic for a biological environment.*

The number of receptors will vary from experiment to experiment and will be stated for each experiment.

Unless otherwise stated, the molecular circuits at the receiver voxel is the circuit given in (1).

In the following experiments, for the partitioned configuration, the total number of \ceeX and \ceeX^ in each receiver voxel is given by the parameter $M$. For the mixed configuration, the value of $M$ should be interpreted as the total number of \ceeX and \ceeX^* in a receiver voxel at the beginning of the simulation. Note that this parameter is common to all receiver voxels.*

The parameter $D_{r}$ is the diffusion coefficient for \ceeX and \ceeX^. In the experiments, we will use the inter-voxel diffusion rate $d_{r}=\frac{D_{r}}{w^{2}}$ (in s) instead. Note that $d_{r}=0$ is the same as the partitioned configuration.*

Since our interest is in the demodulation performance, we assume that there is no inter-symbol interference (ISI). Note that if ISI is present, we can deal with it using the decision feedback algorithm in our earlier work in **[15]**.

We use Stochastic Simulation Algorithm (SSA) **[58]** to simulate the CTMP that models both diffusion and reaction of molecules in the system.

The approximate demodulation filter (1oq) requires the mean $\alpha_{k,p}(t)=E[N_{R,p}(t)|k]$ while the filter (1ot) requires the mean $\beta_{k,p}=E[x_{1}(t)N_{R,p}(t)|k]$. We will use SSA simulation to estimate these means by running SSA simulation 500 times and compute the average.

We numerically integrate the approximate demodulation filters to obtain $Z_{0}(t)$ and $Z_{1}(t)$. We use the initial condition $Z_{k}(0)=0$ for all $k$ which means that all symbols are equally probable in the system. We will use BER as the performance metric. Each BER value is estimated using 300 independent SSA runs.

The BER at time $t$ is determined as follows: If the transmission symbol is 0 (resp. 1), then a bit-error occurs at time $t$ if $Z_{1}(t)>Z_{0}(t)$ (resp. $Z_{0}(t)>Z_{1}(t)$).

VI-B Partitioned case: Comparing the approximate filter (1or) against the optimal filter (1op)

The aim of this section is to verify that the approximate filter (1or) for computing the log-posteriori probability for the partitioned case is a good approximation of the optimal filter (1op).

The optimal demodulation filter (1op) requires the computation of ${\mathbf{E}}[N_{R,p}(t)|k,{\cal X}_{R}^{*}(t)]$ which can be obtained by solving an optimal Bayesian filtering problem. Since filtering problems are computationally expensive to solve, we propose the approximate demodulation filter (1or) which approximates ${\mathbf{E}}[N_{R,p}(t)|k,{\cal X}_{R}^{*}(t)]$ by ${\mathbf{E}}[N_{R,p}(t)|k]$. The aim of this subsection is to compare the accuracy of the optimal and approximate filters.

In this comparison, we consider a transmission medium of 1$\mu$*m $\times$ *\frac{1}{3}$$\mu*m $\times$ *\frac{1}{3}$$\mu*m. With $w=\frac{1}{3}$ *$\mu$m, the medium consists of an array of $3\times 1\times 1$ voxels. We will index these three voxels sequentially by indices 1, 2 and 3. We assume Voxel 1 is the transmitter; and Voxels 2 and 3 are the receiver voxels. A reflecting boundary condition is assumed.

The reason why we have chosen to use such a small number of voxels is because of the dimensionality of the Bayesian filtering problem. The filtering problem requires us to compute the probability of the state vector $N(t)=(N_{1}(t),N_{2}(t),N_{3}(t))$ where $N_{i}(t)$ is the number of signalling molecules in the voxel with index $i$. If each voxel can have a maximum of 100 signalling molecules at a time, then there are approximately 106* possible $N(t)$ vectors and the filtering problem has to estimate the probability ${\rm P}[N(t)|k,{\cal X}_{R}^{*}(t)]$ for each possible $N(t)$ vector. Although there are approximation techniques to solve the Bayesian filtering problem, that would introduce inaccuracies. The use of small number of voxels will allow us to compute ${\rm P}[N(t)|k,{\cal X}_{R}^{*}(t)]$ precisely.*

For this experiment, the $K=2$ symbols are deterministic. Symbol $k$ emits $s_{k}$ signalling molecules at times 0, 0.2 and 0.4, where $s_{0}=8$ and $s_{1}=20$. The number of $M$ per voxel is 10. The simulation time is about 2 seconds.

Figs. 7(a) to 7(d) compare ${\mathbf{E}}[N_{R,p}(t)|k]$ and ${\mathbf{E}}[N_{R,p}(t)|k,{\cal X}_{R}^{*}(t)]$ for $p=1,2$ and $k=0,1$ when Symbol 1 was transmitted. It can be seen that the approximation is quite accurate. The results for Symbol 0 are similar.

Eqs. (1or) and (1op) show that the log-posteriori probability of the receiver can be computed by summing up the log-posteriori probabilities computed from each voxel. Figs. 7(e) to 7(h) compare $L_{0}(t)$ and $L_{1}(t)$ computed from the two voxels when Symbol 1 was transmitted. It can again be seen that the approximation is quite accurate. The results for Symbol 0 are similar.

In the rest of this section, we will use the approximate demodulation filter (1or) because of its lower computational complexity.

VI-C Mixed case: Comparing the approximate filter (1ot) against the optimal filter (1os)

The aim of this subsection is to demonstrate that the approximate filter (1ot) for the mixed configuration is an accurate approximation of the optimal version (1os).

We use the same voxel configuration as subsection VI-B. We allow \ceeX and \ceeX^ to diffuse between the receiver voxels with a $D_{r}=0.2D$. For the filtering problem, we need to compute the probability of the state $N(t)=(N_{1}(t),N_{2}(t),N_{3}(t),X_{R,1}(t),X_{R,2}(t))$ where $N_{i}(t)$ is the number of signalling molecules in the voxel with index $i$ and $X_{R,j}(t)$ is the number of \ceeX molecules in receiver voxel $j$.*

Initially, each receiver voxel has $4$ \ceeX molecules. The transmission symbols are similar to those in subsection VI-B except that $s_{0}=10$ and $s_{1}=15$.

Figs. 8(a) to 8(d) compare ${\mathbf{E}}[X_{p}(t)N_{R,p}(t)|k]$ and ${\mathbf{E}}[X_{p}(t)N_{R,p}(t)|k,{\cal X}_{R}^{*}(t)]$ for $p=1,2$ and $k=0,1$ when Symbol 1 was transmitted. Note that we computed ${\mathbf{E}}[X_{p}(t)N_{R,p}(t)|k]$ by averaging the results from 500 SSA simulations. It can be seen that the approximation is reasonably accurate. The results for Symbol 0 are similar.

Eqs. (1ot) and (1os) show that the log-posteriori probability of the receiver can be computed by summing up the log-posteriori probabilities computed from each voxel. Figs. 8(e) to 8(h) compare $L_{0}(t)$ and $L_{1}(t)$ computed from the two voxels when Symbol 1 was transmitted. It can again be seen that the approximation is quite accurate. The results for Symbol 0 are similar.

VI-D BER for mixed and partitioned configurations

The aim of this experiment is to compare the BER for mixed and partitioned configurations for different receiver locations. We assume that the size of propagation medium is 1$\frac{2}{3}\mu$m$\times$1$\frac{2}{3}\mu$m$\times$1$\frac{2}{3}\mu$*m and the size of voxel is *\frac{1}{3}$$\mum3 (i.e. $w=\frac{1}{3}$ $\mu$ m). This forms a grid of $5\times 5\times 5$ voxels.

We use a 3-tuple $(x,y,z)$ where $1\leq x,y,z\leq 5$ to identify the locations of the voxels. The transmitter voxel is located at (1,1,1). We use two receiver voxels and place them at (4,5,5) and (5,5,5). The transmitter and receiver voxel locations are depicted in the left-most picture in Fig. 6.We use $M=10$ and $P=2$.

We study the impact of $d_{r}$ on the BER. Three different values for $d_{r}$ are used: 0, 0.5 and 1. The SSA simulations is performed up to time 2.5. Fig. 9 shows the BER for Symbols 0 and 1. It can be seen that a lower $d_{r}$ leads to a lower BER. In particular, the partitioned configuration leads to the lowest BER.

The above simulations assume that the receiver voxels are placed at a diagonally opposite corner to the transmitter voxel. For two receiver voxel configuration, We have also performed simulations with the receiver voxels at different locations within the medium. The different receiver voxel placements that we have used are: (3,4,4), (4,4,4); (2,3,3), (3,3,3); and (4,1,1),(5,1,1). Fig. 14 shows BER is lower for smaller $d_{r}$.

VI-E BER for different values of diffusion coefficient

The aim of this experiment is to further study the impact of diffusion of the receiver species on BER. We choose $M=10$ and $P=2$. We vary $d_{r}$ $\mu$m2s*-1** from 0 to 1 with an increment of 0.1. We used the BER at time 2.5 sec for comparison. Figs. 10(a) and 10(b) show how BER varies with $d_{r}$ for, respectively, Symbols 0 and 1. It shows that BER increases monotonically with $d_{r}$.*

We can consider the partitioned configuration as perfect isolation of receptors into clusters where there is a cluster per voxel and the mixed configuration as imperfect isolation where larger values of $d_{r}$ means farther away from perfect isolation. The results in this section show that our approximate demodulators offer a gradual degradation in performance with $d_{r}$.

VI-F Impact of the number of receiver voxels with fixed $M$

This section studies the impact of the number of receiver voxels on BER for partitioned configuration. We maintain $M=10$ and we use three different receivers, with 2, 4 and 6 voxels. The total number of receptors for these receivers are therefore 20, 40 and 60 respectively. The voxel locations are: (4,5,5), (5,5,5) for 2 receiver voxels; (2,5,5), (3,5,5), (4,5,5) and (5,5,5) for 4 receiver voxels; (5,4,5), (1,5,5), (2,5,5), (3,5,5), (4,5,5) and (5,5,5) for 6 receiver voxels. The transmitter and receiver voxel locations are depicted in Fig. 6a-Fig. 6c. Fig. 11 shows that a higher number of voxels leads to a lower BER.

We have also simulated two other configurations for the 6 receiver voxel case by changing the location of receiver voxels in the medium. The different receiver voxel positions that we have used are:

• •

(3,4,5), (4,4,5), (5,4,5), (3,5,5), (4,5,5), (5,5,5) which is depicted in Fig. 6d; and,

• •

(3,4,1), (4,4,1), (5,4,1), (3,5,1), (4,5,1), (5,5,1)

Fig. 15 shows that the results follow the same trend.

VI-G Impact of the number of receiver voxels with a fixed total number of receptors

This section studies the impact of the number of receiver voxels on BER. We use three different number of voxels per receiver, namely 2, 4 and 6 voxels. We maintain the total number of receptors in each receiver at 60. Therefore, the number of receptors per voxel for the three receivers are 30, 15 and 10. Figs. 12(a) and 12(b) show that BER for Symbols 0 and 1 over time. It shows that in general, a higher number of voxels will lead to a lower BER.

VI-H BER using a different molecular receiver circuit

All the above experiments have been carried out with the receiver molecular circuit in (1). In this experiemnt, we use the following molecular circuit which was used in our earlier work **[16]**:

[TABLE]

where $E$ represent an unbound receptor with two binding sites. In forward reaction (1oy), $E$ can bind with $S$ molecule to form the complex $C_{[1]}$ whereas in forward reaction (LABEL:eq:mc1:r2), $C_{[1]}$ can bind with $S$ molecule to form the complex $C_{[2]}$. Furthermore, $\tilde{\lambda}_{1}$, $\mu_{1}$, $\tilde{\lambda}_{2}$ and $\mu_{2}$ are reaction rate constants. The complex $C_{[2]}$ is chosen as the output species.

We assume $M=2$ and $P=2$. All other parameters remain the same. We use three values of $d_{r}$: 0, 0.5 and 1. Figs. 13(a) and 13(b) show the BER for, respectively, Symbols 0 and 1. We witness the same trend as before where the BER increases with $d_{r}$.

VI-I Impact of the voxel size

We know from the literature on RDME that the voxel size has to be chosen correctly in order to produce correct simulation results **[59]**, e.g. the paper **[60]** mentions a lower bound on the voxel size. Here, we present an example to show that, if the voxel size is chosen correctly, then we are able to obtain consistent BER for different voxel sizes.

Unless otherwise stated, we use the same parameter values as Section VI-D. We consider two voxel edge lengths: $w=\frac{1}{3}$ and $w=\frac{1}{6}$ while maintaining the same size for the medium, transmitter and receiver. We make the following adjustments: (i) A receiver consisting of 2 voxels for $w=\frac{1}{3}$ will have 16 voxels for $w=\frac{1}{6}$; (ii) Let $u$ be the rate at which signalling molecules are produced in the transmitter voxel when $w=\frac{1}{3}$, then they are produced at a rate of $\frac{u}{8}$ in each of the 8 transmitter voxels for $w=\frac{1}{6}$. In addition, we also use two diffusion coefficients $D$: 1 $\mu$m2s*-1** and 2 $\mu$m2s**-1**. We assume $M$ = 40. Fig. 16 shows that the voxel dimension $w=\frac{1}{3}$ and $w=\frac{1}{6}$ give consistent BER. We remark that there is some recent work on automating the choice of voxel size for RDME, see [59].*

VI-J Analytical BER formula

In Section V we present a method to analytically compute the BER using LNA. The aim of this section is to study the accuracy of the method. We assume that the the medium is a cube of dimension ($\frac{2}{3}\mu$*m)3. The voxel edge length $w$ is *$\frac{1}{3}\mu$m giving a medium shape of 2-by-2-by-2 voxels. The reason why we have chosen to use a small medium size is that, for the verification of the analytical formula, we will need to simulate many times to obtain an accurate estimation of the BER. The transmitter is located at voxel (1,1,1) and the receiver voxels are at (1,2,2) and (2,2,2). We use the same signalling molecule diffusion coefficient and reaction parameters as before. We assume $K=2$. Since LNA works best when the system is in steady state, we create a steady state in the receiver by using short pulses (0.2s) as the transmission symbol and a much longer simulation time of 20s to allow the system to get to steady state. Also, a reflective boundary condition is used.

We first verify that LNA can be used to accurately compute the mean and the covariance matrix of $Z_{0}(t)$ and $Z_{1}(t)$ for the case where $d_{r}=0.2$. The verification is done by performing SSA simulations 5000 times per input symbol. Fig. 17 shows LNA can accurately estimate the mean and the covariance matrix of $Z_{0}(t)$ and $Z_{1}(t)$ for Symbol 0. Note that there are 3 distinct elements in the covariance matrix, which are the variance of $Z_{0}(t)$ and $Z_{1}(t)$ as well as the covariance of $Z_{0}(t)$ and $Z_{1}(t)$. The results for Symbol 1 is similar, see Fig. 18.

Next we use mean and covariance matrix of $Z_{0}(t)$ and $Z_{1}(t)$ to estimate the BER. Fig. 20 plots the average BER estimated by LNA for the cases $d_{r}=0,0.1,0.2$, against those obtained from SSA. It can be seen the prediction obtained from LNA is fairly accurate. Note also that the BER estimation is more accurate in the later part of the simulation when the system is in steady state. Intuitively, this is because the approximation in LNA is based on approximating the number of reactions in a time interval, see the discussion in the beginning of Section V-B. We remark that the BER estimation using SSA becomes inaccurate for $d_{r}=0$ after time 10 because of the small BER and limited number of SSA simulations.

Remark 4

A practical way of making use of the results in this section for molecular communications is as follows. The transmitter sends a symbol and the receiver runs until steady state. The receiver then makes a decision on the transmission symbol that is sent. Once the decision is made, then the receiver should reset itself to get ready for the next transmission symbol to be received. The resetting of the receiver can be done using chemical reactions, see our recent work [52].

VI-K What if it is possible to compute $A_{p}(t)$?

In the approximate demodulator (1ot), we used $\left[\frac{dX^{*}_{p}(t)}{dt}\right]_{+}$ instead of $A_{p}(t)$ because there does not appear to be a way to distinguish between activation events from diffusion events. In this numerical study, we ask what if it is possible to compute $A_{p}(t)$ and use it in (1ot) instead of $\left[\frac{dX^{*}_{p}(t)}{dt}\right]_{+}$.

We use the same set up in Section VI-J and use SSA simulations to estimate the average BER for $d_{r}=0.1$ and $d_{r}=0.2$ if $A_{p}(t)$ is used. We compare these BERs against that of the partitioned configuration ($d_{r}=0$) in Fig. 20. It shows that if $A_{p}(t)$ can be obtained, then the mixed configuration has the same BER as the partitioned configuration. This shows that partitioning is a useful method to reduce BER because it takes away hard task of identifying the activation events in the mixed configuration.

VII Conclusions and Future Work

This paper considers the demodulation in a molecular communication system where the receiver can be in the partitioned or mixed configuration. We derive the MAP demodulator for both configurations. Our numerical experiment shows that partitioning, or a small diffusion coefficient for the receiver species, can lead to a lower BER. The use of partitioning does not seem to have been studied before so this leads to a new degree of freedom to improve the performance of molecular communication.

This paper models the shape of a receiver by using multiple voxels. This can be considered as using the finite difference method to model the shape of a receiver. The finite difference method is a very basic method to model the shape of 3-dimensional objects. The modern approach is based on finite element or similar methods. We see this as an interesting future direction. The diffusion of molecules in the finite element setting can possibly be handled by the method in **[61]**. However, reaction-diffusion setting is a lot more complicated than the pure diffusion case. This is because RDME is fundamentally a spatially discrete method to approximate the behaviour of the fine grained Smoluchowski equation **[43]**. This approximation is only accurate if the discretisation length scale is chosen correctly **[60]**. Fortunately, there is some recent work on using RDME on an unstructured mesh **[60, 62]**. Since our demodulator design is based on RDME, these recent work will allow us to extend our work to unstructured mesh to better model the 3-dimensional shape of the receiver.

Appendix A Partitioned Configuration

In this appendix we show how the demodulation filter Eq.(1op) can be derived. We first explain how the case for $P=2$ can be derived and then explain how it can be generalised.

We follow the method in **[16]**. The first step of the derivation is determine the probability $\mathbf{P}[X^{*}_{1}(t+\Delta t),X^{*}_{2}(t+\Delta t)|k,{\cal X}^{*}_{1}(t),{\cal X}^{*}_{2}(t)]$. In **[16]** we present an algorithm to write down the expression of this probability by identifying all the reactions that can change the count of the output species, i.e. \ceeX^_1 and \ceeX^_2. By using the algorithm in **[16]**, we can show that:

[TABLE]

where $\delta()$ is an indicator function which takes the value of 1 if all the conditions within $()$ are true, otherwise its value is 0. In addition, we have:

[TABLE]

Note that the term $Q_{1,a}$ in Eq. (1oz) corresponds to the case where the activation reaction (1a) takes place in receiver voxel 1 because $\delta(X^{*}_{1}(t+\Delta t)$ is one greater than $X^{*}_{1}(t)$. The subscripts 1 and $a$ in $Q_{1,a}$ refer to receiver voxel 1 and activation reaction. Similarly, $Q_{1,d}$ refers to deactivation reaction in voxel 1. The terms $Q_{2,a}$ and $Q_{2,d}$ are for voxel 2. Lastly, the term $Q_{0}$ corresponds to no reactions taking place.

The next step is to derive the ODE which shows how the log-posteriori probability $L_{k}(t)$ evolves over time. From **[16]**, we have:

[TABLE]

where $L^{\prime}(t)$ is a term independent of symbol $k$. Since $L_{k}(t)$ does not appear on the RHS of the above equation and $L^{\prime}(t)$ adds the same contribution to all $L_{k}(t)$ for all $k=0,...,K-1$, we can therefore ignore $L^{\prime}(t)$ for the purpose of demodulation since it is the relative (rather than the absolute) magnitude of $L_{k}(t)$ which is needed for demodulation.

By dropping $L^{\prime}(t)$, we can compute the shifted version of the log-posteriori probability $L_{k}(t)$. For conciseness, we use $L_{k}(t)$ to denote the shifted version of the log-posteriori probability. We therefore have:

[TABLE]

The next step is to substitute Eq. (1oz) into Eq. (1oac). After some lengthy manipulations, we arrive at:

[TABLE]

which is Eq. (1op) for the case of $P=2$.

For general $P$, there will be $(2P+1)$ terms in the counterpart of Eq. (1oz). Out of these $(2P+1)$ terms, $2P$ of them are $Q_{p,a}$ and $Q_{p,d}$ for $p=1,\ldots,P$. The last term is $Q_{0}$, which equals to $1-\sum_{p=1}^{P}(Q_{p,a}+Q_{p,d})$. After writing down the counterpart of Eq. (1oz) for $P$ voxels, we can follow the above procedure to derive Eq. (1op).

Appendix B Mixed Configuration

In this appendix we derive the demodulation filter for the mixed configuration for $P=2$. The derivation in this appendix is similar to that for the partitioned configuration in A.

The first step of the derivation is to determine the probability $\mathbf{P}[X^{*}_{1}(t+\Delta t),X^{*}_{2}(t+\Delta t)|k,{\cal X}^{*}_{1}(t),{\cal X}^{*}_{2}(t)]$. As mentioned in Appendix A, we can use the method in **[16]** which is to identify all the reactions that can change the counts of the output species, i.e. \ceeX^_1 and \ceeX^_2. It is important to point out here that the term reactions here takes on a generalized meaning. In the modelling framework of RDME, the diffusion of a species from one voxel to another voxel is considered as a first order chemical reactions **[44]**. Therefore, when we consider the reactions that can change the counts of \ceeX^_1 and \ceeX^_2, we will also need to include the diffusion of the \ceeX^ species between the voxels. Recall that $D_{r}$ is the diffusion coefficient of \ceeX^. Let $d_{r}=\frac{D}{w^{2}}$ where $w$ is the lenght of a voxel edge. By using the algorithm in **[16]**, we have:

[TABLE]

where

[TABLE]

The meanings of the terms $Q_{1,a}$, $Q_{1,d}$ etc. are the same as those in Appendix A. The term $Q_{1\rightarrow 2}$ corresponds to the diffusion of an \ceeX^ molecule from receive voxel 1 to receiver voxel 2. Starting with Eq. (1oae), we can now follow the same procedure mentioned in Appendix A to obtain the demodulation filter (1os) for $P=2$.*

Similarly, we can generalise Eq. (1oae) to general $P$. We need to include terms $Q_{p,a}$ and $Q_{p,b}$, as well as diffusion terms $Q_{p_{1}\rightarrow p_{2}}$ and $Q_{p_{2}\rightarrow p_{1}}$ between any pairs of receiver voxels $p_{1}$ and $p_{2}$ which are neighbours. After forming the counterpart of Eq.(1oae) for $P$ receiver voxels, we can obtain the demodulation filter (1os).

Appendix C Example for Section V

The evolution of the mean concentration of the species in the system is governed by:

[TABLE]

where $u_{k}(t)$ is the concentration of the transmission symbol $k$ at time $t$.