Information and communication technology (ICT) framework provides a novel perspective to fight human diseases [1] - [3] . In this respect, molecular communication could pave the way for a solution to develop therapeutic and diagnostic platforms. Recent Severe Acute Respiratory Syndrome-CoronaVirus 2 (SARS-CoV2) pandemic have resulted in a [4] , [5] . In [4] , the authors model Dengue virus transmission inside the body from its entrance to the host to the transmission to affected organs. The channel considered, which is from skin to the receiver organs, is characterized in terms of noise sources and path loss. Aerosol transmission, in which droplets carry virus, is the another means of virus transport mechanism. In [5] , the authors determine the aerosol channel impulse response and find the response of their system for the sources such as breathing, coughing and sneezing. On the other hand, a study considering the SARS-CoV2 transmission process through the human respiratory tract from molecular communication perspective is is yet to be studied in the literature.
SARS-CoV2 enters the host human through the nose, mouth and eyes. We consider the case that droplets carrying viruses enter the host human from the nose. Viruses travel via mucus flow in the respiratory tract and reach host cells as illustrated in Fig. 1 . SARS-CoV2 virus binds a special receptor on the host cell called angiotensin-converting enzyme or ACE2. Binding is followed by a time delay, τ , which is due to the mechanisms needed for virus replication. In this study, we consider this system and accordingly develop a model for the human respiratory tract by separating the respiratory tract into seven segments. Our aim is to determine the impulse response of the SARS-CoV2-ACE2 binding process to investigate the probability distribution of binding locations. The binding location distribution, which depends on several system parameters including ACE2 density and mucus flow rate offers vital information on the course of disease.
Our contributions can be summarized as follows:
• Proposing a novel model of human respiratory tract that reduces complexities of the original system: We model human respiratory tract by partitioning the tract into seven segments from nasal cavity to alveoli. • Determining impulse response of SARS-CoV2 infection process for the first time in literature • Calculating ACE2 receptor densities in the different regions of the respiratory tract: Based on the available data on surface parameters, we calculate ACE2 receptor density crudely. • Investigating the effects of mucus layer thickness, mucus flow rate and ACE2 density on the virus population reaching the different regions of the respiratory tract: Our results shows that mucus flow rate and ACE2 densities affect the respiratory regions where the viruses reach drastically. The rest of the paper is organized as follows. In Section II, we provide a brief background about SARS-CoV2. In Section III, the developed system model is outlined. In Section III, the diffusion model for viruses diffusing through the mucus layer is derived. Next, in Section V, the impulse response of the system for different receptor and virus concentration is determined. In Section VI, Markov Chain model of the events following the binding process are stated. In Section VII, the simulation results are presented. Finally, conclusions are stated in Section VIII.
Severe Acute Respiratory Syndrome -CoronaVirus 2 (SARS-CoV2), also named novel-coronavirus (2019-n-Cov), has been identified as the causative infectious agent of coronavirus disease- , responsible for the current pandemic. Covid-19 has turned from a local pneumonia outbreak, which originated in Wuhan, China in December 2019, into a global pandemic in a matter of months, which has as of now, October 2020, caused more than a million deaths worldwide and spread to more than 200 countries [6] . Belonging to the family of coronaviruses, SARS-CoV2 is the third and the newest coronavirus in the family to cause an epidemic, just as SARS-CoV in 2003 and MERS-CoV in 2012, and the only one to cause a pandemic. SARS-CoV2 is reported to be a zoonotic viral disease. Bats, snakes, and pangolins have been cited as potential reservoirs based on genome sequencing studies [7] - [9] .
Although it predominantly causes pneumonia and associated comorbidities, Covid-19 is considered to be a syndrome, given that it affects multiple different organs and systems within the human body. Typical clinical symptoms of the patients include fever, dry cough, difficulty of breathing (dyspnea), fatigue, joint pain (arthralgia), muscle pain (myalgia), and loss of sense of smell (anosmia) [10] - [12] . The presence of high variety of pathological events are attributed to different pathophysiological mechanisms involved in SARS-CoV2 and proves that it is more than a respiratory syndrome.
Current epidemiological data suggests that SARS-CoV2 is an airborne viral disease, meaning that it is transmitted through respiratory droplets and droplet nuclei, which are mostly spread during human-to-human contact [13] - [15] . Respiratory droplets (> 5−10µm in diameter) and droplet nuclei (aerosols) (< 5µm in diameter ), are generated and expelled/disseminated from an infected person during speaking, shouting, coughing, or sneezing [16] . Indirect surface transmission, i.e., fomite transmission, and orofecal transmission have also been reported [17] - [19] . Some studies have detected stable SARS-CoV2 viral RNA on solid surfaces such as plastic, aluminum, and stainless steel, yet the significance of fomite transmission is still debated with contradicting views [17] , [18] .
The main pathway of SARS-CoV2 inside the human host is reported to be the respiratory tract. Mucosal openings such as the nose, eyes, or mouth have been identified as the principal sites, where the initial viral entry takes place [20] . Although there are numerous possibilities for viral entry, one pathway a virus particle could take on the macroscopic level is as follows. A virus laden particle enters through the nasal cavity, with the help of the downstream flow of mucosal secretions and gravity, it travels down through the pharynx, larynx, and trachea, enters a bronchi, passes down to bronchioles and finally reaches alveoli. On a microscopic level, once the virus laden droplets reach mucosal membranes, they diffuse through the mucosa (consisting of mucus, periciliary layer, and midlayer) and attach to certain membrane receptors on host cell surfaces, the most prominent one being ACE2, which has been identified as the primary functional receptor for SARS-CoV2, just as for SARS-CoV [21] - [24] .
The current knowledge on SARS-CoV2 infection indicates that the elderly are more susceptible and vulnerable to the infection, while children seem to be the least affected group. Numerous studies report lower rates of SARS-COV2 infection with milder symptoms in children compared to adults [25] - [27] . Some studies attribute these results to their findings that ACE2 expression in children is lower than that of adults.
[28]- [30] . Other possible reasons held responsible for lower rates of morbidity and mortality from SARS-COV2 in children include: the differences in immune responses between children and adults, differences in ACE2 receptor distribution patterns, and lower rates of testing in children due to abundance of asymptomatic cases [31] .
The morphological structure of the virus comes to prominence when discussing viral binding processes. SARS-CoV2 is an enveloped, positive-sense, single-stranded RNA virus and similar to its prior relatives SARS-CoV and MERS-CoV, it belongs to Betacoronavirus genus of the coronavirinae family. SARS-CoV2 viral genome contains four major structural proteins: the S (spike) protein, the M (membrane) protein, the N (nucleocapsid) protein, and the E (envelope) protein [32] . The S protein has a trimeric structure, consisting of an S1 receptor binding subunit and an S2 fusion subunit. During viral infection, S1 and S2 subunits are cleaved by a metalloprotease, TMPRSS-2 (transmembrane protease serine 2), which facilitates viral entry. The S1 subunit functions as the part, which directly binds to the host cell receptor, i.e., ACE2 receptor, creating a Receptor Binding Domain (RBD). The S2 subunit takes role in membrane fusion [33] .
Following viral binding, there are two possible pathways of viral entry for enveloped viruses into host cells: either cytoplasmic fusion in which their envelope fuses with plasma membrane and they release their genome into cytosol, or endosomal membrane fusion (endocytosis) in which they are engulfed by an endosome, and their membrane is fused with the endosomal membrane [34] , [35] . There are multiple mechanisms of endocytic entry suggested by various studies, involving clathrin dependent, caveolae dependent endocytosis [36] , [37] , and clathrin independent, caveolae independent endocytosis [38] , [39] .
In Section II, we presented physio-morphological structure and behavior of the virus, regarding its entry mechanisms into human body and target cells. Here, we will present our system model.
We assumed that the virus carrying aerosol particles enter the human host through the nose, and diffuse through the mucus layer in the nasal cavity, where ACE2 receptors are found most abundantly [40] . The diffusion of the virus takes place in the mucus layer, which renders the shape of the respiratory tract insignificant. Given the fact that the mucus layer is continuous within the respiratory tract [41] , we assume a cylindrical tube with radius r(y) and length l. The change in the radius throughout the tract has limited effect, unless it also modulates the properties of mucus, especially the mucus thickness.
For a large portion of the respiratory tract, mucus layer covers the periciliary layer and a thin surfactant layer separates the two [42] . In our work, we assume that the surfactant layer reduces the surface tension between these two layers to a negligible value, and consequently ignored. Furthermore, we assume that the diffusion coefficient, D, of the virus in periciliary and mucus layers to be the same. In a healthy respiratory system, the mucus inflow to the alveoli is countered by the mucus outflow due to the pericilliary layer. We ignored the mucus outflow mechanism as it may turn the mucus flow into a very complex turbulent fluid model. In other words, we treat it as if it is one single layer.
The existing works studying ACE2 distribution and mucus flow do not comment on differentiations within a region, i.e., ACE2 are homogeneously distributed. Hence, our model assumes cylindrical symmetry.
The virus moves under the influence of the mucus flow from nasal cavity to the alveoli. We partition the respiratory system into seven parts, namely Nasal Cavity, Larynx, Pharynx, Trachea, Bronchi, Bronchiole and Alveoli. Our model is presented in Fig. 2 . Due to the complicated structure of the tracheobronchial tree, we assign transition regions to the closest region. Furthermore, since after each branching, the individual branches become narrower but more numerous, we used the surface area, S i , of each of the seven regions, i ∈ {1, 2, ..., 7}, to calculate its corresponding radii values, r i as
where l i is the length of the i th region. The resultant the respiratory tract is shown in Fig. 3 . Note that Fig. 3 is not to scale, as the corresponding radii for alveol region is two orders larger than the next region, i.e., bronchiole region. Due to the cylindrical symmetry assumption, we can make a longitudinal Upon entering the mucus and periciliary layer, viruses use their viral S-spike proteins to bind to ACE2 receptors on host cell surfaces [43] . We will use the binding rate, λ, to describe the binding process. Due to the spherical shape of the coronavirus, we safely ignore the effect of the orientation of the virus at the time when it makes contact with the ACE2 receptor.
As viruses bind to ACE2 receptors on the host cell's membrane surface, ACE2 receptors are downregulated. Therefore, the number of surface receptors decreases [44] , [45] , making it less likely for other viruses to bind.
We consider two scenarios depending on the ACE2 receptor concentration and the virus population:
• Large virus and large receptor concentration • Only large virus concentration
As stated in Section III we assume a constant mucus flow rate, v from the nasal cavity to pharynx, larynx, trachea, bronchi, bronchiole and finally to the alveoli. Furthermore, the viruses also diffuse with a diffusion coefficient, D, in the mucus layer. The virus concentration is derived using Brownian Motion with drift. We assign y axis for the distance from the entrance of nasal cavity, x axis for the distance from a longitudinal cutting point and z axis as the depth in the mucus layer. Due to the assumption of cylindrical symmetry, the reference point for x coordinate is arbitrary. If a droplet containing N viruses is incident to the mucus level at the location (x 0 , y 0 , z 0 ), the virus concentration at time t is 
The standard deviation for Brownian motion is given as σ = √ 2Dt. Therefore, 95.45% of the population of viruses falls into a sphere of radius 2 √ 2Dt, centred at (x 0 , y 0 +vt, z 0 ), while 99.7% into a sphere with the same centre and radius of 3 √ 2Dt. Hence, for vt 3 √ 2Dt, drift dominates the diffusion and diffusion along the y-axis can be ignored. Similarly, for vt 3 √ 2Dt, drift is dwarfed by diffusion and can be ignored. Fig. 4 shows the dominating trends for Brownian Motion with drift. For v = 50µms −1 [46] and D = 1.5 × 10 −11 m 2 s −1 [47] , the diffusion of the virus in the respiratory tract is shown in Fig. 5 . Clearly, the effects of the diffusion is only visible for large t.
The ACE2-Virus binding can be modelled by obtaining the virus population distribution over the respiratory tract. To achieve this, we start with modelling the kinematics of a single virus incident on the mucus layer. Later, we use our findings as stepping stones to reach impulse response for different scenarios as described in Section III.
We begin our analysis by considering a single virus is moving under the influence of mucus flow. The mucus layer has a thickness of h(y) and a velocity of v, while the respiratory tract radius is r(y), where y lies in the direction of the respiratory tract from nose to lungs. Then, at any segment dy, the concentration of the virus due to a single virus is given by
The time ∆t that the virus spends in a segment of length ∆y is
The probability that it binds to a single ACE2 receptor in the segment with length ∆y becomes
where λ is the molar association constant, N A is the Avogadro's constant and λ 1 = λ/N A is the association constant for a single virus. Note that in the last step, we used first order Taylor series expansion, i.e., e x = 1 + x for small x. Then, p nb , the probability of not binding during ∆t is
If the ACE2 concentration per unit area at y is f (y), then number of ACE2 receptors, n(y), in the patch of length ∆y becomes n(y) = 2πr(y)f (y)∆y, (11) and the probability of the virus evading all ACE2 receptors in the same patch, p e (y), is expressed as p e (y) = p n(y) nb
where from (13) to (14) we use the first order truncation of the binomial expansion, i.e., (1 + x) n = 1 + nx for |nx| 1, which holds due to N A being much larger than any other value in (14) . This assumption is especially effective for ∆y → dy.
From (15), we reach the rate of binding in the patch of length dy as
Then, we find to the number of viruses at y, N (y) using an initial value problem with rate p b (y)
where we used the fact that dy = vdt by definition. An important observation is that N (y) obtained in (21) does not necessarily normalise. DefiningV as
where l is the total length of the respiratory tract,V gives us the rate of the viruses that reach the end of the respiratory tract, i.e., alveoli. Here, since the viruses cannot travel neither forward nor backward, we assume that they will eventually bind to an ACE2 in alveoli.
So far, we only assumed the existence of a single virus to reach (21) . As stated in Section III, there are several scenarios depending on • N (y) = 2C(y)πr(y)h(y)dy, the total number of viruses number on a dy thick strip, • n(y) = 2πr(y)f (y)dy, total number of ACE2 receptors in the same strip, • E b , the expected number of virus bindings in the same strip. E b is loosely calculated by replacing C 1 with C and carrying out (3) to (15) . Hence, (3) becomes
and replacing f (y) with n(y)/2πr(y)dy,
A(y)N (y)n(y),
where A(y) does not depend on ACE2-virus bindings.
Since each ACE2-virus binding destroys both a virus and a receptor, both the virus and the receptor concentrations are affected. Hence, our model must incorporate variations in the concentrations. Each binding causes the number of viruses, N (y) and number of receptorsn(y) to change as
Using (25), (26) and (27) we reach two assumptions: 1) Large N (y): If N (y) is large and E b N (y), the total virus concentration remains constant within the same segment.
2) Large n(y): If n(y) is large and E b n(y), the total ACE2 concentration remains constant within the same segment.
Since A(y) is quite low, i.e., on the order of 10 −10 , for some cases, both of these assumptions hold. Fig. 6 illustrates under which conditions these assumptions hold.
As Fig. 6 shows, when both N (y) and n(y) are large, the assumptions may not hold. Since A(y) depends on r(y) and h(y), the boundaries may change.
Note that changes in n(y) causes a change in the system. As the system parameters change with the input, the system is no longer linear time-invariant. As a result, obtaining the impulse response when (27) does not hold is of no practical use.
1) Large n(y): For large n(y), the virus-ACE2 bindings do not change the receptor number in the same segment. Hence, regardless of how many bindings happen in a given segment, the binding probability of any virus in the same segment is constant. As a result, (17) , multiplied by N (y) gives the expected number of binding on the segment. Furthermore, (21) multiplied by the incident virus count, N 0 , gives us the virus population reaching to any location y of the respiratory tract, i.e.,
Hence, the virus concentration is simply the derivative of (28), i.e.,
.
We proceed to obtain impulse response, I(y, t), by adding the unbounded or free virus distribution,V (y). The free virus population is situated at the location y = vt of the respiratory tract, due to the fact that virus movement on the respiratory tract is solely under the influence of mucus flow. The total number of free viruses is equal to the difference between the initial number of viruses and the total number of bound viruses.V
where δ(.) is the Dirac Delta function. Thus, the impulse response becomes
2) Large N (y) only: In case E b N (y), viruses outnumber the ACE2 receptors. This causes all ACE2 receptors to bind to a virus. Bound virus distribution in the respiratory tract is the same as the ACE2 receptor concentration. Note that, since bound ACE2 receptors downregulate, large N (y) only case cannot be modelled as a linear time invariant system. n(y) = 2πr(y)f (y)dy,
and
where u(.) is the step function and is used to assured that virus distribution is limited to the region 0 − vt, i.e., the range 
We find the impulse response by adding (33) and (34), i.e.,
3) No Limiting Case: In case both assumptions fail, active number of ACE2 receptors constantly changes due to the binding viruses. Therefore, no assumption can be made for this case, and neither (31) nor (34) holds. As a result, there is no closed form expression for this case.
In the host cell, the virus replicate and new virions are released out of the cell via exocytosis. We can model this process with a stationary Markov Chain with six states, namely, binding (B), endocytosis (EN D), release of viral RNA (C), replication (R), degradation (D), exocytosis (EXC) as illustrated in Fig. 7 . The bound virus, can enter the cell via endocytosis, which is mediated by ACE2 receptors. In the state-B, the virus is found bounded to the ACE2 receptor. COVID-19 is an RNA virus, i.e., virus can replicate in the cytoplasm. Thus, in the state-C viral RNA is released to the cytoplasm. In the host cell, the virus can be degraded by lysosomes [48] , which is represented by the state-D.
The transition matrix of the Markov Chain representing the life cycle of the virus in the host cell, Q, is given by
where we set G 1 = −r BE , G 2 = −(r EC + r ED ), G 3 = −(r CR + r CD ), and G 4 = −(r RD + r RE ). The transition rates are provided in Table II . 
represents corresponding state occupancy probabilities. The relation between the states can be expressed as [49] 
which has a solution of (38) is in the form of
Using eigenvalue decomposition, we can express Q as
where v i is an eigenvector of the matrix and λ i the corresponding eigenvalue. As a result, we can express (39) as
Thus, the probability of transition from the state j to the state k in time t is given by
Morphometric measurements of the respiratory tract such as length, diameter, surface area, and mucosal thickness were obtained from estimation studies, journals, databases, and anatomy literature. Given the non-uniform shape and the continuously narrowing nature of the respiratory tract, as in the tracheobronchial tree, we use median values for branching or narrowing structures. Divided into 23 generations of dichotomous sections, the tracheobronchial tree designates a generation for each divided branch starting from trachea, generation 0, and ending at alveoli, (generation 23). The first 16 generations, from generation 0 to generation 16, are defined as the conducting zone, i.e., no gas exchange takes place in this region. From generation 17 to generation 23 is called the transitional and respiratory zone, where gas is exchanged within functional units [50] . The generation 0 directly gives the dimensions for the trachea. Generations 1 to 4 are assumed to be bronchi, 5 to 16, bronchiole and 17 to 23 alveoli respectively. For our parameter data, we mainly use Weibel's "Dimensions of Human Airway Model A" and Gehr's "Annexe A. Anatomy and Morphology of the Respiratory Tract" [51] .
Although there are studies investigating ACE2 receptor gene and protein expressions across different tissues and in specific cell types using single cell RNA sequencing gene expression profiling datasets, mass spectrometry and immune assay techniques, [24] , [52] - [54] , to the best of our knowledge, data on the number of ACE2 receptors on different tissues is not explicitly stated in studies. Most studies provide relative expressions of the receptor in different tissues, shown as proportions, percentages, or plots with no numeric values. Some studies address circulating ACE2 levels, which we cannot directly utilize as we need tissue-specific values. There exist some studies which report ACE2 expression data in animals, which are not compatible with our work either [55] .
The primary challenge of this work is to obtain the ACE2 receptor densities in different tissues of respiratory tract. The lack of studies giving these values is mostly due to the difficulty of measuring ACE2 receptor concentration in a diverse population of all ages. To address this challenge, we exhaustively search among various literature to calculate our estimated values. The specific works that we use are referenced in Table III .
Therefore, we first gather data on the percentage of ACE2 expressing cells for the seven region model described in Section III. Then, we search for the total number of cells in each region. For tissues in which there is no sufficient quantitative data on the percentage of ACE2 expressing cells, the relative proportions of ACE2 expressions of two or more tissues, one of which we have previously calculated are used. Then, we calculate the number of ACE2 expressing cells in each tissue accordingly. Note that these preliminary calculations are the estimates based on the currently available data in the literature.
Due to the lack of data, the effect of age in SARS-CoV2 susceptibility cannot be directly analyzed. However, we investigate the effects of thicker mucus as seen more in elderly and effects of higher ACE2 concentration in nasal cavity as observed in smokers.
In this section, we first present the impulse response simulation and then continue with simulating the effect of mucus flow rate, v, ACE2 receptor density, f (y) and mucus thickness, h(y), on the virus-ACE2 binding.
1) Impulse Response of Unobstructed Viral Progression: In Section. V, we find an analytic expression for the impulse response of unobstructed viral progression through the respiratory tract. Here, we confirm our analytic expression with a Monte Carlo simulation in Fig. 8 . The physiological parameters that we use in the simulations are presented in Table III .
For the Monte Carlo simulation, we divide the respiratory tract to ∆y = 5µm patches. The initial number of viruses are N 0 = 50000. Each virus is independent of each other, i.e., a new number is generated using the Marsenne Twister for each virus in each segment. We see that our analytical solution is in full agreement with the Monte Carlo simulation of the system for large N 0 .
2) Mucus Flow Rate: As it can be seen in Fig. 9 , mucus flow rate has a significant impact on the reach of the virus population. If the patient suffers from another condition causing nasal drip or any other faster mucus flow, the virus spends less time in the upper respiratory system. Therefore, ACE2- virus bindings in the upper respiratory tract is limited, causing the bulk of the virus population to migrate to the lower parts of the respiratory tract, especially bronchioles and alveoli. The virus population in the alveoli is 20-folds more if the mucus drop rate is v = 100µms −1 compared to v = 25µms −1 . This causes the virus to take hold in the alveoli before an immune response can be launched.
3) Nasal ACE2 Receptor Density: Fig. 10 shows us the impact of the ACE2 receptor concentration in the nasal cavity. Assuming distribution of the ACE2 receptors in the other parts of the respiratory tract is the same for different age groups, the difference in the ACE2 levels in nasal cavity has a significant effect on the virus population reaching the lower respiratory tract. The impact of ten-fold increase in ACE2 receptor concentration is six-fold increase in virus concentration in the lower respiratory system. 4) Mucus Thickness: Our model suggests an impact of the mucus thickness. Since we assume that the virus can move freely in the mucus layer via diffusion, thicker mucus implies that there is less chance for the ACE2-virus binding. Fig. 11 shows the effect of the mucus thickness. The virus population in the alveoli is 4.45 times more in the four times thicker mucus compared to the base mucus level.
In this study, we analyze SARS-CoV-2-ACE2 receptor binding event to determine bound virus population at the different regions of respiratory system. To this end, we develop a molecular communication model of SARS-CoV-2 transmission through the human respiratory tract, which reduces the inherent complexity of the respiratory system. We perform an analysis of the developed model using mucus flow rate and ACE2 receptor densities, which are calculated based on the realistic data of respiratory tract surface parameters. Based on the analysis, we reach that higher mucus flow rate results in virus migration to the lower respiratory tract, which is compatible with the experimental results found in the literature. Our model will be useful to describe the travel of the virus through the respiratory tract and to simulate the effect of interventions (e.g. antivirals) to decrease the viral load.