Diffusion-controlled interface kinetics-inclusive system-theoretic propagation models for molecular communication systems

Abstract

Inspired by biological systems, molecular communication has been proposed as a new communication paradigm that uses biochemical signals to transfer information from one nano device to another over a short distance. The biochemical nature of the information transfer process implies that for molecular communication purposes, the development of molecular channel models should take into consideration diffusion phenomenon as well as the physical/biochemical kinetic possibilities of the process. The physical and biochemical kinetics arise at the interfaces between the diffusion channel and the transmitter/receiver units. These interfaces are herein termed molecular antennas. In this paper, we present the deterministic propagation model of the molecular communication between an immobilized nanotransmitter and nanoreceiver, where the emission and reception kinetics are taken into consideration. Specifically, we derived closed-form system-theoretic models and expressions for configurations that represent different communication systems based on the type of molecular antennas used. The antennas considered are the nanopores at the transmitter and the surface receptor proteins/enzymes at the receiver. The developed models are simulated to show the influence of parameters such as the receiver radius, surface receptor protein/enzyme concentration, and various reaction rate constants. Results show that the effective receiver surface area and the rate constants are important to the system’s output performance. Assuming high rate of catalysis, the analysis of the frequency behavior of the developed propagation channels in the form of transfer functions shows significant difference introduce by the inclusion of the molecular antennas into the diffusion-only model. It is also shown that for t > > 0 and with the information molecules’ concentration greater than the Michaelis-Menten kinetic constant of the systems, the inclusion of surface receptors proteins and enzymes in the models makes the system act like a band-stop filter over an infinite frequency range.

1 Introduction

There has been increasing interest in the field of nanotechnology in the recent times, a trend that is evident from the volume of accessible literature on this subject. Nanotechnology emerged as a promising technology with contributions from diverse science fields and with potential applications in wide range of areas. On the basic, this technology dwells on novel concepts, techniques, tool, and devices that can be employed to probe and manipulate materials at the particulate and molecular levels. These novel concepts, techniques, tools, and devices are envisioned to be used to substantially improve existing technologies as well as introduce entirely new technological innovations.

An aspect of nanotechnology that is receiving much attention recently is the subject of nanocommunication. Nanocommunication is a new research area that aims at realizing communications between nano-scale computing devices or nanomachines [1]. Essentially, the communications among nanomachines can be realized by means of molecular, electromagnetic, or quantum information exchange. This paper focuses on nanocommunication by means of molecular communication. Molecular communication (MC) is realized through the transmission and reception of biochemical information encoded in the concentration and type of molecules. The emergence of MC has found potential applications in fields that include bioengineering [2, 3], nanomedicine [4–10], and environmental safety [11]. Readers shall refer to [1, 12–16] and the references therein for reviews of recent work in MC as well as some potential applications.

Just like in every communication among entities, the reliability of information exchange between nanomachines depends on the characteristics of the communication medium. The ability to integrate the knowledge of the medium characteristics into the design and analysis of the nanomachines is dependent on the appropriateness of the model that is used to mimic the propagation behavior. Popular molecular propagation mechanisms include the use of molecular motors [17–19], bacteria [20, 21], and molecular diffusion [22, 23]. The focus of this paper is on the molecular diffusion technique for information exchange among nanomachines. In molecular diffusion technique, information in the form of molecular concentration or type is broadcasted and received by a nanoreceiver, which processes and decodes the information based on predefined biochemical algorithms.

Over the past years, different diffusion models have been developed. Some of these models include those based on the theories of turbulent diffusion [24], electro-diffusion [25], and Fick’s diffusion [26, 27]. In this work, we are concern with Fick’s diffusion technique because of its generality. In the context of molecular communication, limited research efforts have addressed the problem of modeling the MC channel. In [28], Eckford presented a model of molecular communication based on the Brownian motion. And in [29, 30], MC channel models for the scenarios where the diffusion medium is in motion are presented. However, the models presented in [28–30] do not incorporate any biochemical reaction possibilities. It should be noted that the aim of molecular communication is to develop communication systems using underlying biological communication components [31]. For instance, the reception mechanism and information decoding of an MC system can be modeled using some known cell signaling pathways. An example is the case of the popular ligand-receptor binding action on the cell membrane [32]. The mechanism can be employed as the reception mechanism in nanoreceivers and modeled as a biochemical reaction. Molecular receiver reception models that consider the ligand-receptor binding mechanism can be found in [33–37]. Hence, while one can compute the concentration of a molecule at a given space around a receiver, the mechanism of coupling the molecular information into the receiver biological circuit will definitely define how much information it actually transferred. Therefore, the reaction rate of say the ligand-receptor binding action and other factors of influence have to be incorporated into the channel model. Pierobon and Akyldiz [38] developed an interesting mathematical framework that aims at an interpretation of the diffusion-based MC with consideration to the use of ligand-receptor binding action at the receiver. However, this model does not include the actual biochemical reaction kinetics. The model in [38] is used in [34] to model and analyze the noise in a typical MC system. A molecular diffusion model that considers the ligand-receptor binding kinetics in a moving propagation medium is presented in [39]; however, the biochemical kinetics of the reception process is not accounted for in the overall model.

In this paper, we present a deterministic system-theoretic molecular communication diffusion model that incorporates the physical/biochemical kinetics of the interfacing units between communicating nanomachines. These interfacing units such as the surface receptors protein/enzymes and nanopores are termed molecular antennas. The term kinetics is a generic term that is used in this work to describe the biochemical and physical characteristics, attributes, and behaviors of the molecular antennas. The transmitter molecular antenna in this work is the nanopores, and the receiver antennas are modeled as surface receptors proteins and surface receptor enzymes. Hence, the reception mechanism of the receiver is modeled separately as ligand-receptor binding kinetics and as enzyme-catalyzed biochemical reaction kinetics. For the rest of this paper, we shall refer to the surface receptors proteins and surface receptor enzymes simple as surface receptors and surface enzymes, respectively. In summary, the contributions of this work are as follows.

• 1. A nanotransmitter-diffusion channel model is presented, which takes into account the effect of the use of nanopores as the medium through which synthesized molecular signals leave the transmitter space. This model introduces a diffusion coefficient that accounts for the randomness in the density of the emitted molecules.

• 2. Based on the kinetic possibilities that arise due to the use of nanopores to emit molecular signals at the transmitter, and the reception at the receiver using surface receptor and surface enzymes, we derive explicit analytical expressions for the output molecular signal of the system for an impulsive input signal. The use of the impulsive input [30, 39, 40] is analogous to modeling the conventional electronic wireless communication systems by employing a dirac input.

• 3. We also derived system-theoretic propagation models for the different transmit-receive configurations. The arising transfer function-based models typifies the channel models for the scenarios where the influence of the molecular antennas are taken into consideration.

For the derived propagation models, the time and frequency responses of the systems to variation in designs and kinetic parameters and pulsatile signals are evaluated. We note that the idea of the system-theoretic approach is to circumvent the detailed characterization and description of the entire MC system, by making an analogy with linear time invariant systems. In this mode, the different transition states of the input molecular signal can be followed through the communication network in both time and frequency. However, just like many biological systems, the MC system is not linear, but there is a subtle analogy with linear systems which drives analytical tractability. We also point out that an interesting seminal work on system-theoretic model for MC is presented in [32]. However, unlike the approach in this present work, [32] does not include the transmitter antenna model and the use of enzymatic circuit at the receiver interface, and the frequency analysis is restricted to the signal attenuation.

The context and organization of this paper are made in such a way that make it easy for experienced and new researchers in the field of MC to follow and familiarize with the basics that underlay the MC concepts. The rest of this paper is organized as follows. In “Related works,” existing works that are related to the ideas in this paper are discussed. The system model and problem definition are presented in “System model and problem definition.” Preliminary assumptions that are considered in the development of the models in this work are discussed in “Preliminaries to molecular propagation model.” In “Diffusion-only molecular propagation model” and “Interface kinetic-diffusion molecular propagation model,” the diffusion-only and the diffusion-controlled interface kinetic-inclusive propagation models are presented, respectively. Simulation results and discussion are provided in “Illustrative examples, simulation results and discussion”.

2 Related works

To the best of our knowledge, existing works that considered system-theoretic approach to MC analysis are presented in [32, 41]. In [32], this approach is used to model the MC network as a biological circuit. The transmitter is modeled as comprising of an enzyme encoder that synthesizes the acts of some substrates concentration to produce the information molecules. The information molecules diffuse through the intercellular environment and are received by means of ligand-receptor binding process. The transfer function of the entire communication stages are provided, and analytical results of signal attenuation as a function of transmission radius and frequency are presented. Unlike our work, the transfer function and results in [32] do not take the transmitter antenna and receiver radius into consideration. More also, the pattern of the received signal as well as the use of surface enzymatic circuit for reception is not addressed in [32]. In [41], without incorporating any biochemical possibility, the system-theoretic approach is used to develop end-to-end concentration propagation model based on microfluidic channel configurations. Using the developed model, the design of microfluidic channels is performed based on the least-squares finite impulse response filtering to achieve the desired end-to-end transfer function. Other works that support or relate to some of the various concepts used in this work can be grouped into the ones that employ the assumption of linear models, surface receptor models, and diffusion models.

The assumption of linear time invariance model for the MC system as is made in this work has been used in [32, 40]. The transfer function model for diffusion-only and that of ligand-receptor reception model can be found in [32, 39] [42] and are used as reference and comparison to our work. We note that the MC system, just like many biological systems, is a nonlinear system. Hence, to apply linearity assumption, we shall restrict our analysis to behaviors near nominal operating range within which nonlinearity approximates linearity [43]. While the use of surface receptor proteins as the receiver molecular antenna and discussions of the ensuing biochemical kinetics can be found in [33–37], the incorporation of the transmitter antenna is not addressed. A detailed discussion on the ligand-receptor binding, trafficking and signaling is presented in [44], where the explicit expression for ligand-receptor binding action in the case of uniformly distributed ligands around the receptor location is provided.

Our work involves the use of enzymes as surface receptors and the derivation of explicit response expressions for this type of receiver molecular antenna. While enzyme-catalyzed reactions are ubiquitous in biochemical systems [45–50], in the context of MC only few works [35, 51–54] exist. In [52], an enzyme circuit is used to decode frequency coded signal in MC networks. In that regard, enzyme-catalyzed reaction that is modeled after the oscillatory behavior of phosphorylation-dephosphorylation cycle is used to model an MC encoder. But, no explicit expression for the system input and output concentrations were given. In contract to our deterministic model, a stochastic MC model based on the classical reaction–diffusion master equation is proposed in [35] based on the ligand-receptor mechanism. The derived closed-form expression of the mean receiver output is simplified only up to the differentials. In [53], the mitigation of intersymbol interference in MC systems using enzymes that freely diffuse in the propagation environment is presented. While our model considers enzymes that are immobilized at the receiver membrane, in [53], the enzymes diffuse freely along with the transmitted information molecules. More also, [53] does not provide any explicit close-form expressions for the output concentrations of the MC system. In [54], Nakano et al., also considered enzyme-catalyzed kinetics for the application of MC to targeted drug delivery but used the popular Michaelis-Menten equation for evaluation and also did not provide any explicit expressions for the output signals. An MC-based diffusion-controlled enzyme-catalyzed architecture for targeted drug delivery is presented in [51] where enzyme catalysis is considered in deriving analytical expressions for information transfer. However, [51] does not address system-theoretic concern, transfers functions, and does not consider the transmitter antenna in the model.

In relation to molecular diffusion channel models, more recent works such as [55] presented a 3D channel model of a molecular communications and analyzed its attenuation and propagation delay characteristics for an absorbing receiver. An end-to-end analysis of propagation noise and memory for molecular communication over microfluidic channels is presented in [56]. In [57], a realistic channel model for a table-top molecular communication platform that is capable for transmitting short text messages across a room is proposed based on an experimental platform presented in [58]. We note that [55–58] do not consider any biochemical possibilities in the respective models.

3 System model and problem definition

Before proceeding with the system model and definitions, some notation used in the rest of this work are explained below.

X ~ N(ϑ, σ 2)	Implies that the scalar variable X has a normal distribution with mean ϑ and variance σ 2.
X n  ~ N(ϑ n , Σ)	Implies that the vector variable X n has a multivariate normal distribution with mean vector ϑ n and covariance matrix Σ.
h(.)	System impulse responses or time dependent functions are symbolized in small letters
H(.)	The Laplace transform of h(.) written in capital letters.
|H(.)	Magnitude of H(.)
∇2	Laplace operator
Re(.)	Real part of a complex function
Ε[.]	Expectation operator
$A+B\kern1em \underset{\beta_{-1}}{\overset{\beta_1}{\leftrightarrow }}\kern1em C$	Representation of a biochemical reaction between molecular species A and B to for C under the influence of positive constants β and β − 1.

The schematic diagram of a typical molecular communication system is shown in Fig. 1, in which we have differentiated the diffusion channel from the propagation channel. We employ the term molecular propagation channel to include the channel that consists of the diffusion channel and the molecular antenna. Here, the term molecular antenna represents the coupling interface between the transmitter and the diffusion channel and between the diffusion channel and the receiver bio-circuit. The examples of nature’s molecular antenna include the protein receptor on the surface of the biological cell membrane, the cell membrane itself, and the gated ion channels. The focus of this paper is to model the molecular propagation channel illustrated in Fig. 1 (enclosed in red dotted lines). In Fig. 1, a(t) represents the molecular concentration that triggers information processing/encoding at the transmitter. The output of this information/encoding process w(t) is then transmitted using the transmitter molecular antenna. We denote the molecular concentration at the output of the transmit antenna by g(t). The function g(t) is conceptually the concentration of an arbitrary molecule G present at the output of the transmit antenna and the input of the diffusion channel. The molecular concentration observable at the output of the diffusion channel and which serves as the input to the receiver antenna is represented as ψ _C (t). The receiver antenna picks up this signal and outputs p(t). The molecular concentration p(t) is decoded and processed to yield the desired system output y(t).

Fig. 1

Block diagram of a bio-inspired molecular communication system (Best viewed in color)

Let us consider the representation of the molecular propagation channel as a cascade of the responses given as h _i (t,K _i ), i = 1, 2, 3 and shown in Fig. 2, where h _i (t,K _i ), represents the impulse responses of each of the ith unit. The term K _i represents the vector of all other parameters associated with each block. These parameters include the diffusion coefficients and rate constants. Also, let Fig. 2 be a linear system, hence, its transfer function is

Fig. 2

Block diagram of a molecular propagation channel

$$ {H}_{123}\left(s,{\boldsymbol{K}}_1,{\boldsymbol{K}}_2,{\boldsymbol{K}}_3\right) = {H}_1\left(s,{\boldsymbol{K}}_1\right){H}_2\left(s,{\boldsymbol{K}}_2\right){H}_3\left(s,{\boldsymbol{K}}_3\right) $$

(1)

This implies that

$$ {H}_{123}\left(s,{K}_1,{K}_2,{K}_3\right)=\frac{P(s)}{W(s)} $$

(2)

where W (s) and P(s) are the Laplace domain versions of w(t) and p(t), respectively, and s = σ + jω is the Laplace variable.

Depending on the task at hand and the design architecture of the transmitter and receiver, H ₁₂₃ (s, K ₁ , K ₂ , K ₃ ) will differ from one type of molecular antenna to another. For instance, H ₃ (s, K ₃ ) may be a membrane receptor in one design and a cell membrane itself in another. Hence, while H ₃ (s, K ₃ ) for the former depends on the ligand-receptor kinetics, the latter depends on the rate of endocytosis. This difference will significantly impact on the total system response as will be shown in this work.

In this work, we model the transmitter as a spherical nanoporous structure with a relative point source inside it. The nanopores are formed by perforating the surface of the transmitter, and through these nanopores, the molecules emitted from the point source transude as shown in Fig. 3. Existing methods and techniques for forming nanopores can be found in [59–63]. We designate the transfer function of the nanopores as H _N (s, K _N ). The dynamics of the molecules from the transmitter to the receiver space is modeled using Fick’s diffusion theory. In this sense, the transfer function of the classical diffusion medium is represented as H _C (s, K _C ). At the receiver, we consider two different configurations of molecular antenna kinetics, the surface receptor-ligand binding kinetics, and the enzyme-catalyzed reaction kinetics. The transfer functions of the surface receptor and enzyme receptor systems are given as H _R (s, K _R ) and H _E (s, K _E ), respectively. The schematic of molecular communication that involves the two scenarios are shown in Fig. 4, with their corresponding transfer function given in relation to (1) as

Fig. 3

An illustration of diffusion-controlled interface kinetics-inclusive molecular communication between two nano devices

Fig. 4

A block diagram of a diffusion-controlled surface receptor-inclusive propagation model and b diffusion-controlled surface enzyme-inclusive propagation model

$$ {H}_{NCR}\left(s,{\boldsymbol{K}}_N,{\boldsymbol{K}}_C,{\boldsymbol{K}}_R\right) = {H}_N\left(s,{\boldsymbol{K}}_N\right){H}_C\left(s,{\boldsymbol{K}}_C\right){H}_R\left(s,{\boldsymbol{K}}_R\right) $$

(3)

$$ {H}_{NCE}\left(s,{\boldsymbol{K}}_N,{\boldsymbol{K}}_C,{\boldsymbol{K}}_E\right) = {H}_N\left(s,{\boldsymbol{K}}_N\right){H}_C\left(s,{\boldsymbol{K}}_C\right){H}_E\left(s,{\boldsymbol{K}}_E\right) $$

(4)

In the preceding sections, we shall provide the system-theoretic model for the diffusion-only molecular propagation model and models that include the transmitter/receiver molecular antenna kinetics. The derivations of the analytical expressions of the various input and output signals in time and frequency domain are also provided. However, before doing that, we first present the preliminaries and generalize assumptions employed in the rest of this work.

4 Preliminaries to molecular propagation model

Before we proceed with the model development and derivation of expressions, we make the following simplifying assumptions in addition to those made in [51].

• 1. The nanopores and surface receptors/enzymes are assumed to be immobile and evenly distributed on the membrane of a spherical transmitter and receiver, respectively.

• 2. All the surface receptors/enzymes are equally active.

• 3. The process of surface receptors/enzyme immobilization has not altered the structure of the surface receptors/enzymes and the reaction rate coefficients.

• 4. The charge and hydrophobicity of the immobilized surface receptors/enzyme surface are not taken into consideration; hence, the influence of the surface receptors/enzyme surface net charge on the system is neglected.

• 5. Constant temperature and pH, as well as the homogeneity of the diffusion medium, are assumed.

• 6. The influence of medium velocity is not considered. Hence, the molecular motion is purely Brownian.

• 7. Both the transmitter and the receivers are stationary

• 8. The information molecules have spherical shape so that the diffusion coefficient due to rotation is neglected.

• 9. The receptors and enzymes on the receivers’ surface are of the same kind and specificity.

• 10. This model assumes that there are only one transmitter and one receiver, and there are no other competing molecules in the medium.

By making these assumptions, the overall output concentration and system-theoretic models are determined by mass transfer (diffusion) of the information molecules and the surface nanopores/receptors/enzyme kinetics.

5 Diffusion-only molecular propagation model

In relation to Fig. 4, the system-theoretic model for the diffusion-only phenomena is depicted in Fig. 5 and expressed as

Fig. 5

Equivalent system diagram for the diffusion-only propagation model

$$ {H}_C\left(s,{K}_C\right)=\frac{\varPsi_C(s)}{G(s)} $$

(5)

where Ψ _C (s) and G (s) are the Laplace transforms of ψ _C (t) and g(t), respectively.

The expression in (5) does not take into account the influence of the molecular antennas, hence, g(t) = w(t). The molecular signal g(t) can be expressed as the molecular concentration at a reference distance r ₀ about the transmitter and reference time t ₀, hence g(t) ≡ g ₀(t ₀, r ₀). We model h _C (t,K _C ) as the classical diffusion process using Fick’s diffusion law, where K _C = {D, r, θ, ϕ}. In this case, the molecular signal g ₀(t ₀, r ₀) that is initially at the input of the channel with diffusion coefficient D can be evaluated at any coordinate (r, θ, ϕ) and time t using the expression

$$ \frac{\partial {g}_c\left(t,r,\theta, \varphi \right)}{\partial t}=D{\nabla}^2{g}_c\left(t,r,\theta, \varphi \right)+f\left(t,r\right) $$

(6)

where ${\nabla}^2=\frac{1}{r^2}\;\left.\left({\partial}_r\left({r}^2{\partial}_r\right)+\frac{1}{ \sin \theta}\right.\kern0.5em {\partial}_{\theta}\left( \sin \theta \kern0.5em {\partial}_{\theta}\right)+\frac{1}{{ \sin}^2\theta }{\partial}_{\varphi}^2\right)$ and g _c (t, r, θ, φ) are the concentration of the molecule at any spatial distance defined by the coordinate (r, θ, φ) at time t. The function f(t,r) is the source of the information molecules, whose characteristics influence the solution to (6). We shall henceforth assume that D is independent of g(t), and hence of (r, θ, φ), which implies the homogeneity of the diffusion medium.

For g ₀(t ₀, r ₀) = f(t,r) = N _A δ(r-r ₀), the solution to (6) is the fundamental solution

$$ {g}_c\left(t,r\right)=\frac{N_A}{{\left(4\pi Dt\right)}^{1/2}} \exp\;\left(-\frac{{\left|\left(r-{r}_0\right)\right|}^2}{4Dt}\right) $$

(7)

The expression in (7) is a point concentration at a distance r [64] and does not take into consideration the receiver configuration and parameters. To model the portion of g _c (t, r) that is received without recourse to the molecular antenna, but is comparable to what is received using the molecular antennas, let v _ass + be a hypothetical volume termed reception space. This concept of reception space has been employed in [40] and is illustrated in Fig. 6, where the transmitter is located at the distance r∈ R from the center of v _ass +. The sphere of radius r _ass + that encompasses v _ass + is concentric with that of the spherical receiver volume of radius r _R . We consider that any molecule that diffuses into this volume is assumed to have been received.

Fig. 6

Spherical diffusion model of a molecular communication system

Let Y _n denote some arbitrary continuous random variable and suppose that the probability density of Y _n evaluated at y _n is the function f _n (y _n ). The probability of Y _n falling into some arbitrary interval [a b] can be expressed by a multivariate spherically symmetric distribution.

$$ {P}_{\vartheta}\left[a\le {Y}_n\le b\right]={E}_{\vartheta}\left[{\left({Y}_n-\vartheta \right)}^{\prime }{f}_n\left(a\le {Y}_n\le b\right)\right] $$

(8)

where (Y _n  − ϑ _n )′ is essentially the derivative of (Y _n  − ϑ _n ), ϑ _n is the unknown mean number of Y _n having a normal distribution, and f _n (a ≤ Y _n  ≤ b) is a smooth function. In relation to estimating the mean of a multivariate normal distribution, it was shown in [65] that if Y _n  ~ N (ϑ _n , Σ), then

$$ {E}_{\vartheta}\left[{\left({Y}_n-{\vartheta}_n\right)}^{\prime }{f}_n\left(a\le {Y}_n\le b\right)\right]=\varSigma \kern0.1em {E}_{\vartheta}\left[{\nabla}^{\prime }{f}_n\left(a\le {Y}_n\le b\right)\right] $$

(9)

Relating (6) to (9) and bearing in mind that variance is proportional to diffusion coefficient, it can be shown that

$$ r\frac{\partial {g}_c\left(t,r\right)}{\partial t}=rD\left[\right({\nabla}^2{g}_c\left(t,r\right)\Big] $$

(10)

To obtain the radial solution to (10), let us assume that the distance between the transmitter and the receiver is large enough so that we can consider the diffusing information molecules from the transmitter as being equally distributed around the receiver at time t. With this assumption, we can express the initial condition for (10) as the radial source

$$ {g}_c\left(0,r\right)=\frac{N_A\delta \left(r-{r}_0\right)}{4\pi {r}^2} $$

(11)

The solution to (10) that includes the reception space parameter is given by [66]

$$ {g}_c\left(t,r,{r}_{ass+}\right)=\frac{N_A}{4\pi r{r}_{ass+}{\left(4\pi Dt\right)}^{1/2}}\left\{ \exp\;\left(-\frac{{\left|\left({r}_{ass+}-r\right)\right|}^2}{4Dt}\right)\kern0.5em - \exp\;\left(-\frac{{\left|\left({r}_{ass+}+r-2{r}_R\right)\right|}^2}{4Dt}\right)\right\} $$

(12)

The concentration that is received in the space v _ass + = (4π/3)(r _ass + − r _R )³ for r _ass + ≈ r _R is given by

$$ {\psi}_C\left(t,{r}_R\right)={\left.4\pi {r}_R^2D{d}_r{g}_c\left(t,r,{r}_{ass+}\right)\right|}_{r_{ass+}={r}_R} $$

(13)

Hence,

$$ {\psi}_C\left(t,r,{r}_R\right)=\frac{D{N}_A\left({r}_Rr-{r}_R^2\right)}{2{\pi}^{1/2}{(Dt)}^{3/2}r}\kern1em {e}^{-\frac{{\left({r}_R-r\right)}^2}{4Dt}} $$

(14)

Therefore, using Laplace transform, it can be shown that

$$ {H}_C\left(s,{K}_C\right)=\frac{\varPsi_C(s)}{W(s)}=\kern0.5em \frac{4\pi s{X}_v\left({r}^2{r}_R-{r}_R^2r\right)}{\sqrt{D}}{e}^{-{\left({X}_v{s}^{-1/2}\right)}^{-1}} $$

(15)

where ${X}_v=\kern0.5em \sqrt{\frac{D}{{\left(r-{r}_R\right)}^2}}$

The result in (15) is validated by a related Fourier domain result presented in [67] where the transfer function of the diffusion channel is given as

$$ {H}_{CC}\left(\omega, {K}_C\right)=\mathrm{R}\mathrm{e}\left(\frac{\kern0.5em j\omega }{\pi Dr}{e}^{-\left(1+j\right)\sqrt{\frac{\omega }{2D}}}\right) $$

(16)

The graph of (15) and (16) are separately shown in Fig. 7. It can be seen that the diffusion channel behave like a low-pass filter. However, we note that unlike (15), (16) does not take the receiver parameters into consideration.

Fig. 7

Effect of a receiver radius on the concentration of the received signal, b receiver radius on transfer function, c transmit-receive distance on transfer function, and d transmit-receive distance on transfer function in [67]

6 Interface kinetic-diffusion molecular propagation model

While the diffusion-only model provides one with functions for quantifying the molecular information available within a space at some distance away from the transmitter, it does not avail us with the actual expression of the molecular information coupled into the receiver with respect to the molecular antenna used. To address this challenge, we consider the inclusion of the effects of the transmitter molecular antenna and the receiver molecular antenna into the diffusion-only model. We shall consider the effect of incorporating these antennas on the system response, where the transmitter molecular antenna is a nanopore, and the receiver antennas are surface receptors and surface enzymes.

6.1 Diffusion-controlled transmitter antenna kinetic-inclusive model

The system-theoretic model for the scenario, where the effect of the transmitter molecular antenna is taken into consideration, is depicted in Fig. 8 and expressed as

Fig. 8

System diagram of diffusion-controlled transmitter nanopores-only inclusive propagation model

$$ {H}_{NC}\left(s,{K}_{NC}\right)={H}_N\left(s,{K}_N\right){H}_C\left(s,{K}_C\right)=\frac{\varPsi_{NC}(s)}{W(s)} $$

(17)

In this case, G (s) ≠ W(s), hence, H _NC (s, K _NC ) ≠ H _C (s, K _C ). As stated earlier, the transmitter in this work is modeled as a spherical structure with a point source insider it and nanopores perforated on its surface through which molecules transude. The expressions (5)–(15) do not take the presence of the nanopores into account. To appreciate the effect of integrating the influence of the nanopores, we consider the illustrative model shown in Fig. 9. In Fig. 9a, there is one single pore whose kinetic implies that the source cannot be simply modeled as a point source since the emission through the pore is a function of time. When the number of nanopores is increased as depicted in Fig. 9b, and depending on the ration of the nanopores to the initial source concentration, the emission can be modeled as a distributed point source at various locations on the transmitter surface. In Fig. 9c, there are many nanopores, in which case, the emission can be approximated to the diffusion-only case depicted in Fig. 9d. The difference in scenarios Fig. 9c, d can be related to the diffusion coefficient of the medium and that of the interface, that is, the nanopores. In the case of Fig. 9d, the diffusion coefficient D in (5)–(15) can be expressed as that of the homogeneous case for which [68]

Fig. 9

Illustrative diagrams of the influence of number of nanopores on the molecular emission process

$$ D=\frac{K_BT}{6\pi \eta {r}_M} $$

(18)

where K _B is the Boltzmann constant, T is the absolute temperature of the system, η is the viscosity of the diffusion medium, and r _M is the radius of the spherically shaped diffusing particle.

However, in the case of Fig. 9c, certain factors need to be taken into consideration in computing D. These factors include the influence of the nanopores size r _N , the distribution f _N of the nanopores over transmitter surface ∂_Ω, and the ratio of the nanopores footprint to ∂_Ω. For Fig. 9c to be approximated to Fig. 9d, r _M  < < r _N must hold. More also, the ratio of the nanopores footprint ∂ _N to transmitter surface area ∂_Ω should be unity. Since the practicality of ∂_Ω/∂ _N  = 1 is technically unrealistic in this case, we have to include in the model the possibility of some of the diffusing molecules at an arbitrary time t being reflected back at some boundary points on the transmitter surface. The expression for the diffusion coefficient that accounts for the ratio of the nanopore radius to that of diffusing molecules can be found in [69]. In this present work, we consider the case where r _M  < < r _N , and the ratio of ∂_Ω to ∂ _N tends to unity, but not exactly unity, which is the scenario that is related to Fig. 9c. The condition ∂_Ω/∂_N ➔ 1 indirectly addresses the symmetric distribution concern of the pores for isotropic diffusion. Let P _d be the probability that a diffusing molecule at an arbitrary time t is reflected back. A more generalized diffusion coefficient is given by

$$ {D}_{eq}=\frac{K_BT}{6\pi \eta {r}_M}{P}_d\left({n}_N,\frac{\partial_{\varOmega }}{\partial_N},{f}_N\right) $$

(19)

where n _N is the number of the pores of equal dimensions. We define the probability function in (19) by

$$ {P}_d\left({n}_N\right)=1-\frac{1}{n_N}\mathrm{N}\left(\mu, {\sigma}^2\right) $$

(20)

where μ is the mean and σ ² is the variance of n _N . As n _N tends to infinity, D _eq = D. Therefore, for the transmitter interface kinetic-diffusion model, D _eq replaces D in (14), thus

$$ {\psi}_{NC}\left(t,r,{r}_R\right)=\frac{D_{eq}{N}_A\left({r}_Rr-{r}_R^2\right)}{2{\pi}^{1/2}{\left({D}_{eq}t\right)}^{3/2}r}{e}^{-\frac{{\left({r}_R-r\right)}^2}{4{D}_{eq}t}} $$

(21)

Accordingly,

$$ \begin{array}{c}\hfill {H}_{NC}\left(s,{K}_{NC}\right)=\kern0.5em \frac{4\pi s{X}_{v, eq}\left({r}^2{r}_R-{r}_R^2r\right)}{\sqrt{D_{eq}}}{e}^{-{\left({X}_{v, eq}{s}^{-1/2}\right)}^{-1}}\hfill \\ {}\hfill {X}_{v, eq}=\kern0.5em \sqrt{\frac{D_{eq}}{{\left(r-{r}_R\right)}^2}}\hfill \end{array} $$

(22)

where K _NC  = {D _eq, r, r _R }. For zero-mean unit variance, the effect of variation in n _N and invariably ∂_Ω/∂ _N is illustrated in Fig. 10. It can be observed that D _eq produces a noisy version of H _NC (s, K _NC ), and as n _N tends to infinity, Fig. 10d approximates Fig. 10c. This noise can be regarded as a fluctuation in the concentration of the information molecules. Since it is this concentration that is detected by the receiver antennas, implicitly, the random fluctuation (noise) in the information molecule concentration may result in deviations in the number of bound receptors from the mean behavior predicted by deterministic models [44]. Hence, the lesser n _N is the more the effect of molecule random fluctuation is felt and the more the deterministic model assumption pursued in this work is violated.

Fig. 10

Graphical examples of the effect of number of nanopores on the diffusion characteristics where in a n _N  = 5, b n _N  = 10, c n _N  = 100, and d n _N  = 1000

6.2 Diffusion-controlled receiver antenna kinetic-inclusive model

6.2.1 Surface receptor kinetic-diffusion model

We shall now consider the inclusion of the influence of the surface receptor as the receiver molecular antenna to the H _NC (s, K _NC ) model. We consider the scenario where the result of the information reception by the surface receptor initiates a signal transduction cascade in the receiver biological circuit. We term molecular antennas that produce this type of effect enzyme-linked receptors since in actual sense, these surface receptors are not necessarily enzymes but possess intracellular domains that are associated with an enzyme. In this case, after a ligand binds to the extracellular domain of the receptor, the receptor undergoes conformational change, which subsequently activates the intracellular enzymes thereby triggering various activities in the biological circuit of the receiver. An example of this type of molecular antenna is the tyrosine kinase receptor, which is a high-affinity cell surface receptor for many polypeptide growth factors, hormones, and cytokines. This form of ligand-receptor binding reaction kinetics can be modeled as

$$ {\varPsi}_M+Q\kern1em \underset{\xi_{-1}}{\overset{\xi_1}{\leftrightarrow }}\kern1em {P}_R $$

(23)

where Ψ _M is the information carrying molecule or received ligand present inside v _ass+, Q is the surface receptor, and P _R is the ligand-receptor complex. The terms ξ ₁ and ξ _-1 are positive rate reaction constants. The system-theoretic model for the scenario where the effect of the surface receptor is taken into consideration is depicted in Fig. 11 and expressed as

Fig. 11

System diagram of diffusion-controlled transmitter nanopores/surface receptor-inclusive propagation model

$$ {H}_{NCR}\left(s,{K}_N,{K}_C,{K}_R\right)=H\left(s,{K}_{NC}\right)H{}_R\left(s,{K}_R\right)=\frac{P_R(s)}{\varPsi_{NC}(s)} $$

(24)

where P _R (s) is the Laplace transform of p _R (t).

To obtain H _NCR (s, K _N , K _C , K _R ), we apply the law of mass action to obtain the following differential equation model for (23)

$$ \frac{d{p}_R\left(t,r,{r}_R\right)}{dt}={\xi}_1{\psi}_{NC}\left(t,r,{r}_R\right)q\left(t,r,{r}_R\right)-{\xi}_{-1}{p}_R\left(t,r,{r}_R\right) $$

(25)

where p _R (t,r,r _R ), ψ _NC (t,r,r _R ), and q(t) are the concentrations of P _R , Ψ _M , and Q, respectively. If we assume that the concentrations of Ψ _M and Q are independent of r and r _R , that is, of the diffusion process, and q(t,r,r _R ) is constant, then the transfer function can be written as [32, 39]

$$ H{}_{RR}(s)\propto \frac{\xi_1}{s+{\xi}_{-1}} $$

(26)

However, the rate of complex formation, that is the formation of P _R , is proportional to the rate at which Ψ _M and Q encounter themselves, as well as the rate at which the encounter is broken. The rate of encounter is in itself proportional to the joint probability of finding Ψ _M and Q in the same space, v _ass+. The probability of finding a number of Ψ _M in v _ass+ is dependent on r, r _R , and r _ass+. Assuming that once a received molecule dissociated from bonding with a receptor it diffuses away, the probability that a Ψ _M observed in v _ass+ at a time t ₁ is still inside v _ass+ at time t ₂ and have not diffused away also defines the rate of encounter. Obviously, as is illustrated in Fig. 12, the larger x = (r _ass+ −r _R ) is the longer a Ψ _M stays in v _ass+ and the more the probability of rebinding and viz a viz. Hence, to ensure that the no-rebinding assumption holds, it is desirable that x is very small. The rate at which the encounter is broken depends on the affinity constant ξ _a of the binding [70], which can be calculated as ξ _a  = ξ ₁/ξ _− 1. A high ξ _a implies that the received molecules stay long in the bonding; hence, more molecules may have diffused away within the bonding time t _bond and not get the opportunity to bond with the receptors. Recall that we have assumed in the derivation of (14) that r _ass + ≈ r _R ; hence x = (r _ass+ −r _R ) is very small. Let us express t _d = t ₂ − t ₁ and t _bond  = t _1/2 = 0.693/ξ _− 1 [70]. If we assume the same time scale for the binding and diffusion processes, then considering the scenario where t _bond ≤ t _d, we can write that

Fig. 12

Illustration of the ligand-receptor binding inside the conceptual reception volume

$$ {\xi}_{-1}\ge \frac{1.386D}{x^2} $$

(27)

This implies that the larger the value of reception space x, the higher the probability that many molecules may have diffused away from the space without getting the opportunity to bond with the receptors. In this case, a more accurate model must include the probability of escape from the space as a function of time. Such model has been presented in [71, 72]. On the other hand, the smaller the value of reception space x, the lower the probability that some molecules may have diffused away from the space without getting the opportunity to bond with the receptors. In this case, a deterministic model can be pursued results in which the rate of complex formation will now depend only on the concentration of receptors and Ψ _M at each point in time.

Hence, if we assume that ψ _NC (t,r,r _R ) > > q(t) and with no consideration to the effect of receptor clustering [73], then for p _R (t, r, r _R ) = q ₀ − q(t), where q ₀ is the total concentration of the receptors, the rate of complex formation is expressed as

$$ \frac{d{p}_R\left(t,r,{r}_R\right)}{dt}={\xi}_1{\psi}_{NC}\left(t,r,{r}_R\right)\left[{q}_0-{p}_R\left(t,r,{r}_R\right)\right]-{\xi}_{-1}{p}_R\left(t,r,{r}_R\right) $$

(28)

Hence, it can be shown that

$$ {p}_R\left(t,r,{r}_R\right)=\frac{\xi_1{q}_0{\psi}_{NC}\left(t,r,{r}_R\right)}{\xi_1{\psi}_{NC}\left(t,r,{r}_R\right)+{\xi}_{-1}}\left(1-{e}^{-\left({\xi}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{\xi}_{-1}\right)t}\right) $$

(29)

We note that a related result is presented in (2–12), [44]), but without considering the influence of the diffusion channel.

Hence, while the received molecular signal within the space x → 0 is expressed by (14), the actual information concentration detected by the surface receptor antenna is expressed in (29).

Applying Laplace transform method to (28) yield

$$ s{P}_R\left(s,r,{r}_R\right)={\xi}_1{q}_0\;\varPsi \left(s,r,{r}_R\right)-\left({\xi}_1{q}_0\;{\varPsi}_{NC}\left(s,r,{r}_R\right)+{\xi}_{-1}\right){P}_R\left(s,r,{r}_R\right) $$

(30)

Hence,

$$ \frac{P_R\left(s,r,{r}_R\right)}{\varPsi_{NC}\left(s,r,{r}_R\right)}=\frac{\xi_1{q}_0\;}{s+\left({\xi}_1{q}_0\;{\varPsi}_{NC}\left(s,r,{r}_R\right)+{\xi}_{-1}\right)} $$

(31)

It can be seen that (30) is not a linear system; hence, its transfer function cannot be generalized but has to be defined for a specific input. In this case, a specific transfer function for (30) in the case of a delta input can be obtained from the ratio of the Laplace transform of (29) to that of (21). Thus,

$$ \begin{array}{c}\hfill {H}_{NCR}\left(s,{K}_P,{K}_C,{K}_R\right)=\frac{q_0\sqrt{D}\kern0.5em r\kern0.5em {e}^{{\left({X}_{v, eq}{s}^{-1/2}\right)}^{-1}}\;}{N_A{X}_{v, eq}\left({r}_Rr-{r_R}^2\right)}\left(\frac{1}{s}-\frac{G_{0,3}^{3,0}\left({r}_G\left|0,\frac{1}{2},1\right.\right)}{\left({k}_{-1}+s\right)\kern0.5em \sqrt{\pi }}\right)\hfill \\ {}\hfill {r}_G=\frac{{\left({N}_A{k}_1{r}_R\right)}^2{\left(r-{r}_R\right)}^2\left({k}_{-1}+s\right)}{16\pi D{r}^2}\hfill \end{array} $$

(32)

where ${G}_{0,3}^{3,0}\left({r}_G\left|0,\frac{1}{2},1\right.\right)$ is a Meijer-G function. Therefore

$$ H{}_R\left(s,{K}_R\right)=\frac{2{q}_0\kern0.5em {e}^{{\left({X}_{v, eq}{s}^{-1/2}\right)}^{-1}}\;}{N_A4\pi s{r_R}^2}\left(\frac{1}{s}-\frac{G_{0,3}^{3,0}\left({r}_G\left|0,\frac{1}{2},1\right.\right)}{\left({k}_{-1}+s\right)\kern0.5em \sqrt{\pi }}\right) $$

(33)

The transfer function plots H _RR (s), H _NCR (s, r _R ), and H _R (s, r _R ) are shown in Fig. 13 for n _N  = 1000, k ₁ = 50, k ₋₁ = 1, r = 100 nm, r _R  = 30 nm, N _A  = 1 μM, and q ₀ = 0.002 μM. It can be observed that for the range of frequency considered, and bearing in mind the assumptions made in arriving at (26), (32), and (33), the plot of (33) as compared to (26) indicates that the only-receptor system tends towards the system in (26) where the diffusion process is not taken into account. However, with the inclusion of the diffusion process, the response significantly deviates from the no diffusion state.

Fig. 13

Transfer function of the magnitude of a H _NCR (s, r _R) and H _R (s, r _R) and b H _RR (s)

6.2.2 Surface enzyme kinetic-diffusion model

We have considered the scenario where the information reception is achieved using surface receptors. In that case, the reception process initiates a signal transduction cascade in the biological circuit located inside the receiver. In another scenario of concern in this work, instead of surface receptors, we have surface enzymes. These enzymes act as receptors on the surface of the receiver. However, instead of causing a signal transduction, these enzymes catalyzes the production of product molecules from their interaction with the information molecules. In this sense, the enzyme-Ψ _M interaction forms a biological circuit at the surface instead of the inside. This type of receiver has been presented in [51] for the purpose of targeted drug delivery. The system-theoretic model for the scenario where the effect on the surface enzymes is taken into consideration is depicted in Fig. 14 and expressed as

Fig. 14

System diagram of diffusion-controlled transmitter nanopores/surface enzyme-inclusive propagation model

$$ {H}_{NCE}\left(s,{K}_N,{K}_C,{K}_E\right)=H\left(s,{K}_{NC}\right)H\left(s,{K}_E\right)=\frac{P_E(s)}{\varPsi_{NC}(s)} $$

(34)

where P _E (s) is the Laplace transform of p _E(t).

Let E represent the enzyme receptor, we can represent the enzyme-Ψ _M binding action by the popular enzyme-catalyzed kinetic model expressed as

$$ {\varPsi}_M+E\kern1em \underset{k_{-1}}{\overset{k_{\;1}}{\leftrightarrow }}\kern1em C\kern1em \overset{k_2}{\to}\kern1em {P}_E+E $$

(35)

where k ₁, k ₋₁, and k ₂ are positive rate constants, C is the enzyme-Ψ _M complex, and P _E is the product molecule. Let e(t), c(t), and p _E (t) be the concentrations of E, C, and P _E , respectively. Since the number of the immobilized enzymes on the membrane of the receiver are fixed, the total enzyme e ₀ at any instant of time is such that c(t) = e ₀ − e(t). Applying the law of mass action, the following differential equations arise from (35)

$$ \frac{d{\psi}_{NC}\left(t,r,{r}_R\right)}{dt}=-{k}_1\psi \left(t,r,{r}_R\right){e}_0+{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)c(t)+{k}_{-1}c(t) $$

(36)

$$ \frac{dc(t)}{dt}={k}_1{\psi}_{NC}\left(t,r,{r}_R\right){e}_0-{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)c(t)-{k}_{-1}c(t)-{k}_2c(t) $$

(37)

$$ \frac{d{p}_E(t)}{dt}={k}_2c(t) $$

(38)

To derive p _E (t), we have to first of all derive c(t). The expression for c(t) can be obtained from (37) by assuming as previously done in the case of the surface receptor that ψ _NC (t,r,r _R ) > > e(t); hence, (37) simplifies to a pseudo first-order equation. This implies that

$$ {\displaystyle \int \frac{dc(t)}{k_1{\psi}_{NC}\left(t,r,{r}_R\right){e}_0-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]c(t)}}={\displaystyle \int dt} $$

(39)

Integrating both sides of (39) yields

$$ \frac{{ \log}_e\kern0.5em \left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right){e}_0-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]c(t)\right]}{-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]}=t+\vartheta $$

(40)

Since c(t = 0) = 0,

$$ \vartheta =\frac{{ \log}_e\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right){e}_0\right]}{-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]} $$

(42)

Hence

$$ { \log}_e\kern0.5em \left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right){e}_0-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]c(t)\right]=-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]t+{ \log}_e\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right){e}_0\right] $$

(43)

Therefore

$$ c(t)=\frac{\psi_{NC}\left(t,r,{r}_R\right){e}_0}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]}-\frac{e^{-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]t+ \log \left({k}_1{e}_0{\psi}_{NC}\left(t,r,{r}_R\right)\right)}\Big)}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]{k}_1} $$

(44)

where k _M  = (k _− 1 + k ₂)/k ₁ is the classical Michaelis–Menten kinetic constant.

Having obtained c(t), we simply plug it into (38) and integrate, thus

$$ \begin{array}{c}\hfill {p}_E(t)={\displaystyle \int {k}_2c(t)}\kern0.5em dt\hfill \\ {}\hfill ={\displaystyle \int \frac{\psi_{NC}\left(t,r,{r}_R\right){e}_0{k}_2}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]}-\frac{k_2{e}^{-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]t+ \log \left({k}_1{e}_0{\psi}_{NC}\left(t,r,{r}_R\right)\right)}}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]{k}_1}}\kern0.5em dt\hfill \end{array} $$

(45)

Solving (45) as a pseudo first-order equation yields

$$ \begin{array}{c}\hfill {p}_E(t)=\frac{\psi_{NC}\left(t,r,{r}_R\right){e}_0{k}_2}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]}t-{\displaystyle \int \frac{k_2{e}^{-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]t+ \log \left({k}_1{e}_0{\psi}_{NC}\left(t,r,{r}_R\right)\right)}}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]{k}_1}}\kern0.5em dt\hfill \\ {}\hfill =\frac{\psi_{NC}\left(t,r,{r}_R\right){e}_0{k}_2}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]}t-\frac{k_2}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]{k}_1}{\displaystyle \int {e}^{-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]t+ \log \left({k}_1{e}_0{\psi}_{NC}\left(t,r,{r}_R\right)\right)}}\kern0.5em dt\hfill \end{array} $$

(46)

Therefore

$$ {p}_E\left(t,r,{r}_R\right)={k}_2\left(\frac{\psi_{NC}\left(t,r,{r}_R\right){e}_0}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]}t+\frac{e^{-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]t+ \log \left({k}_1{e}_0{\psi}_{NC}\left(t,r,{r}_R\right)\right)}}{{\left(\left[{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_M\right]{k}_1\right)}^2}\right) $$

(47)

Hence, to obtain H _NCE (s, K _N , K _C , K _E ), the Laplace domain analytical expression for p _E (t, r, r _R ) is required. This expression can be obtained by making the simplifying assumption that the rate of catalysis is high so that ψ _NC (t, r, r _R ) > > k _M . Hence,

$$ {p}_E\left(t,r,{r}_R\right)={k}_2\left({e}_0t+\frac{e_0{e}^{-\left[{k}_1{\psi}_{NC}\left(t,r,{r}_R\right)+{k}_{-1}+{k}_2\right]t}}{\left[{\psi}_{NC}\left(t,r,{r}_R\right)\right]{k}_1}\right) $$

(48)

Recall that ${\psi}_{NC}\left(t,r,{r}_R\right)=\frac{D_{eq}{N}_A\left({r}_Rr-{r}_R^2\right)}{2{\pi}^{1/2}{\left({D}_{eq}t\right)}^{3/2}r}{e}^{-\frac{{\left({r}_R-r\right)}^2}{4{D}_{eq}t}}$, if we consider the case where t > 0, then we can approximate ${e}^{-{\left({r}_0-{r}_R\right)}^2/(4Dt)}\approx 1$ in both (21) and (46). Therefore,

$$ \begin{array}{c}\hfill {P}_E\left(s,r,{r}_R\right)=\frac{e_0{k}_2}{s^2}+\frac{2{e}_0{k}_2r\sqrt{D_{eq}}}{N_A{k}_1\left(r{r}_R-{r}_R^2\right){\left({k}_{-1}+{k}_2+s\right)}^{5/2}}{G}_{0,3}^{3,0}\left({r}_G\left|0,\frac{1}{2},\frac{5}{2}\right.\right)\hfill \\ {}\hfill {r}_G=\frac{{\left({N}_A{k}_1{r}_R\right)}^2{\left(r-{r}_R\right)}^2\left({k}_{-1}+{k}_2+s\right)}{16\pi D{r}^2}\hfill \end{array} $$

(49)

where ${G}_{0,3}^{3,0}\left({r}_G\left|0,\frac{1}{2},\frac{5}{2}\right.\right)$ is again a Meijer-G function.

Hence, for a dirac input,

$$ {H}_{NCE}\left(s,{K}_N,{K}_C,{K}_E\right)=\frac{\frac{e_0{k}_2}{s^2}+\frac{2{e}_0{k}_2r\sqrt{D_{eq}}}{N_A{k}_1\left(r{r}_R-{r}_R^2\right){\left({k}_{-1}+{k}_2+s\right)}^{5/2}}{G}_{0,3}^{3,0}\left({r}_G\left|0,\frac{1}{2},\frac{5}{2}\right.\right)}{\frac{\kern0.5em {N}_A{r}_R\left(r{r}_R-{r}_R^2\right){X}_{v, eq}}{r\kern0.5em \sqrt{D}\kern0.5em }\kern0.5em {e}^{-{\left({X}_{v, eq}{s}^{-1/2}\right)}^{-1}}} $$

(50)

Consequently

$$ H{}_E\left(s,{K}_E\right)=\frac{2{q}_0\kern0.5em {e}^{{\left({X}_{v, eq}{s}^{-1/2}\right)}^{-1}}\;}{N_A4\pi s{r_R}^2}\left(\frac{e_0{k}_2}{s^2}+\frac{2{e}_0{k}_2r\sqrt{D_{eq}}}{N_A\left(r{r}_R-{r}_R^2\right){s}^{5/2}}{G}_{0,3}^{3,0}\left({r}_G\left|0,\frac{1}{2},\frac{5}{2}\right.\right)\right) $$

(51)

The plots of H _NCE (s, r _R ) and H _E (s, r _R ) are shown in Fig. 15 for n _N  = 1000, k ₁ = 50, k ₋₁ = 1, k ₂ = 5, r = 100 nm, r _R  = 30 nm, N _A  = 1 μM, and e ₀ = 0.002 μM. And just like the analysis made in the surface receptor case, Fig. 15 indicates that the receptor-only system, which can be considered as a system where the diffusion process is not taken into account, deviates significantly from the no diffusion plus receptor scenario.

Fig. 15

Transfer function of the magnitude of H _NCE (s, r _R) and H _E (s, r _R)

7 Illustrative examples, simulation results and discussion

In this section, we present the time and frequency analysis of the system-theoretic models and expressions presented in this work. Firstly, a comparison between the reception space-only model and the antenna-inclusive model is presented so as to quantify how much information it actually received. Secondly, we evaluate the effect on the system output of the geometric parameters associated with the receiver and the medium and the effect of the reaction parameters associated with the surface receptor-based receiver model and the enzyme receptor-based receiver model. Thirdly, we extend the parameter evaluation to the state of equilibrium when the rate of reactants association is equal to the rate of dissociation. Finally, we provide illustrative examples that are channeled towards evaluating the frequency characteristics of the various propagation models presented to different pulsatile inputs. In all the evaluations and simulations, we use D = 0.85 mPa s⁻¹, which is the average values for a typical cerebrospinal fluid at 37 °C [74]. For clarity, ψ _C (t, r, r _R ) and ψ _CC (t, r) represent the concentration of the molecules received in the reception space with and without consideration to the receiver radius, respectively. The quantities p _R (t, r, r _R ) and p _E (t, r, r _R ) represent the concentration of the information molecules couples into the receiver by the surface receptor and surface enzymes, respectively.

7.1 Comparison between reception space-only model and antenna-inclusive model

Previously, we have defined the reception space as a hypothetical volume around the receiver within which any molecule that diffuses into this space is assumed to have been received [40, 51]. The assumption of a hypothetical reception space is related to the case of perfect absorption of molecules incident upon the surface of the receiver. However, it is most likely that not all molecules that hit the surface receptors are captures since some of the receptors might be busy processing already captured molecules at the time of hit. Hence, a more practical model should include a hitting/capture probability function [71, 72]. In the model that this paper presents, we pursued a deterministic model and consider that the reception space is very small such that the probability that some molecules may have diffused away from the space without getting the opportunity to bond with the receptors is negligible. In this sense, a deterministic model can be pursued in which the rate of complex formation depends only on the concentration of receptors and information molecules at each point in time. Based on this explanation and in pursuit of one of the objectives of this paper, comparison between the reception space-only model and the model that is inclusive of the molecular antennas is presented as shown in Fig. 16. The distance between molecular transmitter and receiver as well as the receiver radius is fixed in all the simulations in Fig. 16. For simplicity, all the values used are kept dimensionless. The values of the transmitter/receiver distance and the receiver radius are kept at unity and 0.3, respectively, and the value of N _A is unity.

Fig. 16

Plots for comparison between the signal concentration of the reception space-only model and antenna-inclusive (a) at varying transmitted signal concentration for the surface receptor-inclusive model (b) at varying kinetic constant values for the surface receptor-inclusive model (c) at varying transmitted signal concentration for the surface enzyme-inclusive model (d) at varying kinetic constant values for the surface receptor-inclusive model.

In Fig. 16a, c, the comparison between the concentrations of the reception space-only model and the models that include the surface receptor and surface enzymes are made, respectively, with respect to variation in the concentrations of the antennas. In obtaining Fig. 16a, ξ ₁ = ξ _− 1 = 10, and for Fig. 16c, k ₁ = k ₂ = 1 and k ₋₁ = 0. It can be observed that the actual information molecule concentration that is captured by the antenna is less than that defined within the reception space. And as anticipated, the received concentration decreases with decrease in the concentration of the receiver antennas. In Fig. 16b, d, the comparison between the concentrations of the reception space-only model and the models that include the surface receptor and surface enzymes are made, respectively, with respect to variation in the reaction parameters. In obtaining Fig. 16b, d, the concentrations of the antennas are kept at 0.9 and 0.5, respectively. The affinity constant is kept at unity such that ξ ₁/ξ _− 1 = 1 and $\frac{k_{-1}+{k}_2}{k_1}=1$. It can be observed that varying the values of the reaction parameters results in change in the antenna output concentration.

7.2 Time evolution of the output functions

Let us consider the communication between a molecular transmitter and a receiver placed less than 400 nm apart. As pointed out in [51], it is important to maintain a constant ratio ρ of ${r}_R^2$ to the square of the enzyme base radius ${r}_E^2$ or receptor base radius ${r}_Q^2$ in other to get a linear response. Hence, we define the effective concentrations of the surface enzymes and receptors as e _eff = e ₀ = bn _E ρ and q _eff = q ₀ = bn _R ρ, respectively, where n _E and n _R are the molar masses of an enzyme and receptor, respectively. The term b is the active enzyme factor used to keep e _eff / q _eff constant as r _R is varied. The following values are used in the simulation, r _E = r _Q = 5 nm, b = 1, n _E = n _R  = 0.002 μM, and n _N  = 1000. For different receiver radius as well as different values of the reaction constants k ₁, k ₋₁, k ₂, ξ ₁, and ξ _− 1, the output signals p _R (t,r,r _R ) and p _E (t,r,r _R ), are plotted in Fig. 17. For good visual, the parameters r _E, r _Q, n _E, and n _R are scaled up by × 10⁶.

Fig. 17

Plots of the output signal variation with (a) receiver radius for surface receptor-inclusive model (b) receiver radius for surface enzyme-inclusive model (c) forward rate constant for surface receptor-inclusive model (d) forward rate constant for surface enzyme-inclusive model (e) reverse rate constant for surface receptor-inclusive model (f) reverse rate constant for surface enzyme-inclusive model

Observation on Fig. 17a, b indicates that as we increase r _R , the output concentration increases as envisaged. Fig. 17c, d indicates that as k ₁ and ξ ₁ decrease, and keeping k ₋₁ (or ξ _− 1) and k ₂ constant, the outputs decrease accordingly. The point it peaks varies depending on the ratio of k ₂ to k ₋₁. Considering Fig. 17e, f, it can be observed that as k ₋₁(or ξ _− 1) increases, the output signal understandably reduces. This is so since the rate of successful complex formation is reduced with increase in the reverse rate constants. In Fig. 17g, h, the variation in the enzyme-catalyzed receiver process is observed. As expected, at k ₋₁ = 50 and k ₋₁ = 0, the output concentration increases with increase in k ₂. Hence, aside increasing the concentrations of the information molecules, increase in output concentration can be obtained by increasing the effective receiver radius, {k ₁, ξ ₁, k ₂}, and decreasing {k ₋₁, ξ _− 1}.

Let us consider the scenarios where the systems are in equilibrium. In this state, the rate of complex formation is expected to be equally to the rate of complex dissociation. This condition arises when in the case of the surface receptor kinetics, ξ ₁ = ξ _− 1. The results for the surface receptor kinetics are shown in Fig. 18. It can be observe that though the ratio ξ _− 1/ξ ₁ = 1 is the same in all cases, significant variation in output can be observed especially at low values of ξ ₁ = ξ _− 1. This variation comes from the exponential term in (29). The value of this term is more significant at lower values of ξ ₁ = ξ _− 1 than at higher values.

Fig. 18

Variation of surface receptor MC system output signal with rate constants at equilibrium for (a) low value of rate constants (b) high value of rate constants

In the case of the enzyme receptor, the state of equilibrium is also defined by the zero state of the rate of change of the enzyme-ligand complex. It is important to point out that the presence of the enzymes does not alter the equilibrium position of the reaction, rather it only lowers the energy required to speed up the reaction. Since, the interest here is not on the concentration of the complexes formed, we assume that k ₋₁ = 0, and consider equilibrium as the case where ${k}_M=\frac{k_{-1}+{k}_2}{k_1}=1$. Figure 19 shows the results for different values of k _M where again significant variation in output can be observed especially at low values of k ₁ = k ₂ due to the exponential term in (47).

Fig. 19

Variation of surface enzyme MC system output signal with rate constants at equilibrium for (a) low value of rate constants (b) high value of rate constants

We note that the values of the reaction rate constants k ₁, k ₋₁, k ₂, ξ ₁, and ξ _− 1 are typically associated with the type of reactions and are influenced by temperature variation. Hence, for applications in which environmental temperature change is a factor, the above analysis should be taken into consideration while keeping in mind that temperature variation also affects the D [75]. The relationship between the rate at which a reaction occurs and the reaction temperature are determined by the Arrhenius equation [76].

7.2.1 Frequency analysis of the propagation models

The ideas of frequency and frequency response are well developed in the context of wireless communication and macroengineering. However, in the context of MC, the notions of frequency and frequency response are not yet well defined; hence, there is knowledge gap that needs to be filled before well-articulated frequency domain analyses of MC systems can be made in the general literature. In this subsection, we shall attempt to provide some verifiable frequency and frequency response concepts in MC systems.

Typically, frequency is related to rate. And “rate” is, of course, a description of how frequently something happens. A unimolecular rate describes how fast one molecule does something, and relatively, a multimolecular rate describes how fast a group of molecules does something. The “something” may be the underlying process/event that defines a phenomenon or an abstraction of the phenomena. In this context, these processes or events include any biological or biophysical process in the MC system. One may be interested in the rate at which molecules traverse a given distance to reach a destination if observed repeatedly. In this case, we may resort to the use of quantum physics concepts [77–79] or Brownian trajectories that follow detailed geometrical structure [80] to quantify the behavior of the motion process and develop a classical model of the diffusing molecule dynamics. From this point of view, the frequency of the system of interest may be given as being inversely proportional to the period it takes signaling molecules to traverse a given distance. This rate is dependent on the square root of the diffusion coefficient of the system, and distance. The diffusion coefficient is defined by the size and mass of the diffusing molecules, as well as the viscosity of the diffusion medium. The period it takes for a signaling molecule to traverse a distance can subsequently be used to define the bandwidth of the medium or channel. Examples of documented signaling periods and their associated bandwidth are given in Table 1. While it may look like it will take a molecule a very long time to cover a distance and effect information interchange/computation, in nanocommunication, the distance of interest is always very small. In fact, if you release a typical protein inside a bacterial cell, within one-hundredth of a second, it is equally likely to be found anywhere in the cell [81].

Table 1 Examples of frequencies of some biophysical signaling

On the other hand, we may be interested in the rate at which biological processes occur. Such processes include enzymatic process, mitosis, meiosis, ligand-receptor association, metabolic/catabolic processes, and the rate of molecular reactions in the interior of a cell. The rate at which these processes occur depends on the rate at which constituent molecules encounter themselves and the rate at which reactions occur, the latter being dependent on temperature and concentration. The frequency of such biological system may be expressed as being inversely proportional to the period it takes the event to occur. The period ranges from few seconds to years and in some cases infinity. The frequencies of some biological events are given in Table 2.

Table 2 Examples of frequencies of some biological processes

In typical MC system, the phenomena of molecule dynamics and biological processes are coupled. This is so since molecules dynamically traverse a distance to deliver information to the receiver, where biological processes take place with some influences again from molecular dynamics in the interior of the receiver. Hence, to obtain the knowledge of the rate at which an entire communication takes place depends on detailed model of the molecular dynamics and the biochemical kinetics of every path in the entire communication sequence. At this juncture, one is compelled to ask the following question. Of what importance is the knowledge of the frequency of occurrence of the biological/biophysical processes and events in the MC system? To answer this question, it is important to point out that due to the nano nature and limited computation capability of nanomachines, a number of them need to work in a cooperative manner to execute a task. Hence, the knowledge of the rate at which communication takes place between any two nodes is crucial to the synchronization of nanodevices in a given MC network. Obtaining such knowledge accurately is a challenging and complex task, which this present work does not address. In some situation, rather than seeking the knowledge of the rate of the entire process, one may be interested in having an understanding of how the rate at which an event or process occurs in one entity influences the behavior of another system. This is the case in where a dynamic analysis of system response is crucial. This brings us to another interpretation of frequency with respect to the MC system. This form on interpretation is in consideration to the representation of a biophysical activity in the form of a mathematically amenable analytical expression. In this sense, the output of a biophysical process/system or a phenomenon is approximately modeled as some composite or elementary signals such as sine wave, pulses and spikes, or even a combination of all. An example is the representation of Ca²⁺ concentrations at the apical and basal membrane as sine waves [82]. Another example is the modeling of insulin secretion by beta-cells and artificial injection as train of Gaussian pulse [83, 84]. The different elementary molecular signals have different properties, which can be in terms of shape, concentration, and rate of repetition. The rate of repetition of these signals can be related to their frequency. However, a more insightful interpretation of frequency in this case will be to relate the characteristics of these signals to the well-explored frequency interpretation in the domain of conventional signal processing. From this perspective, given that we have an analytical model of a biophysical system, we can evaluate the response of the system to some excitations that are abstractions of another system. This approach is often termed frequency response analysis and is a handy tool in the case where the frequency behavior of the exciting system is analytically tractable. We shall pursue this approach in evaluating the frequency behavior of the diffusion-controlled interface kinetics-inclusive system-theoretic propagation models presented in this work. The response can be evaluated in time and frequency domain. However, the ease of computation in the frequency domain endears it to the analysis of complex systems such as biophysical systems.

Consider the overall model of the diffusion-controlled interface kinetics-inclusive system-theoretic propagation system presented in this work and depicted in Fig. 20. Let the system (molecular propagation channel) in Fig. 20 be excitable by some signal w(t), whose Laplace transform is W(s). Let {H _NCR (s, r ,r _R ), H _NCE (s, r ,r _R ), H _R (s, r ,r _R ), H _E(s, r ,r _R ) } ⊂ H(s, r ,r _R ), the comparison between the diffusion-controlled surface receptor kinetics-inclusive system-theoretic propagation system to diffusion-controlled surface enzyme kinetics-inclusive system-theoretic propagation system is shown in Fig. 21. It can be seen that unlike the low-pass filter behavior of the diffusion-only channel shown in Fig. 7, the inclusion of molecular antenna added a band-stop effect to the overall system response. This assertion is made bearing in mind the various simplifying assumptions made in arriving at the expressions of the propagation models. The consequence or practical application of this type of response in MC paradigm is not clear. However, since the frequency range of interest in most biological signaling systems is in the lower frequency range (see Tables 1 and 2), the idea of low-pass filter suffices.

Fig. 20

Diffusion-controlled interface kinetics-inclusive system-theoretic propagation system

Fig. 21

Frequency characteristics of H _NCR (s, r _R), H _NCE (s, r _R), H _R (s, r _R), and H _E (s, r _R) at a–c low frequency range and b–d high frequency range

Let us go further in evaluating the response of the system to an exemplary input signal. The input signal of concern in this analysis is a pulsatile signal [85], which is a signal representation of a molecular flow with periodic variations. Our interest in this type of signaling is due to the findings in the literature that in many biological signaling systems, oscillations or rhythmic pulses can be more efficient in evoking responses than the same input dose given at a constant level [86]. For instance, calcium oscillations can evoke enhanced gene expression compared to fixed level in lymphocytes [87]. It has also been shown that during ovarian cycle in women, a pulsatile release pattern of the gonadotropin-releasing hormone within an appropriate frequency range supports reproductive function, while the release of the equivalent constant level of the hormone is ineffective [88, 89]. In the context of MC, the pulsatile signal can be obtained by transmitter modulation of molecule concentration. Our objective is to model three different pulsatile signals concentration, all having the same total molecular concentration. In relation to (11), we express the pulsatile signal as

$$ {w}_N(t)\propto {\displaystyle \sum_{i=1}^N{\alpha}_i\delta \left(r-i{r}_0\right)} $$

(52)

where α is the amplitude of the pulse and r ₀ is the repetition interval. For practical purpose, instead of the dirac function, a Gaussian function ℘(r) approximation is used. In this present paper, we model the inputs as pulsatile signals of different characteristics. The three pulsatile signals are plotted in Fig. 22 and expressed as

Fig. 22

Input molecular signals a w ₁(t), b w ₂(t), and c w ₃(t)

$$ {w}_1(t)={\displaystyle \sum_{i=1}^1{\wp}_1\left(r-i{r}_{0,1}\right)} $$

(53)

$$ {w}_2(t)={\displaystyle \sum_{i=1}^2{\wp}_2\left(r-i{r}_{0,2}\right)} $$

(54)

$$ {w}_3(t)={\displaystyle \sum_{i=1}^3{\wp}_3\left(r-i{r}_{0,3}\right)} $$

(55)

where ℘ ₁(r), ℘ ₂(r) and ℘ ₃(r) are Gaussian pulses of different widths, and r _0,1,r _0,2, and r _0,3 are their different repetition intervals. The following condition applies

$$ {\displaystyle \sum {w}_1}={\displaystyle \sum {w}_2}={\displaystyle \sum {w}_3} $$

(56)

The spectrograms of the three signals are shown in Fig. 23. It can be observed that frequency of w ₃(t) is the largest followed by w ₂(t) and lastly w ₁(t). In practice, w ₁(t) can be achieved by the stimuli-triggered instantaneously releases of information molecules encapsulated in a liposome, which acts as a molecular transmitter. Such instantaneous release for liposome has been achieved in [90] for laser-triggered partially polymerized liposome, where about 70 % content release was achieved. Timed sequential liposomes content release can be employed to obtain the signals by w ₂(t) and lastly w ₃(t).

Fig. 23

Spectrogram W1(omega), W3(omega) and W3(omega) of the molecular signals (a) w ₁(t) (b) w ₂(t) and (c) w ₃(t), respectively

Based on the analytical representation of these input signals, we excite the system in Fig. 20 with these inputs to obtain the corresponding output responses shown in Fig. 22. It can be observed from Fig. 24 that the diffusion channel behaves like a frequency selective channel especially, which is evident from the effect on the various inputs. More also, the use of higher frequency signals evoke higher response in terms of the peak of the output concentration. This assertion may be related to the result in [91] where it is observed that the deliverable information in an MC system of genetically engineered bacteria increases with clock frequency.

Fig. 24

Effects of a H _CC (s), b H _C (s, r _R), c H _NCR (s, r _R), and d H _NCE (s, r _R) on the oscillating input signals w ₁(t), w ₂(t), and w ₃(t)

It is important to note that the receiver model presented in this work is a non-absorbing receiver model that arises from the inherent characteristics of the particular molecular antenna under consideration. The molecular antennas of interest within the scope addressed in this work are the surface enzymes and surface receptors. Such model can find application in scenarios that are related to the activation of prodrugs by enzymes at targeted drug delivery sites [51] and any reception mechanism that mimics the ligand-receptor binding mechanisms such as those presented in [34, 92, 93]. However, for absorbing receivers [94], which in relation to cellular activities mimic information transfer through membrane diffusion and endocytosis/exocytosis processes, the model has to be modified to account for information molecules detection and removal from the propagation environment.

8 Conclusions

In this work, we have presented and investigated molecular propagation models that incorporate both diffusion and molecular antenna kinetics possibilities for molecular communications. Specifically, we have derived time and frequency system-theoretic analytical expressions for the diffusion-controlled interface (antenna) kinetics-inclusive propagation models. The influence of the geometric and reaction parameters on the derived propagation models was investigated. Our investigation shows that the performance of such system depends crucially on the effective enzyme concentration at the receiver, as well as the rate constants of the reaction. The dependency of the propagation model characteristics on the reaction constants implies that for cases where the temperature of the operating environment significantly varies with time, the influence of temperature on the system may be significant and should be taken care of. We have also provided some insight into the notion of frequency in molecular communication and used such notions to analyze the frequency behavior of the MC propagation models. The frequency analysis indicates that based on the simplifying assumptions made in deriving the analytical expressions, while the diffusion-only channel acts as a low-pass filter, the diffusion-controlled interface (antenna) kinetics-inclusive propagation models is a band-stop filter. We have also shown that the models show signs of frequency selectivity for repetitive pulse signals with varying frequencies.