- Email: [email protected]
Design of a Molecular Communication Channel by Modelling Enzyme Kinetics
8th Vienna International Conference on Mathematical Modelling 8th Vienna18International Conference on Mathematical Modelling February - 20, 2015. Vienna University of Technology, Vienna, 8th Vienna18International Conference on Mathematical Modelling February - 20, 2015. Vienna University of Technology, Austria Available online atVienna, www.sciencedirect.com February 18 - 20, 2015. Vienna University of Technology, Vienna, Austria Austria
ScienceDirect
IFAC-PapersOnLine 48-1 (2015) 035–040
Design of a Molecular Communication Design of a Molecular Communication Design of a Molecular Communication Channel by Modelling Enzyme Kinetics Channel by Modelling Enzyme Channel by Modelling Enzyme Kinetics Kinetics
Rudolf Rabenstein ∗∗ Rudolf Rabenstein ∗∗ Rudolf Rabenstein ∗ at Erlangen-N¨ urnberg (FAU), ∗ Friedrich-Alexander-Universit¨ at Erlangen-N¨ urnberg (FAU), ∗ Friedrich-Alexander-Universit¨ ∗ D-91058 Erlangen, Cauerstr. 7, Germany (e-mail: [email protected]). Friedrich-Alexander-Universit¨ at Erlangen-N¨ urnberg (FAU), D-91058 Erlangen, Cauerstr. 7, Germany (e-mail: [email protected]). D-91058 Erlangen, Cauerstr. 7, Germany (e-mail: [email protected]). Abstract: Molecular communication channels transmit information by the release and detection Abstract: communication channels information by the release and of detection of chemical Molecular substances. Single information bitstransmit can be encoded as individual releases certain Abstract: Molecular communication channels transmit information by the release and detection of chemical substances. Single information bits can be encoded as individual releases of certain amounts of substances. molecules. Intersymbol interference between successive releases can be minimized of chemical Single information bits can be encoded as individual releases of certain amounts of molecules. Intersymbol interference between successive releases can be minimized by enforcing a decay of the substrate concentration through enzyme reactions. A model for amounts of molecules. Intersymbol interference between successive releases can be minimized by enforcing a behaviour decay of the substrate through to enzyme A model for the short term of the enzymeconcentration kinetics is employed derivereactions. the achievable symbol by enforcing a decay of the substrate concentration through enzyme reactions. A model for the short term behaviour of the enzyme kinetics is employed to derive the achievable symbol transmission rate in termsofof the the enzyme parameters of the enzyme reaction andthe theachievable sensitivitysymbol of the the short term behaviour kinetics is employed to derive transmission rate inThe terms of the relation parameters of the the enzyme reaction the sensitivity of the substrate receptor. resulting facilitates design of the and biochemical environment transmission rate inThe terms of the relation parameters of the the enzyme reaction the sensitivity of the substrate receptor. resulting facilitates design of the and biochemical environment for efficient molecular communication. substrate receptor. The resulting relation facilitates the design of the biochemical environment for efficient molecular communication. for efficient molecular communication. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Communication channels, Chemical microsensors, Models, Molecular Keywords: Communication channels, Chemical microsensors, Models, Molecular communication, Enzyme kinetics, Michaelis-Menten model Models, Keywords: Communication channels, Chemical microsensors, Molecular communication, Enzyme kinetics, Michaelis-Menten model communication, Enzyme kinetics, Michaelis-Menten model 1. INTRODUCTION Sec. 2 presents the channel model under investigation and 1. INTRODUCTION Sec. 2 presents the channel model under investigation and thechannel corresponding set of investigation differential equa1. INTRODUCTION Sec. 32 introduces presents the model under and Sec. 3 introduces the corresponding set of differential equations.3 Of special interest are the extreme discussed in Sec. introduces the corresponding set ofcases differential equations.4 Of special interest are the extreme cases discussed Sec. which leadinterest to the are useful derived in in tions. Of special the approximations extreme cases discussed in Sec. 5. 4 which lead to the useful approximations derived in Classical communication channels are based on the prop- Sec. Sec. 6 lead shows application to the above channel 4 which totheir the useful approximations derived in Classical communication channels are based on the propSec. 5. Sec. 6 shows their application to the above channel agation ofcommunication electromagnetic or acoustic waves.onThey span Sec. model the main result to of the thisabove contribution Classical channels are based the prop5. and Sec. 6presents shows their application channel agationfrom of electromagnetic or acoustic waves. They span model conclusions and presents the main result of this contribution ranges centimeters to or astronomical dimensions. Re- model before drawn Sec. 7.of agation of electromagnetic acoustic waves. They span and presentsarethe maininresult this contribution ranges the from centimeters to astronomical dimensions. Re- before conclusions are drawn in Sec. 7. cently, of signalto transduction networks in biologranges fromstudy centimeters astronomical dimensions. Re- before conclusions are drawn in Sec. 7. cently, the study of signal transduction networks in biological systems has led to interest in communication concepts 2. CHANNEL MODEL cently, the study of signal transduction networks in biological systems has ledtransfer to interest in communication concepts 2. CHANNEL MODEL that rely onhas theled of molecules instead concepts of wave ical systems to interest in communication 2. CHANNEL MODEL that rely on Possible the transfer of molecules instead of wave propagation. applications are on the nano-scale that rely on the transfer of molecules instead of wave The channel model under investigation has been proposed propagation. Possible applications are on the nano-scale TheNoel channel model under investigation has been proposed where particlePossible transport is a fast are process while senders The propagation. applications on the nano-scale by et al. (2014b) is illustrated in Fig. 1. A channel model underand investigation has been proposed where particle transport is a fast process while senders by Noel etx(n) al. at (2014b) and is illustrated inemission Fig. 1. of A and receivers electromagnetic difficult to by where particle for transport is a fast waves processare while senders bitstream the input stimulates the Noel et al. (2014b) and is illustrated in Fig. 1. A and receivers for electromagnetic waves are difficult to bitstream x(n) at the input stimulates the emission of build.receivers Furthermore, the transport of waves particles potentially and for electromagnetic areis difficult to substrate at input the source. They the are emission emitted at bitstream molecules x(n) at the stimulates of build. Furthermore, the transportlike of particles isorpotentially substrate of molecules at the t source. They are emitted at supported by natural diffusion is advection substrate build. Furthermore, theprocesses transport of particles potentially , i.e. at t = nt . A logical one multiples the bit interval molecules at the source. They are emitted at 0 0 supported by naturalthe processes like diffusion or advection multiples oftriggers the bit interval t , i.e. at t = nt0 corresponding . A logical one and does not supply like of additional energy as in multiples supported by require natural processes diffusion or advection 1)of release of molecules the bit the interval t000 , i.e. at t = nt00 corresponding . A logical one and does not require the supply of additional energy as in (x(n) = = 1) triggers the release of molecules electrical devices. and does not require the supply of additional energy as in (x(n) to an increase of the total substrate concentration s(t) by (x(n) = 1) triggers the release of molecules corresponding electrical devices. to an increase of. A thelogical total zero substrate concentration s(t) by electrical devices. an amount of s (x(n) = 0) does not trigger to an increase of the total substrate concentration by These concepts have been developed during the last decade an amount of s00 . A logical zero (x(n) = 0) does not s(t) trigger These concepts have been developed during the last decade the release any molecules. an amount of s . A logical zero (x(n) = 0) does not trigger 0 and are often have refered as molecular 0 These concepts beentodeveloped duringcommunication. the last decade the release of any molecules. and are often refered to mostly as molecular communication. release of any between molecules. Practical applications visionarycommunication. and published the and are often referedare to as molecular The environment source and receptor is characPractical applications are mostly visionary and published The environment between source receptor iss(t), characresults are based on more or less simplified models and Practical applications are mostly visionary and published terized by the concentrations of and the substrate the The environment between source and receptor iss(t), characresults areassumptions. based on more or less simplified models and by the concentrations of the substrate the idealizing Summaries the current state and can terized results are based on more or less of simplified models enzyme e(t), and the complex c(t) (see (1)). The receptor terized by the concentrations of the substrate s(t), the idealizing assumptions. Summaries of the current state can e(t), and the complexconcentration c(t) (see (1)).s(t) Theand receptor be found in Akyildiz et Summaries al. (2008); of Pierobon and state Akyildiz idealizing assumptions. the current can enzyme senses the current substrate tries enzyme e(t), and the complex c(t) (see (1)). The receptor be foundNakano in Akyildiz et al. (2008); Pierobon and Akyildiz thethe current s(t) and tries (2010); et al. et (2013). be found in Akyildiz al. (2008); Pierobon and Akyildiz senses to detect input substrate bitstream.concentration senses the current substrate concentration s(t) and tries (2010); Nakano et al. (2013). to detect the input bitstream. (2010); Nakano establish et al. (2013). Recent studies the properties of molecular com- to detect the input bitstream. Recent studies establish the properties ofparticle molecular comreceptor source e(t) munication channels from for of transfer. Recent studies establish the models properties molecular come(t) receptor source s(t) municationmodels channels from the models for particle transfer. x(n) y(n) e(t) receptor Statistical describe probability that the arrival source munication channels from models for particle transfer. s(t) x(n) y(n) Statistical models describe the probability that the arrival c(t) s(t) x(n) y(n) of a molecule is detected a receptor, e.g.that Pierobon and Statistical models describebythe probability the arrival c(t) of a molecule is detected by a receptor, e.g. Pierobon and c(t) Akyildiz (2014); Noel et al. Analogies with elecof a molecule is detected by (2014c). a receptor, e.g. Pierobon and Akyildiz (2014); Noel et al. (2014c). Analogies with electrical circuits andNoel coupled equationswith describe Akyildiz (2014); et al.differential (2014c). Analogies elec- Fig. 1. Channel model for transmission of a bitstream x(n) trical circuits and coupled differential equations describe Channel model foroutput transmission of a bitstream x(n) the dynamics molecule concentrations, see Pierobon and Fig. 1. trical circuits of and coupled differential equations describe at Channel the inputmodel to thefor y(n). of transmission a bitstream x(n) the dynamics of molecule concentrations, see Pierobon and Fig. 1. at the input to the output y(n). Akyildiz (2010); Noel et al. (2014b,a). see the dynamics of molecule concentrations, Pierobon and at the input to the output y(n). Akyildiz (2010); Noel et al. (2014b,a). Akyildiz (2010); Noelupon et al.the (2014b,a). This paper builds channel model proposed The continued release of particles would eventually flood This paper builds upon the channel model proposed The continued release of particles would eventually flood the channel with substrate and overdrive the receptor. in Noel et al.builds (2014b) andthe analyses themodel short term dy- The continued release of particles would eventually flood This paper upon channel proposed in Noel et al. (2014b) and analyses the short term dythe channel with substrate and overdrive the receptor. Therefore, mustoverdrive exist to the decrease the namics model channelsome in Noel ofetthe al. involved (2014b) components. and analyses The the resulting short term dy- the withmechanism substrate and receptor. namicstoofderive the involved components. The resulting model Therefore, some mechanism must exist to decrease the substrate concentration before the exist next to particle release. allows a simplecomponents. relation for the achievable some mechanism must decrease the namics of the involved The resultingsymbol model Therefore, allows to derive a simple relation for the achievable symbol substrate concentration before the next particle release. Possible mechanisms are before diffusion, or chemical rate considering given relation detectorfor sensitivity. concentration theadvection, next particle release. allows to derive aasimple the achievable symbol substrate rate considering a given detector sensitivity. Possible mechanisms are diffusion, advection, or chemical Possible mechanisms are diffusion, advection, or chemical rate considering a given detector sensitivity. Copyright © 2015, 2015,IFAC IFAC (International Federation of Automatic Control) 35 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © Copyright © 2015, IFAC 35 Peer review©under of International Federation of Automatic Copyright 2015,responsibility IFAC 35 Control. 10.1016/j.ifacol.2015.05.054
MATHMOD 2015 36 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
reactions. Only a chemical reaction is considered here, neglecting diffusion and advection by the assumption of a well-mixed environment. This restriction is to be removed in further research.
As an abbreviation for the introduction of different time scales in Sec. 3.3 the functions f µ(t), ν(t) = −µ(t) + ρ ν(t) + µ(t)ν(t) (14) g µ(t), ν(t) = µ(t) − Kν(t) − µ(t)ν(t) (15) are defined. Then a short notation for the amplitude-scaled equations reads as ˙ = f µ(t), ν(t) (16) Te µ(t) Ts ν(t) ˙ = g µ(t), ν(t) (17) with the scaled initial conditions c0 e0 =1− ≤1. (18) µ(0) = 1, ν(0) = ν0 = eT eT
The communication channel is characterized by a chemical reaction in the form of the following enzyme kinetics (see Noel et al. (2014b)) k1
k
S + E C →2 E + P . k−1
(1)
The substrate S and the enzyme E form a complex C at a rate k1 . In turn the complex decomposes again into S and E at a rate k−1 and into E and the product P at a rate k2 . The receptor is sensitive to the substrate S but not to the complex C or the product P. This type of enzyme reaction is well investigated especially for high substrate concentrations, see e.g. Murray (2008).
3.3 Time Scaling Time scaling does not only allow a further simplification of (16),(17). It serves also to highlight the short term behaviour and the long term behaviour of the reaction dynamics. These different kinds of behaviour appear in the extreme cases of high substrate concentration or high enzyme concentration discussed in Sec. 4.
3. MATHEMATICAL MODEL This section gives a derivation of a mathematical model in the form of coupled differential equations. Amplitude and time scaling is used for a representation in dimensionless variables. The presentation is rather concise, details can be found e.g. in Murray (2008) and references therein.
Two different time scales are considered here. They are characterized by the reference times Te and Ts or through (10) by the concentrations eT and s0 . The respective scaled time variables are defined as Ts t t τ eT = τ= and θ = where = . (19) Te Ts θ Te s0 The time-scaled variables with respect to (w.r.t.) Te are γe (τ ) = ν(t) (20) σe (τ ) = µ(t) with d dt d σe (τ ) = σe (τ ) = µ(t) = µ(t)T ˙ (21) e dτ dt dτ and similar for γe (τ ). In the same way the time-scaled variables w.r.t. Ts are introduced as σs (θ) and γs (θ). Their definitions and derivatives w.r.t. the respective time variables are compiled here as σe (τ ) = Te µ(t) ˙ (22) σe (τ ) = µ(t) γe (τ ) = Te ν(t) ˙ (23) γe (τ ) = ν(t) σs (θ) = Ts µ(t) ˙ (24) σs (θ) = µ(t) ˙ . (25) γs (θ) = Ts ν(t) γs (θ) = ν(t) Applying these relations to the amplitude-scaled equations (16),(17) gives two versions of fully amplitude- and time-scaled systems of differential equations with either the reference time Te σe (τ ) = f σe (τ ), γe (τ ) (26) Ts γ (τ ) = g σ (τ ), γ (τ ) , (27) e e Te e
3.1 Differential Equations The differential equations for the enzyme kinetics (1) are (2) s(t) ˙ = −k1 s(t)e(t) +k−1 c(t) e(t) ˙ = −k1 s(t)e(t) +k−1 c(t) +k2 c(t) (3) (4) c(t) ˙ = k1 s(t)e(t) −k−1 c(t) −k2 c(t) with the initial concentrations e(0) = e0 , c(0) = c0 , (5) s(0) = s0 , where s(0) is triggered by logical one, see Fig. 1. In a closed system, the total amount of enzyme eT is constant (6) e(t) + c(t) = e0 + c0 = eT , such that one of the equations (2) − (4) can be eliminated. E.g. (3) is removed by setting e(t) = eT − c(t). The remaining equations are +k−1 c + k1 sc , (7) s˙ = −k1 eT s c˙ = k1 eT s −(k−1 + k2 ) c − k1 sc . (8) 3.2 Amplitude Scaling The concentrations s(t) and c(t) are made dimensionless by introducing the scaled concentrations µ(t) and ν(t) as s(t) c(t) µ(t) = ν(t) = . (9) s0 eT The time constants 1 1 Ts = (10) Te = k1 eT k1 s 0 and the dimensionless constants λ = k2 Ts K =ρ+λ (11) ρ = k−1 Ts lead to the representation in scaled concentrations (12) Te µ˙ = −µ + ρ ν + µν , Ts ν˙ = µ − Kν − µν . (13) 36
or the reference time Ts Te (28) Ts σs (θ) = f σs (θ), γs (θ) (29) γs (θ) = g σs (θ), γs (θ) . For a further simplification, abbreviate the relations between the reference times as eT Te s0 1 Ts (30) = =ε, = =δ= . Te s0 Ts eT ε Then (26,27) can be repesented in two different versions (σe = σe (τ ) and γe = γe (τ )) σe = f σe , γe (31) σe = f σe , γe ε γ e = g σ e , γe γ e = δ g σ e , γe , (32)
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
37
same is true for approaches to molecular communication which are currently envisioned (e.g. Pierobon and Akyildiz (2010), Pierobon and Akyildiz (2014), Noel et al. (2014b)).
and similar for (28,29) with σs = σs (θ) and γs = γs (θ) δ σs = f σs , γs σs = ε f σs , γs (33) γ s = g σs , γs . (34) γ s = g σs , γs Note that the left hand side (l.h.s.) and the r.h.s of (31,32) are identical due to the reciprocal relationsship between ε and δ in (30). In the same way also the l.h.s and the r.h.s of (33,34) are identical. Furthermore, the two relations in (31,32) and in (33,34) represent the same system of differential equations, differing only in the scale of the time axis.
The condition for a high enzyme concentration s0 eT induces from (30) s0 =δ 1 and Te Ts . (36) eT The relation for the reference times shows that the scaled time variables from (19) describe now τ short term behaviour with reference time Te , θ long term behaviour with reference time Ts .
Nevertheless, in the extreme cases of either high substrate concentration or high enzyme concentration, these equivalent formulations give rise to different approximations for either the short term behaviour or the long term behaviour of the enzyme reaction system (1). These extreme cases are discussed in the next section.
The starting point for an approximation are the two pairs of equations from (31-34) which contain the multiplier δ σe = f σe , γe (37) γ e = δ g σe , γ e (38)
and
4. EXTREME CASES
δ σs = f σs , γs (39) γ s = g σs , γs . (40) The behaviour described by these differential equations is now investigated in detail.
The derivation of systems of differential equations above did not make any assumptions on the relation between the amounts of substrate and enzyme in the reaction. However, different application fields are characterized by either a high substrate concentration or a high enzyme concentration.
5. APPROXIMATIONS FOR HIGH ENZYME CONCENTRATION
4.1 High Substrate Concentration Molecular communication on the nano-scale has to work with low numbers of molecules, i.e. with low substrate concentrations. Of main interest is the dynamical behaviour shortly after substrate injections. It is thus desireable to separate short term and long term behaviour by suitable approximations.
The classical application is to convert a raw material (the substrate S) into a chemical product. To produce large quantities the amount of substrate S should be high. The concentration of the enzyme can be as low as possible for a reasonable yield.
Such approximations are well known from the classical case of high substrate concentration. However, it has been observed that they may not be suitable for high enzyme concentrations, see Pedersen et al. (2008). This case has also been investigated e.g. by Tzafriri (2003); Pedersen and Bersani (2010); Bersani and Dell’Acqua (2011). They show that (37-40) can be solved separately for long term behaviour (σe (τ ), γe (τ )) and short term behaviour (σs (θ),γs (θ)) by considering δ as a perturbation parameter. A first order approximation is obtained by simply setting δ = 0. This procedure is briefly outlined below.
The condition for a high substrate concentration eT s0 induces from (30) eT =ε 1 and Ts Te . (35) s0 The latter relation for the reference times shows that the scaled time variables from (19) describe τ long term behaviour with reference time Te , θ short term behaviour with reference time Ts . A suitable starting point for an approximation are those two pairs of equations from (31-34) which contain the multiplier ε σe = f σe , γe σs = ε f σs , γs ε γe = g σe , γe γs = g σs , γs . Setting ε to zero and adopting certain steady state assumptions greatly simplifies these relations and leads to the classical Michaelis-Menten model. References to original work by Henri, Michaelis and Menten, Briggs and Haldane, and reviews thereof are given e.g. in Pedersen and Bersani (2010).
5.1 Long Term Behaviour From (39,40) follows with δ = 0 0 = f σs , γs = −σs + ρ γs + σs γs γs = g σs , γs = σs − Kγs − σs γs Adding (41) and (42) gives (see (11)) γs (θ) = −λ γs (θ) with the solution γs (θ) = γs (0) e−λθ . From (41) follows now γs (θ) γs (0) = ρ λθ . σs (θ) = ρ 1 − γs (θ) e − γs (0)
4.2 High Enzyme Concentration More recently signal transduction in biological systems has been investigated. Here the substrate serves as carrier of information and its concentration may be rather low compared to the amount of enzyme in the environment. The 37
(41) (42) (43) (44)
(45)
MATHMOD 2015 38 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
5.2 Short Term Behaviour
−2
10
concentration
In the same way follows from (37,38) (46) σe = f σe , γe = −σe + ρ γe + σe γe (47) γe = 0 The latter equation states that there is no short term variation of γe (τ ) γe (τ ) = γe (0) . (48) Inserting this constant value into (46) gives σe (τ ) = − 1 − γe (0) σe (τ ) + ργe (0) (49) with the solution γe (0) 1 − e−(1−γe (0))τ +σe (0) e−(1−γe (0))τ . σe (τ ) = ρ 1 − γe (0) (50)
−5
θ→0
τ →∞
lim γs (θ) = lim γe (τ ) τ →∞
2
4 time
6
8
6. APPLICATION TO MOLECULAR COMMUNICATION 6.1 Scenario Consider the communication channel from Fig. 1 with possible injections of substrate at fixed time intervals with equidistant spacing τ0 i.e. at τ = nτ0 . Injecting the concentration s0 represents a logical 1 while no injection at a certain time interval represents a logical 0. In amplitudescaled variables the injection at τ = nτ0 is thus described by a function σ0 (n) ∈ {0, 1}.
τ →∞
θ→0
0
Fig. 2. Substrate concentration s(t) for an initial value of s0 = 10−2 eT . The time axis is labeled in multiples of the reference time Te . Solid line: approximation for the short term behaviour according to (50), dashed line: exact solution computed numerically from (7,8). .
The separate solution for short term and long term behaviour gives two different initial values γe (0) and γs (0). They are related by the requirement that the short term behaviour in the limit τ → ∞ must match the initial value θ → 0 of the long term behaviour. From (45) and (50) follows γs (0) γe (0) , lim σe (τ ) = ρ (51) lim σs (θ) = ρ θ→0 1 − γs (0) τ →∞ 1 − γe (0) and from (44) and (48) lim γs (θ) = γs (0), lim γe (τ ) = γe (0) . (52) θ→0
−4
10
10
5.3 Matching of Short and Long Term Behaviour
The matching conditions lim σs (θ) = lim σe (τ ),
−3
10
The behaviour between two injection points, e.g. for (n − 1)τ0 < τ < nτ0 , is then governed by the model for the short term behaviour as discussed in Sec. 5. The initial condition at the beginning of this interval is σ((n − 1)τ0 ). The initial condition for the next interval σ(nτ0 ) is valid at the second injection point nτ0 . It is the superposition of the substrate concentration decayed during the interval of length τ0 and a new injection σ0 (n) σ(nτ0 ) = σ (n − 1)τ0 e−τ0 + σ0 (n) . (56) Whether substrate has been injected (σ0 (n) = 1) or not (σ0 (n) = 0) can be detected by observing the substrate concentration σ(τ ) during the next interval nτ0 < τ < (n + 1)τ0 . The reliability of this detection is now investigated by considering two limiting cases.
(53)
are thus satisfied for γs (0) = γe (0) = γ(0).
5.4 Short term behaviour of the substrate concentration Since the receiver of the molecular communication channel from Fig. 1 responds to the substrate concentration s(t), only the substrate behaviour is of interest here. In addition, it is the short term behaviour which determines the dynamics after a substrate injection. Fig. 2 shows the behaviour of the substrate concentration σe (τ ) after unscaling w.r.t. amplitude and time in comparison to the numerical solution of the system of differential equations (7,8). The latter has been obtained by the MATLAB routine ode15s for stiff systems. Obviously, the first order approximation from (50) holds quite well for 0 < t < 5Te or for a decay of the substrate concentration by about two decades.
Case 0:
No injection (0) after a series of injections (1) 1 n
After substrate has been injected regularly for n < N , a certain amount of leftover substrate can be expected at n = N although there is no new injection at this point.
The solution (50) can also be written in a simpler way as γ(0) σ ˆe (τ ) = σe (τ ) − ρ =σ ˆe (0)e−τ /ˆτ (54) 1 − γ(0) with τˆ = 1/(1 − γ(0)). Note that σ ˆe (τ ) = σe (τ ) for γ(0) = 0. This special case is investigated here for brevity σe (τ ) = σe (0)e−τ . (55) The more general case σ ˆe (τ ) = σe (τ ) can be elaborated in the same way.
Case 1:
Injection (1) after a series of no injections (0): 0 n
The injection of substrate at n = N is the first one in a long time. The concentration caused by this lone injection may be less than the leftover substrate concentration in Case 0. 38
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
interval τ0 = 1 and N = 9. The substrate concentration at τ = 9 for Case 0 is obviously lower than for Case 1 (see black dots). Thus it is possible to distinguish these two cases. Note also that for Case 0 the initial values settle after a few injections at σ (0) (n) = 1/(e − 1) ≈ 0.58 (see (65)).
The problem to be solved is: Can these two cases be distinguished by observing the substrate concentration σ(τ ) during the observation interval N τ0 < τ < (N +1)τ0 ? 6.2 Detection Reliability A reliable detection requires that the substrate concentration σ (0) (τ ) for case 0 is significantly lower than the substrate concentration σ (1) (τ ) for case 1
scaled concentration
2
σ (0) (τ ) < σ (1) (τ ) for N τ0 < τ < (N + 1)τ0 . (59) From (55) follows the equivalent requirement σ (0) (N τ0 ) < σ (1) (N τ0 ) .
39
(60)
A more practical formulation of this problem has to take also the sensitivity of the measurement into account. The specific measurement method for the substrate concentration depends largely on the chemical or biological properties of the substrate. Nevertheless, all measurement methods have some concentration threshold d below which no reliable detection of concentration differences is possible. The released amount of substrate s0 should be above this threshold, i.e. d or 0< =∆<1. (61) d < s0 s0 Then the requirement (60) can be formulated as
1.5 1 0.5 0 0
1
2
3
4 5 6 scaled time τ
7
8
9
10
Fig. 3. Substrate injection according to Case 0 (dashed line) and Case 1 (solid line) for τ0 = 1 and N = 9. The initial values σ (0) (9) and σ (1) (9) are marked by black dots. A similar oscillating behaviour has been observed in (Noel et al., 2014c, Fig. 6) by simulations in a stochastical framework. Kim et al. (2014) investigated the shape of a single substrate pulse but no repeated injections.
σ (1) (N τ0 ) − σ (0) (N τ0 ) > ∆ , (62) i.e. the difference of the concentrations must be above the detectable treshold.
Inserting the results for the concentration σ (0) (N τ0 ) from (65) and for σ (1) (N τ0 ) from (66) into the requirement (62) allows to determine the maximum injection rate for a reliable detection of logical zeros and ones from measured concentrations.
6.3 Initial concentrations σ (0) (N τ0 ) and σ (1) (N τ0 ) An evaluation of condition (62) requires to determine the initial concentrations σ (0) (N τ0 ) and σ (1) (N τ0 ) for the cases 0 and 1 as introduced in Sec. 6.1.
6.4 Maximum Injection Rate The requirement (62) can be formulated in terms of σ (0) (N τ0 ) from (65) and σ (1) (N τ0 ) from (66) as 2 − e−τ0 e−τ0 = >∆. (67) 1− 1 − e−τ0 1 − e−τ0 Of practical interest is the inverse relation: What is the minimum value for τ0 such that (67) is fullfilled? Solving for τ0 gives t0 2−∆ < τ0 = , (68) ln 1−∆ Te with the bit interval t0 according to (19). Its reciprocal value f0 1 f0 = (69) t0 can be interpreted as an injection rate or symbol frequency with the physical dimension Hertz (Hz). It can be expressed by the nominal frequency fe using (10) as 1 fe = = k1 e T (70) Te which is obtained directly from the total amount of enzyme eT and the rate constant k1 of the enzyme kinetics (7,8).
Case 0: If no substrate is injected at time τ = N τ0 (0) (σ0 (N ) = 0) then the substrate concentration σ (0) (N τ0 ) is determined only by its previous value at τ = (N − 1)τ0 according to (56) σ (0) (N τ0 ) = σ (0) (N − 1)τ0 e−τ0 . (63) (0)
However, since σ0 (n) = 1 for n < N the relation (56) reads for n = N − 1 (64) σ (0) (N − 1)τ0 = σ (0) (N − 2)τ0 e−τ0 + 1 . and similar for n < N −1. Since there is a long series of ones preceeding the zero at n = N , the initial concentrations in each cycle do not change anymore, i.e. σ (0) (N − 1)τ0 = σ (0) (N − 2)τ0 . Then from (64) and (63) follows 1 e−τ0 (0) , σ (N τ ) = . (65) σ (0) (N −1)τ0 = 0 −τ 1− e 0 1− e−τ0
Case 1: Since there is a long series of no injections (0) preceeding the injection (1) at τ = N τ0 , the initial concentration at τ = (N − 1)τ0 is zero and according to (56) (1) (66) σ (1) (N τ0 ) = σ0 (N ) = 1 .
Thus the maximum injection rate or symbol frequency is bound by −1 fe 2−∆ ≈ 1.44 fe . ≤ (71) f0 < fe ln 1−∆ ln 2
Fig. 3 shows the substrate concentration as determined from the model (55) for the short time behaviour for a time 39
MATHMOD 2015 40 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
relative symbol frequency
The upper bound of the symbol frequency in dependency of the threshold ∆ is shown in Fig 4.
The final result is a simple expression for the achievable symbol frequency in terms of the chemical parameters of the molecular channel (rate constant k1 and total amount of enzyme eT ) and the detection threshold of the receiver.
1.5
Nevertheless, the presented analysis can only be regarded as the starting point of more detailed research. At first, a comparison with a second order perturbation approach is desireable. Then the assumption of a well-mixed environment has to be replaced by models for diffusion and advection and the excitation functions are better described by finite-duration pulses. Another direction of research is the intepretation of the relation between symbol frequency and receptor sensitivity in the light of rate-distortion theory.
1
0.5
0
0
0.2
0.4 0.6 scaled threshold ∆
0.8
1
Acknowledgements: I wish to thank Robert Schober for introducing me to the subject, Ian Akyildiz for visions and discussions, and Vandith Ega for programming and editing.
Fig. 4. Relative symbol rate f0 /fe according to (71). The values f0 = fe / ln 2, f0 = fe , f0 = 94 · fe are marked by a white, grey, and black dot respectively. The maximum f0 = fe / ln 2 is possible for an ideal receptor with infinite precision ∆ = 0 (white dot in Fig 4). But also receptors with nonzero thresholds achieve reasonable symbol frequencies. The symbol frequency corresponds to the nominal frequency f0 = fe for a threshold of ∆ = (e − 2)/(e − 1) ≈ 0.42 (grey dot). This means that a detection threshold of 42% of the initial injection s0 is sufficient. Even if the threshold is as high as 88% of s0 , still a symbol rate f0 = 49 · fe is achieved (black dot).
REFERENCES Akyildiz, I.F., Brunetti, F., and Bl´ azquez, C. (2008). Nanonetworks: A new communication paradigm. Computer Networks, 52(12), 2260–2279. Bersani, A.M. and Dell’Acqua, G. (2011). Asymptotic expansions in enzyme reactions with high enzyme concentrations. Mathematical Methods in the Applied Sciences, 34(16), 1954–1960. doi: 10.1002/mma.1495. Kim, N.R., Farsad, N., Chae, C.B., and Eckford, A.W. (2014). A realistic channel model for molecular communication with imperfect receivers. In Communications (ICC), 2014 IEEE International Conference on, 3987–3992. doi:10.1109/ICC.2014.6883944. Murray, J.D. (2008). Mathematical Biology I: An Introduction. Springer, New York, 3. edition. Nakano, T., Eckford, A.W., and Haraguchi, T. (2013). Molecular communication. Cambridge University Press. Noel, A., Cheung, K., and Schober, R. (2014a). Optimal receiver design for diffusive molecular communication with flow and additive noise. IEEE Transactions on NanoBioscience, PP(99), 1–1. doi:10.1109/TNB.2014.2337239. Noel, A., Cheung, K., and Schober, R. (2014b). Improving receiver performance of diffusive molecular communication with enzymes. IEEE Transactions on NanoBioscience, 13(1), 31–43. doi:10.1109/TNB.2013.2295546. Noel, A., Cheung, K., and Schober, R. (2014c). A unifying model for external noise sources and ISI in diffusive molecular communication. IEEE Journal on Selected Areas in Communications, 32(12), 2330–2343. doi:10.1109/JSAC.2014.2367693. Pedersen, M.G., Bersani, A.M., and Bersani, E. (2008). Quasi steadystate approximations in complex intracellular signal transduction networks - a word of caution. Journal of Mathematical Chemistry, 43(4), 1318–1344. doi:10.1007/s10910-007-9248-4. Pedersen, M. and Bersani, A. (2010). Introducing total substrates simplifies theoretical analysis at non-negligible enzyme concentrations: pseudo first-order kinetics and the loss of zero-order ultrasensitivity. Journal of Mathematical Biology, 60(2), 267–283. doi:10.1007/s00285-009-0267-6. Pierobon, M. and Akyildiz, I. (2010). A physical end-to-end model for molecular communication in nanonetworks. IEEE Journal on Selected Areas in Communications, 28(4), 602–611. doi: 10.1109/JSAC.2010.100509. Pierobon, M. and Akyildiz, I. (2014). A statistical-physical model of interference in diffusion-based molecular nanonetworks. IEEE Transactions on Communications, 62(6), 2085–2095. doi: 10.1109/TCOMM.2014.2314650. Tzafriri, A. (2003). Michaelis-Menten kinetics at high enzyme concentrations. Bulletin of Mathematical Biology, 65(6), 1111– 1129. doi:10.1016/S0092-8240(03)00059-4.
It is instructive to compare the respective dynamics of the substrate concentrations. The case τ0 = 1 or f0 = fe is shown in Fig. 3 and the cases for τ0 = ln 2 and τ0 = 94 in 9 ≈ 0.69 the minimum bit interval is Fig. 5. Since ln 2 ≈ 13 t0 ≈ 0.69 · Te . scaled concentration
2 1.5 1 0.5 0 0
1
2
3
4 5 6 scaled time τ
7
8
9
10
Fig. 5. Substrate injection according to Case 0 (dashed line) and Case 1 (solid line) for τ0 = ln 2 and τ0 = 94 . 7. CONCLUSIONS AND FURTHER RESEARCH This contribution has shown that molecular communication channels based on enzyme kinetics can be analyzed with the mathematical tools of classical Michaelis-Menten models. Care has to be taken to observe the correct relation between substrate and enzyme concentration and between short time and long term behaviour as these are different from the classical Michaelis-Menten case. A first order perturbation approach appears to be sufficent for modelling the substrate dynamics. Its range of accuracy (0 < t < 5 Te from Fig. 2) matches favourably with the lenghts of the bit intervals for a reliable detection (Fig. 5). 40
ScienceDirect
IFAC-PapersOnLine 48-1 (2015) 035–040
Design of a Molecular Communication Design of a Molecular Communication Design of a Molecular Communication Channel by Modelling Enzyme Kinetics Channel by Modelling Enzyme Channel by Modelling Enzyme Kinetics Kinetics
Rudolf Rabenstein ∗∗ Rudolf Rabenstein ∗∗ Rudolf Rabenstein ∗ at Erlangen-N¨ urnberg (FAU), ∗ Friedrich-Alexander-Universit¨ at Erlangen-N¨ urnberg (FAU), ∗ Friedrich-Alexander-Universit¨ ∗ D-91058 Erlangen, Cauerstr. 7, Germany (e-mail: [email protected]). Friedrich-Alexander-Universit¨ at Erlangen-N¨ urnberg (FAU), D-91058 Erlangen, Cauerstr. 7, Germany (e-mail: [email protected]). D-91058 Erlangen, Cauerstr. 7, Germany (e-mail: [email protected]). Abstract: Molecular communication channels transmit information by the release and detection Abstract: communication channels information by the release and of detection of chemical Molecular substances. Single information bitstransmit can be encoded as individual releases certain Abstract: Molecular communication channels transmit information by the release and detection of chemical substances. Single information bits can be encoded as individual releases of certain amounts of substances. molecules. Intersymbol interference between successive releases can be minimized of chemical Single information bits can be encoded as individual releases of certain amounts of molecules. Intersymbol interference between successive releases can be minimized by enforcing a decay of the substrate concentration through enzyme reactions. A model for amounts of molecules. Intersymbol interference between successive releases can be minimized by enforcing a behaviour decay of the substrate through to enzyme A model for the short term of the enzymeconcentration kinetics is employed derivereactions. the achievable symbol by enforcing a decay of the substrate concentration through enzyme reactions. A model for the short term behaviour of the enzyme kinetics is employed to derive the achievable symbol transmission rate in termsofof the the enzyme parameters of the enzyme reaction andthe theachievable sensitivitysymbol of the the short term behaviour kinetics is employed to derive transmission rate inThe terms of the relation parameters of the the enzyme reaction the sensitivity of the substrate receptor. resulting facilitates design of the and biochemical environment transmission rate inThe terms of the relation parameters of the the enzyme reaction the sensitivity of the substrate receptor. resulting facilitates design of the and biochemical environment for efficient molecular communication. substrate receptor. The resulting relation facilitates the design of the biochemical environment for efficient molecular communication. for efficient molecular communication. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Communication channels, Chemical microsensors, Models, Molecular Keywords: Communication channels, Chemical microsensors, Models, Molecular communication, Enzyme kinetics, Michaelis-Menten model Models, Keywords: Communication channels, Chemical microsensors, Molecular communication, Enzyme kinetics, Michaelis-Menten model communication, Enzyme kinetics, Michaelis-Menten model 1. INTRODUCTION Sec. 2 presents the channel model under investigation and 1. INTRODUCTION Sec. 2 presents the channel model under investigation and thechannel corresponding set of investigation differential equa1. INTRODUCTION Sec. 32 introduces presents the model under and Sec. 3 introduces the corresponding set of differential equations.3 Of special interest are the extreme discussed in Sec. introduces the corresponding set ofcases differential equations.4 Of special interest are the extreme cases discussed Sec. which leadinterest to the are useful derived in in tions. Of special the approximations extreme cases discussed in Sec. 5. 4 which lead to the useful approximations derived in Classical communication channels are based on the prop- Sec. Sec. 6 lead shows application to the above channel 4 which totheir the useful approximations derived in Classical communication channels are based on the propSec. 5. Sec. 6 shows their application to the above channel agation ofcommunication electromagnetic or acoustic waves.onThey span Sec. model the main result to of the thisabove contribution Classical channels are based the prop5. and Sec. 6presents shows their application channel agationfrom of electromagnetic or acoustic waves. They span model conclusions and presents the main result of this contribution ranges centimeters to or astronomical dimensions. Re- model before drawn Sec. 7.of agation of electromagnetic acoustic waves. They span and presentsarethe maininresult this contribution ranges the from centimeters to astronomical dimensions. Re- before conclusions are drawn in Sec. 7. cently, of signalto transduction networks in biologranges fromstudy centimeters astronomical dimensions. Re- before conclusions are drawn in Sec. 7. cently, the study of signal transduction networks in biological systems has led to interest in communication concepts 2. CHANNEL MODEL cently, the study of signal transduction networks in biological systems has ledtransfer to interest in communication concepts 2. CHANNEL MODEL that rely onhas theled of molecules instead concepts of wave ical systems to interest in communication 2. CHANNEL MODEL that rely on Possible the transfer of molecules instead of wave propagation. applications are on the nano-scale that rely on the transfer of molecules instead of wave The channel model under investigation has been proposed propagation. Possible applications are on the nano-scale TheNoel channel model under investigation has been proposed where particlePossible transport is a fast are process while senders The propagation. applications on the nano-scale by et al. (2014b) is illustrated in Fig. 1. A channel model underand investigation has been proposed where particle transport is a fast process while senders by Noel etx(n) al. at (2014b) and is illustrated inemission Fig. 1. of A and receivers electromagnetic difficult to by where particle for transport is a fast waves processare while senders bitstream the input stimulates the Noel et al. (2014b) and is illustrated in Fig. 1. A and receivers for electromagnetic waves are difficult to bitstream x(n) at the input stimulates the emission of build.receivers Furthermore, the transport of waves particles potentially and for electromagnetic areis difficult to substrate at input the source. They the are emission emitted at bitstream molecules x(n) at the stimulates of build. Furthermore, the transportlike of particles isorpotentially substrate of molecules at the t source. They are emitted at supported by natural diffusion is advection substrate build. Furthermore, theprocesses transport of particles potentially , i.e. at t = nt . A logical one multiples the bit interval molecules at the source. They are emitted at 0 0 supported by naturalthe processes like diffusion or advection multiples oftriggers the bit interval t , i.e. at t = nt0 corresponding . A logical one and does not supply like of additional energy as in multiples supported by require natural processes diffusion or advection 1)of release of molecules the bit the interval t000 , i.e. at t = nt00 corresponding . A logical one and does not require the supply of additional energy as in (x(n) = = 1) triggers the release of molecules electrical devices. and does not require the supply of additional energy as in (x(n) to an increase of the total substrate concentration s(t) by (x(n) = 1) triggers the release of molecules corresponding electrical devices. to an increase of. A thelogical total zero substrate concentration s(t) by electrical devices. an amount of s (x(n) = 0) does not trigger to an increase of the total substrate concentration by These concepts have been developed during the last decade an amount of s00 . A logical zero (x(n) = 0) does not s(t) trigger These concepts have been developed during the last decade the release any molecules. an amount of s . A logical zero (x(n) = 0) does not trigger 0 and are often have refered as molecular 0 These concepts beentodeveloped duringcommunication. the last decade the release of any molecules. and are often refered to mostly as molecular communication. release of any between molecules. Practical applications visionarycommunication. and published the and are often referedare to as molecular The environment source and receptor is characPractical applications are mostly visionary and published The environment between source receptor iss(t), characresults are based on more or less simplified models and Practical applications are mostly visionary and published terized by the concentrations of and the substrate the The environment between source and receptor iss(t), characresults areassumptions. based on more or less simplified models and by the concentrations of the substrate the idealizing Summaries the current state and can terized results are based on more or less of simplified models enzyme e(t), and the complex c(t) (see (1)). The receptor terized by the concentrations of the substrate s(t), the idealizing assumptions. Summaries of the current state can e(t), and the complexconcentration c(t) (see (1)).s(t) Theand receptor be found in Akyildiz et Summaries al. (2008); of Pierobon and state Akyildiz idealizing assumptions. the current can enzyme senses the current substrate tries enzyme e(t), and the complex c(t) (see (1)). The receptor be foundNakano in Akyildiz et al. (2008); Pierobon and Akyildiz thethe current s(t) and tries (2010); et al. et (2013). be found in Akyildiz al. (2008); Pierobon and Akyildiz senses to detect input substrate bitstream.concentration senses the current substrate concentration s(t) and tries (2010); Nakano et al. (2013). to detect the input bitstream. (2010); Nakano establish et al. (2013). Recent studies the properties of molecular com- to detect the input bitstream. Recent studies establish the properties ofparticle molecular comreceptor source e(t) munication channels from for of transfer. Recent studies establish the models properties molecular come(t) receptor source s(t) municationmodels channels from the models for particle transfer. x(n) y(n) e(t) receptor Statistical describe probability that the arrival source munication channels from models for particle transfer. s(t) x(n) y(n) Statistical models describe the probability that the arrival c(t) s(t) x(n) y(n) of a molecule is detected a receptor, e.g.that Pierobon and Statistical models describebythe probability the arrival c(t) of a molecule is detected by a receptor, e.g. Pierobon and c(t) Akyildiz (2014); Noel et al. Analogies with elecof a molecule is detected by (2014c). a receptor, e.g. Pierobon and Akyildiz (2014); Noel et al. (2014c). Analogies with electrical circuits andNoel coupled equationswith describe Akyildiz (2014); et al.differential (2014c). Analogies elec- Fig. 1. Channel model for transmission of a bitstream x(n) trical circuits and coupled differential equations describe Channel model foroutput transmission of a bitstream x(n) the dynamics molecule concentrations, see Pierobon and Fig. 1. trical circuits of and coupled differential equations describe at Channel the inputmodel to thefor y(n). of transmission a bitstream x(n) the dynamics of molecule concentrations, see Pierobon and Fig. 1. at the input to the output y(n). Akyildiz (2010); Noel et al. (2014b,a). see the dynamics of molecule concentrations, Pierobon and at the input to the output y(n). Akyildiz (2010); Noel et al. (2014b,a). Akyildiz (2010); Noelupon et al.the (2014b,a). This paper builds channel model proposed The continued release of particles would eventually flood This paper builds upon the channel model proposed The continued release of particles would eventually flood the channel with substrate and overdrive the receptor. in Noel et al.builds (2014b) andthe analyses themodel short term dy- The continued release of particles would eventually flood This paper upon channel proposed in Noel et al. (2014b) and analyses the short term dythe channel with substrate and overdrive the receptor. Therefore, mustoverdrive exist to the decrease the namics model channelsome in Noel ofetthe al. involved (2014b) components. and analyses The the resulting short term dy- the withmechanism substrate and receptor. namicstoofderive the involved components. The resulting model Therefore, some mechanism must exist to decrease the substrate concentration before the exist next to particle release. allows a simplecomponents. relation for the achievable some mechanism must decrease the namics of the involved The resultingsymbol model Therefore, allows to derive a simple relation for the achievable symbol substrate concentration before the next particle release. Possible mechanisms are before diffusion, or chemical rate considering given relation detectorfor sensitivity. concentration theadvection, next particle release. allows to derive aasimple the achievable symbol substrate rate considering a given detector sensitivity. Possible mechanisms are diffusion, advection, or chemical Possible mechanisms are diffusion, advection, or chemical rate considering a given detector sensitivity. Copyright © 2015, 2015,IFAC IFAC (International Federation of Automatic Control) 35 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © Copyright © 2015, IFAC 35 Peer review©under of International Federation of Automatic Copyright 2015,responsibility IFAC 35 Control. 10.1016/j.ifacol.2015.05.054
MATHMOD 2015 36 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
reactions. Only a chemical reaction is considered here, neglecting diffusion and advection by the assumption of a well-mixed environment. This restriction is to be removed in further research.
As an abbreviation for the introduction of different time scales in Sec. 3.3 the functions f µ(t), ν(t) = −µ(t) + ρ ν(t) + µ(t)ν(t) (14) g µ(t), ν(t) = µ(t) − Kν(t) − µ(t)ν(t) (15) are defined. Then a short notation for the amplitude-scaled equations reads as ˙ = f µ(t), ν(t) (16) Te µ(t) Ts ν(t) ˙ = g µ(t), ν(t) (17) with the scaled initial conditions c0 e0 =1− ≤1. (18) µ(0) = 1, ν(0) = ν0 = eT eT
The communication channel is characterized by a chemical reaction in the form of the following enzyme kinetics (see Noel et al. (2014b)) k1
k
S + E C →2 E + P . k−1
(1)
The substrate S and the enzyme E form a complex C at a rate k1 . In turn the complex decomposes again into S and E at a rate k−1 and into E and the product P at a rate k2 . The receptor is sensitive to the substrate S but not to the complex C or the product P. This type of enzyme reaction is well investigated especially for high substrate concentrations, see e.g. Murray (2008).
3.3 Time Scaling Time scaling does not only allow a further simplification of (16),(17). It serves also to highlight the short term behaviour and the long term behaviour of the reaction dynamics. These different kinds of behaviour appear in the extreme cases of high substrate concentration or high enzyme concentration discussed in Sec. 4.
3. MATHEMATICAL MODEL This section gives a derivation of a mathematical model in the form of coupled differential equations. Amplitude and time scaling is used for a representation in dimensionless variables. The presentation is rather concise, details can be found e.g. in Murray (2008) and references therein.
Two different time scales are considered here. They are characterized by the reference times Te and Ts or through (10) by the concentrations eT and s0 . The respective scaled time variables are defined as Ts t t τ eT = τ= and θ = where = . (19) Te Ts θ Te s0 The time-scaled variables with respect to (w.r.t.) Te are γe (τ ) = ν(t) (20) σe (τ ) = µ(t) with d dt d σe (τ ) = σe (τ ) = µ(t) = µ(t)T ˙ (21) e dτ dt dτ and similar for γe (τ ). In the same way the time-scaled variables w.r.t. Ts are introduced as σs (θ) and γs (θ). Their definitions and derivatives w.r.t. the respective time variables are compiled here as σe (τ ) = Te µ(t) ˙ (22) σe (τ ) = µ(t) γe (τ ) = Te ν(t) ˙ (23) γe (τ ) = ν(t) σs (θ) = Ts µ(t) ˙ (24) σs (θ) = µ(t) ˙ . (25) γs (θ) = Ts ν(t) γs (θ) = ν(t) Applying these relations to the amplitude-scaled equations (16),(17) gives two versions of fully amplitude- and time-scaled systems of differential equations with either the reference time Te σe (τ ) = f σe (τ ), γe (τ ) (26) Ts γ (τ ) = g σ (τ ), γ (τ ) , (27) e e Te e
3.1 Differential Equations The differential equations for the enzyme kinetics (1) are (2) s(t) ˙ = −k1 s(t)e(t) +k−1 c(t) e(t) ˙ = −k1 s(t)e(t) +k−1 c(t) +k2 c(t) (3) (4) c(t) ˙ = k1 s(t)e(t) −k−1 c(t) −k2 c(t) with the initial concentrations e(0) = e0 , c(0) = c0 , (5) s(0) = s0 , where s(0) is triggered by logical one, see Fig. 1. In a closed system, the total amount of enzyme eT is constant (6) e(t) + c(t) = e0 + c0 = eT , such that one of the equations (2) − (4) can be eliminated. E.g. (3) is removed by setting e(t) = eT − c(t). The remaining equations are +k−1 c + k1 sc , (7) s˙ = −k1 eT s c˙ = k1 eT s −(k−1 + k2 ) c − k1 sc . (8) 3.2 Amplitude Scaling The concentrations s(t) and c(t) are made dimensionless by introducing the scaled concentrations µ(t) and ν(t) as s(t) c(t) µ(t) = ν(t) = . (9) s0 eT The time constants 1 1 Ts = (10) Te = k1 eT k1 s 0 and the dimensionless constants λ = k2 Ts K =ρ+λ (11) ρ = k−1 Ts lead to the representation in scaled concentrations (12) Te µ˙ = −µ + ρ ν + µν , Ts ν˙ = µ − Kν − µν . (13) 36
or the reference time Ts Te (28) Ts σs (θ) = f σs (θ), γs (θ) (29) γs (θ) = g σs (θ), γs (θ) . For a further simplification, abbreviate the relations between the reference times as eT Te s0 1 Ts (30) = =ε, = =δ= . Te s0 Ts eT ε Then (26,27) can be repesented in two different versions (σe = σe (τ ) and γe = γe (τ )) σe = f σe , γe (31) σe = f σe , γe ε γ e = g σ e , γe γ e = δ g σ e , γe , (32)
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
37
same is true for approaches to molecular communication which are currently envisioned (e.g. Pierobon and Akyildiz (2010), Pierobon and Akyildiz (2014), Noel et al. (2014b)).
and similar for (28,29) with σs = σs (θ) and γs = γs (θ) δ σs = f σs , γs σs = ε f σs , γs (33) γ s = g σs , γs . (34) γ s = g σs , γs Note that the left hand side (l.h.s.) and the r.h.s of (31,32) are identical due to the reciprocal relationsship between ε and δ in (30). In the same way also the l.h.s and the r.h.s of (33,34) are identical. Furthermore, the two relations in (31,32) and in (33,34) represent the same system of differential equations, differing only in the scale of the time axis.
The condition for a high enzyme concentration s0 eT induces from (30) s0 =δ 1 and Te Ts . (36) eT The relation for the reference times shows that the scaled time variables from (19) describe now τ short term behaviour with reference time Te , θ long term behaviour with reference time Ts .
Nevertheless, in the extreme cases of either high substrate concentration or high enzyme concentration, these equivalent formulations give rise to different approximations for either the short term behaviour or the long term behaviour of the enzyme reaction system (1). These extreme cases are discussed in the next section.
The starting point for an approximation are the two pairs of equations from (31-34) which contain the multiplier δ σe = f σe , γe (37) γ e = δ g σe , γ e (38)
and
4. EXTREME CASES
δ σs = f σs , γs (39) γ s = g σs , γs . (40) The behaviour described by these differential equations is now investigated in detail.
The derivation of systems of differential equations above did not make any assumptions on the relation between the amounts of substrate and enzyme in the reaction. However, different application fields are characterized by either a high substrate concentration or a high enzyme concentration.
5. APPROXIMATIONS FOR HIGH ENZYME CONCENTRATION
4.1 High Substrate Concentration Molecular communication on the nano-scale has to work with low numbers of molecules, i.e. with low substrate concentrations. Of main interest is the dynamical behaviour shortly after substrate injections. It is thus desireable to separate short term and long term behaviour by suitable approximations.
The classical application is to convert a raw material (the substrate S) into a chemical product. To produce large quantities the amount of substrate S should be high. The concentration of the enzyme can be as low as possible for a reasonable yield.
Such approximations are well known from the classical case of high substrate concentration. However, it has been observed that they may not be suitable for high enzyme concentrations, see Pedersen et al. (2008). This case has also been investigated e.g. by Tzafriri (2003); Pedersen and Bersani (2010); Bersani and Dell’Acqua (2011). They show that (37-40) can be solved separately for long term behaviour (σe (τ ), γe (τ )) and short term behaviour (σs (θ),γs (θ)) by considering δ as a perturbation parameter. A first order approximation is obtained by simply setting δ = 0. This procedure is briefly outlined below.
The condition for a high substrate concentration eT s0 induces from (30) eT =ε 1 and Ts Te . (35) s0 The latter relation for the reference times shows that the scaled time variables from (19) describe τ long term behaviour with reference time Te , θ short term behaviour with reference time Ts . A suitable starting point for an approximation are those two pairs of equations from (31-34) which contain the multiplier ε σe = f σe , γe σs = ε f σs , γs ε γe = g σe , γe γs = g σs , γs . Setting ε to zero and adopting certain steady state assumptions greatly simplifies these relations and leads to the classical Michaelis-Menten model. References to original work by Henri, Michaelis and Menten, Briggs and Haldane, and reviews thereof are given e.g. in Pedersen and Bersani (2010).
5.1 Long Term Behaviour From (39,40) follows with δ = 0 0 = f σs , γs = −σs + ρ γs + σs γs γs = g σs , γs = σs − Kγs − σs γs Adding (41) and (42) gives (see (11)) γs (θ) = −λ γs (θ) with the solution γs (θ) = γs (0) e−λθ . From (41) follows now γs (θ) γs (0) = ρ λθ . σs (θ) = ρ 1 − γs (θ) e − γs (0)
4.2 High Enzyme Concentration More recently signal transduction in biological systems has been investigated. Here the substrate serves as carrier of information and its concentration may be rather low compared to the amount of enzyme in the environment. The 37
(41) (42) (43) (44)
(45)
MATHMOD 2015 38 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
5.2 Short Term Behaviour
−2
10
concentration
In the same way follows from (37,38) (46) σe = f σe , γe = −σe + ρ γe + σe γe (47) γe = 0 The latter equation states that there is no short term variation of γe (τ ) γe (τ ) = γe (0) . (48) Inserting this constant value into (46) gives σe (τ ) = − 1 − γe (0) σe (τ ) + ργe (0) (49) with the solution γe (0) 1 − e−(1−γe (0))τ +σe (0) e−(1−γe (0))τ . σe (τ ) = ρ 1 − γe (0) (50)
−5
θ→0
τ →∞
lim γs (θ) = lim γe (τ ) τ →∞
2
4 time
6
8
6. APPLICATION TO MOLECULAR COMMUNICATION 6.1 Scenario Consider the communication channel from Fig. 1 with possible injections of substrate at fixed time intervals with equidistant spacing τ0 i.e. at τ = nτ0 . Injecting the concentration s0 represents a logical 1 while no injection at a certain time interval represents a logical 0. In amplitudescaled variables the injection at τ = nτ0 is thus described by a function σ0 (n) ∈ {0, 1}.
τ →∞
θ→0
0
Fig. 2. Substrate concentration s(t) for an initial value of s0 = 10−2 eT . The time axis is labeled in multiples of the reference time Te . Solid line: approximation for the short term behaviour according to (50), dashed line: exact solution computed numerically from (7,8). .
The separate solution for short term and long term behaviour gives two different initial values γe (0) and γs (0). They are related by the requirement that the short term behaviour in the limit τ → ∞ must match the initial value θ → 0 of the long term behaviour. From (45) and (50) follows γs (0) γe (0) , lim σe (τ ) = ρ (51) lim σs (θ) = ρ θ→0 1 − γs (0) τ →∞ 1 − γe (0) and from (44) and (48) lim γs (θ) = γs (0), lim γe (τ ) = γe (0) . (52) θ→0
−4
10
10
5.3 Matching of Short and Long Term Behaviour
The matching conditions lim σs (θ) = lim σe (τ ),
−3
10
The behaviour between two injection points, e.g. for (n − 1)τ0 < τ < nτ0 , is then governed by the model for the short term behaviour as discussed in Sec. 5. The initial condition at the beginning of this interval is σ((n − 1)τ0 ). The initial condition for the next interval σ(nτ0 ) is valid at the second injection point nτ0 . It is the superposition of the substrate concentration decayed during the interval of length τ0 and a new injection σ0 (n) σ(nτ0 ) = σ (n − 1)τ0 e−τ0 + σ0 (n) . (56) Whether substrate has been injected (σ0 (n) = 1) or not (σ0 (n) = 0) can be detected by observing the substrate concentration σ(τ ) during the next interval nτ0 < τ < (n + 1)τ0 . The reliability of this detection is now investigated by considering two limiting cases.
(53)
are thus satisfied for γs (0) = γe (0) = γ(0).
5.4 Short term behaviour of the substrate concentration Since the receiver of the molecular communication channel from Fig. 1 responds to the substrate concentration s(t), only the substrate behaviour is of interest here. In addition, it is the short term behaviour which determines the dynamics after a substrate injection. Fig. 2 shows the behaviour of the substrate concentration σe (τ ) after unscaling w.r.t. amplitude and time in comparison to the numerical solution of the system of differential equations (7,8). The latter has been obtained by the MATLAB routine ode15s for stiff systems. Obviously, the first order approximation from (50) holds quite well for 0 < t < 5Te or for a decay of the substrate concentration by about two decades.
Case 0:
No injection (0) after a series of injections (1) 1 n
After substrate has been injected regularly for n < N , a certain amount of leftover substrate can be expected at n = N although there is no new injection at this point.
The solution (50) can also be written in a simpler way as γ(0) σ ˆe (τ ) = σe (τ ) − ρ =σ ˆe (0)e−τ /ˆτ (54) 1 − γ(0) with τˆ = 1/(1 − γ(0)). Note that σ ˆe (τ ) = σe (τ ) for γ(0) = 0. This special case is investigated here for brevity σe (τ ) = σe (0)e−τ . (55) The more general case σ ˆe (τ ) = σe (τ ) can be elaborated in the same way.
Case 1:
Injection (1) after a series of no injections (0): 0 n
The injection of substrate at n = N is the first one in a long time. The concentration caused by this lone injection may be less than the leftover substrate concentration in Case 0. 38
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
interval τ0 = 1 and N = 9. The substrate concentration at τ = 9 for Case 0 is obviously lower than for Case 1 (see black dots). Thus it is possible to distinguish these two cases. Note also that for Case 0 the initial values settle after a few injections at σ (0) (n) = 1/(e − 1) ≈ 0.58 (see (65)).
The problem to be solved is: Can these two cases be distinguished by observing the substrate concentration σ(τ ) during the observation interval N τ0 < τ < (N +1)τ0 ? 6.2 Detection Reliability A reliable detection requires that the substrate concentration σ (0) (τ ) for case 0 is significantly lower than the substrate concentration σ (1) (τ ) for case 1
scaled concentration
2
σ (0) (τ ) < σ (1) (τ ) for N τ0 < τ < (N + 1)τ0 . (59) From (55) follows the equivalent requirement σ (0) (N τ0 ) < σ (1) (N τ0 ) .
39
(60)
A more practical formulation of this problem has to take also the sensitivity of the measurement into account. The specific measurement method for the substrate concentration depends largely on the chemical or biological properties of the substrate. Nevertheless, all measurement methods have some concentration threshold d below which no reliable detection of concentration differences is possible. The released amount of substrate s0 should be above this threshold, i.e. d or 0< =∆<1. (61) d < s0 s0 Then the requirement (60) can be formulated as
1.5 1 0.5 0 0
1
2
3
4 5 6 scaled time τ
7
8
9
10
Fig. 3. Substrate injection according to Case 0 (dashed line) and Case 1 (solid line) for τ0 = 1 and N = 9. The initial values σ (0) (9) and σ (1) (9) are marked by black dots. A similar oscillating behaviour has been observed in (Noel et al., 2014c, Fig. 6) by simulations in a stochastical framework. Kim et al. (2014) investigated the shape of a single substrate pulse but no repeated injections.
σ (1) (N τ0 ) − σ (0) (N τ0 ) > ∆ , (62) i.e. the difference of the concentrations must be above the detectable treshold.
Inserting the results for the concentration σ (0) (N τ0 ) from (65) and for σ (1) (N τ0 ) from (66) into the requirement (62) allows to determine the maximum injection rate for a reliable detection of logical zeros and ones from measured concentrations.
6.3 Initial concentrations σ (0) (N τ0 ) and σ (1) (N τ0 ) An evaluation of condition (62) requires to determine the initial concentrations σ (0) (N τ0 ) and σ (1) (N τ0 ) for the cases 0 and 1 as introduced in Sec. 6.1.
6.4 Maximum Injection Rate The requirement (62) can be formulated in terms of σ (0) (N τ0 ) from (65) and σ (1) (N τ0 ) from (66) as 2 − e−τ0 e−τ0 = >∆. (67) 1− 1 − e−τ0 1 − e−τ0 Of practical interest is the inverse relation: What is the minimum value for τ0 such that (67) is fullfilled? Solving for τ0 gives t0 2−∆ < τ0 = , (68) ln 1−∆ Te with the bit interval t0 according to (19). Its reciprocal value f0 1 f0 = (69) t0 can be interpreted as an injection rate or symbol frequency with the physical dimension Hertz (Hz). It can be expressed by the nominal frequency fe using (10) as 1 fe = = k1 e T (70) Te which is obtained directly from the total amount of enzyme eT and the rate constant k1 of the enzyme kinetics (7,8).
Case 0: If no substrate is injected at time τ = N τ0 (0) (σ0 (N ) = 0) then the substrate concentration σ (0) (N τ0 ) is determined only by its previous value at τ = (N − 1)τ0 according to (56) σ (0) (N τ0 ) = σ (0) (N − 1)τ0 e−τ0 . (63) (0)
However, since σ0 (n) = 1 for n < N the relation (56) reads for n = N − 1 (64) σ (0) (N − 1)τ0 = σ (0) (N − 2)τ0 e−τ0 + 1 . and similar for n < N −1. Since there is a long series of ones preceeding the zero at n = N , the initial concentrations in each cycle do not change anymore, i.e. σ (0) (N − 1)τ0 = σ (0) (N − 2)τ0 . Then from (64) and (63) follows 1 e−τ0 (0) , σ (N τ ) = . (65) σ (0) (N −1)τ0 = 0 −τ 1− e 0 1− e−τ0
Case 1: Since there is a long series of no injections (0) preceeding the injection (1) at τ = N τ0 , the initial concentration at τ = (N − 1)τ0 is zero and according to (56) (1) (66) σ (1) (N τ0 ) = σ0 (N ) = 1 .
Thus the maximum injection rate or symbol frequency is bound by −1 fe 2−∆ ≈ 1.44 fe . ≤ (71) f0 < fe ln 1−∆ ln 2
Fig. 3 shows the substrate concentration as determined from the model (55) for the short time behaviour for a time 39
MATHMOD 2015 40 February 18 - 20, 2015. Vienna, Austria Rudolf Rabenstein et al. / IFAC-PapersOnLine 48-1 (2015) 035–040
relative symbol frequency
The upper bound of the symbol frequency in dependency of the threshold ∆ is shown in Fig 4.
The final result is a simple expression for the achievable symbol frequency in terms of the chemical parameters of the molecular channel (rate constant k1 and total amount of enzyme eT ) and the detection threshold of the receiver.
1.5
Nevertheless, the presented analysis can only be regarded as the starting point of more detailed research. At first, a comparison with a second order perturbation approach is desireable. Then the assumption of a well-mixed environment has to be replaced by models for diffusion and advection and the excitation functions are better described by finite-duration pulses. Another direction of research is the intepretation of the relation between symbol frequency and receptor sensitivity in the light of rate-distortion theory.
1
0.5
0
0
0.2
0.4 0.6 scaled threshold ∆
0.8
1
Acknowledgements: I wish to thank Robert Schober for introducing me to the subject, Ian Akyildiz for visions and discussions, and Vandith Ega for programming and editing.
Fig. 4. Relative symbol rate f0 /fe according to (71). The values f0 = fe / ln 2, f0 = fe , f0 = 94 · fe are marked by a white, grey, and black dot respectively. The maximum f0 = fe / ln 2 is possible for an ideal receptor with infinite precision ∆ = 0 (white dot in Fig 4). But also receptors with nonzero thresholds achieve reasonable symbol frequencies. The symbol frequency corresponds to the nominal frequency f0 = fe for a threshold of ∆ = (e − 2)/(e − 1) ≈ 0.42 (grey dot). This means that a detection threshold of 42% of the initial injection s0 is sufficient. Even if the threshold is as high as 88% of s0 , still a symbol rate f0 = 49 · fe is achieved (black dot).
REFERENCES Akyildiz, I.F., Brunetti, F., and Bl´ azquez, C. (2008). Nanonetworks: A new communication paradigm. Computer Networks, 52(12), 2260–2279. Bersani, A.M. and Dell’Acqua, G. (2011). Asymptotic expansions in enzyme reactions with high enzyme concentrations. Mathematical Methods in the Applied Sciences, 34(16), 1954–1960. doi: 10.1002/mma.1495. Kim, N.R., Farsad, N., Chae, C.B., and Eckford, A.W. (2014). A realistic channel model for molecular communication with imperfect receivers. In Communications (ICC), 2014 IEEE International Conference on, 3987–3992. doi:10.1109/ICC.2014.6883944. Murray, J.D. (2008). Mathematical Biology I: An Introduction. Springer, New York, 3. edition. Nakano, T., Eckford, A.W., and Haraguchi, T. (2013). Molecular communication. Cambridge University Press. Noel, A., Cheung, K., and Schober, R. (2014a). Optimal receiver design for diffusive molecular communication with flow and additive noise. IEEE Transactions on NanoBioscience, PP(99), 1–1. doi:10.1109/TNB.2014.2337239. Noel, A., Cheung, K., and Schober, R. (2014b). Improving receiver performance of diffusive molecular communication with enzymes. IEEE Transactions on NanoBioscience, 13(1), 31–43. doi:10.1109/TNB.2013.2295546. Noel, A., Cheung, K., and Schober, R. (2014c). A unifying model for external noise sources and ISI in diffusive molecular communication. IEEE Journal on Selected Areas in Communications, 32(12), 2330–2343. doi:10.1109/JSAC.2014.2367693. Pedersen, M.G., Bersani, A.M., and Bersani, E. (2008). Quasi steadystate approximations in complex intracellular signal transduction networks - a word of caution. Journal of Mathematical Chemistry, 43(4), 1318–1344. doi:10.1007/s10910-007-9248-4. Pedersen, M. and Bersani, A. (2010). Introducing total substrates simplifies theoretical analysis at non-negligible enzyme concentrations: pseudo first-order kinetics and the loss of zero-order ultrasensitivity. Journal of Mathematical Biology, 60(2), 267–283. doi:10.1007/s00285-009-0267-6. Pierobon, M. and Akyildiz, I. (2010). A physical end-to-end model for molecular communication in nanonetworks. IEEE Journal on Selected Areas in Communications, 28(4), 602–611. doi: 10.1109/JSAC.2010.100509. Pierobon, M. and Akyildiz, I. (2014). A statistical-physical model of interference in diffusion-based molecular nanonetworks. IEEE Transactions on Communications, 62(6), 2085–2095. doi: 10.1109/TCOMM.2014.2314650. Tzafriri, A. (2003). Michaelis-Menten kinetics at high enzyme concentrations. Bulletin of Mathematical Biology, 65(6), 1111– 1129. doi:10.1016/S0092-8240(03)00059-4.
It is instructive to compare the respective dynamics of the substrate concentrations. The case τ0 = 1 or f0 = fe is shown in Fig. 3 and the cases for τ0 = ln 2 and τ0 = 94 in 9 ≈ 0.69 the minimum bit interval is Fig. 5. Since ln 2 ≈ 13 t0 ≈ 0.69 · Te . scaled concentration
2 1.5 1 0.5 0 0
1
2
3
4 5 6 scaled time τ
7
8
9
10
Fig. 5. Substrate injection according to Case 0 (dashed line) and Case 1 (solid line) for τ0 = ln 2 and τ0 = 94 . 7. CONCLUSIONS AND FURTHER RESEARCH This contribution has shown that molecular communication channels based on enzyme kinetics can be analyzed with the mathematical tools of classical Michaelis-Menten models. Care has to be taken to observe the correct relation between substrate and enzyme concentration and between short time and long term behaviour as these are different from the classical Michaelis-Menten case. A first order perturbation approach appears to be sufficent for modelling the substrate dynamics. Its range of accuracy (0 < t < 5 Te from Fig. 2) matches favourably with the lenghts of the bit intervals for a reliable detection (Fig. 5). 40