Chain Modeling of Molecular Communications for Body Area Network

Molecular communications provide an attractive opportunity to precisely regulate biological signaling in nano-medicine applications of body area networks. In this paper, we utilize molecular communication tools to interpret how neural signals are generated in response to external stimuli. First, we propose a chain model of molecular communication system by considering three types of biological signaling through different communication media. Second, communication models of hormonal signaling, Ca2+ signaling and neural signaling are developed based on existing knowledge. Third, an amplify-and-forward relaying mechanism is proposed to connect different types of signaling. Simulation results demonstrate that the proposed communication system facilitates the information exchange between the neural system and nano-machines, and suggests that proper adjustment can optimize the communication system performance.


Introduction
Rapid progress of nano-technology enables the manufacturing of nano-machines for medical applications of body area network (BAN) [1]. Thus far, a bio-inspired manufacturing method has attracted considerable attention. It assembles biological molecules or atoms into a micro-to nano-scale unit based on the functional structure of nature cells. To finish a complex task, the cooperation of nano-machines is indispensable due to their limited size and power. Molecular communication (MC) is a promising technology for the cooperation of nano-machines, exploiting biological molecules as intra-body information carriers [2]. Design of MC is inspired by mechanisms of biological activities in nature (e.g., hormone communication between ants and neural signaling inside human body), which can be promoted by the experience from diverse traditional communication systems [3,4]. Promising MC applications include health monitoring systems, drug delivery systems and bio-sensor networks [5].
MC could be divided into wireless and wired types based on different properties of biological mechanisms. Existing research efforts generally focus on a single MC type and lack the comprehensive understanding of the complex biological activities in human bodies. Biologically speaking, the human body is a complex system composed of various biological substances including organs, tissues, cells, molecules, etc. It is widely accepted that, even for a simple biological activity, different substances interact with each other through exchange of power and substance. Moreover, processes of biological activities usually involve multiple types of biological signaling, which should be investigated jointly

•
First, we develop a chain model of molecular communication based on biological signaling. The proposed model considers the biological interactions among hormone, Ca 2+ and neural signals, which is more general than the model of [8]. • Second, we propose an implementable amplify-and-forward relaying mechanism instead of decode-and-forward relaying as in [8]. In addition, multiple astrocytes are utilized instead of a nano-machine, elevating the reliability of relaying from Ca 2+ signaling to neural signaling. • Third, based on the work in [8], we examined the relations between communication performance and more parameters of the proposed model. We also found that source coding is efficient in improving the communication performance, which may provide a guidance for nano-machine design.

Overview of the Communication System
The proposed MC system is illustrated in Figure 1. The proposed communication system is motivated by the experience of multiple traditional communication systems [18][19][20][21][22][23]. Specifically, the communication system is modeled as a chain, which is composed of a sender, a receiver and three types of communication media. The sender and receiver are nano-machines with different functions, respectively, denoted by Tx and Rx. Tx is assumed to release hormones in response to external stimuli. Tx can be, for example, engineered endocrine cells [2], which are able to release hormones under control. Rx is a bio-inspired or electronic nanomachine that is able to detect neural signals. The first type of communication medium is the fluid, where hormones propagate around the astrocytes. There are two types of propagation mechanisms, determined by types and properties of the hormones and fluid medium. First, hormones perform approximative random walk towards any direction with low energy consumption in the passive diffusion mechanism. For example, steroid hormones propagate in the brain via lipid mediator and passive diffusion [24]. Second, hormones drift towards specific directions with high energy consumption in active diffusion. For example, thyroid hormones permeate blood-brain barrier by a slow flux rate in bloodstream via active diffusion [25].
The second type of communication medium is astrocytes, which are specific cells connected with neurons. Astrocytes are closely related to neural activities, such as physically supporting structures of neurons, providing nutrition for neurons, and promoting the growth of neurons [26]. Astrocytes can promote the generation of neural signals through adjustment of Ca 2+ oscillation, thereby serving as nature bridges between another biological system and neural system. Reportedly, astrocytes can respond to hormones (e.g., norepinephrine) in generation of Ca 2+ signals [27]. Based on this knowledge, astrocytes function as the perfect relay in our proposed model, connecting hormonal signaling with neural signaling.
The third type of communication medium is neuron, which processes and transmits information through the human body. Neural communication is an efficient and reliable biological signaling, exploiting both electrical and chemical methods. Different types of neurons respond to different external stimuli (e.g., electricity, chemical and light) [28].
As shown in Figure 1, the communication system is developed to interpret the neural signaling evoked by chemical stimuli. We assume that the proposed biological communication system transmits information by encoding the presence or absence of signals in time scale. Moreover, we assume that time-slotted mechanism works for the proposed communication system, based on the biological clocks guaranteeing the synchronization of biological cells [29].
The communication process is described as follows. Initially, Tx releases hormones in a controlled pattern, responding to external stimuli. Then, hormones diffuse in the fluid medium, some of which are absorbed by the astrocytes, inducing the oscillation of Ca 2+ signaling. Oscillated Ca 2+ ions in astrocytes flow into the neurons and trigger neural signaling. Finally, the neural signals are detected and decoded by Rx.

Hormonal Signaling Model
In Tx, a list of binary information is encoded with the number of released hormone molecules per time slot. One bit is transmitted in a symbol period, denoted by T. It is found that secretion pattern of hormones is usually with burst and not sustaining [30]. Accordingly, we apply a simple coding for MC, i.e., on-off key (OOK) with square pulses. During a time of T, a certain number of hormone molecules Q 0 are released into fluid medium, if Tx transmits bit "1"; no hormone molecules are released into fluid medium, if Tx transmits bit "0". Let R H (t) be the released molecules during the ith slot, the expression of OOK coding is given by, where δ(t) is the Dirac's delta function of t. Hormone molecules propagate with random or directed walks, motivated by the passive or active diffusion mentioned in Section 3. Plenty of factors impact the distance and velocity of diffusion processes, including molecules diameter, medium type, temperature and pressure. For simplicity, we only consider passive diffusion of hormones, i.e., hormones performing random walks. Concentration of the hormone molecules presents a decreasing gradient with the distance, which could be expressed by the well-known Fick's second law, where ∇ 2 is the Laplace operator, indicating the sum of the second derivatives of C H (x, t). R H (t) is the releasing rate of molecules at Tx in time scale. D H is the diffusion coefficient of hormones, indicating the diffusion rate. Let x S be the position of Tx, and x A denote the position of an astrocyte. By solving Equation (2), the concentration accumulated at the position of x A is given by, where g(x, t) is the Green function of hormones [9], expressing the concentration distribution of hormone molecules in space and time. d S,A denotes the distance between Tx and an astrocyte, given by d S,A = |x A − x S |. A single hormonal molecule could be released by Tx in the n S th slot and absorbed by an astrocyte in the n A th slot. The corresponding probability is given by [31], where n = n A − n S , and m H (v) denotes the probability function distribution of the molecular life expectancy, given by m H (v) = ιe −ιv with mean of 1/ι.

Calcium Signaling Model
Astrocytes are special star-shaped glial cells, physically close to neurons, and play a crucial role in regulation of neural activities. It is observed that astrocytes propagate Ca 2+ waves in response to external stimuli, which further induce and impact the behaviors of neural signaling [32]. In this section, we explain how Ca 2+ signals are generated and propagate in astrocytes.
We introduce a mathematical model of Ca 2+ signaling in astrocytes according to experimental observations [27]. The dynamic of Ca 2+ signaling in an astrocyte is shown in Figure 2. One major organelle that stores Ca 2+ ions is endoplasmic reticulum (ER). With an external stimulus in astrocyte, Inositol 1,4,5-trisphosphate (IP 3 ) molecules are generated by G-protein coupled receptors (GR), changing the permeability of IP 3 -sensitive Ca 2+ receptors (IP 3 R) of the ER and induce the Ca 2+ fluxes. Three major Ca 2+ fluxes are considered, namely channel fluxes released from the ER to the cytoplasm via receptors of IP 3 (U chan ), uptake fluxes from the cytoplasm to the ER by Ca 2+ pumps (U pump ), and leak fluxes from from the ER to the cytoplasm (U leak ). Induced Ca 2+ fluxes incur the oscillation of Ca 2+ , and patterns of Ca 2+ waves are positively related to the concentration of external stimuli. Let C CY and C ER denote the Ca 2+ concentrations in cytoplasm and ER. Variable h ∞ is employed to denote the Ca 2+ channel inactivation gate. On the time scale, the introduced model of Ca 2+ signaling is expressed as a set of equations, where, In the above equations, a 1 denotes the ratio of ER to cytoplasm, a 2 denotes the time constant of the inactivation gate, and a 3 denotes the steady value of the inactivation gate. v chan , v leak and v pump denote the maximum flux rates of U chan , U leak and U pump , respectively. b 1 -b 6 are constants of the model, and C IP3 denotes the concentration of IP 3 .
On the spatial scale, the propagation of Ca 2+ ions in astrocytes is based on the diffusion of Ca 2+ [12]. ER is modeled as a point source of releasing Ca 2+ . Ca 2+ ions around the boundary of astrocytes can flow into the connected neurons. For any position x inside an astrocyte, the dynamic of the Ca 2+ concentration via diffusion can be expressed by, where D Ca denotes diffusion coefficient of Ca 2+ smf R Ca (t) denotes the variation rate of Ca 2+ concentration, determined by Equations (5)- (13). Equation (14) can be solved similarly to Equation (3).

Neural Signaling Model
Neural communication is a fast biological signaling, and its velocity varies due to different neuron types. Neural communication is also biologically reliable, as the neurons are able to eliminate errors in communications. In this section, we introduce a signaling model across two neighboring neurons, including the processes of neural firing, axonal transmission, gap junctional transmission, and postsynaptic response.
One neuron is normally composed of different number of dendrites, axons, somas, and axonal terminals. In a communication, neural signals are first generated by ions through neural firing process in soma, transmitted along the axon fiber in the pattern of electrical impulses. Then, the electrical impulses trigger the release of neurotransmitters from the vesicles located in axon terminals. The neurotransmitters propagate to the neighboring neurons, until they are absorbed by the dendrites. Finally, the postsynaptic response is induced that might trigger firing of neighboring neurons.

Amplify-and-Forward Neural Firing
The generation of neural signals is called neural firing, which is an all-or-none process. Various types of ions including Ca 2+ are released by astrocytes and flow into neurons. During the process, local potentials of neuron membrane are gradually elevated and aggregated into a global potential. If the global potential exceeds a threshold, neural signals are generated. Otherwise, it fails to generate neural signals since the global potential falls quickly.
Neural firing is a relaying process, in which information carrier varies from Ca 2+ to neural signals. In our model, we apply amplify-and-forward (AF) relaying in neural firing process. AF relaying is based on the concentration-mediated mechanism in astrocyte-neuron communications [33]. We assume that local neural signals are positively elevated according to the concentration of Ca 2+ without decoding.
To generate neural signals, there must be adequate Ca 2+ ions flowing into neurons. The number of astrocytes should be sufficient to generate the required intensities of Ca 2+ flows. The local potentials are elevated by Ca 2+ from connected astrocytes, which cooperate for achieving neural firing. This process is expressed by the well-known Hodgkin-Huxley model [34], given by, where M A is the total number of astrocytes contributing to the neural firing, V N (t) denotes the global potential on the neural membrane during the firing process, V r is the resting global potential without any stimuli, ν is the linear response of the neural membrane to an input pulse of Ca 2+ flow, A N,i (t) denotes the intensity of Ca 2+ flow released by the ith astrocyte, injected into the neuron at time t i , and A N,i (t) is proportional to Ca 2+ concentration of astrocyte, expressed by A N,i (t) ∝ C CY (t).
The action potential spikes are typical neural signal S(t), generated when V N (t) is strong enough to exceed firing threshold θ 1 , where ϕ(t) denotes the waveform of action potential spike, and t j denotes the time that the spike occurs. Based on experimental observations, S(t) follows or approaches Poisson distribution in response to random intensities of Ca 2+ flows [35]. We employ λ(t) to denote the firing rate, which is determined by not only the Ca 2+ flows, but also the inherent capability of neurons. The distribution of S(t) could be impacted by the noise or some other stochastic factors.

Axonal Transmission
In the process of axonal transmission, action potential spikes transmit along the axon until they arrive in the axon terminals. During the process, various ions (e.g., K + , Na + and Cl − ) exchange between the inside and outside of the axon fiber. The axonal transmission process could be modeled by the cable theory, expressed by the quantitative equation in [36], where γ m , γ l and γ c denote the membrane resistance, longitudinal resistance and capacitances, respectively. V N (x, t) denotes the potential of axonal fiber at time t and location x. The boundary condition of the axonal transmission is given by, The above equations can be solved using the Fourier method [37], namely V N (x, t) is calculated as, where L N is the length of the axon fiber. ζ i denotes the Fourier coefficient, following the boundary conditions in Equation (18), and can be obtained by, 6.3. Gap Junctional Transmission

Vesicle Release
We show the process of gap junctional transmission in Figure 3. Action potential spikes arrive in the axonal terminals, triggering vesicle movements and inducing vesicle-membrane fusion. As a result, neurotransmitters stored inside vesicles are released into the gap junction. The release process of neurotransmitters is stochastic, related to the vesicle-membrane distance [38]. Therefore, we exploit a pool model to classify vesicles into readily releasable pool (RRP) and reserve pool (RP), respectively. The vesicles of RP are difficult to release, since they are far from the membrane, whereas the vesicles of RRP are easy to release, since they are close to the membrane.
Let N RRP and N RP , respectively, denote the sizes of RRP and RP. The vesicles of RP are assumed to fill the RP with rate of τ f , if N RRP is reduced due to the release process. The average firing rate of vesicles is denoted byλ, which is calculated as the integral of firing rate λ(t) during a time slot, The release probability of a vesicle is related toλ and N RRP , given by, In local scale, P rele can be approximated as the linear relation withλ, i.e., P rele =λN RRP . Due to the stochastic behavior of neurotransmitters, it is possible to cause errors during the vesicle release process. For binary information, the error probabilities of bit "1" and "0" are, respectively, calculated by, where k i is the number of action potential spikes during the ith slot. Equations (23) and (24) indicate that error transmission of binary information is caused by miss or redundant release of vesicles.

Queueing Model of Neurotransmitters
The neurotransmitters diffuse in the gap junctions until absorbed by the dendrite receptors of neighboring neurons. Due to the stochastic behavior of the vesicle release process, the absorbed time of neurotransmitters is random. The diffusion of neurotransmitters is slow; this process is normally ignored since gap junctions are extremely short (several nanometers). However, neurotransmitters can be absorbed with delay due to the limited absorbing ability of dendrites [39]. Generally, neurotransmitters gather around the receptors of dendrites and wait to be absorbed. To study this phenomenon, we use the queuing theory to calculate the delay of gap junctional transmission.
In our model, the neurotransmitters are defined as the batched customers, which are independent with identical distributions. We adopt the M x |G|1|K model, where M x is the absorbing time of these batched neurotransmitters. The time intervals of vesicle release follow the negative exponential distribution, i.e., f (t i ) = σλ b exp(−σλ b t i )(0 < σ < 1). µ b and τ b , respectively, denote the batch size and batch service rate, satisfying λ b = µ b τ b . K denotes the capacity of the receptor. G denotes absorbing time t s of the receptor with a distribution G(t), which is given by 0 < t s = 1 µ = ∞ The distribution of the batch sizes is expressed as (x) = (1 − σ) x−1 σ. Based on the result of [40], the blocking probability is given by, where p T (k) is the probability that neurotransmitters are absorbed by k customers. The blocking phenomenon is more serious for bigger λ b . Based on the previous derivations, the average delay is calculated as, where g 0 (x, t) is the probability that customers are served at time t.

Postsynaptic Response
The process of postsynaptic response includes the formation of postsynaptic potential and the decoding of neural signals.

Postsynaptic Potential
The neurotransmitters are absorbed by dendrites of neighboring neurons via various receptors, such as the AMPA and NMDA receptors. Neurotransmitters accumulate on the dendrite membrane, namely postsynaptic terminals. During this process, the postsynaptic potential of dendrite membrane is elevated to generate the waveform of output signals, which could be expressed by an alpha function [41], where G max is the maximum amplitude of the potential, and t p is the corresponding time to peak.

Neural Decoding
Decoding process of the neural signal includes two steps: (1) estimation of the action potential spikes; and (2) extraction of the binary information.
In the first step, a waveform of output signal is denoted by y(t), which could be calculated as a convolution of the action potential spikes and (t), where q denotes a stochastic variable in the signal amplification, following the Gamma distribution [42].
To extract the binary information of y(t), we first sample y(t) into discrete signal y(n) with interval τ s , where N sam is the number of samplings in waveform of y(t), which depends on G max and t p . Let t p be the origin point; N sam 2 sampling points are determined on both sides of the origin point, where kτ s < min(t j − t j−1 ) is satisfied. The sampling interval τ s should be properly regulated to alleviate the distortion of sampling signals. We adopt threshold ϑ 2 to determine whether a spike exists or not at the position of sampling. If ∑ k n=1 y(n) k < ϑ 2 , a spike exists. Accordingly, we obtain a pattern of estimated spikes, denoted by y (t).
In the second step, we extract the binary information by examining the distribution of y (t). The Kolmogorov-Smirnov test [43] is employed to verify whether the distribution of y (t) follows the Poisson distribution in time slot T. We say y (t) follows the Poisson distribution if the following condition is satisfied, where variable K α meets condition P(K K α ) = 1 − α, K denotes the Kolmogorov distribution, and κ i is the Kolmogorov-Smirnov statistic of the Poisson process during the ith slot. We obtain λ using the Maximum Likelihood Estimation (MLE). Threshold ϑ 3 is employed for decoding, which is smaller than λ. For the ith bit, if the condition in Equation (30) is satisfied and estimated λ is higher than ϑ 3 , the decoded information is bit "1". Otherwise, the decoded information is bit "0". Further, we deduce the transmission error probabilities of decoding for both bit "1" and "0", According to Equations (23), (24), (31) and (32), the error probability of transmitting one bit neural signal is expressed as,

Channel Capacity and Transmission Delay
In this section, we deduce the expressions of channel capacity and transmission delay for the proposed communication system.

Channel Capacity
Let X i , Y i and Z i , respectively, represent the input, relay and output signal, indicating the diffusive concentration of hormone, oscillating concentration of Ca 2+ , and decoded action potential spikes. To deduce the channel capacity, we define the following equations, In the above equations, P (1,1) i and P (1,2) i are related to the probability of transmitting bit "1", denoted by p. According to the conclusions in [31], the expressions of P (1,1) i and P (1,2) i are, respectively, given by, where u i is the probability of Equation (4), and N b is the total number of time slots. The expressions of P[Z i = 1|Y i = 0] and P[Z i = 0|Y i = 1] can be calculated according to Equations (33) and (34). Based on information theory, the channel capacity is expressed as, We assume that, under the condition of Y, X is uncorrelated to Z. During the ith time slot, the mutual information is given by, The mutual information between X and Y is further calculated as, . The second item of Equation (44) is further expanded as,

Delay
Transmission delay of the proposed communication system is denoted by Ω, which is calculated as the delay summation of hormonal signaling Ω H , Ca 2+ signaling Ω Ca and neural signaling Ω N , The estimation of Ω H is related to the emission, propagation and reception processes of hormones. For simplicity, we only consider the propagation process in channel for the transmission delay calculation. The transmission delay is related to two factors on the scales of time and space. The first factor is the transmission distance on the spatial scale. One feasible method is to estimate based on the probability distribution function along propagation distance d S,A , given by, where er f (x) is the error function, expressed as er f (x) = 2 √ π x 0 e −u 2 du. Another factor determining transmission delay is the variation rate of molecules at the sender on the time scale. Due to the fixed number of molecules Q 0 in one time slot, the variation rate of molecules is determined by the length of each time slot. The average delay of hormonal diffusion is calculated as, The estimation of Ω Ca is based on the time difference, which starts with the generation of Ca 2+ waves in connected astrocytes, ends with the absorbing of Ca 2+ by a neuron. Similar with the estimation of Taynnan et al. [13], the variational concentrations of Ca 2+ in the astrocytes are assumed to be equal with the absorbed concentrations of Ca 2+ during an increasing delay of Ω Ca [13], i.e., where ∆(x) denotes the variation of x during a time period, C i CY (t) denotes the Ca 2+ concentration in cytoplasm of the ith connected astrocyte with the neuron, and C N (t) is the variational concentration of Ca 2+ absorbed by the neuron.
In the calculation of the neural signaling delay, we mainly consider the contribution of gap junctional transmission. Based on the proposed queuing model of gap junctional transmission, the average delay of the vesicle release process is approximated as the average waiting time ω of the queueing model. The delay of axonal transmission Ω Ax is calculated based on Equation (19). Let N N be the number of neurons, and the average delay between two neighboring neurons is expressed as,

Performance Evaluation
A computer simulation was employed to examine the performance of proposed communication system. The control variate method was exploited since there are too many parameters in the proposed system. We verified the relations between some important parameters and the performance indicators (i.e., channel capacity and transmission delay).

Simulation Design
Simulation parameters were grouped according to hormonal, Ca 2+ and neural signaling, respectively.
Parameters of the hormonal signaling are listed in [9,31,44]: diffusion coefficient D H was fixed at 4.8 µm 2 /ms, number of released molecules per bit Q 0 varied from 1000 to 10,000, the distance between Tx and astrocytes were randomly generated within 1-20 µm, the length of time slot T varied from 1 s to 20 s, molecule life expectancy ι was set to 0.2, and probability p of transmitting bit "1" varied between 0 and 1.
Parameters of the Ca 2+ signaling are listed in [27]: ratio of ER to cytoplasm a 1 was 0.185, the maximum released flux rate from ER to cytoplasm v chan was 6 s −1 , the maximum leaked flux rate from ER to cytoplasm v leak was 0.11 s −1 , the maximum pumped flux rate from cytoplasm to ER v pump was 0.9 µMs −1 , b 1 was 0.13 µM, b 2 was 0.08234 µM, b 3 was 0.1 µM, b 4 was 0.9434 µM, b 5 was 1.049 µM, b 6 was 0.2 µM −1 s −1 , diffusion coefficient D Ca was 20 µ 2 s −1 , the average propagation distance of Ca 2+ signaling d ER,B was 10 µm, and the length of time slot T varied from 1 to 20 s.
Parameters of neural signaling is listed in [15,36,38]: firing rateλ of neural signaling varied from 15 to 30 Hz, the average release probability of vesicles p v was set to 0.4, pool size of readily release pool N RRP varied from 3 to 10, the pool size of reserve pool N RP was set to 60, and the average filling time τ f varied from 10 to 40 ms. In axonal transmission, membrane resistance γ m was set as 64.1 MΩ, longitudinal resistance γ l was set as 8 MΩ, and capacitances γ c was set to 1 mF/cm 2 , the average length of axons L N was set to 43.2 µm, and the length of time slot T varied from 1 to 20 s.

Behaviors of Signals
In the simulation, we present how various signals behave in 40 s to transmit 4 bits of binary information. Figure 4a,d,g presents hormonal signals, Figure 4b,e,h presents Ca 2+ signals, and Figure 4c,f,i presents neural signals. The results reveal that the hormonal and Ca 2+ signals propagate with low frequency, and the neural signals propagate with high frequency.  We analyzed the behaviors of various signals by sending different bits. Comparing Figure 4a,d,g, we see that concentrations of hormones are elevated by the percentage of sending bit "1". Similarly, the frequency of Ca 2+ oscillation is also related to percentage of sending bit "1". In Figure 4h, there are totally four peaks, and the intervals between different peaks are decreasing, since the frequency of Ca 2+ oscillation is regulated. Moreover, Ca 2+ oscillation stops after 30 s, when hormonal concentration is too high, and the concentration of Ca 2+ stabilizes at a high level. In addition, the number of action potential spikes seems positively related to the oscillation frequency of Ca 2+ .

Channel Capacity
Channel capacity of the proposed communication system is jointly determined by the three types of signaling. Our target was to check the relations between channel parameters and the channel capacity. Figure 5 shows how channel capacity C changes with different Q 0 , p and λ. Apparently, C is elevated by increasing number of released molecules Q 0 , because that astrocytes can easily absorb hormones due to the increase of Q 0 and the error probability of hormonal signaling is decreased. In addition, C reaches its peak if p is closer to 0.5. By comparing Figure 5a-c, we note that increase of firing rate λ of neural signaling restrains C, because a larger λ increases the error probability of vesicle release process.  Figure 6 shows how channel capacity C changes with different M A , p, N RRP and τ f . Obviously, the optimal p of maximizing C is not 0.5, which is different with traditional wireless communication. Thereby, the source coding of MC should be particularly designed. Increasing the number of astrocytes M A also elevates C, since the relaying capability of Ca 2+ signaling is enhanced. Moreover, Figure 6b reveals that C is larger with a larger readily release pool N RRP and a shorter filling time τ f .

Transmission Delay
We checked the relations between various channel parameters and the transmission delay. Figure 7a presents how transmission delay of hormonal signaling Ω H changes with different T and d S,A . Apparently, increasing the length of time slot T results in the growing Ω H , because releasing rate is decreased by increasing T and a fixed number of molecules Q 0 is released for each bit. The increasing rate of Ω H degrades fast when T is large. Therefore, we deduce that it is important to choose a proper T in the proposed model. Besides, we note that Ω H increases with growing distance d S,A between astrocyte and Tx, since the molecular concentration degrades significantly on the spatial scale. Figure 7b presents how Ω Ca changes with different T and M A . Similar with Ω H , the increased length of time slot T signifies the increased Ω Ca . We guess that the variation of T impacts the hormonal signaling and Ca 2+ signaling. Moreover, a larger M A results in a smaller Ω Ca , because the proposed system could transmit more Ca 2+ in a unit time, with the stronger relaying ability. Figure 7c presents how the average waiting time ω in neural signaling changes with different N RRP , τ f andλ. The values of ω is much smaller than Ω Ca and Ω H . We see that ω increases with firing rate λ, because firing rate promotes the vesicle release process, more neuro-transmitters are blocked in gap junctions, and ω increases. In addition, ω grows with the increasing size of readily release pool N RRP and the decreasing filling time τ f , since vesicle release process is promoted by those parameters.

Conclusions
In this paper, we design a chain model for molecular communication to understand and interpret neural signaling. The proposed model contains three types of signals, where Ca 2+ signals of astrocytes function as the relay of communication. Due to this advantage, the proposed model becomes more implementable. Behaviors of the different signals demonstrate their close relations with each other. Simulation results validate the effectiveness of the proposed chain model, and also reveal that proper source encoding helps improve the system performance.
Our work provides a philosophy of exploiting complex biological reactions for communication engineering in body area network. Future work will address the expansion of chain modeling and molecular communications in more complex networks, and explore more practical and complex biological activities.

Conflicts of Interest:
The authors declare no conflict of interest.