arXiv:1204.4253v5 [] 2 Nov 2013

Extended master equation models for molecular
communication networks
Chun Tung Chou
School of Computer Science and Engineering,
University of New South Wales,
Sydney, Australia. E-mail: ctchou@cse.unsw.edu.au
January 27, 2018
Abstract
We consider molecular communication networks consisting of transmitters and receivers
distributed in a fluidic medium. In such networks, a transmitter sends one or more signalling
molecules, which are diffused over the medium, to the receiver to realise the communication.
In order to be able to engineer synthetic molecular communication networks, mathematical
models for these networks are required. This paper proposes a new stochastic model for
molecular communication networks called reaction-diffusion master equation with exogenous
input (RDMEX). The key idea behind RDMEX is to model the transmitters as time series
of signalling molecule counts, while diffusion in the medium and chemical reactions at the
receivers are modelled as Markov processes using master equation. An advantage of RDMEX
is that it can readily be used to model molecular communication networks with multiple
transmitters and receivers. For the case where the reaction kinetics at the receivers is linear,
we show how RDMEX can be used to determine the mean and covariance of the receiver
output signals, and derive closed-form expressions for the mean receiver output signal of the
RDMEX model. These closed-form expressions reveal that the output signal of a receiver
can be affected by the presence of other receivers. Numerical examples are provided to
demonstrate the properties of the model.

Keywords: Molecular communication networks, nano communication networks, synthetic molecular communication networks, master equations, stochastic models, synthetic biology

1

Introduction

We consider molecular communication networks consisting of transmitters and receivers distributed in a fluidic medium. In such networks, a transmitter sends one or more signalling
molecules, which are diffused over the medium, to the receiver to realise the communication.
The study of molecular communication has its origin in biology and biophysics. Molecular communication is a vital mechanism in multi-cellular organisms. The human body, which has an
estimated 1014 cells, uses molecular communication to keep the body in a healthy state. In fact,
cells in the human body constantly communicate with other cells using molecular communication.
There are a couple of reasons why synthetic molecular communication networks, which are
inspired by molecular communication in living organisms, should be studied [1, 2, 3]. Firstly,
synthetic molecular communication networks can be combined with nano-sensors and molecular
computing [4] to form nano-sensor networks [5] for health monitoring, medical diagnosis and
cancer therapy. Secondly, the study of synthetic molecular communication networks can be used
to enhance our understanding of their biological counterparts.
In order to be able to engineer synthetic molecular communication networks, we need mathematical models which can be used to predict the performance of these networks. For example,
if a transmitter in a molecular communication network emits a number of signalling molecules
1

to communicate with a receiver, we would like to be able to determine the receiver output signal in order to determine the probability of correct reception at the receiver. Such evaluations
can be realised if a mathematical model is available to determine receiver output signal based
on the transmitter’s input signal. The main contribution of this paper is that we propose a new
stochastic model for molecular communication networks. Our model is based on reaction-diffusion
master equation (RDME) [6] which is a well known model in physics and chemistry for modelling
systems with both diffusion and chemical reactions. In this paper, we propose an extension to
RDME, which we call reaction-diffusion master equation with exogenous input (RDMEX). The
key idea behind RDMEX is to model the transmitters as time series of signalling molecule counts,
while diffusion in the medium and chemical reactions at the receivers are modelled by RDME.
An advantage of RDMEX is that it can readily be used to model molecular communication networks with many transmitters and many receivers. For the case where the reaction kinetics at
the receivers is linear, we show how RDMEX can be used to determine the mean and covariance
of receiver output signal of molecular communication networks. These results allow us to derive
closed-form expressions showing how the receiver outputs relate to the transmitter signals when
there are multiple transmitters and receivers. These expressions show that the output of a receiver
can be influenced by the presence of other receivers in a molecular communication network. They
also reveal the coupling between diffusion and chemical reactions at the receivers.
This paper is organised as follows. In Section 2, we present background materials on master
equations. In Section 3, we present the RDMEX and show how it can be used to model molecular
communication networks with multiple transmitters and receivers. We also show in this section,
for the case where the reaction kinetics at receivers is linear, how we can determine the mean
and covariance of the receiver output signals in molecular communication networks. The rest
of the paper is focused on determining the mean receiver output signal. We approach this by
using two methods, which will be discussed in sections 4 and 5. In Section 4, we determine the
continuum limit (i.e. infinite spatial resolution) of the RDMEX and show that it results in a
reaction-diffusion partial differential equation (RDPDE). We derive a closed-form solution to this
RDPDE and interpret the results. In Section 5, we determine the mean receiver output signal of
RDMEX with finite spatial resolution and derive another closed-form solution. Numerical results
are then presented in section 6 to show the accuracy of our solutions. Finally, section 7 describes
the related work and section 8 gives the conclusions.

2

Background on master equations

The aim of this section is to provide the necessary background on master equations. The treatment
here is brief and includes only the results needed for this article. The reader can refer to the texts
[6, 7] or tutorial article [8] for a more complete treatment of this subject. This section is divided
into two parts. We first introduce the general master equation and give a simple example on how
to use master equation to model a chemical reaction. We then quote some results on mean and
covariance of the Markov processes modelled by master equations.

2.1

General master equation

Consider a continuous-time integer-value vector Markov process Q(t) ∈ Zp , where p is the number
of vector components, Z is the set of all integers and t is time. When the Markov process Q(t) is
in state q ∈ Zp , it jumps to the state q + rj (where rj ∈ Zp with j = 1, 2, ...J and J is the total
number of possible jumps) at a transition rate of Wj (q). Let P (q, t|q0 , t0 ) denote the conditional
probability that Q(t) = q given that Q(t0 ) = q0 .
We are interested to determine how P (q, t|q0 , t0 ) evolves over time. We can do this by using a
coupled set of ordinary differential equations (ODEs) known as the master equation:
dP (q, t|q0 , t0 )
dt

=

J
X
j=1

Wj (q − rj )P (q − rj , t|q0 , t0 ) −
2

J
X
j=1

Wj (q)P (q, t|q0 , t0 )

(1)

where one equation of the form (1) is needed for each valid state q. Note that the first and second
terms on the right-hand side of (1) can be interpreted, respectively, as the rates of entering and
leaving the state q. In order to simplify notation, we will write P (q, t) instead of P (q, t|q0 , t0 ) from
now on.
A common application of master equation is to model the dynamics of chemical reactions [9].
We will give a simple example to illustrate that.
Example 1 Consider the chemical reaction:
k+

−−
*
L+R)
−
−C
k−

where the chemical species are L, R and C, and the forward and reverse reaction constants are k+
and k− respectively. This chemical reaction can be described by a Markov process with state space
q = [nL , nR , nC ]T where nL is the number of molecules of chemical L etc., and T denotes matrix
transpose.
There are two possible types of jumps (i.e. J = 2) in this system. The forward reaction is
modelled by the jump r1 = [−1, −1, 1]T where the entries of r1 reflects the fact that one molecule
of L and one molecule of R react to form a molecule of C. The rate of jump, which according
to standard result in chemical kinetics, is W1 (q) = k+ q(1)q(2) = k+ nL nR where q(1) is the
first component of the vector q, etc. Similarly, the reverse reaction is modelled by the jump
r2 = [1, 1, −1]T with W2 (q) = k− q(3) = k− nC . The master equation for this chemical reaction is:
dP (q, t)
dt

=

2
X
j=1

2
X

Wj (q − rj )P (q − rj , t) −

Wj (q)P (q, t)

(2)

j=1

Equations of the type (2) is also known as chemical master equations.

2.2

Results on mean and covariance

The master equation (1) shows the time evolution of the state probability of the Markov process
Q(t). However, (1) can be difficult to work with because we need one equation for each valid state,
hence the number of equations can be very large. Therefore, it is easier if we can determine the
mean and covariance of the Markov process Q(t), which is defined as follows:
X
X
hQ(t)i =
qProb(Q(t) = q) =
qP (q, t)
(3)
q

Σ(t)

=

q

X
(q − hQ(t)i)(q − hQ(t)i)T P (q, t)

(4)

q

where h•i will be used in this paper to denote the mean operator.
It is possible to use (1) to derive how the mean and covariance of the state of the Markov
process Q(t) evolve over time. The following proposition is taken from [10].
Proposition 1 For the general master equation (1), we have
dhQ(t)i
dt

=

J
X
j=1

rj hWj (Q(t))i

In particular, if Wj (q) is a linear function of q. Let
dhQ(t)i
dt

=

AhQ(t)i

dΣ(t)
dt

=

AΣ(t) + Σ(t)AT +

PJ

j=1 rj Wj (q)

(5)
= Aq, then
(6)

J
X
j=1

3

rj rjT Wj (hQ(t)i)

(7)

3

Modelling molecular communication network using RDMEX

We consider a molecular communication network with multiple transmitters and multiple receivers
in an isotropic fluidic medium, see Figure 1. The communication takes place by the transmitter
emitting one or more signalling molecules over time. Once the molecules leave the transmitter,
they diffuse in the medium according to Brownian motion. The receivers are assumed to consist of
one or more receptors. When a signalling molecule L reaches a receptor R, they may bind together
to form a complex C. We will consider the number of complexes at the receiver as the output of
the receiver. For example, a receiver may infer that a bit has been sent by the transmitter when
the number of complexes exceeds a threshold. Note that use of this definition of output is also
used in chemotaxis models in biophysics [11, 12] and molecular communication network models in
engineering [13].
Based on the above description, we see that a model for molecular communication networks
must at least capture the diffusion of signalling molecules and the reactions at the receivers. Both
reactions and diffusion can be modelled as Markov processes, so the master equation (1) is a natural
choice. The remaining issue is how we can model the transmitter. In this paper, we model the
transmitter by a time sequence which specifies the number of molecules emitted by the transmitter
at a particular time. This approach is fairly general and can be used to model encoding methods
that have been considered in the literature, such as molecular coding [1] and concentration coding
[14]. We will consider other modelling approaches, e.g. modelling the internal mechanism of the
transmitters, in future work.
In section 3.1, we introduce the RDMEX model by way of an example and then we prove some
results on mean and covariance of the RDMEX model in section 3.2.

3.1

The RDMEX model

In this section, we will introduce the RDMEX model for a molecular communication network
with 2 transmitters, 2 receivers in a 1-dimensional medium. The reason for that is to simplify the
presentation. It is fairly simple to generalise the model to multiple transmitters, multiple receivers
in a 3-dimensional medium, which we will discuss at the end of the section.
Another simplification is that we will assume, at the receiver, the rate at which the complexes
are formed is a linear function of the number of local signalling molecules and is independent
of the number of receptors. (This will be made precise below.) This simplification allows us to
produce closed form expressions in the continuum limit. Note that, it is straightforward to model
non-linear reaction rates or incorporate more complex receivers in our model.
The following is a list of model assumptions, parameters and notation.
1. We assume that both transmitters use one and the same type of signalling molecule L.
2. For transmitter 1, we assume that it emits k1,1 signalling molecules of L at time t1,1 , k1,2
molecules of L at time t1,2 , ..., k1,b molecules of L at time t1,b , where k1,b (b = 1, 2, 3, ...) are
positive integers and t1,b ∈ R. Similarly, transmitter 2 emits k2,b signalling molecules of L at
time t2,b where b = 1, 2, 3, .... The number of molecules ka,b emitted at time ta,b is assumed
to be independent of the state of the system at or before ta,b .
3. The medium is assumed to be a 1-dimensional space of length X. The medium is partitioned
into N equal width voxels of width ∆x such that N ∆x = X. We index the voxels by using
1, ..., N . See Figure 2 for an illustration.
4. The medium is assumed to be isotropic. The rate at which a signalling molecule L diffuses
from one voxel to a neighbouring voxel is d̆ per molecule per unit time. The rate of diffusion
from a voxel to a non-neighbouring voxel is zero. Also, the molecule cannot leave the medium
and we assume the boundary is reflective. The parameter d̆ is related to one-dimensional
D̆
macroscopic diffusion constant D̆ by d̆ = ∆x
2.
4

5. We assume that each transmitter or receiver occupies only one voxel. Transmitters 1 and
2 are located respectively in the voxel indexed by T1 and T2 . Similarly, we assume that
receivers 1 and 2 are located in the voxels indexed by R1 and R2 . The indices T1 , T2 , R1
and R2 are integers in the interval [1, N ] and are assumed to be distinct. (Note that it is
simple to modify the model so that a transmitter or a receiver occupies multiple voxels.)
6. We assume both receivers 1 and 2 use the same type of receptors R and these receptors are
fixed in space, i.e. do not diffuse. With the simplification mentioned in the introductory
part of this section, we assume that the reaction between the signalling molecule L and the
complex C is:
k+

−
*
L−
)
−
−C
k−

where k+ and k− are the macroscopic reaction rate constants. With this assumption, the
rate of formation of complexes at receiver 1, at any time, is proportional to the number of
signalling molecules L in the R1 -th voxel. The situation at receiver 2 is similar.
Remark 1 The reaction kinetics assumed above can be viewed as a linearisation of the
−
*
second order reaction L + R )
− C. Similar assumption is also used in [11, 12] to model
receptor kinetics in chemotaxis.
7. The state vector q consists of (N + 2) elements where

nL,1 nL,2 .... nL,N
q =

nC,1

nC,2



(8)

where nL,j represents the number of signalling molecules at the j-th voxel and nC,u represents
the number of complexes at the u-th receiver.
8. We define two indicator vectors 1T1 , 1T2 ∈ ZN +2 . The T1 -th element of 1T1 is 1 and are
otherwise zero. 1T2 is similarly defined.
9. In order to write down the diffusion and reaction within the molecular communication network, we define the following state transition vectors rj and transition rates Wj . The total
number of possible jumps in this system is J = 2N + 4 where 2N of them model diffusion
and the rest models reactions at the receivers. We will state these jumps below in four
categories. Note that all rj ∈ ZN +2 and only the non-zero elements of rj are stated, and q
is the state vector defined in (8).
(a) The diffusion of L from voxel j to j + 1, where 1 ≤ j ≤ N − 1, is modelled by rj and
Wj (q). Specifically, rj (j) = −1, rj (j + 1) = 1 and Wj (q) = d̆q(j) = d̆nL,j .
Explanation: If a signalling molecule diffuses from voxel j to j + 1, it means the number
of signalling molecules in voxel j is decreased by one (hence rj (j) = −1) and that in
voxel j + 1 is increased by one (hence rj (j + 1) = 1). The rate at which this particular
type of jumps takes place is proportional to the number of molecules in the j-th voxel,
which is given by the j-th element of the state vector q.
For convenience, we define rN to be a zero vector and WN (q) = 0.
(b) For the diffusion of L from voxel j to j − 1 where 2 ≤ j ≤ N , rN +j (j) = −1, rN +j (j −
1) = 1 and WN +j (q) = d̆q(j) = d̆nL,j . For convenience, we define rN +1 to be a zero
vector and WN +1 (q) = 0.
(c) For receiver 1, the vector r2N +1 and the rate W2N +1 (q) are used to model the forward
reaction of the conversion of a signalling molecule L to a complex C. Specifically,
k+
k+
q(R1 ) = ∆x
nL,R1 .
r2N +1 (R1 ) = −1, r2N +1 (N + 1) = 1 and W2N +1 (q) = ∆x
Explanation: In the forward reaction, a signalling molecule is removed in the R1 -th
voxel, hence r2N +1 (R1 ) = −1 and a complex is formed, hence r2N +1 (N + 1) = 1
5

because the number of complexes at receiver 1 is the (N + 1)-th element in the state
vector (8). The rate W2N +1 (q) is proportional to the number of signalling molecules in
the R1 -th voxel where the receiver is located.
For the reverse reaction, r2N +2 (N +1) = −1, r2N +2 (R1 ) = 1 and W2N +2 (q) = k− q(N +
1) = k− nC,1 .
(d) For receiver 2, the forward reaction: r2N +3 (R2 ) = −1, r2N +3 (N +2) = 1 and W2N +3 (q) =
k+
k+
∆x q(R2 ) = ∆x nL,R2 . For the reverse reaction, r2N +4 (N + 2) = −1, r2N +4 (R2 ) = 1
and W2N +4 (q) = k− q(N + 2) = k− nC,2 .
The RDMEX model for the 2-transmitter 2-receiver molecular communication network is
dP (q, t)
dt

=

2 X
∞
X
a=1 b=1

+

J
X
j=1

{P (q − ka,b 1Ta ) − P (q, t)}δ(t − ta,b )

Wj (q − rj )P (q − rj , t) −

J
X

Wj (q)P (q, t)

(9)

j=1

where δ(t) denotes the Dirac delta function.
Let us, for the time being, assume that the first term on the right-hand side of (9) is not there.
If this is the case, then (9) is of the same form as the master equation (1) and the equation models
a Markov process. Given this model includes both reaction and diffusion, equation (9) without
the first term is known in the literature as reaction-diffusion master equation (RDME).
The novelty of the RDMEX model is the introduction of the first term on the right-hand side
of (9). This term can be viewed as a deterministic input because molecules are emitted by the
transmitters at pre-determined times. Let us look at this term more closely. At time ta,b , the
a-th transmitter emits ka,b signalling molecules into the Ta -th voxel (where the a-th transmitter
is located). This means that if the system is in the state q just before the time ta,b (denoted as
+
t−
a,b ), then it will be in state q + ka,b 1Ta just after ta,b (denoted as ta,b ). In addition, we have
+
P (q, t−
a,b ) = P (q + ka,b 1Ta , ta,b ), which is modelled by the first term in (9). Note that it is possible
to give a stochastic interpretation of ka,b , see Remark 2.
The deterministic input in (9) can be thought of as an external arrival of the system. We will
refer to (9) as reaction-diffusion master equation with exogenous input, or RDMEX for short. The
name is inspired by time series models such as ARX and ARMAX [15].
Note that the RDMEX model is no longer a Markov process due to the deterministic arrivals.
However, RDMEX is piecewise Markovian in the sense that, it is Markovian in between two
consecutive deterministic arrivals.
We have given an example of RDMEX for a 2-transmitter 2-receiver model in 1-dimension.
The model can be readily generalised to include more transmitters and receivers. In order to
generalise the model to 3-dimensional space, we will need to divide the space into 3-dimensional
cubic voxels of equal volume. (The use of more complicated geometry is possible, see [16].) The
molecules in a voxel are only allowed to diffuse to any of its neighbouring voxels. This can also be
readily be done. Lastly, we remark that it is also possible to use more complex receiver structure
or to consider non-isotropic medium.

3.2

Mean and covariance of receiver output in the RDMEX model

We will now generalise the result of Proposition 1 to the case of RDMEX model.
Proposition 2 For the RDMEX model in (9), assuming that Wj (q) is a linear function of q. Let

6

PJ

j=1 rj Wj (q)

= Aq, then
dhQ(t)i
dt

=

dΣ(t)
dt

=

AhQ(t)i +

2 X
∞
X
a=1 b=1

ka,b 1Ta δ(t − ta,b )

AΣ(t) + Σ(t)AT +

J
X

rj rjT Wj (hQ(t)i)

(10)

(11)

j=1

Proof: For the time evolution on the mean hQ(t)i, we can start with the derivative of (3):
P
(q,t)
(q,t)
dhQ(t)i
= q q dPdt
and then use (9) for dPdt
. This is fairly straightforward and uses exactly
dt
the same argument as the proof in [10].
Alternatively, one can argue the correctness of (10) as follows. Given that the difference
between (1) and (9) is the deterministic arrivals modelled by impulses, this means that between
two consecutive deterministic arrivals, the evolution of the state in RDMEX can be described by
a standard master equation. Hence, (6) holds between two consecutive deterministic arrivals. It
can be readily shown that the effect of ka,b molecules arriving at time ta,b is to add ka,b 1Ta to the
state vector. Hence the form of (10).
For the evolution of covariance, one can follow the derivation in [10] provided that the impulses are handled correctly because the multiplication of Dirac deltas (or distributions) is not
well defined. However, one can argue the correctness of (11) using the same argument in the
last paragraph. We know that between two consecutive deterministic arrivals, (7) holds for the
RDMEX model. It remains to show that the covariance matrix just before a deterministic arrival
is equal to that just after the deterministic arrival.
−
Let Q(t−
a,b ) and hQ(ta,b )i be the state and mean state just before the deterministic arrival at
+
time ta,b . At time ta,b , just after ta,b , the state of the system will become Q(t−
a,b ) + ka,b 1Ta . Also,
+
the mean state at ta,b is:
hQ(t+
a,b )i =
=

X

qP (q, t+
a,b ) =

q

X
q

qP (q − ka,b 1Ta , t−
a,b )

X
−
(q + ka,b 1Ta )P (q, t−
a,b ) = hQ(ta,b )i + ka,b 1Ta

(12)

q

−
Note that we have used the fact that P (q, t+
a,b ) = P (q − ka,b 1Ta , ta,b ) in the above derivation. The
overall result is that, at time ta,b , the state Q(t) and mean state hQ(t)i are incremented by the
same vector ka,b 1Ta .
Given that, at any deterministic arrival, both the state and mean state change by exactly the
same amount, therefore, deterministic arrivals do not cause discontinuity in covariance. Hence
(11).

For the rest of the paper, we will focus on studying the properties of equation (10), though we
will present a numerical example in Section 6 to demonstrate the accuracy of (11). A detail study
on (11) is also important and will be done in a future paper.

Remark 2 We will now briefly discuss a generalisation of the RDMEX model. Instead of assuming a deterministic emission of exactly ka,b signalling molecules by the a-th transmitter at time
ta,b , one may assume that the number of molecules emitted is a random variable Ka,b with mean
hKa,b i and covariance cov(Ka,b ). Provided that the random variable Ka,b is independent of the
state q at time ta,b or earlier, similar results to Proposition 2 can be derived. For equation (10), we
need to replace ka,b by hKa,b i, and we need to add cov(Ka,b ) to the right-hand side of (11). This
generalisation says that one can interpret ka,b in (10) as the mean number of molecules emitted
at time ta,b by the a-th transmitter. With this stochastic interpretation of ka,b , one can consider
the signalling molecules are generated by an irreversible chemical reaction.

7

4

Continuum limit of RDMEX

In section 3.1, we present an example of the RDMEX model for a 2-transmitter, 2-receiver molecular communication network in an isotropic 1-dimensional medium. We also show that if the
reaction kinetics at the receiver is linear, then the mean number of molecules in the network
evolves according to the ODE (10). In section 4.1, we will determine the continuum limit of (10)
as ∆x → 0. In order to simplify the presentation, we have so far limited our study to 1-dimensional
but given most molecular communication networks are 3-dimensional, we will generalise the continuum limit to 3-dimensional case as well.
The continuum limit of the RDMEX is in fact a RDPDE. A nice property of the resulting
RDPDE is that a closed form solution is available. This closed form solution shows that the
output signal of a receiver can be affected by the presence of other receivers in the network. This
will be discussed in section 4.2.

4.1

Continuum limit and generalisation to 3-dimensional space

In this section, we will study the continuum limit of equation (10) and show that in the limit,
when the interval ∆x goes to zero, (10) converges to a RDPDE and a number of chemical kinetics
ODEs. In order that we will be able to solve the RDPDE analytically later on, we will assume
from now onwards that the 1-dimensional medium is infinite, which means that the molecules in
each voxel can diffuse to either of its neighbouring voxel and the state vector q is


.... nL,−2 nL,−1 nL,0 nL,1 nL,2 .... nC,1 nC,2
q =
(13)
where, as before, nL,j is the number of L in the j-th voxel where j ∈ Z, and nC,u is the number
of complexes formed at the u-th receiver. Given this state vector, we can write equation (10) as
dhnL,j (t)i
dt

=

d̆(hnL,j−1 (t)i − 2hnL,j (t)i + hnL,j+1 (t)i) +
−

dhnC,u (t)i
dt

=

2
X
u=1

δK (j − Ru )(

∞
2 X
X
a=1 b=1

δK (j − Ta )ka,b δ(t − ta,b )

k+
hnL,Ru (t)i − k− hnC,u (t)i) ∀j ∈ Z
∆x

k+
hnL,Ru (t)i − k− hnC,u (t)i for u = 1, 2
∆x

(14)
(15)

where δK (j) is the Kronecker delta1 .
Suppose the centre of the voxel j is at position xj , we replace hnL,j (t)i by `(xj , t)∆x where
`(xj , t) is the mean concentration in voxel j at time t. By dividing both sides of (14) by ∆x,
taking the limit ∆x → 0 and noting that d̆(∆x)2 = D̆, we have
∂`
∂t

= D̆

2
∞
2
X
X
∂2` X
dhnC,u (t)i
+
δ(x
−
x
)
k
δ(t
−
t
)
−
δ(x − xR,u )
(16)
T,a
a,b
a,b
∂x2 a=1
dt
u=1
b=1
|
{z
}
=ka (t)

dhnC,u (t)i
dt

= k+ `(xR,u , t) − k− hnC,u (t)i for u = 1, 2

(17)

where xT,a (resp. xR,u ) is the centre of the voxel Ta (Ru ) where the a-th transmitter (u-th receiver)
is located. Note that we have also used the following conversion between the Kronecker and Dirac
(j)
deltas: lim∆x→0 δK∆x
= δ(x).
This shows that in the continuum, the RDMEX converges to a RDPDE (16) and a number of
ODEs describing the kinetics at the receivers (17). The RDPDE (16) has a simple interpretation.
1 Note: We use both Konecker delta and Dirac delta in this paper. They are denoted, respectively, as δ (•) and
K
δ(•).

8

The second term in (16) says that the transmitter at Ta adds signalling molecules to the system
according to time sequence ka (t), while the third terms says the signalling molecules are absorbed
from the system if they form complexes at the receivers. Given that we assume that one signalling
molecule reacts to form one complex, therefore the rate of absorption of signalling molecules is
equal to the rate of complex formation, which is given by (17).
Note that it is well known in literature, see [6, 16], that a RDME with linear reaction rates
converges to a RDPDE in continuum. In the above, we show analogues result holds for the
RDMEX model.
4.1.1

Generalisation to 3-dimensional space

In order to simplify the presentation, we have so far limited to the 1-dimensional case. Given
that most molecular communication networks are 3-dimensional, we state that in a 3-dimensional
infinite medium the RDMEX will converge to the following RDPDE in the continuum. Here v
denotes a point in the 3-dimensional space (i.e. v is a 3-dimensional vector) and `(v, t) is the mean
concentration of the signalling molecule at the location v at time t.
∂`
∂t

=

dcu (t)
dt

=

D∇2 ` +

2
X
a=1

δ(v − vT,a )ka (t) −

2
X
u=1

δ(v − vR,u )

dcu (t)
dt

k+ `(xR,u , t) − k− cu (t) for u = 1, 2

(18)
(19)

where ∇2 is the Laplacian in 3-dimensional space, and vT,a (resp. vR,u ) is a 3-dimensional vector
specifying the location of a-th transmitter (u-th receiver). Note that we have also introduced a
new notation cu (t) to denote the mean number of complexes hnC,u (t)i at receiver u; this is to
simplify the notation later.
The derivation assumes that the molecule in a voxel diffuses to a neighbouring voxel at a rate
of d per molecule per unit time and each voxel is a cube of size ∆3 . The parameter d is related
to the 3-dimensional macroscopic diffusion constant D by d = ∆D2 . Also, the rate of formation
k+
of complexes at the receiver voxel is given by ∆
3 times the number of signalling molecules L in
the receiver voxel. Given that the derivation for the 3-dimensional is essentially the same as the
1-dimensional case, we do not present it here.

4.2

Solution to RDPDE

In this section, we present a solution to the RDPDE (18) and (19), which is the continuum
limit of the 3-dimensional RDMEX model. The key result is a closed-form expression of the
multivariate transfer function from the transmitter signals k1 (t) and k2 (t) (which models the
number of molecules injected by the transmitters into the medium at time t and can be viewed
as the inputs to the system) to the mean number of complexes formed at the receivers c1 (t) and
c2 (t) (which can be viewed as the outputs). We will divide this section into two parts. We will
first derive the transfer function and then provide an interpretation of the transfer function.
4.2.1

Derivation of transfer function

The aim of this part is to derive a multivariate transfer function from k1 (t) and k2 (t) to c1 (t) and
c2 (t) using (18) and (19).
√
We first define a few notation. Let ι = −1, and Cu (ω), Ka (ω) and L̃(v, ω) be the temporal
Fourier transform of, respectively, cu (t), ka (t) and `(v, t) where ω is the transform variable. It is
shown in Appendix A that
L̃(v, ω)

=

2
X
a=1

φ(v − vT,a , ω)Ka (ω) −
9

2
X
u=1

φ(v − vR,u , ω)ιωCu (ω)

(20)

where
φ(v, ω)

=

1
exp(−
4πDkvk

r

ιω
kvk)
D

(21)

is the temporal Fourier transform of the 3-dimensional diffusion kernel

2

1
3
(4πDt) 2

exp(− kvk
4Dt ). Given

(20) holds for any location v, we can use it to determine the concentration of the signalling
molecules at the two receivers. By substituting v = vR,1 and then v = vR,2 in (20), we have:
L̃(vR,1 , ω)
L̃(vR,2 , ω)

= φ11 (ω)K1 (ω) + φ12 (ω)K2 (ω) − φ0 (ω)ιωC1 (ω) − φ∆R (ω)ιωC2 (ω)

(22)

= φ21 (ω)K1 (ω) + φ22 (ω)K2 (ω) − φ∆R (ω)ιωC1 (ω) − φ0 (ω)ιωC2 (ω)

(23)

where
φua (ω)

=

φ∆R (ω)

=

φ0 (ω)

=

φ(vR,u − vT,a , ω) for u, a = 1, 2

(24)

φ(0, ω) where 0 = zero vector

(26)

φ(vR,1 − vR,2 , ω)

(25)

One can interpret φua (ω) as the transfer function which models the dynamics of the diffusion
of molecules from the a-th transmitter located at vT,a (where the molecules are injected into the
medium) to the u-th receiver located at vR,u .
When signalling molecules are absorbed to form complexes, it creates a concentration gradient.
The diffusion dynamics between the two receivers are modelled by φ∆R (ω). Given that we assume
that the locations of the transmitters and receivers are distinct, both φua (ω) and φ∆R (ω) are well
defined.
The transfer function φ0 (ω) models the local impact of absorption of signalling molecule at each
receiver. This transfer function is unfortunately not well defined as can be seen from substituting
v = 0 in the definition of φ(v, ω) in (21). The transfer function φ0 (ω) also appears in the modelling
of receptor noise in chemotaxis in the biophysics literature [11, 12]. In fact [11, Equation (18)]
and [12, Equation (6)] are special cases of (18) where the input terms ka (t) are absent. Both
[11, 12] deal with the indefiniteness of φ0 (ω) by cutting off an integral to evaluate φ0 (ω) at a finite
frequency. However, this requires us to make an assumption on the size of the receptor molecule.
Instead, in section 5, we derive a new method to approximate φ0 (ω), and we will show using
numerical examples in section 6 that our approximation gives accurate results. For the rest of this
section, we will continue to use equations (22) and (23) with the understanding that φ0 (ω) is not
well defined but can be well approximated.
Both equations (22) and (23) are obtained from (18). We still need to work on (19). By taking
the Fourier transform of (19), we have
ρ1 (ω)

=

ρ2 (ω)

=

k+
C1 (ω)
=
ιω + k−
L̃(vR,1 , ω)
C2 (ω)
k+
=
ιω
+ k−
L̃(vR,2 , ω)

(27)
(28)

where ρ1 (ω) and ρ2 (ω) are transfer functions that model the reaction kinetics at the receivers.
Given that we have assumed that the binding and unbinding rate constants at both receivers are
identical, it is not surprising that ρ1 (ω) and ρ2 (ω) are the same here. It is straightforward to
generalise to the case where the receivers have different reaction kinetics.
By using equations (22), (23), (27) and (28), we can eliminate L̃(vR,1 , ω) and L̃(vR,2 , ω) to
obtain the transfer function from the inputs k1 (t) and k2 (t) to the outputs c1 (t) and c2 (t):

−1


C1
C2


=





ρ1 0
φ0 φ∆r 

I + ιω

0 ρ2
φ∆r φ0 

|
{z
}|
{z
}
R

Φ0

10



ρ1
0

0
ρ2


|

φ11
φ21

φ12
φ22
{z
Φ


}

K1
K2


(29)

where, for conciseness, we have dropped the dependence on transform variable ω. Equation (29)
is the key result of this section. It is the solution to the RDPDE (18) and (19), which are in turn
the continuum limit of the mean concentration in the stochastic RDMEX model. We will provide
some physical interpretation of (29) in a moment. Before that, we want to point out that (29)
can be used to compute c1 (t) and c2 (t) given k1 (t), k2 (t) and the system parameters by numerical
Laplace transform. This will be done in Section 6.
Note that it is numerically more efficient to solve for the receiver outputs using (29) rather
than (10). This is because one also needs to solve for the number of signalling molecules in the
voxels in (10) but this is not needed when (29) is used. Since the number of voxels is much larger
than the number of transmitters and receivers, numerical solution via (29) is more efficient.
4.2.2

Interpretation of transfer function (29)

Equation (29) may not look easy to interpret in the first instance, so we will first specialise it to
the case of 1-transmitter and 1-receiver. In this case, we have
C1

=

ρ1 φ11
K1
1 + ιωφ0 ρ1

(30)

One can readily show that the input-output transfer function in (30) corresponds to the block
diagram representation of figure 3. The negative feedback loop occurs because the net number
of signalling molecules at the receiver `(vR,1 , t) is given by the difference between those that
arrive via diffusion (modelled by the feedforward block of φ11 ) minus those reacted to form the
complexes (modelled by the feedback block of ιωφ0 ). The forward loop consists of φ11 which
models the diffusion dynamics of signalling molecule from the location of the transmitter to that
of the receiver, and ρ1 which models the conversion of the signalling molecules to complexes.
We will now take a closer look at the denominator of (30) with the aim to determine the
1
at low
strength of the feedback term ιωφ0 ρ1 . It has been shown in [11, 12] that φ0 ∝ 2πD
k+
frequency. Consider the case where D  k+ or D ≈ 0, which corresponds to the situation where
the chemical kinetics is not diffusion-limited. We have ιωφ0 ρ1 ≈ 0 base on the expression of ρ1
in (27) and consequently C1 ≈ ρ1 φ11 K1 . Since the chemical kinetics is not diffusion-limited, the
chemical kinetics and diffusion are basically ”decoupled”, so the transfer function from K1 to C1
is the multiplication of the transfer function φ11 from K1 to L̃(vR1 , ω) (which models diffusion)
and the transfer function ρ1 from L̃(vR1 , ω) to C1 (which models reaction kinetics). We will show
ρ1 φ11
in section 6 using numerical examples to show that the transfer function 1+ιωφ
≈ ρ1 φ11 holds
0 ρ1
when D  k+ .
Given the interpretation of the 1-transmitter 1-receiver case as a system with feedback in figure
3, it can be shown that equation (29) corresponds to multivariate feedback system with 2 inputs
and 2 outputs. The block structure of the multivariate feedback system is the same as that in
figure 3 but we need to replace the single-input single-output transfer functions φ11 , ρ1 and φ0 by
their multivariate counterparts of Φ1 , R and Φ0 in (29).
For the 2-transmitter 2-receiver case, we see from (29) that the response at each receiver is
affected by both transmitters, as well as by the other receiver. Let us assume that the transmitters
and receivers form two communication pairs where transmitters 1 and 2 intend to communicate
with, respectively, receivers 1 and 2. We want to determine the unintended signal at the receivers.
Let us assume for the time being that the two receivers are sufficiently far apart so that φ∆r (ω)
is negligible compared with φ0 (ω). In addition, we assume the two receivers are identical, so
ρ1 (ω) = ρ2 (ω). In this case, we can simplify (29) to





ρ1
φ11 φ12
K1
C1
(31)
≈
C2
K2
1 + ιωφ0 ρ1 φ21 φ22
Comparing with the single-transmitter single-receiver transfer function in (29), we can see that
ρ1 φ12
the unintended signal due to transmitter 2 on receiver 1 is 1+ιωφ
K2 . One can readily see that
0 ρ1
the magnitude of this unintended signal can be reduced if transmitter 2 is well separated in space
11

from receiver 1 because φ is a decreasing function of distance. Similar conclusion can be drawn
for transmitter 1 and receiver 2. Note that the above argument requires that the receivers are well
separated. If this is not the case, the matrix inverse in (29) will create a complicated mixture of
signal at both receivers.
In general, spatial separation is a strategy to reduce the magnitude of unintended signal that
one communication pair has on the others. It is interesting to point out that (29) allows one to
explore other methods to reduce the magnitude of the unintended signal. An interesting case to
study is if the transmitters emit molecules at different frequencies and the receivers are frequency
sensitive. We will not explore this further here and leave this for future work.

5

Discrete solution for mean concentration in RDMEX

In the last section, we presented a closed-form solution to the RDPDE (18). A problem that
we face is that the frequency response φ0 (ω) in (26) is not well defined. This problem arises
because the size of voxel, in the continuum limit, becomes zero. Therefore, a solution to overcome
this problem is to consider finite voxel size instead. This means that we need to work with the
3-dimensional analogue of equations (14) and (15).
Consider an isotropic 3-dimensional space. We divide the space into identical cubic voxels of
volume ∆3 each. We index the voxel using a 3-dimensional vector ξ = [i j k] where i, j, k ∈ Z;
note that we will use ξ and [i j k] interchangeably in this section. We use ξT,a (ξR,u ) to index the
voxel that the a-th transmitter (u-th receiver) is located. We will also use `R,u (t) to denote the
mean concentration of signalling molecules at the voxel where the u-th receiver is located. Let
c̃u (t) denote the mean number of complexes at the u-th receiver for this discrete model and C̃u (ω)
be its continuous Fourier transform. (Just to avoid any possible confusion. We discretise only
space, not time. So, t remains continuous.)
Let `i,j,k (t) denotes the mean concentration of the signalling molecule in voxel i, j, k. The
mean concentration in a voxel is given by the mean number of molecules divided by the volume
of a voxel which is ∆3 . One can show that the generalisation of (14) and (15) to 3-dimensional
space — with mean concentration of signalling molecules, rather than mean number, per voxel —
is:
d`i,j,k (t)
dt

=

d(`i−1,j,k (t) − 2`i,j,k (t) + `i+1,j,k (t)) +
d(`i,j−1,k (t) − 2`i,j,k (t) + `i,j+1,k (t)) +
d(`i,j,k−1 (t) − 2`i,j,k (t) + `i,j,k+1 (t)) +

dc̃u (t)
dt

=

2
2
X
X
dc̄u (t)
1
1
δ
(ξ
−
ξ
)k
(t)
−
δ (ξ − ξR,u )
T,a a
3 K
3 K
∆
∆
dt
a=1
u=1

(32)

k+ `R,u (t) − k− c̃u (t) for u = 1, 2

(33)

where D = ∆d2 . One can readily show that the continuum limit of (32) and (33) is (18) and (19)
Let L̃d (ξ, ω) denote the temporal Fourier transform of `i,j,k (t). (Recall that ξ = [i, j, k].) It is
shown in Appendix B that

L̃d (ξ, ω)

=

2
X
a=1

ψ(ξ − ξT,a , ω)Ka (ω) −

2
X
u=1

ψ(ξ − ξR,u , ω)ιω C̃u (ω)

(34)

where
ψ(ξ = [i, j, k], ω)

=

1
2
4π D̃∆

I I
C

12

C

|k|+1

Wz∗
Wxi−1 Wyj−1 dWx dWy
2 −1
Wz∗

(35)

where Wx and Wy are complex variables and the contour C is the unit circle on the complex plane;
also, Wz∗ is the solution of the following quadratic equation in Wz with modulus less than unity:
Wz2 − (2 + (Wx − Wx−1 )2 + (Wy − Wy−1 )2 + ιω

∆2
)Wz + 1 = 0
D

(36)

Note that (34) is the discrete space analogue of (20). Thus one can identify φ(0, ω) with
ψ(0, ω). In addition, one can show that:


C̄1
C̄2




=


I + ιω

ρ1
0

0
ρ2



ψ0
ψ∆r

ψ∆r
ψ0

−1 

ρ1
0

0
ρ2



ψ11
ψ21

ψ12
ψ22



K1
K2


(37)

where for conciseness we have dropped the dependence on ω, and ψua (ω) = ψ(ξR,u − ξT,a , ω),
ψ∆r (ω) = ψ(ξR,1 − ξR,2 , ω) and ψ0 (ω) = ψ(0, ω).
Numerical integration can be used to compute ψ(ξ, ω) in (35). This will be used in Section 6.

6

Numerical examples

6.1

Overview

In this section, we will give a number of numerical examples to show that equations (29) in section
4.2 and (37) in section 5 can be used to accurately predict the mean output of the receivers
in the stochastic model RDMEX. We will also use these numerical examples to illustrate the
issues of using molecular signalling for communication. We will present three sets of results:
single-transmitter single-receiver in section 6.2, single-transmitter two-receiver in section 6.4 and
two-transmitter two-receiver in 6.5.
In order to verify the accuracy of (29) and (37), we will use simulation to compute the mean
receiver output. One method is to simulate RDMEX many times and compute the mean. Alternatively, one can use the fact that if the number of molecules is large, then the behaviour of one
simulation run is fairly close to the mean. We will mainly use the latter method in this paper.
We simulate the RDMEX model using the τ -leaping method [17]. The τ -leaping algorithm uses a
constant time step to advance the simulation and is a faster alternative to the Gillespie’s algorithm
[17]. We will refer to the simulation result as RDMEX-τ .
The fluidic medium is assume to have a D of 0.05. (Since the parameters in the diffusionreaction system can be scaled to some dimensionless quantities [18, Section 8.2], we do not specify
the units for the parameters here.) The locations of the transmitters and receivers, as well as the
reaction rate constants, vary between experiments and will be specified later. Different transmitter
signals will be used to demonstrate the accuracy of the RDMEX model. Our goal is to compare
the output c1 (t) and c2 (t) from RDMEX-τ with that from the following analytical models:
1. Equation (37). We will refer to this as RDMEX-M.
2. Equation (29) of the continuum model with φ(0, ω) replaced by ψ(0, ω). We will refer to
this as RDPDE.
1
. This approxi3. Equation (29) of the continuum model with φ(0, ω) approximated by 2πD∆
mation is inspired by the one used in [11, 12] where we have replaced the size of receptors
by the voxel dimension parameter ∆. We will refer to this as RDPDE-X.

4. For the single-transmitter and single-receiver case, we consider the “decoupled” model C1 =
ρ1 φ11 K1 . This will be referred to as DE.
For all these analytical models, we first determine the Laplace transforms C1 and C2 by using the Laplace transforms K1 and K2 , and the transfer function. We then invert the Laplace
transform numerically using the matlab function invlap.m [19].
13

6.2

Single-transmitter single-receiver case

The system consists of a transmitter at the voxel [0, 0, 0] and a receiver at voxel [3, 0, 0].
6.2.1

Model accuracy and the effect of k+

In this set of experiments, the transmitter emits 10 molecules every 10−4 time units for a duration
of 0.2 time units and then stops emitting for 0.3 time units. The signal k1 (t) is obtained by
concatenating this emission pattern 3 times. The value of k− is 0.05. We determine the output
signal c1 (t) for t ∈ [0, 2].
We use two different values for k+ . Figure 4 shows the results for k+ = 2.5 × 10−3 . The figure
compares the mean number of complexes formed in the time interval [0, 2]. Results are obtained
from RDMEX-τ (simulation) and RDPDE, RDMEX-M, RDPDE-X and DE (analytical models).
It can be seen that both RDPDE and RDMEX-M (our analytical solutions) match RDMEX-τ
well. The model RDPDE-X does not give good approximation because the voxel size is not a good
approximation for the receptor size. Henceforth, we will not consider RDPDE-X further. The
decoupled model DE does not match RDMEX-τ for this value of k+ .
We then change k+ to 2.5 × 10−4 . The results are plotted in Figure 5. Both RDMEX-M and
RDPDE again match RDMEX-τ well. We also see that the decoupled model DE gives better
prediction than before. This validates the discussion in Section 4.2.2 that the decoupled model
holds when k+ is sufficiently small.
6.2.2

Accuracy of mean and standard deviation computation

In this experiment, we use a different transmitter signal to show the accuracy of RDMEX-M and
RDPDE. We define two transmitter symbols s0 and s1 . The symbol duration is 2 time units.
When a transmitter sends s1 , it emits 10 molecules every 10−4 time units for a duration of 0.2
time units and then stops emitting for 1.8 time units. When a transmitter sends s0 , it does not
emit any molecules for 2 time units. The transmitter signal k1 (t) in this experiment is simply s1 .
The receiver parameters are k+ = 2.5 × 10−3 and k− = 8. We determine the output signal c1 (t)
for t ∈ [0, 2].
We first verify the accuracy of using RDMEX-M and RDPDE to compute the mean of receiver
output. For this experiment, we simulate the system using RDMEX-τ 125 times and compute the
mean receiver output as the reference. Figure 6 shows the mean receiver output from RDMEX-M,
RDPDE and RDMEX-τ . It can be readily seen that RDMEX-M and RDPDE are accurate also
for a different transmitter signal.
Our next goal is to verify the accuracy of using equation (11) to determine the standard
deviation of the receiver output. We extend the RDMEX-M to solve for both the receiver output
as well as the system states, which are the number of signalling molecules in the voxels. The
system states are used as the input to (11) and numerical integration is used to solve for (11). For
verification, we simulate the system using RDMEX-τ 125 times to compute the standard deviation
of the receiver output. Figure 7 plots the standard deviation of the receiver output from the two
methods. It can be seen that (11) is accurate.
One can envisage using the symbols s1 and s0 to encode the communication between the
transmitter and the receiver. The communication scheme is similar to ON-OFF keying. For
decoding, the receiver can use, say, the peak number of complexes to detect the symbol transmitted.
If the peak number of complexes is above a threshold, then s1 has been transmitted; otherwise,
s0 has been transmitted.
Remark 3 The above method of calculating the variance using (11) is computational intensive.
A future work is to derive efficient algorithm to solve (11).

14

6.2.3

Assumption of linear receiver model

In this paper, we assume that the rate of formation of complex C is a linear function of the number
of signalling molecules in the receiver voxel, according to the following chemical reaction kinetics:
k1+

−−
*
L)
−−
−
−C
k1−

(38)

Under certain assumptions, the above reaction kinetics can be used to approximate more
complex reactions. Consider the following chemical reactions:
g1+

g2

−−
*
L+E)
−−
−
− I −→ C + E
g1−

g3

C −→ L

(39)
(40)

Reaction (39) is of Michaelis-Menten [20] type where a molecule L reacts with an enzyme E to form
an intermediate product I, following by the decomposition I into a product C and the enzyme
E. It can be shown that, for suitable choice of reaction constants in (39), the reaction kinetics of
(39) can be approximated by the forward reaction in (38). Similarly, reaction (40) can be made
identical to the reverse reaction in (38) by choosing g3 = k− .
We want to show that a receiver with reactions (39) and (40) gives similar output signal
compared to one with reactions (38). We implement the receiver kinetics (39) and (40) in τ -leaping
simulation and refer to this method as RDMEX-τ -MM. We compare this against RDMEX-M with
identical parameters to those used in section 6.2.2. The parameters in reactions (39) and (40)
have been chosen to approximate those in (38). The results are plotted in Figure 8 and it can be
seen that RDMEX-M is able to approximate more complicated reactions.

6.3

Effect of ∆

In this section, we study the effect of the size of voxels. We use the same parameter setting for
the single-transmitter single-receiver case in Section 6.2.1. We assume both the transmitter and
receiver is a cube whose length of an edge is χ = 0.01. We use four different voxel sizes with
∆ = χ, χ2 , χ3 and χ4 . Since the size of the transmitter (or receiver) is a constant, this means that
the transmitter (or receiver) occupies, respectively, 1, 8, 27 and 64 voxels for these 4 different voxel
sizes. We assume that the emissions from the transmitter is uniformly distributed across all the
voxels that it occupies. Figure 9 show the mean number of complexes at the receiver, which is the
sum of the number of complexes at all receiver voxels. It can be seem that the predicted output
for ∆ = χ2 , χ3 and χ4 are almost the same. This shows that as long as ∆ is sufficiently small, the
prediction is independent of ∆.

6.4

Single-transmitter two-receiver case

Equation (29) shows that when there are multiple receivers, it is possible for a receiver to affect
the output of another receiver. We will illustrate this phenomenon. We consider three different
networks. Network 0 consists of a transmitter and two receivers, called 1 and 2. The transmitter,
receivers 1 and 2 are located, respectively, at voxel [0, 0, 0], [1, 0, 0] and [2, 0, 0]. Network 1 is
composed of the transmitter and receiver 1 of network 0. Network 2 consists of the transmitter
and receiver 2 of network 0. For all the three networks, the transmitted signal is s1 and the
receiver parameters are k+ = 2.5 × 10−3 and k− = 8.
Figure 10 shows the output for receiver 1 for networks 0 and 1, while Figure 11 shows the output
for receiver 2 for networks 0 and 2. (The curve with labelled 1t15r in Figure 11 will be explained
later.) It can be seen that the output of receiver 1 is almost the same for both network 0 (receiver
2 present) and network 1 (receiver 2 absent). However the output of receiver 2 for network 0
(receiver 1 present) is very different from that in network 2 (receiver 1 absent). Specifically, when
receiver 1 is present, the output of receiver 2 has a lower peak number of complexes and a higher
number of complexes at the tail.
15

An explanation of how receiver 1 affects the output of receiver 2 is as follows. Note that
receiver 1 is situated in between the transmitter and receiver 2. Some signalling molecules that
reach receiver 2 have to pass through the voxel containing receiver 1. When these signalling
molecules are in the receiver 1 voxel, some of them react to form complexes and are held up in
the voxel. This means less signalling molecules reach receiver 2 in the early part of the symbol
duration, thus resulting in a lower peak number of complexes in receiver 2. During the later part
of the symbol duration, the complexes in receiver 1 dissociate to release the signalling molecules.
Some of these signalling molecules, which are held up earlier in receiver 1, reach receiver 2 later
on. This means more signalling molecules reach receiver 2 in the later part of the symbol duration
if receiver 1 is present. This explains the behaviour at the tail of the output of receiver 2.
When the number of receivers in a network is high, the output of some receivers in a network
can be highly affected by the presence of the other receivers. We create a network with a large
number of receivers by adding an additional 14 receivers to network 2, making a total of 15 receivers
in the network. We again focus on the output of receiver 2. The curve labelled with 1t15r in
Figure 11 shows the output of receiver 2 for this 1-transmitter 15-receiver network. It can be seen
that the receiver output in this network is very different from that receiver 2 in network 2. The
reason is that the other 14 receivers are affecting the output of receiver 2.
The above results mean that the design of molecular communication need to take all receivers
in the network into consideration. For example, if receiver 2 uses a threshold on peak number of
complexes to detect s1 , we can see from Figure 11 that a threshold that works for the 1-transmitter
2-receiver network may not necessary work for the 1-transmitter 15-receiver network.

6.5

Two-transmitter Two-receiver case

We consider a molecular communication network with 2 transmitters located at [0, 0, 0] (transmitter 1) and [4, 0, 0] (transmitter 2), and 2 receivers located at [1, 1, 0] (receiver 1) and [2, 0, 1]
(receiver 2). Both transmitters use s1 and s0 defined earlier as their transmission symbols. The
signals for transmitters 1 and 2 are, respectively, s0 s1 and s1 s0 .
We first verify the accuracy of our proposed models. Figures 12 and 13 show the output signals
for, respectively, receivers 1 and 2. The three curves in each figure are obtained from RDMEX-τ
(simulation) and RDPDE and RDMEX-M (our analytical solutions). It can be seen from both
figures that the prediction from both RDPDE and RDMEX-M match that of RDMEX-τ . In fact,
the curves in the figures match so well that they overlap.
We now assume the transmitters and receivers form two unicast communication pairs: transmitter 1 communicates with receiver 1 while transmitter 2 communicates with receiver 2. If
transmitter 1 were the only transmitter, then receiver 1 should have a zero signal in the first
symbol duration (between 0 and 2 time units). However, Figure 12 shows that receiver 1 has a
non-zero signal during the first symbol duration due to the transmitter 2 sending an s1 during
this time. During the second symbol duration, the output signal of receiver 1 is due entirely to
transmitter 1. If receiver 1 uses a threshold based detector, then a suitable choice of threshold
will enable receiver 1 to correctly decode the two symbols sent by transmitter 1.
Let us now consider the output signal of receiver 2 shown in Figure 13. The signal in the first
symbol duration is due to transmitter 2 (the intended signal) while that in the second symbol is
due to transmitter 1 (the unintended signal). If receiver 2 uses a threshold based detector, then due
to the unintended signal, a bit error will occur in the second symbol duration. This is a typical
example of bit error when multiple transmitters and receivers share a common communication
channel.

7

Related work

Molecular communication networks can be divided into two categories, according to whether
they are natural or synthetic. Natural molecular communication networks are prevalent in living
organisms. Their synthetic counterparts, though still rare, do exist. For example, [21] presents

16

a system with multiple genetically engineered cells that use cell signalling to coordinate their
behaviour.
The modelling of natural and synthetic molecular communication networks is studied in different disciplines. The former is mainly studied in biophysics and mathematical physiology, while
the latter in synthetic biology. There is also a recent interest in the engineering community to
study molecular communication networks from a communication theory point of view [13, 22, 23].
This gives rise to a new research area called nano communication networks [1].
Despite the fact that molecular communication networks are studied in diverse disciplines, the
set of mathematical models that are being used are similar. This is not surprising given that the
primary goal is to model diffusion and reaction kinetics. The classes of mathematical models being
used include molecular dynamics, master equation, partial differential equation (PDE), FokkerPlanck equation, Langevin equation and others [6]. We will focus on the first three classes of
models in this discussion.
Molecular dynamics is commonly used in simulation of molecular communication networks.
Many examples of simulators exist, especially for natural molecular communication networks, see
[24] for a recent overview. For synthetic networks, a recent example is [25]. By analysing the
molecular dynamics of transmitters and receivers, [26] characterises the noise in transmitters and
receivers as, respectively, sampling and counting noise.
There are ample examples in using PDE — in particular diffusion PDE, telegraph equation
and RDPDE — to model molecular communication. For natural networks, [11, 12] use RDPDE
to study the noise in receptor binding in chemotaxis, and [27] uses RDPDE to study signalling
cascades. However, these papers do not consider the transmitters. For synthetic networks, telegraph or diffusion PDEs (or their kernels) have been used to characterise the diffusion of signalling
molecules in [28, 14, 29, 22, 23] and others. However, these papers do not consider the coupling
effect between diffusion and receiver reaction kinetics. In our earlier work in [30], we use a RDPDE
in the form of (18), as a deterministic model for molecular communication network. The RDPDE
in [30] is solved numerically and no analytic solution is provided. In this paper, we derived a RDPDE model (18) for molecular communication and provide an interpretation of the model as the
mean receiver output of molecular communication networks. In addition, we present an analytical
solution to this RDPDE and show that it can be used to accurately predict mean receiver output
in molecular communication networks.
For some time, RDME has been considered to be a phenomenological model because it diverges
in certain cases [31]. Fortunately, the problem has been resolved in [32] and there is now a firm
theoretic basis for RDME. There are many examples of work that use RDME to model natural
molecular communication networks, see [16, 33]. However, these papers do not consider the
transmitters. The use of RDME in studying synthetic molecular communication networks appear
to be novel. To the best of our knowledge, our RDMEX model, which is formed by coupling time
sequences of signalling molecule emission pattern with RDME, has not been proposed before. The
proposed RDMEX model is one of the novel contributions of this paper.

8

Conclusions and future work

In this paper, we have proposed a new stochastic model called reaction-diffusion master equation
with exogenous input (RDMEX) for modelling molecular communication networks with multiple
transmitters and multiple receivers. We show that we can readily derive the mean and covariance
of receiver output of RDMEX model for the case where reaction kinetics at the receiver is linear.
We then solve the mean receiver output of RDMEX model using two different methods. In the first
method, we derive the continuum limit of RDMEX and present a closed-form expression of the
solution to the resulting reaction-diffusion partial differential equation. In the second method, we
solve the mean receiver output of RDMEX explicitly by using transform techniques. We present
numerical examples comparing the accuracy of our analytical solutions against simulation. We
find our analytical solutions give accurate prediction of mean receiver output of molecular communication networks. This paper has focused on studying the mean receiver output of molecular

17

communication networks. However, one needs to understand the properties of noise in these networks in order to evaluate their performance, such as bit error rate or capacity. Future work
includes the study of the properties of noise in the RDMEX model and using the RDMEX model
to evaluate the performance of molecular communication networks.

A

Derivation of (20)

Equation (20) can be derived by using a couple of different methods, e.g. Fourier transform or
Green’s function. We will use Green’s function here.
Note that the RDPDE (18) is an inhomogeneous partial differential equation where the last
two terms act as the forcing function. Let G(v, t) denote the three dimensional kernel of the linear
2
diffusion equation ∂`
∂t = D∇ `. Let also ∗S,T denote convolution in space and time, and ∗T denote
convolution in time only. By using the theory of Green’s function, we have
(
`(v, t)

=
=

G(v, t) ∗S,T
2
X
a=1

2
X

2
X

dcu (t)
δ(v − vT,a )ka (t) −
δ(v − vR,u )
dt
a=1
u=1

G(v − vT,a , t) ∗T ka (t) −

2
X
u=1

G(v − vR,u , t) ∗T

)

dcu (t)
dt

(41)

Now we can obtain (20) from applying temporal Fourier transform to (41), noting that φ(v, ω)
is the temporal Fourier transform of G(v, t).

B

Derivation of (34)

Note that (32) is a linear difference equation in space and a continuous differential equation in
time. We can draw a parallel with the theory of Green’s function and note that (34) holds if
ψ(ξ, ω) is the temporal Fourier transform of the solution of the following equation:
dhi,j,k (t)
dt

=

D(hi−1,j,k (t) − 2hi,j,k (t) + hi+1,j,k (t)) + D(hi,j−1,k (t) − 2hi,j,k (t) + hi,j+1,k (t)) +
D(hi,j,k−1 (t) − 2hi,j,k (t) + hi,j,k+1 (t)) +

1
δK (ξ)δ(t)
∆3

(42)

This equation can be solved by using Fourier transform on the continuous variable t and ztransform on the discrete variables i, j and k. Let Hi,j,k (ω) denote the Fourier transform of
hi,j,k (t). The multi-dimensional z-transform of Hi,j,k (ω) is defined as:
Ψ(Wx , Wy , Wz , ω)

=

∞
X

∞
X

∞
X

Hi,j,k (ω)Wx−i Wy−j Wz−k

(43)

i=−∞ j=−∞ k=−∞

By using (42), it can be shown that
Ψ(Wx , Wy , Wz , ω)

=

D∆{(Wx −

Wx−1 )2

−1
+ (Wy − Wy−1 )2 + (Wz − Wz−1 )2 −

ιω∆2
}
D̃

(44)

We first compute inverse z-transform with respect to Wz . We re-write (44) as
Ψ(Wx , Wy , Wz , ω)

=

Wz
−1
−1
D∆ (Wz − Wz∗ )(Wz − Wz∗
)

(45)

where Wz∗ , as defined in section 5, is the root of (36) with modulus less than unity. Note we have
also used the fact that the roots of (36) are reciprocal of each other. It can now be shown that
18

the inverse z-transform of Ψ(Wx , Wy , Wz , ω) with respect to Wz is
|k|+1

−1 Wz∗
2 −1
D∆ Wz∗

(46)

If we apply the standard inverse z-transform contour integrals, with respect to Wx and Wy , to
this expression, then we obtain the formula for ψ(v, ω) given in (35).

References
[1] I. Akyildiz, F. Brunetti, and C. Blázquez, “Nanonetworks: A new communication paradigm,”
Computer Networks, vol. 52, pp. 2260–2279, 2008.
[2] M. Moore, T. Suda, and K. Oiwa, “Molecular Communication: Modeling Noise Effects on
Information Rate,” IEEE transactions on nanobioscience, vol. 8, no. 2, pp. 169–180, 2009.
[3] S. Hiyama and Y. Moritani, “Molecular communication: Harnessing biochemical materials to
engineer biomimetic communication systems,” Nano Communication Networks, vol. 1, no. 1,
pp. 20–30, May 2010.
[4] Z. Xie, L. Wroblewska, L. Prochazka, R. Weiss, and Y. Benenson, “Multi-Input RNAi-Based
Logic Circuit for Identification of Specific Cancer Cells,” Science, vol. 333, no. 6047, pp.
1307–1311, Sep. 2011.
[5] B. Atakan, O. Akan, and S. Balasubramaniam, “Body area nanonetworks with molecular
communications in nanomedicine,” Communications Magazine, IEEE, vol. 50, no. 1, pp. 28–
34, 2012.
[6] C. Gardiner, Stochastic methods.

Springer, 2010.

[7] N. van Kampen, Ed., Stochastic Processes in Physics and Chemistry, 3rd ed. North-Holland,
2007.
[8] R. Erban, J. Chapman, and P. Maini, “A practical guide to stochastic simulations of reactiondiffusion processes,” arXiv, vol. q-bio.SC, Apr. 2007.
[9] D. Gillespie, “A rigorous derivation of the chemical master equation,” Physica A: Statistical
Mechanics and its Applications, 1992.
[10] I. Lestas, J. Paulsson, N. Ross, and G. Vinnicombe, “Noise in gene regulatory networks,”
IEEE Transactions on Automatic Control, vol. 53, no. 1, pp. 198–200, 2008.
[11] W. Bialek and S. Setayeshgar, “Physical limits to biochemical signaling,” Proceedings of the
National Academy of Sciences, vol. 102, no. 29, pp. 10 040–10 045, 2005.
[12] R. Endres and N. Wingreen, “Accuracy of direct gradient sensing by cell-surface receptors,”
Progress in biophysics and molecular biology, vol. 100, no. 1-3, pp. 33–39, 2009.
[13] M. Pierobon and I. F. Akyildiz, “A physical end-to-end model for molecular communication
in nanonetworks,” IEEE JOURNAL ON SELECTED AREAS IN COMMUNICATIONS,
vol. 28, no. 4, pp. 602–611, 2010.
[14] M. Mahfuz, D. Makrakis, and H. Mouftah, “On the characterization of binary concentrationencoded molecular communication in nanonetworks,” Nano Communication Networks, vol. 1,
pp. 289–300, 2010.
[15] L. Ljung, System identification - Theory for the user, 2nd ed.

19

Prentice Hall, 1999.

[16] S. Isaacson and C. Peskin, “Incorporating diffusion in complex geometries into stochastic
chemical kinetics simulations,” SIAM Journal on Scientific Computing, vol. 28, no. 1, pp.
47–74, 2007.
[17] D. Gillespie, “Stochastic simulation of chemical kinetics,” Annual review of physical chemistry,
2007.
[18] J. Crank, The mathematics of diffusion.

OUP, 1980.

[19] K. J. Hollenbeck, INVLAP.M: A matlab function for numerical inversion of Laplace transforms by the de Hoog algorithm, 1998, http://www.isva.dtu.dk/staff/karl/invlap.htm.
[20] M. Jackson, Ed., Molecular and cellular biophysics.

CPU, 2006.

[21] S. Basu, Y. Gerchman, C. H. Collins, F. H. Arnold, and R. Weiss, “A synthetic multicellular
system for programmed pattern formation,” Nature, vol. 434, no. 7037, pp. 1130–1134, Apr.
2005.
[22] B. Atakan and O. B. Akan, “Deterministic capacity of information flow in molecular nanonetworks,” Nano Communication Networks, vol. 1, no. 1, pp. 31–42, May 2010.
[23] K. V. Srinivas, R. S. Adve, and A. W. Eckford, “Molecular communication in fluid media:
The additive inverse Gaussian noise channel,” arXiv.org, vol. , Dec. 2010.
[24] K. Burrage, P. Burrage, A. Leier, T. Marquez-Lago, and D. Nicolau Jr, “Stochastic simulation
for spatial modelling of dynamic processes in a living cell,” Design and Analysis of Biomolecular Circuits, p. 43, 2011.
[25] E. Gul, B. Atakan, and O. B. Akan, “NanoNS: A nanoscale network simulator framework for
molecular communications,” Nano Communication Networks, vol. 1, no. 2, pp. 138–156, Sep.
2010.
[26] M. Pierobon and I. F. Akyildiz, “Diffusion-based Noise Analysis for Molecular Communication
in Nanonetworks,” IEEE TRANSACTIONS ON SIGNAL PROCESSING, vol. 59, no. 6, pp.
2532–2547, 2011.
[27] Z. Neufeld and B. Kholodenko, “Positional information generated by spatially distributed
signaling cascades,” PLoS Computational Biology, vol. 5, no. 3, 2009.
[28] M. Pierobon and I. F. Akyildiz, “Noise Analysis in Ligand-Binding Reception for Molecular
Communication in Nanonetworks,” IEEE TRANSACTIONS ON SIGNAL PROCESSING,
vol. 59, no. 9, pp. 4168–4182, 2011.
[29] D. Arifler, “Capacity analysis of a diffusion-based short-range molecular nano-communication
channel,” Computer Networks, vol. 55, pp. 1426–1434, 2010.
[30] C. T. Chou, “Molecular circuits for decoding frequency coded signals in nano-communication
networks,” Nano Communication Networks, 2012, doi:10.1016/j.nancom.2011.11.001.
[31] S. Isaacson, “Relationship between the reaction–diffusion master equation and particle tracking models,” Journal of Physics A: Mathematical and Theoretical, vol. 41, p. 065003, Feb.
2008.
[32] D. Fange, O. Berg, P. Sjöberg, and J. Elf, “Stochastic reaction-diffusion kinetics in the microscopic limit,” Proceedings of the National Academy of Sciences, vol. 107, no. 46, pp. 19 820–
19 825, 2010.
[33] J. Elf, “Spontaneous separation of bi-stable biochemical systems into spatial domains of
opposite phases,” Systems Biology, vol. 1, no. 2, pp. 230–236, 2004.

20

List of Figures
1

2
3
4
5
6
7
8
9
10

11

12
13

An illustration of a molecular communication network with two transmitters and
two receivers. The signalling molecules diffuse in a fluidic medium and may bind
with the receptors at the receviers. . . . . . . . . . . . . . . . . . . . . . . . . . . .
The voxels in the 1-dimensional RDMEX model. . . . . . . . . . . . . . . . . . . .
Block diagram showing the transfer function from the input k1 (t) to the output
c1 (t) of a molecular communication network with one transmitter and one receiver.
The mean number of complexes in the 1-transmitter 1-receiver network for K =
2.5 × 10−3 . (Section 6.2.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The mean number of complexes in the 1-transmitter 1-receiver network for K =
2.5 × 10−4 .(Section 6.2.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The mean number of complexes in the 1-transmitter 1-receiver network with transmitter symbol s1 . (Section 6.2.2) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The standard deviation of the number of complexes in the 1-transmitter 1-receiver
network with transmitter symbol s1 . (Section 6.2.2) . . . . . . . . . . . . . . . . .
The mean number of complexes in the 1-transmitter 1-receiver network for section
6.2.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
This figure shows the effect of the size of a voxel on the mean number of complexes.
The ∆ values used are χ, χ2 , χ3 and χ4 . . . . . . . . . . . . . . . . . . . . . . . . . .
This figure compares the output signal of two different networks. Network 0 consists
of a transmitter and 2 receivers (receivers 1 and 2). Network 1 consists of the
transmitter and receiver 1 of Network 0. The figure shows the output signal of
receiver 1 for these two networks. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
This figure compares the output signal of three different networks. Network 0
consists of a transmitter and 2 receivers (receivers 1 and 2). Network 2 consists of
the transmitter and receiver 2 of Network 0. The third network (label as 1t15r)
consists of network 2 plus another 14 receivers. The figure shows the output signal
of receiver 2 for these networks. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The output signal of receiver 1 in the 2-transmitter 2-receiver network. . . . . . . .
The output signal of receiver 2 in the 2-transmitter 2-receiver network. . . . . . . .

21

22
22
22
23
23
24
24
25
25

26

26
27
27

Transmi*ers	
  

T	
  

T	
  

Signalling	
  
molecule	
  
R	
  

Receptor	
  
R	
  
Receivers	
  

Figure 1: An illustration of a molecular communication network with two transmitters and two
receivers. The signalling molecules diffuse in a fluidic medium and may bind with the receptors
at the receviers.

Δx	
  

…	
  
Voxel	
  1	
  

Voxel	
  2	
  

…	
  
Voxel	
  i-­‐1	
  

Voxel	
  i	
  

Voxel	
  i+1	
  

Voxel	
  N-­‐1	
   Voxel	
  N	
  

Figure 2: The voxels in the 1-dimensional RDMEX model.

k1 (t)

�(vR,1 , t)

φ11 (ω)

ρ1 (ω)

c1 (t)

−
ιωφ0 (ω)
Figure 3: Block diagram showing the transfer function from the input k1 (t) to the output c1 (t) of
a molecular communication network with one transmitter and one receiver.

22

7000

6000

Mean number of complexes

5000

4000

3000

2000
RDMEX−τ
RDMEX−M
RDPDE
RDPDE−X
DE

1000

0
−0.5

0

0.5

1

1.5

2

time

Figure 4: The mean number of complexes in the 1-transmitter 1-receiver network for K = 2.5 ×
10−3 . (Section 6.2.1)

700

600

Mean number of complexes

500
RDMEX−τ
RDMEX−M
RDPDE
DE

400

300

200

100

0
−0.5

0

0.5

1

1.5

2

time

Figure 5: The mean number of complexes in the 1-transmitter 1-receiver network for K = 2.5 ×
10−4 .(Section 6.2.1)

23

mean
700
RDPDE
RDMEX−X
RDMEX−τ

600

mean number of complexes

500

400

300

200

100

0
−0.5

0

0.5

1

1.5

2

time

Figure 6: The mean number of complexes in the 1-transmitter 1-receiver network with transmitter
symbol s1 . (Section 6.2.2)

30
RDPDE
RDMEX−τ

standard deviation of number of complexes

25

20

15

10

5

0
−0.5

0

0.5

1

1.5

2

time

Figure 7: The standard deviation of the number of complexes in the 1-transmitter 1-receiver
network with transmitter symbol s1 . (Section 6.2.2)

24

100
RDMEX−M
RDMEX−τ −MM

90

80

Mean number of complexes

70

60

50

40

30

20

10

0
−0.5

0

0.5

1

1.5

2

time

Figure 8: The mean number of complexes in the 1-transmitter 1-receiver network for section 6.2.3.

1400

1200

Number of complexes

1000

800

600

400
∆=χ
∆ = χ/2
∆ = χ/3
∆ = χ/4

200

0

0

0.2

0.4

0.6

0.8

1
time

1.2

1.4

1.6

1.8

2

Figure 9: This figure shows the effect of the size of a voxel on the mean number of complexes.
The ∆ values used are χ, χ2 , χ3 and χ4 .

25

Mean number of complexes at receiver 1

2500
network 0: Both receivers 1 and 2 are present
network 1: Only receiver 1 is present
2000

1500

1000

500

0

−500

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

time

Figure 10: This figure compares the output signal of two different networks. Network 0 consists
of a transmitter and 2 receivers (receivers 1 and 2). Network 1 consists of the transmitter and
receiver 1 of Network 0. The figure shows the output signal of receiver 1 for these two networks.

Mean number of complexes at receiver 2

1200
network 0: Both receivers 1 and 2 are present
network 2: Only receiver 2 is present
1t15r

1000

800

600

400

200

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

time

Figure 11: This figure compares the output signal of three different networks. Network 0 consists
of a transmitter and 2 receivers (receivers 1 and 2). Network 2 consists of the transmitter and
receiver 2 of Network 0. The third network (label as 1t15r) consists of network 2 plus another 14
receivers. The figure shows the output signal of receiver 2 for these networks.

26

Receiver 1
1600
RDMEX−M
RDPDE
RDMEX−τ

mean number of complexes

1400

1200

1000

800

600

Response
due to s1 of
transmitter 1

Response due to s1
of transmitter 2

400

200

0
0

0.5

1

1.5

2

2.5

3

3.5

4

time

Figure 12: The output signal of receiver 1 in the 2-transmitter 2-receiver network.

Receiver 2
1600
RDMEX−M
RDPDE
RDMEX−τ

mean number of complexes

1400

1200

1000

800

600

Response
due to s1 of
transmitter 1

Response due to s1
of transmitter 2

400

200

0
0

0.5

1

1.5

2

2.5

3

3.5

4

time

Figure 13: The output signal of receiver 2 in the 2-transmitter 2-receiver network.

27

