Analysis of Receiver Covered by Heterogeneous Receptors in Molecular Communications

Xinyu Huang1, Yuting Fang2, Stuart T. Johnston2, Matthew Faria2, Nan Yang1, and Robert Schober3 1School of Engineering, Australian National University, Canberra, ACT 2600, Australia 2University of Melbourne, Parkville, VIC 3010, Australia 3Friedrich-Alexander University Erlangen-Nürnberg, 91058 Erlangen, Germany

Abstract

This paper analyzes the channel impulse response of an absorbing receiver (RX) covered by multiple non-overlapping heterogeneous receptors with different sizes and arbitrary locations in a molecular communication system. In this system, a point transmitter (TX) is assumed to be uniformly located on a virtual sphere at a fixed distance from the RX. Considering molecule degradation during the propagation from the TX to the RX, the expected molecule hitting rate at the RX over varying locations of the TX is analyzed as a function of the size and location of each receptor. Notably, this analytical result is applicable for different numbers, sizes, and locations of receptors, and its accuracy is demonstrated via particle-based simulations. Numerical results show that (i) the expected number of absorbed molecules at the RX increases with an increasing number of receptors, when the total area of receptors on the RX surface is fixed, and (ii) evenly distributed receptors lead to the largest expected number of absorbed molecules.

Index Terms:

Molecular communication, channel modeling, size of receptors, location of receptors.

I Introduction

Molecular communication (MC) has emerged as a promising technology to facilitate micro-scale or nano-scale communications [1]. In MC, information is encoded into molecules that are released from a transmitter (TX). After being released, the molecules propagate in a fluid medium until they arrive at a receiver (RX). The RX then detects and decodes the information encoded in the molecules. Therefore, an accurate modeling of practical RXs is of great significance for the design and development of MC systems.

In biology, living cells recognize specific molecules via receptor-ligand interactions, which are fundamental for a cell to communicate with its neighbours and the entire organism [2]. Specifically, a signal is delivered to the cell by binding of ligands to their complementary receptors, which results in a cascade of chemical reactions. Due to the complexity of modeling the entire process of signal delivery, the RX models considered in the MC literature are relatively simple, compared to the actual process. For example, the most widely-adopted RX models in MC are the passive RX, fully absorbing RX, and reactive RX [3], where the passive RX ignores receptor-ligand interactions, while the fully absorbing RX and reactive RX ignore the receptor size and assume an infinite number of receptors covering the entire RX surface. To enhance the practicality of RX models, the authors of [4] considered a partially absorbing RX and some recent studies, e.g., [5, 6, 7, 8], assumed the uniform distribution of a finite number of receptors on the RX surface, where all receptors have the same size. Specifically, the authors of [5] assumed that molecules are absorbed once they hit the receptors, and the authors of [6, 7] assumed that a reversible reaction occurs once molecules hit the receptors. Very recently, [8] investigated the receptor occupancy induced by the competition of molecules when binding to the receptors in synaptic MC.

While interesting, the previous studies may not be accurate in practical scenarios when the receptors on the RX surface have different sizes and arbitrary locations. For example, receptors tend to form clusters on the RX surface at specific locations [9]. Since a given ligand can activate all receptors within a cluster, the cluster can be regarded as a single larger receptor. In [10], the authors investigated a RX covered by receptors with different sizes and arbitrary locations under steady-state conditions and only derived the capacitance of the RX, but did not investigate the time-varying number of molecules absorbed by the RX. In this paper, we model receptors as absorbing patches (APs) and assume that there are multiple non-overlapping APs on the RX surface, where the APs have different sizes and arbitrary but fixed locations. When molecules hit the APs, they are absorbed by the RX. We further assume a point TX uniformly located on a virtual sphere at a fixed distance from the RX. By taking into account that molecules may degrade when they propagate from the TX to the RX, we investigate the expected channel impulse response (CIR) between the TX and the RX averaged over the varying locations of the TX when only the distance between the TX and RX is fixed. Here, the expected CIR is defined as the expected molecule hitting rate at the RX [3].

Our major contributions can be summarized as follows. We derive the expected molecule hitting rate, the expected fraction of absorbed molecules, and the expected asymptotic fraction of absorbed molecules as time approaches infinity at a RX with multiple APs, where all expressions are functions of the size and location of each AP. The desired expressions allow us to investigate the impact of different numbers, sizes, and locations of APs on molecular absorption. Furthermore, we compare three types of APs distributions, namely, APs evenly distributed over the RX surface (i.e., the APs are equally spaced), APs randomly distributed over the RX surface, and APs located within a region on the RX surface. Particle-based simulations (PBSs) are used to verify the accuracy of our expressions. Our numerical results reveal that the expected number of absorbed molecules increases with an increasing number of APs, when the total area of the RX surface covered by APs is fixed. We also show that evenly distributed APs lead to a larger number of absorbed molecules than the other two types of AP distributions considered.

II System Model

Refer to caption — Figure 1: Illustration of the MC system model where a point TX communicates with a spherical RX covered by multiple APs.

In this paper, we consider an unbounded three-dimensional (3D) environment where a point TX communicates with a spherical RX, as depicted in Fig. 1. We choose the center of the RX as the origin of the environment and denote the radius of the RX by $r_{\scriptscriptstyle\mathrm{R}}$ . There are $N_{\mathrm{p}}$ non-overlapping APs on the RX surface. We denote the $i$ th AP by $\mathrm{AP}_{i}$ . We assume the shape of each AP as a circle and denote $a_{i}$ as the radius of $\mathrm{AP}_{i}$ . We define $\mathcal{A}$ as the ratio of the total area of APs to the RX surface, i.e., $\mathcal{A}=\sum_{i=1}^{N_{\mathrm{p}}}a_{i}^{2}/4r_{\scriptscriptstyle\textnormal{R}}^{2}$ . In spherical coordinates, we denote $\vec{l}_{i}=[r_{\scriptscriptstyle\textnormal{R}},\theta_{\mathrm{p},i},\varphi_{\mathrm{p},i}]$ as the location of the center of $\mathrm{AP}_{i}$ , where $\theta_{\mathrm{p},i}$ and $\varphi_{\mathrm{p},i}$ are the polar angle and azimuthal angle, respectively. Once molecules hit any AP, they are absorbed by the RX immediately. As is customary, for analytical tractability, we ignore the occupancy of molecules to the APs such that several molecules can be absorbed by an AP at the same time. Furthermore, we model the part of the RX surface that is not covered by APs as perfectly reflective, which means that molecules are reflected back once they hit this part of the RX surface.

In the considered system, the point TX is uniformly located on a virtual sphere with a distance $r_{0}$ from the center of the RX. For given $r_{0}$ , the CIR of the RX is different for different locations of the TX. This is due to the fact that our considered RX has a heterogeneous boundary condition where the locations and sizes of the APs are determined by $\vec{l}_{i}$ and $a_{i}$ , respectively. Accordingly, this work focuses on the expected CIR averaged over all possible locations of the TX when only the distance between the TX and RX is given. Thus, our results provide valuable insights for the practical case where only the distance between the TX and RX is known, while the accurate angular position of the TX relative to the RX is difficult to obtain. We note that the measurement of the distance between two cells is much easier than determining the relative angle between two cells by using the concentration gradient of molecules released by one of the cells [11].

We assume that the propagation environment between TX and RX is a fluid medium with uniform temperature and viscosity. At time $t=0$ , an impulse of $N_{\sigma}$ molecules of type $\sigma$ is released from the TX. Once molecules are released, they diffuse randomly with a constant diffusion coefficient $D_{\sigma}$ . Moreover, we consider unimolecular degradation in the propagation environment, where type $\sigma$ molecules can degrade to type $\hat{\sigma}$ molecules that cannot be recognized by the RX, i.e., $\sigma\stackrel{{\scriptstyle k_{\mathrm{d}}}}{{\longrightarrow}}\hat{\sigma}$ [12, Ch. 9], where $k_{\mathrm{d}}\;[\mathrm{s}^{-1}]$ is the degradation rate.

III Analysis of Channel Impulse Response

In this section, we derive i) the expected molecule hitting rate, ii) the expected fraction of absorbed molecules, and iii) the expected asymptotic fraction of absorbed molecules as $t\rightarrow\infty$ at a RX with multiple APs by applying boundary homogenization [13]. To this end, we first derive the expected CIR of a RX with a uniform surface reaction rate. A uniform surface reaction rate implies that the reactivity of the molecules is identical for all points on the RX surface. We then assume the system to be in steady state and derive the diffusion current of molecules across the RX surface. Based on the diffusion current in the steady state, we finally determine the effective reaction rate and apply it to derive the expected CIR of a RX with multiple APs.

III-A Problem Formulation

In spherical coordinates, we denote the molecule concentration at time $t$ at location $\vec{r}$ by $C(\vec{r},t)$ , $|\vec{r}|\geq r_{\scriptscriptstyle\textnormal{R}}$ . When an impulse of molecules is released from the TX at time $t=0$ , the initial condition can be expressed as $C(\vec{r},t\rightarrow 0)=\delta(r-r_{0})/(4\pi r_{0}^{2})$ [14, Eq. (3.61)], where $\delta(\cdot)$ is the Dirac delta function. After the release, the diffusion of molecules in the propagation environment is described by Fick’s second law as follows [15]

\displaystyle\frac{\partial C(\vec{r},t)}{\partial t}=D_{\sigma}\nabla^{2}C(\vec{r},t)-k_{\mathrm{d}}C(\vec{r},t),

(1)

where $\nabla^{2}$ is the 3D spherical Laplacian. When molecules hit the RX surface, the reaction between the molecules and the RX surface is described by the radiation boundary condition [14, Eq. (3.64)], which is given by

\displaystyle D_{\sigma}\frac{\partial C(|\vec{r}|,t)}{\partial|\vec{r}|}\bigg{|}_{|\vec{r}|=r_{\scriptscriptstyle\textnormal{R}}}=wC(r_{\scriptscriptstyle\textnormal{R}},t),

(2)

where $w$ denotes the reaction rate. The unit of $w$ is $\mu\mathrm{m}/\mathrm{s}$ , which is also validated by (2). We note that $w\rightarrow\infty$ when $\vec{r}\in\Omega_{\mathrm{AP}}$ while $w=0$ when $\vec{r}\in\Omega_{\mathrm{R}}$ , where $\Omega_{\mathrm{AP}}$ and $\Omega_{\mathrm{R}}$ represent the parts of the RX surface that are fully covered and not covered by APs, respectively. Due to the heterogeneous boundary condition in (2), it is difficult to directly solve (1) to obtain $C(\vec{r},t)$ . In this paper, we apply boundary homogenization to derive the expected CIR. The main idea of boundary homogenization is replacing the heterogeneous boundary condition in (2) by a uniform boundary condition with an appropriately chosen effective surface reaction rate, denoted by $w_{\mathrm{e}}$ , which means that we replace the RX surface in Fig. 1 with an equivalent uniform surface with reaction rate $w_{\mathrm{e}}$ . Hence, (2) can be rewritten by replacing $w$ with $w_{\mathrm{e}}$ . We derive the expected CIR of the RX and $w_{\mathrm{e}}$ in the following subsections.

III-B Expected CIR of RX with Uniform Surface Reaction Rate

In this subsection, we analyze the expected CIR of a RX with uniform surface reaction rate $w$ . Based on the initial condition in $C(\vec{r},t\rightarrow 0)$ and the boundary condition in (2), the authors of [14] derived $C(|\vec{r}|,t)$ by solving (1) when $k_{\mathrm{d}}=0$ . We denote $h_{\mathrm{u}}(t,w)$ as the expected molecule hitting rate at a RX with a uniform surface reaction rate. We specify $h_{\mathrm{u}}(t,w)$ including the effect of molecule degradation in the following lemma.

Lemma 1

The expected molecule hitting rate at a RX with uniform surface reaction rate $w$ at time $t$ is given by

		$\displaystyle h_{\mathrm{u}}(t,w)=\frac{r_{\scriptscriptstyle\mathrm{R}}w}{r_{0}}\left[\frac{1}{\sqrt{\pi D_{\sigma}t}}\exp\left(-\frac{\varepsilon^{2}}{t}-k_{\mathrm{d}}t\right)-\gamma(w)\right.$
		$\displaystyle\left.\times\exp\left[\gamma(w)(r_{0}-r_{\scriptscriptstyle\mathrm{R}})+\zeta(w)t\right]\mathrm{erfc}\left(\frac{\varepsilon}{\sqrt{t}}+\gamma(w)\sqrt{D_{\sigma}t}\right)\right],$		(3)

where $\varepsilon\!=\!\frac{r_{0}-r_{\scriptscriptstyle\mathrm{R}}}{\sqrt{4D_{\sigma}}}$ , $\gamma(w)\!=\!\frac{wr_{\scriptscriptstyle\mathrm{R}}+D_{\sigma}}{D_{\sigma}r_{\scriptscriptstyle\mathrm{R}}}$ , $\zeta(w)\!=\!\gamma(w)^{2}D_{\sigma}-k_{\mathrm{d}}$ , and $\mathrm{erfc}(\cdot)$ is the complementary error function.

Proof:

According to [16], $h_{\mathrm{u}}(t,w)$ can be obtained via $h_{\mathrm{u}}(t,w)=h_{\mathrm{u}}(t,w)\big{|}_{k_{\mathrm{d}}=0}\times\exp\left(-k_{\mathrm{d}}t\right)$ . We also find that $h_{\mathrm{u}}(t,w)\big{|}_{k_{\mathrm{d}}=0}=4\pi r_{\scriptscriptstyle\textnormal{R}}^{2}wC(r_{\scriptscriptstyle\textnormal{R}},t)\big{|}_{k_{\mathrm{d}}=0}$ based on [14, Eq. (3.107)]. Combining these two results, we obtain (1). ∎

We further denote $H_{\mathrm{u}}(t,w)$ as the fraction of molecules captured by the RX by time $t$ and present it in the following corollary.

Corollary 1

The fraction of molecules captured by a RX with uniform surface reaction rate $w$ by time $t$ is given by

\displaystyle H_{\mathrm{u}}(t,w)=\frac{r_{\scriptscriptstyle\mathrm{R}}w}{r_{0}}\left[\alpha_{1}(t)-\alpha_{2}(t,w)\right],

(4)

where

	$\displaystyle\alpha_{1}(t)=$	$\displaystyle\frac{1}{2\sqrt{k_{\mathrm{d}}D_{\sigma}}}\left[\exp(-\beta)\mathrm{erfc}\left(\frac{\varepsilon}{\sqrt{t}}-\sqrt{k_{\mathrm{d}}t}\right)\right.$
		$\displaystyle\left.-\exp\left(\beta\right)\mathrm{erfc}\left(\frac{\varepsilon}{\sqrt{t}}+\sqrt{k_{\mathrm{d}}t}\right)\right]$		(5)

and

\displaystyle\alpha_{2}(t,w)\!=\!\frac{1}{2\zeta(w)}\left[\psi_{1}(t,w)-\psi_{2}(t,w)\right]-\frac{\gamma(w)\exp(-\beta)}{\zeta(w)}

(6)

with

	$\displaystyle\psi_{1}(t,w)=$	$\displaystyle 2\gamma(w)\exp\left(\gamma(w)(r_{0}-r_{\scriptscriptstyle\mathrm{R}})+\zeta(w)t\right)$
		$\displaystyle\times\mathrm{erfc}\left(\frac{\varepsilon}{\sqrt{t}}+\gamma(w)\sqrt{D_{\sigma}t}\right),$		(7)

		$\displaystyle\psi_{2}(t,w)=\left(\gamma(w)^{2}\sqrt{\frac{D_{\sigma}}{k_{\mathrm{d}}}}-\gamma(w)\right)\exp(-\beta)$
		$\displaystyle\times\mathrm{erf}\left(\frac{\varepsilon}{\sqrt{t}}-\sqrt{k_{\mathrm{d}}t}\right)-\left(\gamma(w)^{2}\sqrt{\frac{D_{\sigma}}{k_{\mathrm{d}}}}+\gamma(w)\right)$
		$\displaystyle\times\left[\exp(-\beta)-\exp(\beta)\mathrm{erfc}\left(\frac{\varepsilon}{\sqrt{t}}+\sqrt{k_{\mathrm{d}}t}\right)\right],$		(8)

$\mathrm{erf}(\cdot)=1-\mathrm{erfc}(\cdot)$ is the error function, and $\beta=(r_{0}-r_{\scriptscriptstyle\mathrm{R}})\sqrt{k_{\mathrm{d}}/D_{\sigma}}$ .

Proof:

$H_{\mathrm{u}}(t,w)$ can be obtained as $H_{\mathrm{u}}(t,w)=\int_{0}^{t}h_{\mathrm{u}}(u,w)\mathrm{d}u$ . By substituting (1) into this equation, we obtain (4). ∎

We further denote $H_{\mathrm{u},\infty}(w)$ as the asymptotic fraction of molecules captured by the RX as $t\rightarrow\infty$ . We present $H_{\mathrm{u},\infty}(w)$ in the following corollary.

Corollary 2

As $t\rightarrow\infty$ , the asymptotic fraction of molecules captured by a RX with uniform surface reaction rate $w$ is given by

\displaystyle H_{\mathrm{u},\infty}(w)=\frac{r_{\scriptscriptstyle\mathrm{R}}w\left(\gamma(w)-\sqrt{\frac{k_{\mathrm{d}}}{D_{\sigma}}}\right)}{r_{0}\zeta(w)}\exp(-\beta).

(9)

Proof:

Please see Appendix A. ∎

III-C Determination of Effective Reaction Rate

In this section, we determine the effective reaction rate $w_{\mathrm{e}}$ . First, we investigate the diffusion flux of molecules across the RX surface, denoted by $J$ , in the steady state. The diffusion flux $[\textrm{molecule}\cdot\textrm{m}^{-2}\cdot\textrm{s}^{-1}]$ is the rate at which molecules move across a unit area in a unit time [15], and is given by [15, Eq. (2.6)]

\displaystyle J=-D_{\sigma}\frac{\partial C(|\vec{r}|,t)}{\partial|\vec{r}|}\bigg{|}_{|\vec{r}|=r_{\scriptscriptstyle\textnormal{R}}}.

(10)

We next define the rate of molecule movement across the RX surface in a unit time as the diffusion current, denoted by $I$ , which is given by $I=-4\pi r_{\scriptscriptstyle\mathrm{R}}^{2}J$ . In the steady state, we have $\partial C(\vec{r},t)/\partial t=0$ . If we set $k_{\mathrm{d}}=0$ in (1), we obtain

\displaystyle\nabla^{2}C(\vec{r},t)=0.

(11)

As explained in [17], since (11) is analogous to the Laplace’s equation for the electrostatic potential in charge-free space, the diffusion current to an isolated absorbing RX of any size and shape can be expressed as [15, Eq. (2.24)]

\displaystyle I=4\pi D_{\sigma}GC_{0},

(12)

where we define $G$ as the “capacitance” of the RX and $C_{0}$ is the molecule concentration at $|\vec{r}|\rightarrow\infty$ in the steady state. We note that $C_{0}=1$ was adopted in some previous studies, e.g., [10, 6]. We further note that $G$ measures the ability of RX to absorb molecules, which is different from the conventional electrical capacitance of a conductor, denoted by $\hat{G}$ . According to [15], if the RX and conductor have the same size and shape, $G$ can be obtained as $G=\hat{G}/4\pi\epsilon_{0}$ , where $\epsilon_{0}$ is the vacuum permittivity. Since the expression for $\hat{G}$ has been derived for a variety of conductors, we can obtain $G$ directly if the RX and conductor have the same size and shape. For example, we have $\hat{G}=4\pi\epsilon_{0}r_{\scriptscriptstyle\textnormal{R}}$ for a spherical conductor with radius $r_{\scriptscriptstyle\textnormal{R}}$ . According to the relationship between $G$ and $\hat{G}$ , we can obtain the capacitance of a fully absorbing RX with the same radius, denoted by $G_{\mathrm{a}}$ , as $G_{\mathrm{a}}=r_{\scriptscriptstyle\mathrm{R}}$ . Therefore, the diffusion current of a fully absorbing RX, denoted by $I_{\mathrm{a}}$ , is given by $I_{\mathrm{a}}=4\pi D_{\sigma}r_{\scriptscriptstyle\mathrm{R}}C_{0}$ , which aligns with the derivation in [17, Eq. (1)]. For the RX with multiple APs in Fig. 1, the capacitance of the RX, denoted by $G_{\mathrm{p}}$ , was derived in [10] by using the method of matched asymptotic expansions. Specifically, $G_{\mathrm{p}}$ is a function of the size and location of each AP, and it is given by [10, Eq. (3.37a)]

		$\displaystyle\frac{1}{G_{\mathrm{p}}}=\frac{2}{N_{\mathrm{p}}\overline{m}\kappa r_{\scriptscriptstyle\mathrm{R}}}\Bigg{[}1+\frac{\kappa}{2N_{\mathrm{p}}\overline{m}}\ln\left(\frac{\kappa}{2}\right)\sum_{i=1}^{N_{\mathrm{p}}}m_{i}^{2}+\frac{\kappa}{N_{\mathrm{p}}\overline{m}}$
		$\displaystyle\times\bigg{(}\sum_{i=1}^{N_{\mathrm{p}}}m_{i}s_{i}+2\sum_{i=1}^{N_{\mathrm{p}}}\sum_{j=i+1}^{N_{\mathrm{p}}}m_{i}m_{j}\mathcal{F}(\vec{l}_{i}^{\prime},\vec{l}_{j}^{\prime})\bigg{)}\!+\!\left(\kappa\ln\left(\frac{\kappa}{2}\right)\!\!\right)^{2}$
		$\displaystyle\times\frac{\vartheta}{4N_{\mathrm{p}}\overline{m}}+\mathcal{O}\left(\kappa^{2}\ln\left(\frac{\kappa}{2}\right)\right)\Bigg{]},$		(13)

where $\kappa=\frac{a_{1}}{r_{\scriptscriptstyle\mathrm{R}}}$ , $m_{i}=\frac{2a_{i}}{r_{\scriptscriptstyle\mathrm{R}}\kappa\pi}$ , $\overline{m}=\frac{1}{N_{\mathrm{p}}}\sum_{i=1}^{N_{\mathrm{p}}}m_{i}$ , $s_{i}=\frac{m_{i}}{2}\left(\ln\left(\frac{4a_{i}}{r_{\scriptscriptstyle\mathrm{R}}\kappa}\right)-\frac{3}{2}\right)$ , $\vartheta=\frac{\left(\sum_{i=1}^{N_{\mathrm{p}}}m_{i}^{2}\right)^{2}}{N_{\mathrm{p}}\overline{m}}-\sum_{i=1}^{N_{\mathrm{p}}}m_{i}^{3}$ , and

\displaystyle\mathcal{F}(\vec{l}_{i}^{\prime},\vec{l}_{j}^{\prime})=\left[\frac{1}{|\vec{l}_{i}^{\prime}-\vec{l}_{j}^{\prime}|}+\frac{1}{2}\ln|\vec{l}_{i}^{\prime}-\vec{l}_{j}^{\prime}|-\frac{1}{2}\ln\left(2+|\vec{l}_{i}^{\prime}-\vec{l}_{j}^{\prime}|\right)\right]

(14)

with $\vec{l}_{i}^{\prime}=\vec{l}_{i}/r_{\scriptscriptstyle\mathrm{R}}$ . In (III-C), $\mathcal{O}(\cdot)$ represents the infinitesimal of higher order, which is omitted during calculation.

When all APs have the same size, (III-C) can be further simplied as presented in the following corollary.

Corollary 3

When all APs have identical sizes, i.e., $a_{1}=a_{2}=\cdots=a_{\scriptscriptstyle{N_{\mathrm{p}}}}$ , $G_{\mathrm{p}}$ can be simplified as follows

	$\displaystyle\frac{1}{G_{\mathrm{p}}}=$	$\displaystyle\frac{\pi}{N_{\mathrm{p}}\kappa r_{\scriptscriptstyle\mathrm{R}}}\Bigg{[}1+\frac{\kappa}{\pi}\bigg{(}\!\ln(2\kappa)-\frac{3}{2}+\frac{4}{N_{\mathrm{p}}}\sum_{i=1}^{N_{\mathrm{p}}}\sum_{j=i+1}^{N_{\mathrm{p}}}\mathcal{F}(\vec{l}_{i}^{\prime},\vec{l}_{j}^{\prime})\bigg{)}$
		$\displaystyle+\mathcal{O}\left(\kappa^{2}\ln\left(\frac{\kappa}{2}\right)\right)\Bigg{]}.$		(15)

Proof:

When $a_{1}=a_{2}=\cdots=a_{\scriptscriptstyle{N_{\mathrm{p}}}}$ , we have $m_{i}=\frac{2}{\pi}$ and $\vartheta=0$ . By substituting $m_{i}$ and $\vartheta$ into (III-C), we obtain (3). ∎

The capacitance of a RX with a single AP can be obtained by setting $N_{\mathrm{p}}=1$ in (3). In [10], the authors applied a high order asymptotic expansion to derive a more accurate expression, which is given by [10, Eq. (6.31)]

	$\displaystyle\frac{1}{G_{\mathrm{p}}}\bigg{\|}_{N_{\mathrm{p}}=1}=$	$\displaystyle\frac{\pi}{\kappa r_{\mathrm{R}}}\bigg{[}1+\frac{\kappa}{\pi}\left(\ln(2\kappa)-\frac{3}{2}\right)-\frac{\kappa^{2}}{\pi^{2}}\left(\frac{\pi^{2}+21}{36}\right)$
		$\displaystyle+\mathcal{O}(\kappa^{3}\ln\kappa)\bigg{]}.$		(16)

The diffusion current of molecules across a RX with multiple APs, denoted by $I_{\mathrm{p}}$ , can be obtained by replacing $G$ with $G_{\mathrm{p}}$ in (12). We denote $h_{\mathrm{p}}(t)$ , $H_{\mathrm{p}}(t)$ , and $H_{\mathrm{p},\infty}$ as the expected molecule hitting rate, the expected fraction of absorbed molecules, and the expected asymptotic fraction of absorbed molecules at a RX with multiple APs, respectively. Exploiting the fact that the ratio between the expected asymptotic fraction of absorbed molecules at a RX with multiple APs and that at a fully absorbing RX is equal to the ratio between the diffusion current across a RX with multiple APs and that across a fully absorbing RX [6], we derive $w_{\mathrm{e}}$ and present $h_{\mathrm{p}}(t)$ , $H_{\mathrm{p}}(t)$ , and $H_{\mathrm{p},\infty}$ along with $w_{\mathrm{e}}$ in the following theorem.

Theorem 1

The expected molecule hitting rate, the expected fraction of absorbed molecules, and the expected asymptotic fraction of absorbed molecules as $t\rightarrow\infty$ at a RX with multiple APs are given by

\displaystyle h_{\mathrm{p}}(t)=h_{\mathrm{u}}(t,w_{\mathrm{e}}),~H_{\mathrm{p}}(t)=H_{\mathrm{u}}(t,w_{\mathrm{e}}),~H_{\mathrm{p},\infty}=H_{\mathrm{u},\infty}(w_{\mathrm{e}}),

(17)

respectively, where

\displaystyle w_{\mathrm{e}}=\frac{D_{\sigma}G_{\mathrm{p}}}{r_{\scriptscriptstyle\mathrm{R}}(r_{\scriptscriptstyle\mathrm{R}}-G_{\mathrm{p}})},

(18)

with $G_{\mathrm{p}}$ given in (III-C). When all APs have the same size or there is only a single AP, $G_{\mathrm{p}}$ is given in (3) and (III-C), respectively.

Proof:

According to the relationship between the expected asymptotic fraction of absorbed molecules and the diffusion current, we have

\displaystyle\frac{H_{\mathrm{u},\infty}(w_{\mathrm{e}})|_{k_{\mathrm{d}}=0}}{H_{\mathrm{a},\infty}}=\frac{I_{\mathrm{p}}}{I_{\mathrm{a}}},

(19)

where $H_{\mathrm{u},\infty}(w_{\mathrm{e}})\big{|}_{k_{\mathrm{d}}=0}=r_{\scriptscriptstyle\mathrm{R}}^{2}w_{\mathrm{e}}/(r_{0}(w_{\mathrm{e}}r_{\scriptscriptstyle\mathrm{R}}+D_{\sigma}))$ and $H_{\mathrm{a},\infty}$ is the asymptotic fraction of absorbed molecules as $t\rightarrow\infty$ at a fully absorbing RX, given by $H_{\mathrm{a},\infty}=r_{\scriptscriptstyle\mathrm{R}}/r_{0}$ [14, Eq. (3.116)]. By solving (19), we obtain (18). By substituting $w_{\mathrm{e}}$ into (1), (4), and (9), we obtain (17). ∎

IV Numerical Results

In this section, we present numerical results to validate our theoretical analysis and offer insights for MC system design. Specifically, we use PBSs to simulate the random diffusion of the molecules, where we calculate the coordinates of the locations of molecules to be absorbed at the RX surface by using [18, Eqs. (13)-(15)]. If the locations of to-be-absorbed molecules are within the APs, we treat these molecule as molecules which have been absorbed. Otherwise, these molecules are reflected back to their positions at the start of the current simulation step [6]. To model the location of the TX in the simulation, we fix the distance between the TX and the center of the RX as $r_{0}$ and change the location of the TX for each realization. The TX is located at any point on the virtual sphere with the same probability. The simulation time step is $\Delta t=10^{-6}\;\mathrm{s}$ and all results are averaged over 1000 realizations. Throughout this section, we set $\mathcal{A}=\{0.05,0.1\}$ [10], $r_{\mathrm{\scriptscriptstyle\textnormal{R}}}=10\;\mu\mathrm{m}$ , $r_{0}=20\;\mu\mathrm{m}$ , $D_{\sigma}=79.4\;\mu\mathrm{m}^{2}/\mathrm{s}$ [19], $N_{\sigma}=1000$ [16], $k_{\mathrm{d}}=0.8\;\mathrm{s}^{-1}$ , and specify the value of $N_{\mathrm{p}}$ and the locations of the APs in each figure. In Fig. 2 and Fig. 3, we observe that the simulation results (denoted by markers) match well with the analytical curves (denoted by solid and dashed lines) generated based on Section III, which demonstrates the accuracy of our analysis.

In Fig. 2, we plot the expected molecule hitting rate at the RX at time $t$ in Fig. 2(a) and the expected number of absorbed molecules at the RX by time $t$ in Fig. 2(b). We set $\mathcal{A}=0.05$ and increase the number of APs from 1 to 11. In this figure, all APs have the same size and are evenly distributed over the RX surface. Here, we apply the Fibonacci lattice [20] to determine the locations of the evenly distributed APs. Specifically, the location of $\mathrm{AP}_{i}$ is given by $\theta_{\mathrm{p},i}=\pi/2-\arcsin(2(i-\mathcal{B}-1)/N_{\mathrm{p}})$ and $\varphi_{\mathrm{p},i}=(4\pi(i-\mathcal{B}-1))/(1+\sqrt{5})$ , where an integer $\mathcal{B}$ is given by $\mathcal{B}=(N_{\mathrm{p}}-1)/2$ and $N_{\mathrm{p}}$ is an odd number. In this figure, we observe that the expected hitting rate in Fig. 2(a) and the expected number of absorbed molecules in Fig. 2(b) increase with increasing $N_{\mathrm{p}}$ . This is because APs can absorb molecules from any direction to the maximum when there are a larger number of APs on the RX surface, leading to an increased expected number of absorbed molecules.

In Fig. 3, we vary the distributions and sizes of the APs on the RX surface. We set $\mathcal{A}=0.1$ and $N_{\mathrm{p}}=13$ for different AP distributions. Here, we consider three types of AP distributions, namely, APs evenly distributed on the RX surface, APs randomly distributed on the RX surface, and APs located within a region of the RX surface. For APs located within a region, we assume that APs are evenly distributed within $\theta_{\mathrm{p},i}\in[2.812,3.471]$ and $\varphi_{\mathrm{p},i}\in[0,2\pi]$ . From this figure, we observe that evenly distributed APs achieve the highest expected hitting rate and the largest expected number of absorbed molecules, compared to the other two types of distributions. This is because evenly distributed APs ensure that the APs are equally spaced in all directions of the RX. Thus, for most locations of the TX, the probability of molecules hitting APs is the highest when the APs are evenly distributed, compared to the other two types of distributions. Moreover, we also consider $N_{\mathrm{p}}=4$ to validate our analytical expressions when the APs have different sizes. We denote $\mathcal{A}_{\mathrm{p},i}$ as the ratio of the area of $\mathrm{AP}_{i}$ to the RX surface area and set $\mathcal{A}_{\mathrm{p},1}=0.01$ , $\mathcal{A}_{\mathrm{p},2}=0.02$ , $\mathcal{A}_{\mathrm{p},3}=0.03$ , and $\mathcal{A}_{\mathrm{p},4}=0.04$ . We also set the locations of the four APs to $\vec{l}_{1}=[10\;\mu\mathrm{m},\pi/2,\pi]$ , $\vec{l}_{2}=[10\;\mu\mathrm{m},\pi/2,\pi/2]$ , $\vec{l}_{3}=[10\;\mu\mathrm{m},\pi/2,0]$ , and $\vec{l}_{4}=[10\;\mu\mathrm{m},\pi/2,3\pi/2]$ .

In Fig. 4, we plot $w_{\mathrm{e}}$ versus $N_{\mathrm{p}}$ for different values of $r_{\scriptscriptstyle\textnormal{R}}$ and $D_{\sigma}$ to investigate the impact of $N_{\mathrm{p}}$ , $r_{\scriptscriptstyle\textnormal{R}}$ , and $D_{\sigma}$ on the effective reaction rate $w_{\mathrm{e}}$ . We note that a larger $w_{\mathrm{e}}$ means a higher probability of molecules hitting APs. We set $\mathcal{A}=0.05$ and consider evenly distributed APs on the RX surface. First, we observe that $w_{\mathrm{e}}$ increases when $N_{\mathrm{p}}$ increases, which indicates that the RX absorbs more molecules if there are more APs for a given $\mathcal{A}$ . This observation aligns with the observations in Fig. 2. Second, we observe that $w_{\mathrm{e}}$ increases when $r_{\scriptscriptstyle\textnormal{R}}$ decreases. This is because the AP-TX distances of different APs increase and become more similar when $r_{\scriptscriptstyle\textnormal{R}}$ decreases and $r_{0}$ is kept fixed. Given that $\mathcal{A}$ is relatively small, the probability of molecules hitting APs becomes higher when the AP-TX distances increase and are more similar. Third, we observe that $w_{\mathrm{e}}$ increases when $D_{\sigma}$ increases. This is because the RX can absorb more molecules when the diffusion coefficient of molecules increases, which leads to higher $w_{\mathrm{e}}$ .

V Conclusion

In this paper, we proposed a RX model that is covered by multiple APs with different sizes and arbitrary locations. By considering a point TX uniformly located on a virtual sphere at a fixed distance from the RX, we derived closed-form expressions for the expected CIR between the TX and RX in the presence of molecule degradation. Simulations verified the accuracy of our analytical expressions. Our numerical results showed that when the proportion of the total area of APs to the RX surface is kept fixed, the expected number of absorbed molecules increases with the number of APs. They also showed that evenly distributed APs yield the largest expected number of absorbed molecules. Future directions for research include 1) investigating the optimal spatial and area arrangements of the APs for maximization of the expected number of absorbed molecules and 2) replacing the point TX with a practical TX model to characterize cell-to-cell communication.

Appendix A Proof of Corollary 2

According to (4), we have

\displaystyle H_{\mathrm{u,\infty}}(w)=\frac{r_{\scriptscriptstyle\mathrm{R}}w}{r_{0}}\left[\alpha_{1}(\infty)-\alpha_{2}(\infty,w)\right].

(20)

Based on (1), we have $\alpha_{1}(\infty)=\exp(-\beta)/\sqrt{k_{\mathrm{d}}D_{\sigma}}$ . To obtain $\alpha_{2}(\infty,w)$ , we need to derive $\psi_{1}(\infty,w)$ and $\psi_{2}(\infty,w)$ in (6). In (1), if $\zeta(w)\leq 0$ , we have $\psi_{1}(\infty,w)=0$ . If $\zeta(w)>0$ , we obtain

		$\displaystyle\psi_{1}\!(\!\infty,\!w\!)\!\!=\!\!2\gamma(w)\!\exp(\!\gamma(w)(r_{0}\!-\!r_{\scriptscriptstyle\mathrm{R}}\!)\!)\!\!\lim_{t\rightarrow\infty}\!\!\frac{\mathrm{erfc}\!\left(\!\!\frac{\varepsilon}{\sqrt{t}}\!+\!\gamma(\!w\!)\sqrt{D_{\sigma}t}\!\right)}{\exp(-\zeta(w)t)}$
		$\displaystyle\overset{(a)}{=}\!\!2\gamma(w)\!\!\lim_{t\rightarrow\infty}\!\!\frac{\exp\!\!\left(\!-\!\!\left(\!\frac{\varepsilon^{2}}{t}\!+\!k_{\mathrm{d}}t\!\right)\!\!\right)}{-\zeta(w)}\!\!\left(\!\!\frac{\varepsilon}{\sqrt{t^{3}}}\!-\!\gamma(w)\sqrt{\frac{D_{\sigma}}{t}}\right)\!\!=\!\!0,$		(21)

where step $(a)$ is obtained by applying L’Hôpital’s rule. Therefore, we have $\psi_{1}(\infty,w)=0$ for any value of $\zeta(w)$ . According to (1), we obtain $\psi_{2}(\infty,w)=-2\gamma(w)^{2}\sqrt{{D_{\sigma}}/{k_{\mathrm{d}}}}\exp\left(-\beta\right)$ . By substituting $\psi_{1}(\infty,w)$ and $\psi_{2}(\infty,w)$ into (6), we obtain $\alpha_{2}(\infty,w)$ . Then, by substituting $\alpha_{1}(\infty)$ and $\alpha_{2}(\infty,w)$ into (20), we obtain (9).

References

[1] N. Farsad, H. B. Yilmaz, A. Eckford, C.-B. Chae, and W. Guo, “A comprehensive survey of recent advancements in molecular communication,” IEEE Commun. Surv. Tutor., vol. 18, no. 3, pp. 1887–1919, 3rd Quarter, 2016.
[2] I. Guryanov, S. Fiorucci, and T. Tennikova, “Receptor-ligand interactions: Advanced biomedical applications,” Mater. Sci. Eng. C, vol. 68, pp. 890–903, Nov. 2016.
[3] V. Jamali, A. Ahmadzadeh, W. Wicke, A. Noel, and R. Schober, “Channel modeling for diffusive molecular communication—a tutorial review,” Proc. IEEE, vol. 107, no. 7, pp. 1256–1301, Jul. 2019.
[4] B. C. Akdeniz, N. A. Turgut, H. B. Yilmaz, C.-B. Chae, T. Tugcu, and A. E. Pusane, “Molecular signal modeling of a partially counting absorbing spherical receiver,” IEEE Trans. Commun., vol. 66, no. 12, pp. 6237–6246, Dec. 2018.
[5] A. Akkaya, H. B. Yilmaz, C.-B. Chae, and T. Tugcu, “Effect of receptor density and size on signal reception in molecular communication via diffusion with an absorbing receiver,” IEEE Commun. Lett., vol. 19, no. 2, pp. 155–158, Feb. 2015.
[6] A. Ahmadzadeh, H. Arjmandi, A. Burkovski, and R. Schober, “Comprehensive reactive receiver modeling for diffusive molecular communication systems: Reversible binding, molecule degradation, and finite number of receptors,” IEEE Trans. Nanobiosci., vol. 15, no. 7, pp. 713–727, Oct. 2016.
[7] J. Sun and H. Li, “Expected received signal in diffusive molecular communication with finite binding receptors,” IEEE Commun. Lett., vol. 24, no. 12, pp. 2829–2833, Dec. 2020.
[8] S. Lotter, M. Schäfer, J. Zeitler, and R. Schober, “Saturating receiver and receptor competition in synaptic DMC: Deterministic and statistical signal models,” IEEE Trans. Nanobiosci., Jun. 2021.
[9] T. Duke and I. Graham, “Equilibrium mechanisms of receptor clustering,” Prog. Biophys. Mol. Biol., vol. 100, no. 1-3, pp. 18–24, Sep. 2009.
[10] A. E. Lindsay, A. J. Bernoff, and M. J. Ward, “First passage statistics for the capture of a Brownian particle by a structured spherical target with multiple surface traps,” Multiscale Model. Simul., vol. 15, no. 1, pp. 74–109, Jan. 2017.
[11] A. D. Lander, “How cells know where they are,” Science, vol. 339, no. 6122, pp. 923–927, Feb. 2013.
[12] R. Chang, Physical Chemistry for the Biosciences. Sausalito, CA, USA: Univ. Science Books, 2005.
[13] L. Dagdug, M.-V. Vázquez, A. M. Berezhkovskii, and V. Y. Zitserman, “Boundary homogenization for a sphere with an absorbing cap of arbitrary size,” J. Chem. Phys., vol. 145, no. 21, p. 214101, Dec. 2016.
[14] K. Schulten and I. Kosztin, Lectures in Theoretical Biophysics. Urbana, IL, USA: Univ. of Illinois Press, 2000, vol. 117.
[15] H. C. Berg, Random Walks in Biology. Princeton, NJ, USA: Princeton Univ. Press, 1993.
[16] A. C. Heren, H. B. Yilmaz, C.-B. Chae, and T. Tugcu, “Effect of degradation in molecular communication: Impairment or enhancement?” IEEE Trans. Mol. Biol. Multi-Scale Commun., vol. 1, no. 2, pp. 217–229, Jun. 2015.
[17] H. C. Berg and E. M. Purcell, “Physics of chemoreception,” Biophys. J., vol. 20, no. 2, pp. 193–219, 1977.
[18] X. Huang, Y. Fang, A. Noel, and N. Yang, “Membrane fusion-based transmitter design for molecular communication systems,” in Proc. IEEE ICC 2021, Montreal, QC, Canada, Jun. 2021, pp. 1–6.
[19] H. B. Yilmaz, A. C. Heren, T. Tugcu, and C.-B. Chae, “Three-dimensional channel characteristics for molecular communications with an absorbing receiver,” IEEE Commun. Lett., vol. 18, no. 6, pp. 929–932, Jun. 2014.
[20] Á. González, “Measurement of areas on a sphere using Fibonacci and latitude–longitude lattices,” Math. Geosci., vol. 42, no. 1, pp. 49–64, Nov. 2010.