Extended master equation models for molecular communication networks Chun Tung Chou School of Computer Science and Engineering, 3 University of New South Wales, 1 Sydney, Australia. E-mail: [email protected] 0 2 January 4, 2013 n a J Abstract 3 We consider molecular communication networks consisting of transmitters and receivers ] distributedinafluidicmedium. Insuchnetworks,atransmittersendsoneormoresignalling I N 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 . s models for these networks are required. This paper proposes a new stochastic model for c 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 4 of signalling molecule counts, while diffusion in the medium and chemical reactions at the v receiversaremodelledasMarkovprocessesusingmasterequation. AnadvantageofRDMEX 3 is that it can readily be used to model molecular communication networks with multiple 5 transmitters and receivers. For the case where the reaction kinetics at the receivers is linear, 2 we show how RDMEX can be used to determine the mean and covariance of the receiver 4 output signals, and derive closed-form expressions for the mean receiver output signal of the . 4 RDMEX model. These closed-form expressions reveal that the output signal of a receiver 0 can be affected by the presence of other receivers. Numerical examples are provided to 2 demonstrate the properties of the model. 1 : Keywords: Molecularcommunicationnetworks,nanocommunicationnetworks,syntheticmolec- v ular communication networks, master equations, stochastic models, synthetic biology i X r a 1 Introduction We consider molecular communication networks consisting of transmitters and receivers dis- tributed 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 com- munication 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, cellsinthehumanbodyconstantlycommunicatewithothercellsusingmolecularcommunication. 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 math- ematical 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 sig- nal 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 stochasticmodelformolecularcommunicationnetworks. Ourmodelisbasedonreaction-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 keyideabehindRDMEXistomodelthetransmittersastimeseriesofsignallingmoleculecounts, 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 net- works 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 therearemultipletransmittersandreceivers. Theseexpressionsshowthattheoutputofareceiver 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. InSection3, wepresenttheRDMEXandshowhowitcanbeusedtomodelmolecular 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 Theaimofthissectionistoprovidethenecessarybackgroundonmasterequations. Thetreatment 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 Consideracontinuous-timeinteger-valuevectorMarkovprocessQ(t) Zp,wherepisthenumber ∈ 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+r (where r Zp with j = 1,2,...J and J is the total j j ∈ ∈ number of possible jumps) at a transition rate of W (q). Let P(q,tq ,t ) denote the conditional j 0 0 | probability that Q(t)=q given that Q(t )=q . 0 0 We are interested to determine how P(q,tq ,t ) evolves over time. We can do this by using a 0 0 | coupled set of ordinary differential equations (ODEs) known as the master equation: J J dP(q,tq0,t0) (cid:88) (cid:88) | = W (q r )P(q r ,tq ,t ) W (q)P(q,tq ,t ) (1) j j j 0 0 j 0 0 dt − − | − | j=1 j=1 2 whereone equationof theform(1) isneeded foreachvalid state q. Notethat thefirstand second terms on the right-hand side of (1) can be interpreted, respectively, as the rates of entering and leavingthestateq. Inordertosimplifynotation,wewillwriteP(q,t)insteadofP(q,tq ,t )from 0 0 | 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: L + R (cid:41)k(cid:42)+ 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 =[n ,n ,n ]T where n is the number of molecules of chemical L etc., and T denotes matrix L R C L transpose. There are two possible types of jumps (i.e. J = 2) in this system. The forward reaction is modelled by the jump r = [ 1, 1,1]T where the entries of r reflects the fact that one molecule 1 1 − − 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 W (q) = k q(1)q(2) = k n n where q(1) is the 1 + + L R first component of the vector q, etc. Similarly, the reverse reaction is modelled by the jump r =[1,1, 1]T with W (q)=k q(3)=k n . The master equation for this chemical reaction is: 2 2 − − C − 2 2 dP(q,t) (cid:88) (cid:88) = W (q r )P(q r ,t) W (q)P(q,t) (2) j j j j dt − − − j=1 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)canbedifficulttoworkwithbecauseweneedoneequationforeachvalidstate, 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: (cid:88) (cid:88) Q(t) = qProb(Q(t)=q)= qP(q,t) (3) (cid:104) (cid:105) q q (cid:88) Σ(t) = (q Q(t) )(q Q(t) )TP(q,t) (4) −(cid:104) (cid:105) −(cid:104) (cid:105) q where will be used in this paper to denote the mean operator. (cid:104)•(cid:105) 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 J d Q(t) (cid:88) (cid:104) (cid:105) = r W (Q(t)) (5) j j dt (cid:104) (cid:105) j=1 In particular, if W (q) is a linear function of q. Let (cid:80)J r W (q)=Aq, then j j=1 j j d Q(t) (cid:104) (cid:105) = A Q(t) (6) dt (cid:104) (cid:105) J dΣ(t) (cid:88) = AΣ(t)+Σ(t)AT + r rTW ( Q(t) ) (7) dt j j j (cid:104) (cid:105) j=1 3 3 Modelling molecular communication network using RD- MEX Weconsideramolecularcommunicationnetworkwithmultipletransmittersandmultiplereceivers 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, theydiffuseinthemediumaccordingtoBrownianmotion. Thereceiversareassumedtoconsistof oneormorereceptors. WhenasignallingmoleculeLreachesareceptorR,theymaybindtogether 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 usedinchemotaxismodelsinbiophysics[11,12]andmolecularcommunicationnetworkmodelsin engineering [13]. Based on the above description, we see that a model for molecular communication networks mustatleastcapturethediffusionofsignallingmoleculesandthereactionsatthereceivers. Both reactionsanddiffusioncanbemodelledasMarkovprocesses,sothemasterequation(1)isanatural choice. The remaining issue is how we can model the transmitter. In this paper, we model the transmitterbyatimesequencewhichspecifiesthenumberofmoleculesemittedbythetransmitter at a particular time. This approach is fairly general and can be used to model encoding methods thathavebeenconsideredintheliterature, suchasmolecularcoding[1]andconcentrationcoding [14]. We will consider other modelling approaches, e.g. modelling the internal mechanism of the transmitters, in future work. Insection3.1,weintroducetheRDMEXmodelbywayofanexampleandthenweprovesome 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. Itisfairlysimpletogeneralisethemodeltomultipletransmitters,multiplereceivers 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 k signalling molecules of L at time t , k 1,1 1,1 1,2 molecules of L at time t , ..., k molecules of L at time t , where k (b=1,2,3,...) are 1,2 1,b 1,b 1,b positiveintegersandt R. Similarly,transmitter2emitsk signallingmoleculesofLat 1,b 2,b ∈ time t where b=1,2,3,.... The number of molecules k emitted at time t is assumed 2,b a,b a,b to be independent of the state of the system at or before t . a,b 3. Themediumisassumedtobea1-dimensionalspaceoflengthX. Themediumispartitioned 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 fromonevoxeltoaneighbouringvoxelis˘dpermoleculeperunittime. Therateofdiffusion fromavoxeltoanon-neighbouringvoxeliszero. Also,themoleculecannotleavethemedium and we assume the boundary is reflective. The parameter ˘d is related to one-dimensional macroscopic diffusion constant D˘ by ˘d= D˘ . ∆x2 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 T and T . Similarly, we assume that 1 2 receivers 1 and 2 are located in the voxels indexed by R and R . The indices T , T , R 1 2 1 2 1 and R are integers in the interval [1,N] and are assumed to be distinct. (Note that it is 2 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: L (cid:41)k(cid:42)+ 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 R -th voxel. The situation at receiver 2 is similar. 1 Remark 1 The reaction kinetics assumed above can be viewed as a linearisation of the second order reaction L+R (cid:41)(cid:42) 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 (cid:2) (cid:3) q = n n .... n n n (8) L,1 L,2 L,N C,1 C,2 wheren representsthenumberofsignallingmoleculesatthej-thvoxelandn represents L,j C,u the number of complexes at the u-th receiver. 8. We define two indicator vectors 1 ,1 ZN+2. The T -th element of 1 is 1 and are T1 T2 ∈ 1 T1 otherwise zero. 1 is similarly defined. T2 9. In order to write down the diffusion and reaction within the molecular communication net- work, we define the following state transition vectors r and transition rates W . The total j j 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 r ZN+2 and only the non-zero elements of r are stated, and q j j ∈ 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 r and j ≤ ≤ − W (q). Specifically, r (j)= 1, r (j+1)=1 and W (q)=˘dq(j)=˘dn . j j j j L,j − Explanation: Ifasignallingmoleculediffusesfromvoxelj toj+1,itmeansthenumber of signalling molecules in voxel j is decreased by one (hence r (j) = 1) and that in j − voxel j+1 is increased by one (hence r (j+1)=1). The rate at which this particular j 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 r to be a zero vector and W (q)=0. N N (b) For the diffusion of L from voxel j to j 1 where 2 j N, r (j)= 1, r (j N+j N+j − ≤ ≤ − − 1) = 1 and W (q) = ˘dq(j) = ˘dn . For convenience, we define r to be a zero N+j L,j N+1 vector and W (q)=0. N+1 (c) For receiver 1, the vector r and the rate W (q) are used to model the forward 2N+1 2N+1 reaction of the conversion of a signalling molecule L to a complex C. Specifically, r (R )= 1, r (N +1)=1 and W (q)= k+q(R )= k+n . 2N+1 1 − 2N+1 2N+1 ∆x 1 ∆x L,R1 Explanation: In the forward reaction, a signalling molecule is removed in the R -th 1 voxel, hence r (R ) = 1 and a complex is formed, hence r (N + 1) = 1 2N+1 1 2N+1 − 5 because the number of complexes at receiver 1 is the (N +1)-th element in the state vector (8). The rate W (q) is proportional to the number of signalling molecules in 2N+1 the R -th voxel where the receiver is located. 1 Forthereversereaction,r (N+1)= 1,r (R )=1andW (q)=k q(N+ 2N+2 2N+2 1 2N+2 − − 1)=k n . − C,1 (d) Forreceiver2,theforwardreaction: r (R )= 1,r (N+2)=1andW (q)= 2N+3 2 2N+3 2N+3 − k+q(R ) = k+n . For the reverse reaction, r (N +2) = 1, r (R ) = 1 ∆x 2 ∆x L,R2 2N+4 − 2N+4 2 and W (q)=k q(N +2)=k n . 2N+4 − − C,2 The RDMEX model for the 2-transmitter 2-receiver molecular communication network is 2 ∞ dP(q,t) (cid:88)(cid:88) = P(q k 1 ) P(q,t) δ(t t ) dt { − a,b Ta − } − a,b a=1b=1 J J (cid:88) (cid:88) + W (q r )P(q r ,t) W (q)P(q,t) (9) j j j j − − − j=1 j=1 where δ(t) denotes the Dirac delta function. Letus,forthetimebeing,assumethatthefirsttermontheright-handsideof (9)isnotthere. Ifthisisthecase,then(9)isofthesameformasthemasterequation(1)andtheequationmodels 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 t , the a,b a-th transmitter emits k signalling molecules into the T -th voxel (where the a-th transmitter a,b a is located). This means that if the system is in the state q just before the time t (denoted as a,b t− ), then it will be in state q +k 1 just after t (denoted as t+ ). In addition, we have a,b a,b Ta a,b a,b P(q,t− )=P(q+k 1 ,t+ ), which is modelled by the first term in (9). Note that it is possible a,b a,b Ta a,b to give a stochastic interpretation of k , see Remark 2. a,b The deterministic input in (9) can be thought of as an external arrival of the system. We will referto(9)asreaction-diffusionmasterequationwithexogenousinput,orRDMEXforshort. 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 W (q) is a linear function of q. Let j 6 (cid:80)J r W (q)=Aq, then j=1 j j 2 ∞ d Q(t) (cid:88)(cid:88) (cid:104) (cid:105) = A Q(t) + k 1 δ(t t ) (10) dt (cid:104) (cid:105) a,b Ta − a,b a=1b=1 J dΣ(t) (cid:88) = AΣ(t)+Σ(t)AT + r rTW ( Q(t) ) (11) dt j j j (cid:104) (cid:105) j=1 Proof: For the time evolution on the mean Q(t) , we can start with the derivative of (3): d(cid:104)Q(t)(cid:105) =(cid:80) qdP(q,t) and then use (9) for dP(q,t(cid:104)). Th(cid:105)is is fairly straightforward and uses exactly dt q dt 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 k molecules arriving at time t is to add k 1 to the a,b a,b a,b Ta state vector. Hence the form of (10). For the evolution of covariance, one can follow the derivation in [10] provided that the im- pulses 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− ) and Q(t− ) be the state and mean state just before the deterministic arrival at a,b (cid:104) a,b (cid:105) time t . At time t+ , just after t , the state of the system will become Q(t− )+k 1 . Also, a,b a,b a,b a,b a,b Ta the mean state at t+ is: a,b (cid:88) (cid:88) Q(t+ ) = qP(q,t+ )= qP(q k 1 ,t− ) (cid:104) a,b (cid:105) a,b − a,b Ta a,b q q (cid:88) = (q+k 1 )P(q,t− )= Q(t− ) +k 1 (12) a,b Ta a,b (cid:104) a,b (cid:105) a,b Ta q Note that we have used the fact that P(q,t+ )=P(q k 1 ,t− ) in the above derivation. The a,b − a,b Ta a,b overall result is that, at time t , the state Q(t) and mean state Q(t) are incremented by the a,b (cid:104) (cid:105) same vector k 1 . a,b Ta 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). (cid:3) For the rest of the paper, we will focus on studying the properties of equation (10), though we willpresentanumericalexampleinSection6todemonstratetheaccuracyof (11). Adetailstudy 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 assum- ing a deterministic emission of exactly k signalling molecules by the a-th transmitter at time a,b t , one may assume that the number of molecules emitted is a random variable K with mean a,b a,b K and covariance cov(K ). Provided that the random variable K is independent of the a,b a,b a,b (cid:104) (cid:105) stateq attimet orearlier, similarresultstoProposition2canbederived. Forequation (10), we a,b need to replace k by K , and we need to add cov(K ) to the right-hand side of (11). This a,b a,b a,b (cid:104) (cid:105) generalisation says that one can interpret k in (10) as the mean number of molecules emitted a,b at time t by the a-th transmitter. With this stochastic interpretation of k , one can consider a,b a,b 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 molec- ular 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. Inordertosimplifythepresentation,wehavesofarlimitedourstudyto1-dimensional → but given most molecular communication networks are 3-dimensional, we will generalise the con- tinuum 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, whentheinterval∆xgoestozero, (10)convergestoaRDPDEandanumberofchemicalkinetics 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 (cid:2) (cid:3) q = .... n n n n n .... n n (13) L,−2 L,−1 L,0 L,1 L,2 C,1 C,2 where, as before, n is the number of L in the j-th voxel where j Z, and n is the number L,j C,u ∈ of complexes formed at the u-th receiver. Given this state vector, we can write equation (10) as 2 ∞ d(cid:104)nL,j(t)(cid:105) = ˘d( n (t) 2 n (t) + n (t) )+(cid:88)(cid:88)δ (j T )k δ(t t ) L,j−1 L,j L,j+1 K a a,b a,b dt (cid:104) (cid:105)− (cid:104) (cid:105) (cid:104) (cid:105) − − a=1b=1 2 (cid:88)δ (j R )(k+ n (t) k n (t) ) j Z (14) − K − u ∆x(cid:104) L,Ru (cid:105)− −(cid:104) C,u (cid:105) ∀ ∈ u=1 d n (t) k C,u + (cid:104) (cid:105) = n (t) k n (t) for u=1,2 (15) dt ∆x(cid:104) L,Ru (cid:105)− −(cid:104) C,u (cid:105) where δ (j) is the Kronecker delta1 . K Suppose the centre of the voxel j is at position x , we replace n (t) by (cid:96)(x ,t)∆x where j L,j j (cid:104) (cid:105) (cid:96)(x ,t) is the mean concentration in voxel j at time t. By dividing both sides of (14) by ∆x, j taking the limit ∆x 0 and noting that ˘d(∆x)2 =D˘, we have → ∂(cid:96) = D˘ ∂2(cid:96) +(cid:88)2 δ(x x )(cid:88)∞ k δ(t t ) (cid:88)2 δ(x x )d(cid:104)nC,u(t)(cid:105) (16) ∂t ∂x2 − T,a a,b − a,b − − R,u dt a=1 b=1 u=1 (cid:124) (cid:123)(cid:122) (cid:125) =ka(t) d n (t) C,u (cid:104) (cid:105) = k (cid:96)(x ,t) k n (t) for u=1,2 (17) + R,u − C,u dt − (cid:104) (cid:105) wherex (resp. x )isthecentreofthevoxelT (R )wherethea-thtransmitter(u-threceiver) T,a R,u a u islocated. NotethatwehavealsousedthefollowingconversionbetweentheKroneckerandDirac deltas: lim δK(j) =δ(x). ∆x→0 ∆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. 1Note: WeusebothKoneckerdeltaandDiracdeltainthispaper. Theyaredenoted,respectively,asδK(•)and δ(•). 8 The second term in (16) says that the transmitter at T adds signalling molecules to the system a accordingtotimesequencek (t), whilethethirdtermssaysthesignallingmoleculesareabsorbed a fromthesystemiftheyformcomplexesatthereceivers. Giventhatweassumethatonesignalling 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 denotesapointinthe3-dimensionalspace(i.e.v isa3-dimensionalvector)and(cid:96)(v,t)isthemean concentration of the signalling molecule at the location v at time t. 2 2 ∂(cid:96) = D 2(cid:96)+(cid:88)δ(v v )k (t) (cid:88)δ(v v )dcu(t) (18) T,a a R,u ∂t ∇ − − − dt a=1 u=1 dc (t) u = k (cid:96)(x ,t) k c (t) for u=1,2 (19) + R,u − u dt − where 2 is the Laplacian in 3-dimensional space, and v (resp. v ) is a 3-dimensional vector T,a R,u ∇ specifying the location of a-th transmitter (u-th receiver). Note that we have also introduced a new notation c (t) to denote the mean number of complexes n (t) at receiver u; this is to u C,u (cid:104) (cid:105) 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 = D . Also, the rate of formation ∆2 of complexes at the receiver voxel is given by k+ times the number of signalling molecules L in ∆3 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 k (t) and k (t) (which models the 1 2 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 c (t) and 1 c (t) (which can be viewed as the outputs). We will divide this section into two parts. We will 2 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 k (t) and k (t) to c (t) and 1 2 1 c (t) using (18) and (19). 2 We first define a few notation. Let ι = √ 1, and C (ω), K (ω) and L˜(v,ω) be the temporal u a − Fourier transform of, respectively, c (t), k (t) and (cid:96)(v,t) where ω is the transform variable. It is u a shown in Appendix A that 2 2 (cid:88) (cid:88) L˜(v,ω) = φ(v v ,ω)K (ω) φ(v v ,ω)ιωC (ω) (20) T,a a R,u u − − − a=1 u=1 9 where (cid:114) 1 ιω φ(v,ω) = exp( v ) (21) 4πD v − D(cid:107) (cid:107) (cid:107) (cid:107) isthetemporalFouriertransformofthe3-dimensionaldiffusionkernel 1 exp( (cid:107)v(cid:107)2). Given (4πDt)32 −4Dt (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 =v and then v =v in (20), we have: R,1 R,2 L˜(v ,ω) = φ (ω)K (ω)+φ (ω)K (ω) φ (ω)ιωC (ω) φ (ω)ιωC (ω) (22) R,1 11 1 12 2 0 1 ∆R 2 − − L˜(v ,ω) = φ (ω)K (ω)+φ (ω)K (ω) φ (ω)ιωC (ω) φ (ω)ιωC (ω) (23) R,2 21 1 22 2 ∆R 1 0 2 − − where φ (ω) = φ(v v ,ω) for u,a=1,2 (24) ua R,u T,a − φ (ω) = φ(v v ,ω) (25) ∆R R,1 R,2 − φ (ω) = φ(0,ω) where 0= zero vector (26) 0 One can interpret φ (ω) as the transfer function which models the dynamics of the diffusion ua of molecules from the a-th transmitter located at v (where the molecules are injected into the T,a medium) to the u-th receiver located at v . R,u Whensignallingmoleculesareabsorbedtoformcomplexes,itcreatesaconcentrationgradient. Thediffusiondynamicsbetweenthetworeceiversaremodelledbyφ (ω). Giventhatweassume ∆R that the locations of the transmitters and receivers are distinct, both φ (ω) and φ (ω) are well ua ∆R defined. Thetransferfunctionφ (ω)modelsthelocalimpactofabsorptionofsignallingmoleculeateach 0 receiver. This transfer function is unfortunately not well defined as can be seen from substituting v =0inthedefinitionofφ(v,ω)in(21). Thetransferfunctionφ (ω)alsoappearsinthemodelling 0 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 k (t) are absent. Both a [11,12]dealwiththeindefinitenessofφ (ω)bycuttingoffanintegraltoevaluateφ (ω)atafinite 0 0 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 φ (ω), and we will show using 0 numericalexamplesinsection6thatourapproximationgivesaccurateresults. Fortherestofthis section, we will continue to use equations (22) and (23) with the understanding that φ (ω) is not 0 well defined but can be well approximated. Bothequations(22)and(23)areobtainedfrom(18). Westillneedtoworkon(19). Bytaking the Fourier transform of (19), we have C (ω) k 1 + ρ (ω) = = (27) 1 L˜(vR,1,ω) ιω+k− C (ω) k 2 + ρ (ω) = = (28) 2 L˜(vR,2,ω) ιω+k− where ρ (ω) and ρ (ω) are transfer functions that model the reaction kinetics at the receivers. 1 2 Given that we have assumed that the binding and unbinding rate constants at both receivers are identical, it is not surprising that ρ (ω) and ρ (ω) are the same here. It is straightforward to 1 2 generalise to the case where the receivers have different reaction kinetics. By using equations (22), (23), (27) and (28), we can eliminate L˜(v ,ω) and L˜(v ,ω) to R,1 R,2 obtain the transfer function from the inputs k (t) and k (t) to the outputs c (t) and c (t): 1 2 1 2  −1 (cid:20) C (cid:21)  (cid:20) ρ 0 (cid:21)(cid:20) φ φ (cid:21) (cid:20) ρ 0 (cid:21)(cid:20) φ φ (cid:21)(cid:20) K (cid:21) 1 = I+ιω 1 0 ∆r  1 11 12 1 (29) C  0 ρ φ φ  0 ρ φ φ K 2  2 ∆r 0  2 21 22 2 (cid:124) (cid:123)(cid:122) (cid:125)(cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125) R Φ0 Φ 10