1

Improving the capacity of molecular
communication using enzymatic reaction cycles

arXiv:1707.05912v2 [cs.IT] 20 Jul 2017

Hamdan Awan, Student Member, IEEE, and Chun Tung Chou Member, IEEE,

Abstract
This paper considers the capacity of a diffusion-based molecular communication link assuming
the receiver uses chemical reactions. The key contribution is we show that enzymatic reaction cycles,
which is a class of chemical reactions commonly found in cells consisting of a forward and a backward
enzymatic reaction, can improve the capacity of the communication link. The technical difficulty in
analysing enzymatic reaction cycles is that their reaction rates are nonlinear. We deal with this by
assuming that the amount of certain chemicals in the enzymatic reaction cycle is large. In order to
simplify the problem further, we use singular perturbation to study a particular operating regime of the
enzymatic reaction cycles. This allows us to derive a closed-form expression of the channel gain. This
expression suggests that we can improve the channel gain by increasing the total amount of substrate in
the enzymatic reaction cycle. By using numerical calculations, we show that the effect of the enzymatic
reaction cycle is to increase the channel gain and to reduce the noise, which results in a better signalto-noise ratio and in turn a higher communication capacity. Furthermore, we show that we can increase
the capacity by increasing the total amount of substrate in the enzymatic reaction cycle.

Index Terms
Molecular communication; Communication capacity; Chemical reaction based receiver; Enzymatic
reaction cycles; Signal-to-noise ratio.

H. Awan and C. T. Chou are with the School of Computer Science and Engineering, The University of New South Wales,
Sydney, NSW 2052, Australia. e-mail: (hawan,ctchou@cse.unsw.edu.au)

March 6, 2018

DRAFT

2

I. I NTRODUCTION
Molecular communication is a promising approach to realize communication between nanoscale devices [1], [2], [3] and especially the internet of bio-nano things [4]. Molecular communication has many useful applications such as detection of harmful pathogens in environment and
detection of tumour cells in human bodies [5]. A key characteristic of molecular communication
is the use of molecules as the information or signal carrier. The transmission of signalling
molecules can be carried out by diffusion [6] or active transport [7]. In this paper, we consider
diffusion-based molecular communication and in particular its information transmission capacity.
The information capacity of a communication link sets a fundamental limit on its communication performance [8]. The information capacity of molecular communication has been studied
in a number of papers. We can divide these papers into two categories depending on whether
they consider reactions at the receiver. The first category of articles, e.g. [9], [10], [11], assumes
a receiver that can count the number of molecules within the receiver volume. These articles do
not consider the reactions at the receiver. The second category of articles considers various types
of chemical reactions at the receiver. The typical reaction considered is ligand-receptor binding,
e.g. in [12], [13], [14]. Our earlier work [15] considers a few different types of reactions at the
receiver, including linearised form of ligand-receptor binding, catalysis and regulated catalysis.
This paper will focus on the capacity of chemical reaction based receivers.
Researchers in biology have found that certain networks of chemical reactions, which are
referred to as motifs, appear more often than the others in gene regulatory networks and protein
reaction networks [16]. These motifs can be considered to be the basic modules (or building
blocks) to realise cell functions. Recently, there is a growing interest in the engineering and
synthetic biology communities to interconnect these modules to create artificial molecular circuits
[17], [18]. To be best of our knowledge, there appears to be few work on studying how the
interconnection of modules will impact on the communication performance. In our earlier work
[19], we study the impact on the capacity of a molecular communication link by a module
consisting of a forward and a backward linearised catalytic reaction. In this paper, we consider

DRAFT

March 6, 2018

3

enzymatic reaction cycles (ERCs) of the form:

K`ZÝ
âÝ
Ýá
Ý KZ ÝÝÑ K ` Z˚

(1)

P ` Z˚ Ý
âÝ
Ýá
Ý PZ˚ ÝÝÑ P ` Z

(2)

where K and P are the enzymes catalysing the conversion between substrates Z and Z˚ , and
KZ (resp. P Z˚ ) is a complex1 formed by the binding of K and Z (resp. P and Z˚ ) molecules.
These ERCs are commonly found in cells and perform covalent modifications to proteins; they
cover reactions such as phorsphorylation, methylation and acetylation [20]. The key contribution
of this paper is to show that it is possible to use ERCs to improve the communication capacity
of diffusion-based molecular communication. We do this by combining ERCs with receiver
molecular circuits that we have studied earlier in [15] and show that the combinations with
ERCs have a higher capacity than those without. A technical difficulty with studying ERCs is
that their reaction rates are nonlinear. We assume the amount of certain chemical species in the
ERC is large and appeal to the tool of singular perturbation and consider a particular operation
regime. This allows us to derive a closed-form expression of the channel gain. This expression
suggests that we can improve the channel gain by increasing the total amount of substrate in
the ERC. By using numerical calculations, we show that the effect of the ERC is to increase
the channel gain and to reduce the noise, which results in a better signal-to-noise ratio and in
turn a higher communication capacity. Furthermore, we show that we can improve the system
capacity by increasing the total amount of substrate.
We organize this paper as follows. Section II presents related work. Section III summarizes
the modelling and analysis framework. In Section IV, we present the ERC and how singular
perturbation can be used to derive the channel gain. Section V presents numerical results on
channel gain, noise and system capacity. Finally, Section VI concludes the paper.

1

The reader may refer to the online entry of Compendium of Chemical Terminology (published by the International Union

of Pure and Applied Chemistry) for a detailed definition of the term complex http://goldbook.iupac.org/html/C/C01203.html.

March 6, 2018

DRAFT

4

II. R ELATED WORK
The interest of research community in molecular communication is on the rise as shown by
recent surveys [1], [2], [5], [3], [21].
The main components of a molecular communication system are transmitter, propagation
medium and receiver. On the transmitter side different modulation schemes have been proposed
in literature such as molecule shift keying, frequency shift keying, pulse position modulation,
concentration shift keying and reaction shift keying [22], [23], [9], [24], [25], [26], [27].
For the propagation medium different models have been used in molecular communication
literature. For example the papers [28], [29] assume that medium is continuous while in this
paper, as well as in our previous works [25], [30], [31], we assume that the medium is divided
into voxels. The use of voxels provides a convenient way to integrate diffusion and reactions
into one mathematical model, see [32] for a tutorial introduction and our earlier work for more
details. An alternative end-to-end model appears in [33], [34] which is based on particle tracking.
For the receiver side, different receiver designs have been proposed in literature for molecular
communication systems, e.g. [35], [36], [37], [38], [39]. Similarly different demodulation techniques for molecular communication systems are presented in [36], [38], [40], [41]. An alternative
way of designing receivers for molecular communication is by using molecular circuits, see [15],
[42], [43] for example. The effect of different receiver molecular circuits on the communication
performance has been studied in the literature. For a molecular communication system the key
performance parameters are noise and capacity. The noise properties of ligand-receptor binding
type of receivers are studied in [44], [33]. The information theoretic analysis on capacity of
molecular communication system is discussed in [45]. The information transmission capacity
of different types of receiver molecular circuits is compared in [15]. The capacity analysis for
molecular communication based on ligand receptor binding has been presented in [12], [10],
[13].
One significant area of research is to study and improve the capacity of molecular communication system. For this we start with understanding the behaviour of connecting modules

DRAFT

March 6, 2018

5

(a) The OM-only configuration.

(b) The ERC-OM configuration.

Fig. 1: Two different link configurations considered in this paper.

in molecular communication system similar to that in synthetic biological systems [46]. This
concept of modularity in cell biology and signal transduction networks is discussed in [47], [48],
[49]. However, there appears to be few work in molecular communication where the receiver
is realised by interconnecting modules. In our previous work [19], we presented an idea to
improve the capacity of the molecular communication link by introducing a module consisting
of simplified (i.e. linearised) ERCs. In this paper, we consider ERCs with nonlinear reaction
rates which are harder to analyse due to nonlinearities.
In the existing literature, there are a number of papers which study how enzyme can be used
to improve the performance of molecular communication networks. For instance, the papers
[35] [50], [51] use enzymes to reduce the amount of inter-symbol interference and the paper
[52] studies a molecular communication link whose receiver uses an enzyme to produce output
molecules from signalling molecules. They key difference between our work and these papers
is that we use an ERC which consists of both forward and backward catalytic reactions while
the other papers use only forward catalytic reactions.

III. M ODELLING AND A NALYSIS F RAMEWORK
We consider a molecular communication link which consists of a transmitter and a receiver.

March 6, 2018

DRAFT

6

Fig. 2: The propagation medium

The transmitter emits signalling molecules which diffuse freely in the propagation medium.
When the signalling molecules reach the receiver, they react with the chemical reactions within
the receiver to produce output molecules. The number of output molecules in the receiver over
time is the output signal of the communication link. Our aim is to maximise the communication
capacity of the link.
In this paper, we will consider two different receiver configurations, see Figure 1, where each
receiver configuration is composed of different modules and each module is composed of a
number chemical reactions. In the configuration in Figure 1a, the receiver consists only of an
output module; we will refer to this as the OM-only link where OM is short for output module.
In the configuration in Figure 1b, the receiver consists of an ERC followed by an output module.
We will refer to this as the ERC-OM link. A key contribution of this paper is to show that the
ERC-OM link has a higher communication capacity than the OM-only system.
This section is organised as follows. We describe how the transmitter and propagation medium
are modelled in Subsections III-A and III-B. This is followed by the modelling of the output
module in Subsection III-C. The models for transmitter, propagation medium and output module
are then combined in Subsection III-D to give a model for the OM-only system. Subsection
III-E explains how the capacity the OM-only system can be computed.

A. Propagation Medium and Transmitter
1) Propagation Medium: We assume the medium of propagation is a three dimensional space
of dimension `X ˆ `Y ˆ `Z where each dimension is an integral multiple of length ∆ i.e. there
DRAFT

March 6, 2018

7

exist positive integers Mx , My and Mz such that `X “ Mx ∆, `Y “ My ∆ and `Z “ Mz ∆. The
medium is divided into Mx ˆ My ˆ Mz cubic voxels where the volume of each voxel is ∆3 .
Figure 2 shows an example with Mx = 5 and My = Mz = 1. We assume that each voxel is
given a unique index. The indices of the voxels are given in the top-right corner of the voxels
in Figure 2.
Diffusion is modelled by molecules moving from one voxel to a neighbouring voxel. The
arrows in Figure 2 show the directions of the movement of the molecules. We assume that the
medium is homogeneous with the diffusion coefficient for the signalling molecule in the medium
is D. Define d “

D
.
∆2

The diffusion of molecules from a voxel to a neighbouring voxel takes

place at a rate of d, i.e. within an infinitesimal time δt, the probability that a molecule diffuses
to a neighbouring voxel is dδt.
The set up in Figure 2 can be used to realise both reflecting and absorbing boundary conditions.
A reflecting boundary condition means the molecules do not leave the medium. An absorbing
boundary condition means once a signalling molecule leaves the medium, it will never return.
An example of absorbing boundary condition is shown in Figure 2 where signalling molecules
can leave a surface of Voxel 3 at an escape rate of e.
We assume the transmitter and the receiver each occupies a voxel. However, it is straightforward to generalise to the case where a transmitter or a receiver occupies multiple voxels. The
transmitter and receiver are assumed to be located, respectively, at the voxels with indices T
and R. For example, in Figure 2, Voxel 2 (dark grey) contains the transmitter and Voxel 4 (light
grey) contains the receiver. Hence T “ 2 and R “ 4 for this example.
2) Transmitter: The transmitter emits signalling molecules (denoted by L) at a rate of uptq
at time t. This means that in the time interval [t, t ` δt), the transmitter emits uptqδt signalling
molecules into the transmitter voxel. We consider uptq as the input signal of the molecular
communication link. We further assume that uptq “ c ` wptq , where c is a positive constant
and wptq is a zero-mean stationary Gaussian random process. Later on, we will use the spectral
property of the signal wptq to maximise the mutual information through water filling [8].

March 6, 2018

DRAFT

8

B. Diffusion-only Subsystem
In this section we consider the molecular communication link in Figure 2 assuming that
the receiver reaction mechanism has been removed. This means the system consists of the
transmitter which injects signalling molecules into the transmitter voxel and the diffusion of
signalling molecules in the propagation medium. We will refer to this as the diffusion-only
subsystem and it is the same for both the OM-only and ERC-OM links. Our aim is to illustrate
how this subsystem can be modelled. Although the example is based on the network in Figure
2, generalisation to other networks is straightforward.
The state of the diffusion-only subsystem is the number of signalling molecules in each voxel.
Let nL,i ptq denote the number of signalling molecules in Voxel i, then the state of the diffusiononly subsystem in Figure 2 is given by:
„
T
nL ptq “ nL,1 ptq nL,2 ptq nL,3 ptq nL,4 ptq nL,5 ptq
where the superscript

T

(3)

denotes matrix transpose.

Our basic modelling framework is to consider each movement of a signalling molecule from
a voxel to another as an event. For example let us consider the movement of a molecule from
Voxel 2 to Voxel 3. This movement occurs at a rate of dnL,2 and after this movement has taken
place, nL,2 will be decreased by 1 and nL,3 increased by 1. We can describe the change in the
number of molecules in the voxels due to this movement by a jump vector. For this example,
the jump vector qd “ r0, ´1, 1, 0, 0sT and the state changes from nL ptq to nL ptq ` qd due to this
movement. The rate of this movement can be represented by a jump rate Wd pnL ptqq “ dnL,2 .
We can use exactly the same concept to model the boundary condition.
In general, we write the jump vectors qd,j and jump rates Wd,j pnL ptqq where the subscript d is
used to indicate an event due to diffusion and j is used to index the events. We use Jd to denote
the total number of diffusion events. For Figure 2, Jd “ 9 consisting of 8 inter-voxel diffusion
events and 1 escape event. With the 9 jump vectors and jump rates, we can use stochastic
differential equation (SDE) [53] to model the dynamics of the diffusion-only subsystem. This
modelling framework is based on the two facts: first, diffusion between voxels can be modelled

DRAFT

March 6, 2018

9

as abstract chemical reactions [32]; second, chemical reactions can be modelled by SDE [54].
The diffusion of signalling molecules in the configuration in Figure 2 can be modelled by the
following SDE:
n9 L ptq “

Jd
ÿ

qd,j Wd,j pnL ptqq `

Jd
ÿ

b
qd,j Wd,j pnL ptqqγj ` 1T uptq

(4)

j“1

j“1

where γj is continuous Gaussian white noise with unit variance and 1T is a unit vector with
1 at the T -th element with the subscript T being the index of the transmitter. There are three
terms on the right-hand side of Eq. (4) and we will discuss them one by one. The first term
describes the deterministic dynamics. Since all the jump rates of all the diffusion events are
linear, this term can be written as a product of a matrix H and the state vector nptq. The matrix
H is defined by the following equality:
HnL ptq “

Jd
ÿ

qd,j Wd,j pnL ptqq

(5)

j“1

For Figure 2, the H matrix is given by:
»
0
0
0
— ´d d
—
— d ´2d
d
0
0
—
—
—
H“— 0
d ´2d ´ e d
0
—
—
— 0
0
d
´2d d
—
–
0
0
0
d ´d

fi
ffi
ffi
ffi
ffi
ffi
ffi
ffi
ffi
ffi
ffi
ffi
fl

(6)

The second term of Eq. (4) describes the stochastic dynamics. An intuitive way to understand
the first two terms in Eq. (4) is as follows. Over a finite time interval ∆t, the number of times
that the j-th type of jump occurs can be approximated by a Poisson random variable with mean
Wd,j pnL ptqq ∆t. If Wd,j pnL ptqq ∆t is large, then we know from probability theory that a Poisson
random variable with mean Wd,j pnL ptqq ∆t can be approximated by a Gaussian variable with
both mean and variance given by Wd,j pnL ptqq ∆t. We can therefore approximate the number
a
of times that the j-th type of jumps occurs by Wd,j pnL ptqq ∆t ` Wd,j pnL ptqq ∆tγj where the
first term is the mean number of jumps and the second term is the deviation from the mean;

March 6, 2018

DRAFT

10

these two terms give rise to the first two terms in Eq. (4). For a more detailed explanation, the
reader can refer to [54].
The third term models the transmitter. Since the transmitter emits uptqδt molecules at time t,
we add this number of molecules to voxel T (the index of the transmitter voxel) at time t.
It is important to point out that the elements in nL ptq, which have the interpretation of the
number of molecules, is strictly speaking a discrete random variable. The SDE is an approximation which holds when the order of the number of molecules is Op100q [55]. However, as
far as the first and second order moments are concerned, the SDE (4) gives the same result as
a master equation formulation that assumes the number molecules is discrete [56].

C. Output Module

The output signal of both OM-only and ERC-OM links is the count of the number of output
molecules over time. The aim of the output module is to produce the output molecules, and
hence the output signal. In our earlier work [15], we studied the impact of a number of molecular
circuits on the capacity of the communication link. We will use two of the molecular circuits in
[15] as output modules in this paper. Our aim is to show that the improvement of capacity by
ERCs is general and applies to multiple types of output modules.
Each output module has two chemical species B and X where X is the output molecule.
The identity of B depends on whether we are using OM-only or ERC-OM links. For the case
of OM-only links, the chemical species B is the signalling molecules in the receiver voxel. We
will explain what the identity of B is for the ERC-OM case later on. We assume that the output
molecules X do not diffuse and stay inside the receiver voxel.
Following the terminology of [15], the output modules are referred to as reversible conversion
(RC) and catalysis plus regulation (CARTEG).

DRAFT

March 6, 2018

11

1) The RC output module: The RC output module consists of two linear chemical reactions
and they are described by their chemical reaction equation, jump vector and jump rate.
T

„
BÑX

´1 1

(7)

, k´ nX

(8)

T

„
XÑB

, k` nB

1 ´1

where nB and nX are respectively the number of B and X molecules, and k` and k´ are reaction
rate constants. In reaction (7), the jump vector r´1 1s indicates that a B molecule is converted
to an X molecule in this reaction while the reaction (8) is the reverse of this, hence the name
reversible conversion.
2) The CATREG output module: The CATREG output module consists of three chemical
reactions mentioned below with their respective respective jump functions and jump rates.
T

„
B ÑB`X

0 1
0 ´1

, k´ nX

(10)

, k0 nX

(11)

T

„
B ÑX φ

(9)

T

„
XÑφ

, k` nB

´1 0

In reaction (9) the chemical B acts as a catalyst to produce the output molecule X at a rate of
k` nB . Note that the number of B molecules remains unchanged before and after the reaction.
This is indicated by the jump vector which shows that every time this reaction occurs, the number
of B remains unchanged and the number of output molecules is increased by 1. The reaction
(10) is a degradation reaction where X molecules are degraded to a species φ that we are not
interested to keep track of. Lastly in reaction (11) the degradation of B molecules in the receiver
voxel is driven by the presence of the output molecules X at a rate of k0 nX and we can view
this as using X to regulate the amount of B.
3) Modelling the output module: This section presents an SDE model of the output module.
This model will be used for the modelling of the OM-only and ERC-OM links later on. The

March 6, 2018

DRAFT

12

TABLE I: R matrix for different output modules
R Matrix

Output module
»

fi
´k
k
`
´
—
ffi
–
fl
k`
´k´
»
fi
0
´k
0 ffi
—
–
fl
k` ´k´

RC

CATREG

output module consists of two chemicals. The state vector of the output module is:
T
„
ñR ptq “ nB ptq nX ptq

(12)

We will use qr,j and Wr,j pñR ptq to denote, respectively, the jump vector and jump rates in
the output module. Note that the subscript r indicates the jump vector and jump rate belong to
the output module, which is situated inside the receiver. Let Jr be the number of reactions in
the output module; the values of Jr for RC and CATREG are, respectively, 2 and 3. We will
index j from Jd ` 1 to Jd ` Jr so that we can keep the expression simple when we combine
the diffusion-only subsystem and the output module later. By using the jump vectors and jump
rates of the output module, the SDE that describes the evolution of the number of molecules in
the output module is:
n9̃ R ptq “ RñR ptq `

Jdÿ
`Jr

b
qr,j Wr,j pñR ptqqγj

(13)

j“Jd `1

where γj is continuous Gaussian white noise of unit variance and the R matrix is defined by the
ř d `Jr
relation: RñR ptq “ Jj“J
qr,j Wr,j pñR ptqq. The R matrices for the RC and CATREG output
d `1
modules are shown in Table I. The derivation of Eq. 13 follows from the fact that the dynamics
of chemical reactions can be approximately modelled by SDE [54]. The modelling technique
being used is similar to that in deriving Eq. 4.

D. The OM-only link
In this section, we will combine the SDE models for the diffusion-only subsystem in Section
III-B and the output module in Section III-C to form the complete model for the OM-only

DRAFT

March 6, 2018

13

link. It is important to note that in the OM-only link, the B molecule in the OM is the
signalling molecules in the receiver voxel. Therefore, the interconnection between the diffusiononly subsystem and the output module is the number of signalling molecules in the receiver
voxel, which is common to both of them.
We will use the example in Fig. 2 to explain how the diffusion-only subsystem and the output
module can be combined together. The dynamics of the diffusion-only subsystem for Fig. 2 is
given by Eq. 4. For this example, the receiver voxel has the index R “ 4, so the evolution of the
number of signalling molecules in the receiver voxel nL,R ptq is given by the R-th (i.e. fourth)
row of Eq. 4, which is:
n9 L,R ptq “ dnL,3 ptq ´ 2dnL,R ptq ` dnL,5 ptq `

Jd
ÿ

b
rqd,j sR Wd,j pnL ptqqγj
j“1
looooooooooooooomooooooooooooooon

(14)

ξd ptq

where rqd,j sR denote the R-th element of the vector qd,j .
In the output module, the signalling molecules in the receiver voxel are the B molecules in
Section III-C. The dynamics of the number of signalling molecules in the receiver voxel due to
the reactions in the output module is given by the first element of Eq. 13, which is:
Jdÿ
`Jr

b
n9 L,R ptq “ R11 nL,R ptq ` R12 nX ptq `
rqr,j s1 Wr,j pñR ptqqγj
j“J
`1
d
loooooooooooooooomoooooooooooooooon

(15)

ξr ptq

where rqr,j s1 denotes the first element of the vector qr,j .
For the OM-only link, the dynamics of nL,R ptq is obtained by combining Eq. (14) and (15)
as follows:
n9 L,R ptq “dnL,3 ptq ´ 2dnL,R ptq ` dnL,5 ptq ` R11 nL,R ptq ` R12 nX ptq ` ξtotal ptq

(16)

where ξtotal ptq “ ξd ptq ` ξr ptq.
We are now ready to describe the complete model for the OM-only link. Let nptq be the state
of the OM-only link and it is given by:
T

„
nptq “

March 6, 2018

nL ptq

T

nX ptq

(17)

DRAFT

14

We will also need to modify the jump vectors from the diffusion-only subsystem and the
output module to obtain the jump vectors for the complete model; this will be explained in a
moment. We use qj and Wj pnptqq to denote the jump vectors and jump rates of the combined
model. The SDE for the complete system is:
J
b
ÿ
nptq
9
“ Anptq `
qj Wj pnptqqγj ` 1T uptq

(18)

i“1

where J “ Jd ` Jr , and the matrix A is defined by Anptq “

řJ

i“1 qj Wj pnptqq.

The matrix A

has the block structure:
fi

»
1TR 1R R11

— H`
A “–
R21 1TR

1R R12 ffi
fl
R22

(19)

where H comes from the diffusion only subsystem (Note: an example of H for Figure 2 is in
(6).) and R11 , R12 etc come from the output module. The vector 1R is a unit vector with an 1 at
the R-th position; in particular, note that 1TR nL ptq “ nL,R ptq which is the number of signalling
molecules in the receiver voxel. Note that, the coupling between the diffusion-only subsystem
and the output module, as exemplified by (16), takes place at the R-th row of A.
We now explain how the jump vectors for the OM-only link are formed. Let md denote the
dimension of the vector nL ptq. Note that md is in fact the number of voxels. The dimension of
the jump vectors qj in the complete system is md ` 1. Given jump vector qd,j (j “ 1, ..., Jd )
from the diffusion only subsystem with dimension md , we append a zero to qd,j to obtain qj .
The jump vectors qr,j (j “ Jd ` 1, ..., Jd ` Jr ) from the output module has dimension mr ` 1.
To obtain qj from qr,j , we do the following: (1) take the first element of qr,j and put it in the
R-th element of qj ; (2) take the last element of qr,j and put it in the last element of qj . Note
that jump rates are unchanged when combining the subsystems.

E. Capacity of the OM-only link
The input and output signals for the OM-only link are, respectively, the production rate uptq
of the signalling molecules in the transmitter voxel and the number of output molecules nX ptq
DRAFT

March 6, 2018

15

in the receiver voxel. In this section, we will derive an expression for the mutual information
between the input uptq and output nX ptq.
We begin by stating a result in [57] which states that, for two Gaussian distribution random
processes aptq and bptq, their mutual information Ipa, bq is given by:
ˆ
˙
ż
´1 8
|Φab pωq|2
Ipa, bq “
log 1 ´
dω
4π ´8
Φaa pωqΦbb pωq

(20)

where Φaa pωq (resp. Φbb pωq) is the power spectral density of aptq (bptq), and Φab pωq is the cross
spectral density of aptq and bptq.
In order to apply the above results to the communication link given in Eq. (18), we need a
result from [56] on the power spectral density of systems consisting only of chemical reactions
with linear reaction rates. Following from [56], if all the jump rates Wj pnptqq in (18) are linear
in nptq, then the power spectral density of nptq is given by the following SDE:
J
b
ÿ
9
nptq
“ Anptq `
qj Wj pxnp8qyqγj ` 1T uptq

(21)

i“1

where xnptqyq denotes the mean of nptq and is the solution to the following ordinary differential
equation:
9 “ Axnptqy ` 1T c
xnptqy

(22)

where c, which was defined before, is the mean of input uptq.
As a result, the dynamics of the OM-only link in Eq. (21) are described by a set of linear SDE
with uptq as the input and nX ptq (which is the last element of the state vector nptq) as the output.
The input uptq has the form uptq “ c ` wptq where c is a constant to set the operating point of
the system and wptq is a zero-mean Gaussian random process. The noise in the output nX ptq is
caused by the Gaussian white noise γj ’s in Eq. (21). Therefore, Eq. (21) models a continuoustime linear time-invariant (LTI) stochastic system subject to Gaussian input and Gaussian noise.

The power spectral density ΦX pωq of the signal nX ptq can be obtained from standard results
on the output response of a LTI system to a stationary input [58] and is given by:
ΦX pωq “ |Ψpωq|2 Φu pωq ` Φη pωq

March 6, 2018

(23)

DRAFT

16

where Φu pωq is the power spectral density of uptq and |Ψpωq|2 is the channel gain with Ψpωq “
Ψpsq|s“iω defined by:
xNX psqy “ 1X xN psqy “ looooooooomooooooooon
1X psI ´ Aq´1 1T U psq

(24)

Ψpsq

Note that Eq. (24) can be obtained from Eq. (21) after taking the mean and applying Laplace
transform. The term Φη pωq denotes the stationary noise spectrum and is given by:
Φη pωq “

Jdÿ
`Jr

|1X piωI ´ Aq´1 qj |2 Wj pxnp8qyq

(25)

j“1

where nptq denotes the state of the complete system in (17) and xnp8qy is the mean state of
system at time 8 due to constant input c.
Similarly, by using standard results on the LTI system, the cross spectral density Ψxu pωq has
the following property:
|Ψxu pωq|2 “ |Ψpωq|2 Φu pωq2

(26)

By substituting Eq. (23) and Eq. (26) in the mutual information expression in Eq. (20), we
arrive at the mutual information IpnX , uq between uptq and nX ptq is:
ˆ
˙
ż
|Ψpωq|2
1
log 1 `
Φu pωq dω
IpnX , uq “
2
Φη pωq

(27)

The capacity of the link can be determined by applying the water-filling solution to (27) subject
to power constraint on input uptq [8]. The capacity of the link depends on the channel gain,
noise power spectral density and the input power spectral density.
IV. T HE ERC-OM LINK
The aim of this section is to use analytical methods to study the property of the ERC-OM
link. We learn from Section III-E that the channel gain can be used to influence the capacity of
the communication link. However, channel gain is only defined for LTI systems but the reaction
rates of the ERC are nonlinear functions of the concentration of the reactants. In order to study
the channel gain of the ERC-OM link, we use singular perturbation and assume the amount of

DRAFT

March 6, 2018

17

certain chemical species is large to obtain a linear approximation of the input-output response
of a ERC-OM link. We derive closed-form expressions for the approximate channel gain for
the ERC-OM link with RC or CATREG as the output module. An insight from these closedform expressions is that we can increase the channel gain by increasing the amount of certain
chemical species in the ERC module. We will show using numerical studies in Section V that
these chemical species can be used to increase the capacity of the ERC-OM link.

A. Enzymatic Reaction Cycles (ERC)
The ERC consists of two sets of chemical reactions which facilitate the conversion of a pair
of chemical species. We will refer to the pair of chemicals as Z and Z˚ . We will refer to the
conversion of Z to Z˚ as forward and the conversion of Z˚ to Z as backward. We have already
shown the chemical equations for ERC in Section I but we will rewrite them in a slightly different
form so that we can simplify the notation of some equations later on. The chemical equations
for the ERC are:
β1

k

1
K`ZÝ
âÝ
Ýá
Ý C1 ÝÝÑ K ` Z˚

β2

α

1
k2
á
P ` Z˚ Ý
âÝ
Ý C2 ÝÝÑ P ` Z
αÝ
2

(28)
(29)

where the enzymes K and P catalyze, respectively, the forward and reverse conversions. The
chemical species C1 and C2 are intermediate complexes form by the enzymes and the substrates.
The symbols k1 , k2 , β1 , β2 , α1 and α2 are reaction rate constants. The enzyme P can exist on its
own or as part of the complex C2 . We use PT to denote the total amount of P and we assume
PT is a constant.

B. ERC with the RC output module
The aim of this section is to derive the channel gain for the ERC-OM link where the output
module is RC. In this setup, we assume: (1) The chemical K in the ERC in reaction (28) is
the signalling molecules in the receiver voxel; and, (2) The chemical B in the RC reaction is

March 6, 2018

DRAFT

18

Z˚ . We also assume that all the chemicals of the ERC and output module, with the exception
of signalling molecules, stay within the receiver voxel.
The output module of this link consists of the following reactions:
k`

Z˚ Ý
âÝ
Ýá
ÝX

(30)

k´

where X, as before, represents the output molecule. Note that these reactions are identical to
those in Section III-C1 except that we have replaced B by Z˚ . We will use ZT to denote the
total amount of Z, Z˚ , C1 , C2 and X; and we will assume that ZT is a constant.
Since our aim is to derive the channel gain, which relates the mean input and output signals,
we will therefore use only the deterministic part of the SDEs. In order to simplify the notation,
we will drop the angle brackets x y and all the chemical species counts are assumed to be their
mean value. In addition, for the chemical species in the ERC, we will simply use their chemical
names as their count.
The dynamics of the chemical species in the ERC and output modules are:
C9 1 ptq “ ´pβ2 ` k1 qC1 ptq ` β1 nL,R ptqpZT ´ Z˚ ptq ´ C1 ptq ´ C2 ptq ´ nX ptqq

(31)

C9 2 ptq “ ´pα2 ` k2 qC2 ptq ` α1 Z˚ ptqpPT ´ C2 ptqq

(32)

Z9 ˚ ptq “ k1 C1 ptq ` α2 C2 ptq ´ α1 Z˚ ptqpPT ´ C2 ptqq ´ k` Z˚ ptq ` k´ nX ptq

(33)

n9 X ptq “ k` Z˚ ptq ´ k´ nX ptq

(34)

Note that Eq. (31), (32) and (33) are nonlinear. We first assume that ZT " (Z˚ ptq - C1 ptq - C2 ptq
- nX ptqq, so we can simplify Eq. (31) to:
C9 1 ptq “ ´pβ2 ` k1 qC1 ptq ` β1 nL,R ptqZT

(35)

We next assume that PT " C2 ptq, so we can replace the term pPT ´ C2 ptqq in Eq. (32) and (33)
by PT . As a result of these two assumptions, we have got rid of the nonlinearities in the model.
Even after getting rid of nonlinearities, the resulting equations are still difficult to analyse.
The next step is to use time scale separation between the rate of diffusion and chemical reactions
to simplify the equations. We define G1 “
DRAFT

β1 ZT
d

and G2 “

k´
.
k1

We assume 1 =

1
G1

and 2 =

1
G2

March 6, 2018

19

are small; we will discuss these assumptions further at the end of this sub-section. This allows
us to define the slow variable W ptq:
W ptq “ Z˚ ptq ` nX ptq

(36)

and apply singular perturbation to separate the fast and slow dynamics [47]. The essential idea
behind time scale separation is that the fast dynamics reaches equilibrium before the slow
dynamics. This means we can replace a differential equation describing fast dynamics by an
algebraic equation and we only have to keep track of the slow dynamics. Note that in singular
perturbation, the approximation error in replacing a differential equation describing the fast
dynamics by an algebraic equation is of the order Opmaxp1 , 2 qq if 1 and 2 are sufficiently
small, see [59].
After applying singular perturbation and going through many steps we obtain one algebraic
equation and one ordinary differential equation (ODE), as follows:
nX ptq « rZ˚ ptq

(37)

Z9 ˚ ptq “ rG1 k11 C1 ptq ` G1 a2 C2 ptq ´ G1 a1 Z˚ ptqsr
where k11 =

k1
G1

, a1 =

α1 P T
G1

, a2 =

α2
G1

and r “

k`
.
k´

1
s
p1 ` rq

(38)

By using Laplace transform, we can combine

these two equations as:
NX psq “ rr

pk1 C1 psq ` α2 C2 psq{p1 ` rq
s
s ` α1 PT {p1 ` rq

(39)

By using the ODEs that describe the diffusion of the signalling molecules, we can derive the
expressions of C1 psq and C2 psq in terms of the input signal U psq.
The expressions for C1 psq and C2 psq can be found in Eqs. (85) and (87). Furthermore the
complete derivation for NX psq is presented in Appendix A.
The final result is:
k1 β1 ZT
NX psq “ Qpsq
U psq
ps ` β2 ` k1 q
loooooooooomoooooooooon

(40)

Ψ̃psq

March 6, 2018

DRAFT

20

where
Qpsq “ rr

p1 ` α1 PT {rps ` α2 ` k2 qps ` α1 PT q ´ pα1 α2 PT qsq{p1 ` rq T
s1R psI ´ Hq´1 1T (41)
s ` α1 PT {p1 ` rq

Note that the the channel gain is given by |Ψ̃psq|2 where Ψ̃psq in given in Eq. (40). It can be
seen from the expression of Ψ̃psq that channel gain is an increasing function of ZT .
Discussions: (1) In order to derive the results in this section, we assume that ZT and PT are
large, and limit ourselves to the operating regime where 1 “

d
β1 ZT

and 2 “

k1
k´

are small. If

ZT is large, then 1 can be made small. Therefore, the requirements we need are large ZT , PT
and small k1 relative to k´ . If the output receiver is to be constructed by synthetically, then the
quantities ZT and YT can be chosen by the designer. The condition on k1 can be met by choosing
an appropriate ERC. Similar methods have been used in [60] to realise a synthetic biomolecular
circuit with certain properties. The paper [60] assumes that the quantities of certain chemical
species are large and requires certain reaction rate constant to be large. (2) In this paper, we
have limited ourselves to a particular operating regimes. It may be possible to arrive at the same
result by consider other operating regimes. We will leave this for further work.

C. ERC with the CATREG output module
The aim of this section is to derive the channel gain for the ERC-OM link where the output
module is CATREG. In this setup, we assume: (1) The chemical K in the ERC in Reaction
(28) is the signalling molecules in the receiver voxel; and, (2) The chemical B in the CATREG
reaction is Z˚ . The latter assumption means that the output module of this link consists of the
following reactions:
k`

Z˚ ÝÝÑ Z˚ ` X
k´

(42)

X ÝÝÑ φ

(43)

k

(44)

Z˚ ÝÝ0Ñ Xφ

Note that these reactions are identical to those in Section III-C2 except that we have replaced
B by Z˚ . We will use ZT to denote the total amount of Z, Z˚ , C1 and C2 . Note that due to the
DRAFT

March 6, 2018

21

degradation of Z˚ by X in reaction (44), the value of ZT decreases with time. In this paper, we
will assume that the degradation rate k0 is small so ZT can be considered to be a constant over
a small finite time interval.
The dynamics of the chemical species in the ERC and output modules are:
C9 1 ptq “ ´pβ2 ` k1 qC1 ptq ` β1 nL,R ptqpZT ´ Z˚ ptq ´ C1 ptq ´ C2 ptqq

(45)

C9 2 ptq “ ´pα2 ` k2 qC2 ptq ` α1 Z˚ ptqpPT ´ C2 ptqq

(46)

Z9 ˚ ptq “ k1 C1 ptq ` α2 C2 ptq ´ α1 Z˚ ptqpPT ´ C2 ptqq ´ k0 nX ptq

(47)

n9 X ptq “ k` Z˚ ptq ´ k´ nX ptq

(48)

We will assume ZT " (Z˚ ptq - C1 ptq - C2 ptq) and PT " C2 ptq to remove the nonlinearities in the
model. We further assume that 1 =

1
G1

and 2 =

1
G2

are small, and this allows us to define a slow

variable W ptq in Eq. (36). After applying singular perturbation and going through numerous
steps, we obtain an algebraic equation and an ODE:
nX ptq « rZ˚ ptq

(49)

Z9 ˚ ptq “ rG1 k11 C1 ptq ` G1 a2 C2 ptq ´ G1 a1 Z˚ ptq ´ rpk0 ` k´ qsr1{p1 ` rqs

(50)

where k11 =

k1
G1

, a1 =

α1 P T
G1

, a2 =

α2
,
G1

r “

k`
,
k´

and the expressions of C1 psq and C2 psq are

given in Eq. (85) and (92) in Appendix A. By combining these two equations with those that
describe the diffusion dynamics, we can show that NX psq again has the form as Eq. (40) but
the expression of Q(s) is:

Qpsq “ rr

p1 ` α1 PT {rps ` α2 ` k2 qps ` α1 PT ` k` ´ rk0 q ´ pα1 α2 PT qsq{p1 ` rq T
s1R psI ´ Hq´1 1T
s ` α1 PT {p1 ` rq
(51)

This shows that the channel gain is again an increasing function of ZT . We will show in Section
V that we can use ZT to improve the capacity of the ERC-OM link when CATREG is the output
module.

March 6, 2018

DRAFT

22

D. Computation of capacity
This section explains how the capacity of an ERC-OM link can be computed. The dynamics
of the link are described by the SDE:
9
nptq
“

J
ÿ
i“1

J
b
ÿ
qj Wj pnptqq `
qj Wj pxnptqyqγj ` 1T uptq

(52)

i“1

where qj ’s are jump vector, Wj pnptqq’s are the jump rates, and J is the total number of diffusion
and reaction events. Note that although we have made use of the same notation as before, we
assume qj , Wj pnptqq and J have been adapted for the ERC-OM link. Since the ERC contains
nonlinear reactions, it means that some of the jump rates Wj pnptqq are nonlinear functions of the
state vector nptq. This means we cannot write the first term on the right-hand side of Eq. (52)
as the product of a matrix and the state vector nptq.
We argue in Sections IV-B and IV-C that if we assume ZT " (Z˚ ptq - C1 ptq - C2 ptq - nX ptqq
and PT " C2 ptq, then we can approximate the nonlinear dynamics with linear ones. This is the
same as replacing nonlinear jump rates Wj pnptqq by linear ones. Once this replacement is done,
we can proceed as before to compute the mutual information and capacity.
V. N UMERICAL EVALUATION
This section presents numerical examples to illustrate the use of the ERC to improve the
communication performance of a molecular communication link.
We assume an array of 5 ˆ 2 ˆ 2 voxels with a voxel size of ( 13 µm)3 (i.e. ∆ “

1
3

µm). We

have also used larger arrays of voxels and we obtain similar results. The transmitter and receiver
are located at voxels (2,1,1) and (4,2,2) respectively. We assume the diffusion coefficient D of
the medium is 1 µm2 s´1 . The deterministic emission rate c is chosen to be 10 molecules per
second. We assume an absorbing boundary with an escape rate e equal to

d
.
10

For the CATREG

receiver circuit in the output module, we fix k0 = 0.01 whereas k` and k´ (both or any one
of them) can be varied to obtain different values of association constant r “

k`
.
k´

For the ERC

module, we choose β1 “ β2 “ 1, k1 “ 0.05, α1 “ α2 “ 1 and k2 “ 0.5. Furthermore, we choose
ZT “ 500 and PT “ 200.
DRAFT

March 6, 2018

23

Mean Number of Output Molecules vs Time

100

90

Linear Approximation
SSA Simulation

Mean Number of Output Molecules

80

70

60

50

40

30

20

10

0
0

50

100

150

200

250

300

350

400

450

500

Time

Fig. 3: This figure shows that SSA simulation of ERC-OM with nonlinear reaction rates gives
almost the same output as the linear approximation.

The ERC-OM link is nonlinear and we derive linear approximation for the link in Sections
IV-B and IV-C. We first verify that the linear approximation is sufficiently accurate. We do this
by simulating the ERC-OM with nonlinear reaction rates and compare the results against those
given by the linear approximation. For the simulation of ERC-OM with nonlinear reaction rates,
we use Stochastic Simulation Algorithm (SSA) [61] which is a method to simulate systems with
both chemical reactions and diffusion. The SSA algorithm simulates the Markov chain which
describes the evolution of the number of molecules in a system due to diffusion and reactions
[32]. The reason why we choose SSA is because it can be applied to reactions with non-linear
reaction rates. We use SSA to simulate an ERC-OM link with nonlinear reaction rates and
compute the mean number of output molecules. We compare this against the mean number of
output molecules given by the linear approximation. Figure 3 compares the mean number of
output molecules from SSA simulation and the linear approximation. It can be seen from the
figure that the linear approximation is accurate.
Our next step is to show the improvement made by the ERC-OM link over that of OM-only.

March 6, 2018

DRAFT

24

Channel Gain For OM-only Link with CATREG

10-2

10-4

Channel Gain

10-6

10-8

10-10

10-12

r=4
r=0.5
r=0.1
r=0.02

10-14
10-4

10-3

10-2

10-1

100

Frequency (Hz)

Fig. 4: The channel gain for the OM-only link with CATREG as the output module.

Noise PSD for the OM-Only Link with CATREG

10 4

Noise PSD

10 3

10 2

10 1

10 0

10 -1
10 -4

r=4
r=0.5
r=0.1
r=0.02

10 -3

10 -2

10 -1

10 0

Frequency (Hz)

Fig. 5: The noise power spectral density for the OM-only link with CATREG as the output
module.

We have performed numerical experiments on using both RC and CATREG as the output module.
We will present the results for CATREG only because the results for RC are similar. For these
numerical experiments, we vary the association constant r “

k`
k´

of the output module.

We first compare the channel gain and noise power spectral density of the OM-only and ERCOM links. Figures 4 and 5 show, respectively, the channel gain and noise power spectral density
of the OM-only link for different values of r. The corresponding results for ERC-OM links are
shown in Figures 6 and 7. These figures are best viewed in colour because we have used the
same coloured line for the same value of r. By comparing Figures 4 and 6, we can see that
the effect of the ERC is to increase the channel gain. By comparing Figures 5 and 7, we can

DRAFT

March 6, 2018

25

Channel Gain for ERC-OM Link with CATREG

10 0

10 -2

Channel Gain

10 -4

10 -6

10 -8

10 -10
r=4
r=0.5
r=0.1
r=0.02

10 -12

10 -14
10 -4

10 -3

10 -2

10 -1

10 0

Frequency (Hz)

Fig. 6: The channel gain for the ERC-OM link with CATREG as the output module.

Noise PSD for the ERC-OM Link with CATREG

10 4

Noise PSD

10 3

10 2

10 1

10

0

10 -1
10 -4

r=4
r=0.5
r=0.1
r=0.02

10 -3

10 -2

10 -1

10 0

Frequency (Hz)

Fig. 7: The noise power spectral density for the ERC-OM link with CATREG as the output
module.

see that the effect of the ERC is to decrease the noise power spectral density. This means the
overall effect of ERC is to increase the signal-to-noise ratio of the link, and therefore capacity.
We now compare the capacity of the OM-only and ERC-OM links. In these numerical
experiments, we vary k` from 0 to 10 and keep other reaction rate constants unchanged. We
use two different pairs of ZT and PT . We compute the capacity of the OM-only and ERC-OM
links. Figure 8 shows the results for ZT “ 500 and PT “ 200 and Figure 9 is for ZT “ 2000
and PT “ 500. We see in both Figures 8 and 9 that the ERC can improve the link capacity for
all values of r. The figures also show that the capacity increases with r initially but plateaus off
later on.

March 6, 2018

DRAFT

26

Capacity of the ERC-OM and OM-only Links with CATREG

1

0.9

For ERC-OM Link: Z

T

= 500 , P T = 200

For OM-Only Link
0.8

0.7

Capacity

0.6

0.5

0.4

0.3

0.2

0.1

0
0

.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

6.5

7

7.5

8

8.5

9

9.5

10

Parameter r = k + / k - for CATREG with Constant k - and k 0

Fig. 8: Capacity of the ERC-OM and OM-only links for ZT “ 500 and PT “ 200.
We next study how the value of ZT impacts on the capacity. The derivation in Sections IV-B
and IV-C shows that ZT can increase the channel gain. We vary ZT from 500 to 5000. We use
two different values of k` and keep all the other reaction rate constants unchanged. Figure 10
shows that we can increase the capacity by increasing ZT and this increase is observed for both
values of k` being used. This is a welcoming news because it gives us a method to increase
the capacity of the link by adjusting the amount of chemical being used. If one thinks about the
degrees of freedom that one can use to influence the chemical reactions in a receiver, one can
change the amount of chemical species or the reaction rate constants. Unfortunately, it may not
always be possible to change reaction rate constants because they depend on temperature and
pressure of the operating environment which may be out of our control. However, the amount of
chemical species is a parameter that can readily be controlled. Therefore, our research provides
a practical method to tune the capacity of a communication link.
VI. C ONCLUSIONS
This paper considers the capacity of a communication link in diffusion-based molecular
communication. We consider the case that the receiver uses chemical reactions. Our key contri-

DRAFT

March 6, 2018

27

Capacity of the ERC-OM and OM-only Links with CATREG

1

0.9

0.8

0.7

Capacity

0.6

0.5

0.4

For ERC-OM Link: Z

0.3

T

= 2000 , P T = 500

For OM-Only Link
0.2

0.1

0
0

.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

6.5

7

7.5

8

8.5

9

9.5

10

Parameter r = k + / k - for Catreg receiver with Constant k - and k 0

Fig. 9: Capacity of the ERC-OM and OM-only links for ZT “ 2000 and PT “ 500.

Capacity of the ERC-OM Link For Varying Z T

1

0.9

0.8

0.7

Capacity

0.6

0.5

0.4

0.3
Parameter r = k + / k - = 2
Parameter r = k + / k -= 5

0.2

0.1

0
500

1000

1500

2000

2500

3000

3500

4000

4500

5000

Z T = 500-5000 , P = 200

Fig. 10: Effect on capacity with varying ZT .

March 6, 2018

DRAFT

28

bution is that we show that enzymatic reaction cycles can be used to increase the capacity of
communication link. We further show that we can increase the capacity of the link by increasing
the amount of a certain chemical species in the enzymatic reaction cycles. This provides a
practical way to adjust the link capacity.

DRAFT

March 6, 2018

29

A PPENDIX A
E XPRESSIONS
Note that in this section we present the derivation of Equations (31)-(33) and (40) for ERCOM link with RC receiver circuit only. The equations for the ERC-OM link with CATREG
receiver circuit can be derived using the same process.
For ease of reference, we first re-write the reactions in the ERC as mentioned in Section IV-A.

β

1
K`ZÝ
âÝ
Ýá
Ý C1

(53)

β2

k

C1 ÝÝ1Ñ K ` Z˚

(54)

α

1
ÝÝ
á
P ` Z˚ â
Ý C2
αÝ
2

(55)

k

C2 ÝÝ2Ñ P ` Z

(56)

Furthermore we present the reactions in the output Module (i.e. RC receiver circuit) as:

k`

Z˚ Ý
âÝ
Ýá
ÝX
k´

(57)

We obtain equation for Z˚ ptq from the above equations as follows:

Z9 ˚ ptq “ k1 C1 ptq ` α2 C2 ptq ´ α1 Z˚ ptqpPT ´ C2 ptqq ´ k` Z˚ ptq ` k´ nX ptq

(58)

The first term k1 C1 ptq is obtained as a result of Reaction (54). The second term α2 C2 ptq is
obtained as a result of backward Reaction in (55). Similarly the third term is obtained as a
result of forward Reaction in (55) where the number of P molecules is given by PT ´ C2 ptq.
Finally the last two terms are obtained as a result of forward and backward reactions in (57)
respectively.

March 6, 2018

DRAFT

30

Similarly, we find the expressions for C9 1 ptq and C9 2 ptq from the reactions in the ERC as
follows:
C9 1 ptq “ ´β2 C1 ptq ´ k1 C1 ptq ` β1 nL,R ptqpZT ´ Z˚ ptq ´ C1 ptq ´ C2 ptq ´ nX ptqq

(59)

C9 2 ptq “ ´α2 C2 ptq ´ k2 C2 ptq ` α1 Z˚ ptqpPT ´ C2 ptqq

(60)

For C9 1 ptq the first term comes from the backward Reaction in (53) whereas the second term
comes from the forward reaction in (55). The last term comes from the forward reaction in (53)
where the number of signalling molecules K in receiver is given by nL,R ptq and the number of
Z molecules is given by ZT ´ Z˚ ptq ´ C1 ptq ´ C2 ptq ´ nX ptq.
Similarly for C9 2 ptq the first term comes from the backward Reaction in (55) whereas the
second term comes from the forward reaction in (56). The last term comes from the forward
reaction in (55) where the number of signalling molecules Z˚ in receiver is given by Z˚ ptq and
the number of P molecules is given by PT ´ C2 ptq.
Finally the forward and backward reactions in (57) give us following equation for the number
of output molecules nX ptq:

n9 X ptq “ k` Z˚ ptq ´ k´ nX ptq

(61)

Note that Eq. (58), (59) and (60) and are nonlinear. To remove this non-linearity we first
assume that the constant ZT " (Z˚ ptq - C1 ptq - C2 ptq - nX ptqq. This enables us to simplify
Eq. (59) as:
C9 1 ptq “ ´pβ2 ` k1 qC1 ptq ` β1 nL,R ptqZT

(62)

We next assume that constant PT " C2 ptq, so we can replace the term PT - C2 ptq in Eq. (60)
by PT .

C9 2 ptq “ ´α2 C2 ptq ´ k2 C2 ptq ` α1 Z˚ ptqPT

(63)

Using the same assumption PT " C2 ptq in Eq. (58) we get:
DRAFT

March 6, 2018

31

Z9 ˚ ptq “ k1 C1 ptq ` α2 C2 ptq ´ α1 Z˚ ptqPT ´ k` Z˚ ptq ` k´ nX ptq

(64)

The resulting equations are still difficult to analyse so we use time scale separation between
the rate of diffusion and chemical reactions for further simplification. We define G1 “
G2 “

k´
.
k1

We assume 1 =

1
G1

and 2 =

1
G2

β1 ZT
d

and

are small and define a slow variable W ptq:

W ptq “ Z˚ ptq ` nX ptq

(65)

Using this value in Eq.(61) we get:

n9 X ptq “ k` Z˚ ptq ´ k´ nX ptq

(66)

Next we apply singular perturbation to separate the fast and slow dynamics. To do this we
multiply both sides by a small 2 “

k1
k´

and get :

2 n9 X ptq “ 2 k` Z˚ ptq ´ 2 k´ nX ptq

2 n9 X ptq “

k1
k1
k` Z˚ ptq ´
k´ nX ptq
k´
k´

(67)

(68)

Next we set 2 equal to zero and simplify the equation as follows:

0 ˆ n9 X ptq « k1

k`
Z˚ ptq ´ k1 nX ptq
k´

nX ptq “ rZ˚ ptq
where r “

(69)

(70)

k
k´

This means that we can simplify Equation (65) as:

March 6, 2018

DRAFT

32

Next we define G1 “

β1 ZT
d

, G2 “

W ptq “ Z˚ ptq ` nX ptq

(71)

W ptq “ Z˚ ptqr1 ` rs

(72)

k´
,
d

k11 “

k1
,
G1

kd “

k´
,
k`

a1 =

α1 PT
G1

, a2 =

α2
G1

and write Z˚ ptq

and nX ptq from Eq. (58) and (61) respectively using these notations:

d
Z9 ˚ ptq “ G1 k11 C1 ptq ` G1 a2 C2 ptq ´ G1 a1 Z˚ ptq ` G2 dnX ptq ´ G2 Z˚ ptq
kd

n9 X ptq “ ´G2 dnX ptq ` G2

d
Z˚ ptq
kd

(73)

(74)

9 ptq as follows :
Next we solve for W

9 ptq “ Z9 ˚ ptq ` n9 X ptq
W

(75)

Putting the values from Eqs. (73) and (74) we get:

9 ptq “ G1 k 1 C1 ptq ` G1 a2 C2 ptq ´ G1 a1 Z˚ ptq
W
1

(76)

9 ptq “ Z9 ˚ ptqr1 ` rs
W

(77)

We also know that:

Comparing both these equations we get following equation:

Z9 ˚ ptqr1 ` rs “ G1 k11 C1 ptq ` G1 a2 C2 ptq ´ G1 a1 Z˚ ptq

DRAFT

(78)

March 6, 2018

33

Z9 ˚ ptq “ rG1 k11 C1 ptq ` G1 a2 C2 ptq ´ G1 a1 Z˚ ptqsr

1
s
p1 ` rq

(79)

Next by using Laplace transform, we can obtain following equation:

Z˚ psq “

rG1 k11 C1 psq ` G1 a2 C2 psqs{p1 ` rq
s ` G1 a1 {p1 ` rq

(80)
(81)

Which reduces to following by using the values of G1 , k11 , a1 and a2 .
Z˚ psq “ r

pk1 C1 psq ` α2 C2 psq{p1 ` rq
s
s ` α1 PT {p1 ` rq

(82)

Furthermore we take Laplace of Eq. (70) and put the value of Z˚ psq from Eq. (82) to obtain:

NX psq “ rr

NX psq “ rZ˚ psq

(83)

pk1 C1 psq ` α2 C2 psq{p1 ` rq
s
s ` α1 PT {p1 ` rq

(84)

The next aim is to derive expressions of C1 psq and C2 psq respectively. We perform the Laplace
transform of (62) and go through some steps to obtain C1 psq in terms of the input signal U psq
as follows:

C1 psq “

β1 ZT
β1 ZT
NL,R psq “
1T psI ´ Hq´1 1T U psq
ps ` β2 ` k1 q
ps ` β2 ` k1 q R

(85)

where NL,R psq “ 1TR psI ´ Hq´1 1T U psq.
Similarly to obtain C2 psq we perform the Laplace transform of (63) and go through some
steps to simplify it as:
C2 psq “

March 6, 2018

α1 P
Z˚ psq
ps ` α2 ` k2 q

(86)

DRAFT

34

Putting value of Z˚ psq from Eq. (82) in Eq. (86) we get:

pk1 α1 P qβ1 ZT
NL,R psq
rps ` α2 ` k2 qps ` α1 P q ´ pα1 α2 P qsps ` β2 ` k1 q
pk1 α1 P qβ1 ZT
1T psI ´ Hq´1 1T U psq
C2 psq “
rps ` α2 ` k2 qps ` α1 P q ´ pα1 α2 P qsps ` β2 ` k1 q R

C2 psq “

(87)

The final result for ERC-OM link with RC receiver circuit is obtained by putting the values
of C1 psq and C2 psq Eq. (84) as follows:

pk1 C1 psq ` α2 C2 psq{p1 ` rq
s
(88)
s ` α1 PT {p1 ` rq
k1 β1 ZT
p1 ` α1 PT {rps ` α2 ` k2 qps ` α1 PT q ´ pα1 α2 PT qsq{p1 ` rq T
NX psq “ r
r
s1R psI ´ Hq´1 1T U psq
ps ` β2 ` k1 q
s ` α1 PT {p1 ` rq
NX psq “ rr

(89)
Taking
Qpsq “ rr

p1 ` α1 PT {rps ` α2 ` k2 qps ` α1 PT q ´ pα1 α2 PT qsq{p1 ` rq T
s1R psI ´ Hq´1 1T (90)
s ` α1 PT {p1 ` rq

We obtain final result as:

k1 β1 ZT
NX psq “ Qpsq
U psq
ps
`
β
`
k
q
2
1
loooooooooomoooooooooon

(91)

Ψ̃psq

Which is the same result as Equation (40).
Note that we can obtain the final result for the ERC-OM link with CATREG circuit in the
output module in similar way, however we are not presenting the complete derivation for this
case. An important similarity is that the expression for C1 psq is same for both RC and CATREG
circuits in the output module. However for CATREG circuit in the output module the expression
for C2 psq is different which is obtained as:
C2 psq “

DRAFT

pk1 α1 PT qβ1 ZT
NL,R psq
rps ` α2 ` k2 qps ` α1 PT ` k` ´ rk0 q ´ pα1 α2 PT qsps ` β2 ` k1 q

(92)

March 6, 2018

35

R EFERENCES
[1] I. Akyildiz, F. Brunetti, and C. Blázquez, “Nanonetworks: A new communication paradigm,” Computer Networks, vol. 52,
pp. 2260–2279, 2008.
[2] 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.
[3] T. Nakano, T. Suda, Y. Okaie, M. J. Moore, and A. V. Vasilakos, “Molecular Communication Among Biological
Nanomachines: A Layered Architecture and Research Issues,” NanoBioscience, IEEE Transactions on, vol. 13, no. 3,
pp. 169–197, 2014.
[4] I. Akyildiz, M. Pierobon, S. Balasubramaniam, and Y. Koucheryavy, “The internet of bio-nano things,” IEEE Communications Magazine, vol. 53, no. 3, pp. 32–40, 2015.
[5] T. Nakano, M. J. Moore, F. Wei, A. V. Vasilakos, and J. Shuai, “Molecular communication and networking: Opportunities
and challenges,” NanoBioscience, IEEE Transactions on, vol. 11, no. 2, pp. 135–148, 2012.
[6] M. Pierobon and I. 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.
[7] N. Farsad, A. W. Eckford, S. Hiyama, and Y. Moritani, “A simple mathematical model for information rate of active transport
molecular communication,” in Computer Communications Workshops (INFOCOM WKSHPS), 2011 IEEE Conference on.
IEEE, 2011, pp. 473–478.
[8] R. G. Gallager, Information theory and reliable communication.

Springer, 1968, vol. 2.

[9] 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.
[10] A. Einolghozati, M. Sardari, A. Beirami, and F. Fekri, “Capacity of discrete molecular diffusion channels.” ISIT, pp.
723–727, 2011.
[11] M. Pierobon and I. Akyildiz, “Capacity of a Diffusion-Based Molecular Communication System With Channel Memory
and Molecular Noise,” Information Theory, IEEE Transactions on, vol. 59, no. 2, pp. 942–954, 2013.
[12] A. Einolghozati, M. Sardari, and F. Fekri, “Capacity of diffusion-based molecular communication with ligand receptors.”
ITW, pp. 85–89, 2011.
[13] P. J. Thomas and A. W. Eckford, “Capacity of a simple intercellular signal transduction channel,” IEEE Transactions on
information Theory, vol. 62, no. 12, pp. 7358–7382, 2016.
[14] G. Aminian, M. Farahnak-Ghazani, M. Mirmohseni, M. Nasiri-Kenari, and F. Fekri, “On the capacity of point-to-point
and multiple-access molecular communications with ligand-receptors,” IEEE Transactions on Molecular, Biological and
Multi-Scale Communications, vol. 1, no. 4, pp. 331–346, 2015.
[15] C. T. Chou, “Molecular communication networks with general molecular circuit receivers,” in ACM The First Annual
International Conference on Nanoscale Computing and Communication. New York, New York, USA: ACM Press, 2014,
pp. 1–9.
[16] U. Alon, An Introduction to Systems Biology: Design Principles of Biological Circuits.

March 6, 2018

Chapman & Hall, 2006.

DRAFT

36

[17] K. S. Nilgiriwala, J. Jiménez, P. M. Rivera, and D. Del Vecchio, “Synthetic Tunable Amplifying Buffer Circuit in E. coli,”
ACS Synthetic Biology, vol. 4, no. 5, pp. 577–584, May 2015.
[18] D. Del Vecchio, A. J. Ninfa, and E. D. Sontag, “Modular cell biology: retroactivity and insulation,” Molecular systems
biology, vol. 4, Feb. 2008.
[19] H. Awan, “Reducing the effect of reaction rate constants on the performance of molecular communication networks,” in
Proceedings of the 3rd ACM International Conference on Nanoscale Computing and Communication, ser. NANOCOM’16.
New York, NY, USA: ACM, 2016, pp. 8:1–8:6. [Online]. Available: http://doi.acm.org/10.1145/2967446.2967455
[20] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular Biology of the Cell, 5th ed.

Garland

Science, 2007.
[21] N. Farsad, H. B. Yilmaz, A. Eckford, C.-B. Chae, and W. Guo, “A comprehensive survey of recent advancements in
molecular communication,” IEEE Communications Surveys and Tutorials, vol. 18, no. 3, pp. 1887–1919, 2016.
[22] H. ShahMohammadian, G. G. Messier, and S. Magierowski, “Optimum receiver for molecule shift keying modulation in
diffusion-based molecular communication channels,” Nano Communication Networks, vol. 3, no. 3, pp. 183–195, Sep.
2012.
[23] M. Kuran, H. Yilmaz, T. Tugcu, and I. Akyildiz, “Modulation Techniques for Communication via Diffusion in
Nanonetworks,” in Communications (ICC), 2011 IEEE International Conference on, 2011, pp. 1–5.
[24] M. Mahfuz, D. Makrakis, and H. Mouftah, “On the characterization of binary concentration-encoded molecular communication in nanonetworks,” Nano Communication Networks, vol. 1, pp. 289–300, 2010.
[25] C. T. Chou, “A Markovian approach to the optimal demodulation of diffusion-based molecular communication networks,”
Communications, IEEE Transactions on, vol. PP, no. 99, pp. 1–1, 2015.
[26] H. Awan and C. T. Chou, “Impact of receiver molecular circuits on the performance of reaction shift keying,” in Proceedings
of the Second Annual International Conference on Nanoscale Computing and Communication, ser. NANOCOM’ 15.
New York, NY, USA: ACM, 2015, pp. 2:1–2:6. [Online]. Available: http://doi.acm.org/10.1145/2800795.2800798
[27] ——, “Generalized solution for the demodulation of reaction shift keying signals in molecular communication networks,”
IEEE Transactions on Communications, vol. 2, pp. 715–727, 2017, doi:10.1109/TCOMM.2016.2638899.
[28] M. Mahfuz, D. Makrakis, and H. Mouftah, “On the characteristics of concentration-encoded multi-level amplitude
modulated unicast molecular communication,” in Electrical and Computer Engineering (CCECE), 2011 24th Canadian
Conference on, 2011, pp. 000 312–000 316.
[29] R. Mosayebi, H. Arjmandi, A. Gohari, M. Nasiri-Kenari, and U. Mitra, “Receivers for diffusion-based molecular
communication: Exploiting memory and sampling rate,” Selected Areas in Communications, IEEE Journal on, vol. 32,
no. 12, pp. 2368–2380, 2014.
[30] C. T. Chou, “Extended master equation models for molecular communication networks,” IEEE Transactions on Nanobioscience, vol. 12, no. 2, pp. 79–92, 2013, doi:10.1109/TNB.2013.2237785.
[31] ——, “Noise properties of linear molecular communication networks,” Nano Communication Networks, vol. 4, pp. 87–97,
2013.

DRAFT

March 6, 2018

37

[32] R. Erban, J. Chapman, and P. Maini, “A practical guide to stochastic simulations of reaction-diffusion processes,” arXiv
preprint arXiv:0704.1908, 2007.
[33] 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.
[34] ——, “Diffusion-based Noise Analysis for Molecular Communication in Nanonetworks,” IEEE Transactions on Signal
Processing, vol. 59, no. 6, pp. 2532–2547, 2011.
[35] A. Noel, K. C. Cheung, and R. Schober, “Improving Receiver Performance of Diffusive Molecular Communication With
Enzymes,” NanoBioscience, IEEE Transactions on, vol. 13, no. 1, pp. 31–43, 2014.
[36] ——, “Optimal Receiver Design for Diffusive Molecular Communication With Flow and Additive Noise,” IEEE
Transactions on Nanobioscience, vol. 13, no. 3, pp. 350–362, Sep. 2014.
[37] D. Kilinc and O. B. Akan, “Receiver Design for Molecular Communication,” Selected Areas in Communications, IEEE
Journal on, vol. 31, no. 12, pp. 705–714, Dec. 2013.
[38] M. U. Mahfuz, D. Makrakis, and H. T. Mouftah, “Strength-based optimum signal detection in concentration-encoded
pulse-transmitted OOK molecular communication with stochastic ligand-receptor binding,” Simulation Modelling Practice
and Theory, vol. 42, pp. 189–209, 2014.
[39] L.-S. Meng, P.-C. Yeh, K.-C. Chen, and I. F. Akyildiz, “On Receiver Design for Diffusion-Based Molecular Communication,” IEEE Transactions on Signal Processing, vol. 62, no. 22, pp. 6032–6044, Nov. 2014.
[40] H. Awan and C. T. Chou, “Demodulation of reaction shift keying signals in molecular communication network with protein
kinase receiver circuit,” in Wireless Communications and Networking Conference (WCNC), 2016 IEEE.

IEEE, 2016.

[41] C. T. Chou, “Maximum a-posteriori decoding for diffusion-based molecular communication using analog filters,”
Nanotechnology, IEEE Transactions on, vol. 14, no. 6, pp. 1054–1067, 2015.
[42] ——, “Impact of Receiver Reaction Mechanisms on the Performance of Molecular Communication Networks,” IEEE
Transactions on Nanotechnology, vol. 14, no. 2, pp. 304–317, Mar. 2015.
[43] ——, “Molecular circuits for decoding frequency coded signals in nano-communication networks,” Nano Communication
Networks, vol. 3, no. 1, pp. 46–56, 2012.
[44] M. Pierobon, “A Molecular Communication System Model Based on Biological Circuits,” in ACM The First Annual
International Conference on Nanoscale Computing and Communication. New York, New York, USA: ACM Press, 2014,
pp. 1–8.
[45] M. Mahdavifar and A. Beirami, “Diffusion channel with poisson reception process: capacity results and applications,” in
In Proceedings of IEEE International Symposium on Information Theory (ISIT), 2015.
[46] S. Cardinale and A. P. Arkin, “Contextualizing context for synthetic biology–identifying causes of failure of synthetic
biological systems,” Biotechnology journal, vol. 7, no. 7, pp. 856–866, 2012.
[47] D. Del Vecchio, “A control theoretic framework for modular analysis and design of biomolecular networks,” Annual
Reviews in Control, vol. 37, no. 2, pp. 333–345, 2013.
[48] D. Del Vecchio, A. J. Ninfa, and E. D. Sontag, “Modular cell biology: retroactivity and insulation,” Molecular systems
biology, vol. 4, no. 1, p. 161, 2008.

March 6, 2018

DRAFT

38

[49] J. Saez-Rodriguez, A. Kremling, H. Conzelmann, K. Bettenbrock, and E. D. Gilles, “Modular analysis of signal transduction
networks,” Control Systems, IEEE, vol. 24, no. 4, pp. 35–52, 2004.
[50] H. B. Yilmaz, Y.-J. Cho, W. Guo, and C.-B. Chae, “Interference reduction via enzyme deployment for molecular
communication,” Electronics Letters, vol. 52, no. 13, pp. 1094–1096, 2016.
[51] Y. J. Cho, H. B. Yilmaz, W. Guo, and C.-B. Chae, “Effective enzyme deployment for degradation of interference molecules
in molecular communication,” in Wireless Communications and Networking Conference (WCNC), 2017 IEEE. IEEE, 2017,
pp. 1–6.
[52] U. A. Chude-Okonkwo, “Diffusion-controlled enzyme-catalyzed molecular communication system for targeted drug
delivery,” in Global Communications Conference (GLOBECOM), 2014 IEEE.
[53] C. Gardiner, Stochastic methods.

IEEE, 2014, pp. 2826–2831.

Springer Berlin, 2009.

[54] D. J. Higham, “Modeling and Simulating Chemical Reactions,” SIAM Review, vol. 50, no. 2, p. 347, 2008.
[55] W. H. de Ronde, F. Tostevin, and P. R. ten Wolde, “Feed-forward loops and diamond motifs lead to tunable transmission
of information in the frequency domain,” Physical Review E, vol. 86, no. 2, p. 021913, Aug. 2012.
[56] P. B. Warren, S. Tănase-Nicola, and P. R. ten Wolde, “Exact results for noise power spectra in linear biochemical reaction
networks,” The Journal of chemical physics, vol. 125, no. 14, p. 144904, 2006.
[57] F. Tostevin and P. ten Wolde, “Mutual information in time-varying biochemical systems,” Physical Review E, vol. 81,
no. 6, Jun. 2010.
[58] A. Papoulis and S. U. Pillai, Probability, Random Variables and Stochastic Processes.

McGraw Hill, 2002.

[59] H. K. Khalil and J. Grizzle, “Nonlinear systems, vol. 3,” Prentice hall Upper Saddle River, 2002.
[60] D. Mishra, P. M. Rivera, A. Lin, D. Del Vecchio, and R. Weiss, “A load driver device for engineering modularity in
biological networks,” Nature Biotechnology, vol. 32, no. 12, pp. 1268–1275, Dec. 2014.
[61] D. Gillespie, “Exact numerical simulation of the Ornstein-Uhlenbeck process and its integral,” Physical Review E, 1996.

DRAFT

March 6, 2018

