key: cord-187700-716af719
authors: Lee, Duan-Shin; Zhu, Miao
title: Epidemic Spreading in a Social Network with Facial Masks wearing Individuals
date: 2020-10-31
journal: nan
DOI: nan
sha: 
doc_id: 187700
cord_uid: 716af719

In this paper, we present a susceptible-infected-recovered (SIR) model with individuals wearing facial masks and individuals who do not. The disease transmission rates, the recovering rates and the fraction of individuals who wear masks are all time dependent in the model. We develop a progressive estimation of the disease transmission rates and the recovering rates based on the COVID-19 data published by John Hopkins University. We determine the fraction of individual who wear masks by a maximum likelihood estimation, which maximizes the transition probability of a stochastic susceptible-infected-recovered model. The transition probability is numerically difficult to compute if the number of infected individuals is large. We develop an approximation for the transition probability based on central limit theorem and mean field approximation. We show through numerical study that our approximation works well. We develop a bond percolation analysis to predict the eventual fraction of population who are infected, assuming that parameters of the SIR model do not change anymore. We predict the outcome of COVID-19 pandemic using our theory.

In December of 2019 a few patients of a new infectious respiratory disease were detected in Wuhan, China. This disease has been called coronavirus disease 2019 (COVID-19) and the virus that causes COVID-19 has been named SARS-CoV-2 by the World Health Organization (WHO). WHO declared the outbreak a public health emergency of international concern at the end of January 2020, and a pandemic on March 11, 2020 . Since the outbreak, most countries have adopted various measures in an attempt to contain the pandemic. These measures include restriction of traveling, shutting down schools, restaurants and businesses, canceling large gatherings such as concerts, sports and religious activities, and even city lockdowns where residents are not allowed to leave home unless emergencies. Clearly these measures seriously affect daily lives and are devastating to the economics. The purpose of this report is to show that wearing facial masks is a simple and inexpensive measure to contain the spread of COVID-19. In fact, we shall show that if a relatively small fraction of population wears facial masks, the disease can be contained.

Facial masks have been shown in labs to be effective to limit the spread of droplets or aerosols if a wearer coughs [1] - [4] . This ability is measured by a quantity called the outward mask filter efficiency of masks. Facial masks also protect their wearers from inhaling droplets or aerosols from a nearby cougher, if the cougher does not wear a mask. This ability is measured by a quantity called the inward mask filter efficiency of masks. Thus, facial masks can be particularly useful to confine the spreading of diseases that transmit through droplets or aerosols. However, the effect of wearing facial masks to the epidemic spreading has never been studied in a network level. In this report we present an epidemic network study to justify this argument. Specifically, we propose a time dependent susceptible-infected-recovered (SIR) model with two types of individuals. Type 1 individuals wear a facial mask and type 2 individuals do not. A randomly selected individual from a population is a type 1 individual with probability p, and is of type 2 with probability 1 − p. There are four types of contacts between two individuals depending on whether the two individuals wear a facial mask or not. These four types of contacts have four different disease transmission rates. From the data published by John Hopkins University [5] we progressively estimate the time dependent disease transmission rates and the recovery rates of the SIR model.

For parameter p, we propose a stochastic version of the SIR model. We derive the transition probability of the number of infected individuals from one time slot to the next. We propose a maximum likelihood estimation of p that maximizes the transition probability. The transition probability is expressed in terms of binomial distributions. Parameters of the binomial distributions corresponding to real data published in [5] are typically very large. That makes the transition probability numerically difficult to compute. We propose an approximation of the transition probability based on central limit theorem and mean field approximation. Through numerical studies, we show that the approximation works well. We derive a percolation analysis of the maximum number of individuals that can eventually be infected. We incorporate the maximum likelihood estimation and the percolation analysis into the progressive estimation. That is, based on the data published by John Hopkins University, we progressively estimate the disease transmission rates and the recovery rates. We then find p from the maximum likelihood estimation. Finally, using percolation analysis, we predict the maximum number of individuals that can eventually be infected.

The outline of this report is as follows. In Section II, we present a time dependent SIR model. We present a progressive estimation of the disease transmission rates and the recovery rates of the model. In Section III, we present a maximum likelihood method to estimate p. In Section IV we present a percolation analysis. In Section V we present results of numerical study and simulation. We present the conclusions in Section VI.

In this section we present a discrete-time susceptible-infected-recovered (SIR) model. Time is divided into periods of equal lengths. There are two types of individuals. Type 1 individuals wear a facial mask and type 2 individuals do not. A randomly selected individual from a population is a type 1 individual with probability p(t) in period t. A randomly selected individual is of type 2 with probability 1 − p(t). Let s(t) and r(t) be the number of susceptible and recovered individuals, respectively, at time t. Similarly, let x i (t) be the number of infected type i individuals in period t for i = 1, 2. In this model, the disease transmission rates are time dependent. Let β ij (t) be the expected number of type j susceptible individuals who receive the disease from one type i infected individual per unit time in period t. Let γ(t) be the recovering rate of the disease in period t. We assume that both types of infected individuals have the same recovering rate.

In this report we assume that the length of a time unit in the discrete-time SIR model is τ days. The dynamics of this discrete-time SIR model is as follows. The existing infected individuals at time t− transmit the disease to newly infected individuals. Those infected individuals existed at time t− become recovered at time (t + 1)−. It is easy to derive the following set of difference equations for the SIR model.

where

is the size of the population. We assume that the epidemics is in the early stage. That is, we assume that s(t) ≈ n. Under this assumption, the difference equations for x 1 (t) and x 2 (t) reduce to

We now determine parameters β ij (t) and p(t). We reduce the number of parameters. Previous study [2] suggested that

where η 1 and η 2 are outward and inward efficiencies of masks, respectively. Study [4] suggested that η 1 > η 2 . In addition, we assume that

Thus, we only need to determine β 22 (t). The other three parameters are determined according to Eqs. (4), (5) and (6) . Substituting (4), (5) and (6) into (2) and (3), we obtain

Now we determine the values of parameters γ(t), β 22 (t) and p(t) based on the data published by John Hopkins University [5] . Note that John Hopkins University publishes total number of daily newly infected individuals and the number of recovered individuals. From the published data, one can easily compute r(t) for each t. The website does not distinguish between infected individuals who wear masks or who do not. Thus, we must determine β 22 (t) and p(t) based on the total number of infected individuals

at time t. From the data published by John Hopkins University, we estimate γ(t), β 22 (t) and p(t). We develop in section IV a percolation analysis to determine the ultimate size of infected population, if these parameters do not change.

Recovering rate is simple to determine. From (1) we have

Next, we determine β 22 (t). To determine the value of β 22 (t) in period t, we let x(t) be the total number of infected individuals in period t, i.e. Adding Eqs. (7) and (8) and solving for β 22 (t), we obtain

To determine the value of β 22 (t) at time t, we use x(t + 1) published in [5] . We assume that x 1 (t) and x 2 (t) and p(t) are available for t. We use (10) to determine β 22 (t). The value of p(t) is determined by a maximum likelihood estimation method. We shall present this method in Section III. After β 22 (t) at time t is determined, we use (7) and (8) to determine x 1 (t + 1) and x 2 (t + 1). Note that the sum of x 1 (t + 1) and x 2 (t + 1) determined in this way agrees with x(t + 1).

We summarize the algorithm that determines parameters of the time-dependent SIR model in Algorithm 1. Compute γ(t) using Eq. (9); 7:

Find β 22 (t) using Eq. (10) and p(t), x 1 (t), x 2 (t) and y(t + 1);

Compute x 1 (t + 1) and x 2 (t + 1) by Eqs. (7) and (8); 9: Find predicted giant component size S(t) using Eq. (V); 10: end for

In this section we present a maximum likelihood estimation method to determine the value p(t) based on x(t), x(t + 1) and β 22 (t − 1). This estimation will be progressively used by Algorithm 1 in periods t = 2, 3, . . .. Thus, in the rest of this report we simplify the notation by eliminating the dependency of time from p(t) and γ(t) and simply use p and γ, respectively.

Recall that x(t) and x(t + 1) are the total number of infected individuals in periods t and t + 1. We would like to determine the value of p such that the likelihood of this sample path is maximized. To determine the likelihood function, we propose a probabilistic version of the SIR model. This probabilistic SIR model is based on independent cascade models. Independent cascade models are popular not only in the study of epidemic spreading but also in the influence maximization problems of viral marketing [6] , [7] . In an independent cascade model, each infected node has exactly one opportunity to transmit the disease to its neighbors. Whether the transmissions are successful or not depend on independent events. In our model, we distinguish between nodes that have used their opportunity to transmit the disease from those who have not. Our model is a discrete time model. Let X t be the total number of currently infected individuals in period t. Let Y t be the number of individual who contract the disease in period t. In our model, we assume that those who contract the disease in period t have ability to transmit the disease to their neighbors, and lose the ability in periods t + 1, t + 2, . . . and so on. Thus, in period t there are X t − Y t infected individuals who can not transmit the disease to others. Those who are infected but cannot transmit the disease to others in period t can remain infected in subsequent periods or become recovered. A graphical illustration of the model is shown in Figure 1 . It is clear that this model is a discrete time Markov chain. We shall derive the transition probability

and findp(t) such thatp

We now specify more details of the probabilistic SIR model. Each individual can be of type 1, or of type 2. Each infectious individual has k contacts, to whom he or she can transmit the disease. A type i infectious individual can transmit the disease to a type j susceptible individual with probability φ ij , where i, j = 1, 2. Parameter φ ij is related to β 22 (t) and γ through

Let C be the number of type 1 individuals among the Y t infectious individuals in period t.

where {U j , j = 1, 2, . . .} and {V j , j = 1, 2, . . .} are two independent and identically distributed (i.i.d.) sequences of Bernoulli random variables with success probabilities p 1 and p 2 , respectively. The two sequences are independent to anything else. Event {U j = 1} indicates a type 1 infectious individual successfully transmits the disease to a neighboring node. This occurs with probability p 1 , where

Similarly, event {V j = 1} indicates a type 2 infectious individual successfully transmits the disease to a neighboring node with probability p 2 , where

An infected individual becomes recovered with probability γ, and remains infected otherwise. Hence,

where {W j , j = 1, 2, . . .} is an i.i.d. sequence of Bernoulli random variable independent of anything else. The success probability of W j is 1 − γ.

We now derive an expression for the transition probability in (11) as a function of p. Let We can rewrite (14) and (17) in terms of I 1 , I 2 and I 3 , i.e.

Conditioning of event {C = i}, we have

Conditioning on {C = i}, random variables I 1 , I 2 and I 3 all have binomial distributions. Define binomial probability mass function

The conditional distributions of I 1 , I 2 and I 3 are

Substituting (22), (23) and (24) into (21) and noticing that I 1 , I 2 and I 3 are independent, we obtain

In addition,

Taking average to (25) with respect to event {C = i} using (26), we have

Eq. (25) is very complicated to evaluate and it is difficult to be directly used in the optimization problem (12) when y(t), y(t + 1) or both are large. We now propose an approximation method to simplify (27). First, we take logarithm on (27) to obtain

From the De Moivre-Laplace central limit theorem we approximate binomial distributions by normal distributions. Let N (z, µ, σ 2 ) denote the probability density function (pdf) of a normal random variable with mean µ and variance σ 2 , i.e.,

2πσ .

First, we approximate probability mass function b(i, y(t), p) in (25) by a normal distribution. Specifically,

where

Next, recall that conditioning on event {C ∈ (z, z+dz)}, I 1 and I 2 are independent binomial random variables. We approximate them by normal pdf's with means

and variances

respectively. Since I 1 and I 2 are conditionally independent, we have P(I 1 + I 2 ∈ (y(t + 1) − 1/2, y(t + 1) + 1/2) | C ∈ (z, z + dz)) ≈ N (y(t + 1), zkp 1 + (y(t) − z)kp 2 ,

We propose to further simplify Eq. (32). We apply a mean-field approximation [8] , [9] to replace z, which is the sample value of C, in (30) and (31) with the expectation E[C] = y(t)p. Specifically, let

where

is the basic reproduction number of the branching process. We further approximateσ 1 defined in (33) bỹ

= y(t)R 0 assuming that p 1 and p 2 are small. With these approximations, the logarithmic transition probability reduces to

Since the second term on the right side of the preceding equation is independent of p, maximizing the logarithmic transition probability is equivalent to maximizing the first term on the right side. Define

We approximate the optimization problem in (12) by the following optimization problem

To study the extrema of f we differentiate f , i.e.

Note that f (R 0 ) has a unique positive root, which is

The last expression gives the most probable basic reproduction number that produces the sample path (y(t), y(t + 1)). At this root, the second derivative is negative, i.e.

achieves maximum when Eq. (39) holds. We solve p from (34) and (39). Recall that we assume (4), (5) and (6) . Under these assumptions, (39) can have a unique root in [0, 1], in which case, the unique root is the solution of the optimization problem (37). Eq. (39) can have no roots in [0, 1], in which case the solution of (37) isp(t) = 0 orp(t) = 1. We summarize the solution of the optimization problem (37) in the following proposition. The proof of the proposition is presented in the appendix at the end of this report. Proposition 1. Eq. (39) has two real roots. Either Eq. (39) has exactly one root in interval [0, 1], or it has no root in this interval. If Eq. (39) has exactly one root in [0, 1], it is the smaller root (denoted by p ) of the two roots. In this case, the optimal solution of (37) is p . If Eq. (39) has no roots in [0, 1], either p < 0 or p > 1. In the former case, the optimal solution isp(t) = 0. In the latter case, the optimal solution isp(t) = 1.

In this section we consider a random contact network, in which a disease transmitted by droplets or aerosols spreads according to an independent cascade model. That is, the disease transmits from a node at one end of an edge to the node at the other end with a probability. In addition, the transmissions along all edges are independent. We shall present a percolation analysis of this model and obtain percolation thresholds and sizes of giant components.

We now describe our model. Consider a random graph (G, V, E). Randomly select a node from the graph. Let Z be the degree of this node. Let g 0 (z) denote the probability generating function of Z, i.e.,

Now randomly select an edge. Let Y be the excess degree of a node reached along the randomly selected edge. Let g 1 (z) be the probability generating function of Y , i.e.

Every node in this graph can be one of two types. A type 1 node denotes an individual who wears a facial mask and a type 2 node denotes an individual who does not. A randomly selected node is of type 1 with probability p and is of type 2 with probability 1 − p. Assume that this event is independent of anything else. As mentioned before, an infectious disease spreads in this network according to an independent cascade model. Consider a randomly selected edge connecting two nodes, say node V 1 and node V 2 . Let φ ij be the conditional probability that the disease transmits from node V 1 to node V 2 , given that the types of nodes V 1 and V 2 are i and j, respectively, where i, j = 1, 2.

We now present a percolation analysis of the random network model described above. Percolation analysis has been a useful tool to study resilience of communication networks [10] , [11] and epidemic networks [12] - [15] . There are two types of percolation models, bond percolation and site percolation, depending on whether edges or nodes, respectively, are removed randomly [8] . Our model is a form of bond percolation, where edges are randomly removed. Removed edges imply that diseases can not be transmitted from one end of these edges to the other end. The size of the largest component in the percolated network is the maximum fraction of possibly infected population in the network.

Randomly select an edge. Let V 1 and V 2 be the two nodes at the two ends of the edge. Suppose that the type of node V 1 is i. Let E i be the event that along the selected edge from V 1 , one can not reach a giant component. Let u i = P(E i ). Now we condition on the event that the type of node V 1 is j. Event E i will occur, if the randomly selected edge is removed. This occurs with probability 1 − φ ij . With probability φ ij the randomly selected edge is present. Let Y be the number of neighbors of node V 2 , not including V 1 . Event E i will occur, if one can not reach a giant component along any one these Y edges. Combining these arguments, we have

where C is the number of type 1 nodes connected to node V 2 not including V 1 . The distribution of C conditioning on Y = k is binomial. Thus, Eq. (41) becomes

for i = 1, 2. Eq. (42) is a system of nonlinear equations, from which we can solve for u 1 and u 2 .

Once we have u 1 and u 2 we can compute the giant component size of the percolated network. Randomly select a node from the network. Let S i , where i = 1, 2, be the conditional probability that the randomly selected node is connected with a giant component, given that the selected node is of type i. Let X be the degree of the randomly selected node. The node is connected with a giant component if along at least one of its edges one can reach a giant component. Conditioning on X = k, let I be the number of type 1 nodes among the k neighbors. Combining all these arguments, we have

for i = 1, 2. By taking average on the conditional probabilities, a randomly selected node is connected to a giant component with probability

This is also the expected size of the giant component.

Let u be a 2 × 1 vector over the set of real numbers R. Specifically, let

where symbol T denotes transposition of vectors. Let f be a vector-valued function mapping from R 2 to R 2 , i.e.,

where x = (x 1 , x 2 ). Eq. (42) implies that u is a root of

The roots of equations of the form (45) are also called the fixed points of function f . It is clear that function f always has fixed point 1 = (1, 1) T . Lee et al. [16] established that f has an additional fixed point if the dominant eigenvalue of the Jacobian matrix evaluated at 1 is greater than one. In addition, this fixed point is attractive. The Jacobian matrix for function f evaluated at x = a is defined as

The Jacobian matrix for the function f defined in (44) evaluated at (1, 1) T is

where E[Y ] = g 1 (1) is the expected excess degree of a node reached by a randomly selected edge. It is easy to derive the eigenvalues of J . Denote the two eigenvalues of J by λ 1 and λ 2 , Then,

Note that both eigenvalues are real, and λ 1 > λ 2 . Thus, spectral radius or the dominant eigenvalue of J is λ 1 . It is also this eigenvalue that controls the percolation threshold of the epidemic network. 

In this section we present numerical and simulation results. We first verify the accuracy of (35) to approximate the transition probability in (27). We select an arbitrary set of parameters:

According to [2] , the efficiency of a typical mask is in the range from 0.58% to 85%. Since masks are typically more efficient to stop viral transmission if sources wear masks [4] , we set η 1 = 0.8 and η 2 = 0.65. We select a value for p between 0 and 1. We randomly generate ten thousand values for y(t). We perform simulation to generate ten thousand values for y(t + 1) based on y(t) and p. We solve (12) based on (27) and (37) for the ten thousand values of y(t) and y(t + 1). We calculate the average distance between the two solutions. The result is shown in Figure 2 . This shows that the approximation method works quite well.

Next, we use the data published by John Hopkins University [5] to predict the spread of the COVID-19 pandemic using our model. Due to large infection rates of some countries, we assume that k = E[Y ] = 500 in order to have φ 22 ≤ 1. We choose τ = 5 days as the width of a time slot in the time dependent SIR model. We first execute Algorithm 1 on the number of infected and recovered individuals in the mainland of China. The time functions of β 22 (t) and p(t) are shown in Figure  3 . The maximum fraction of population who will be infected predicted at time t is also shown in the figure. The epidemic started from the People's Republic of China from January, 2020. The authority managed to contain the epidemic very well by the end of February. However, there are several significant events between January and August of 2020. In mid April, the Chinese authority revised the way that death toll is calculated [17] . In mid June, Beijing faced a second wave of infections [18] In mid July, there was a surge of infected number of individuals in Xinjiang [19] . At the end of July, there was a surge in Dalian [20] . From Figure 3 , we see that β 22 (t) is large initially in January, and is controlled to reach a small value by the end of February. The value of β 22 (t) again rises and drops before and after the events described above. From Figure 3 , note that p(t) also rises and drops before and after the events. However, function p(t) lags behind function β 22 (t). This correlation has a nice interpretation. Rise in β 22 (t) usually manifests itself in the rise of infected number of individuals. As the population sees a rise in the infected number of individuals, more people wear masks to protect themselves. On the other hand, decreases of β 22 (t) result in less number of newly infected individuals. People typically see this as sign of a safe community. Thus, less people wear masks in public places. As β 22 (t) rises, the predicted size of giant component S(t) also increases. However, as p(t) catches up and rises, S(t) decreases. Notice that the time sequence β 22 (t) reflects a joint efforts of many measures to contain the epidemic. Wearing facial masks is a measure in the personal level. Shutting down schools, businesses, and keeping people at home is a measure in the government level. In this paper, we consider only the measure of wearing masks and ignore other measures. Thus, p(t) in Figure ? ? may be higher than the actual fraction of population who wear masks, as it reflects a joint effect of many measures to contain the epidemic.

Next we study the epidemics of the United States of America, India and France. We execute Algorithm 1 to compute β 22 (t), p(t) and S(t). The time functions of β 22 (t), p(t) and S(t) for the three countries are shown in Figure 4 , Figure 5 and Figure  6 , respectively. We present p(t) and S(t) of the three countries in Figure 5 and Figure 6 . Since these functions fluctuate a lot, we use a built-in Matlab function to compute the polynomial regression of these functions in order to show a general trend. Note that India has a smaller infection rate β 22 (t) than that of the U.S.A. However, India has a much smaller predicted size of giant component than that of the U.S.A. This is because India has a very large recovery rate [21] , [22] . We also note that around July France has a β 22 (t) comparable with that of the US. However, France has a high predicted S(t) after August. It is worth noting that p(t) of France in mid August is not close to 1. By raising p(t), France may have a better control of the epidemics.

In this report, we presented a time dependent SIR model, in which some individuals wear facial masks and some do not. Based on the number of infected individuals and the number of recovered individuals published by J. H. University, we estimate the disease infection rates and recovery rates. We proposed a probabilistic version of the SIR model. We derived the transition probability of this random SIR model. By maximizing the transition probability, we estimate the most probable value of the fraction of population who wear masks. This transition probability numerically difficult to compute, if the states of model are large. Based on central limit theorem, we proposed an approximation. Through numerical and simulation study, we show that the approximation works well. Finally we carried out a percolation analysis to predict the eventual fraction of population who will be infected with the disease. We proposed a progress analysis of the epidemics. Using results from the progressive analysis, we analyzed the epidemics of four countries.

Effectiveness of facemasks to reduce exposure hazards for airborne infections among general populations

The effect of mask use on the spread of influenza during a pandemic

Testing the efficacy of homemade masks: Would they protect in an influenza pandemic

Respiratory source control using a surgical mask: An in vitro study

Maximizing the spread of influence through a social network

Networks, crowds and markets reasoning about a highly connected world

Networks: An Introduction

Statistical approach to quantum field theory : An introduction

Resilence of the internet to random breakdowns

Percolation in interdependent and interconnected networks: Abrupt change from second-to first-order transitions

Percolation on heterogeneous networks as a model for epidemics

Epidemics and percolation in small-world networks

Epidemic models and percolation

Percolation and epidemic thresholds in clustered networks

A generalized configuration model with degree correlations and its percolation analysis

The death toll in Wuhan was revised from 2579 to

COVID-19 pandemic in Beijing

Officials rule out domestic transmission as origin of Dalian COVID-19 cluster

Coronavirus recovery: What do India's high COVID recovery numbers really mean? We explain

India's Covid-19 recovery rate surges past 78%, one of the highest globally: Govt

In this appendix, we prove Proposition 1. Proof of Proposition 1. We rewrite R 0 defined in (34) to show dependency with p explicitly, i.e.R 0 (p) = kφ 22 (η 1 η 2 p 2 − (η 1 + η 2 )p + 1).In (49) we have assumed (4), (5) and (6) . Note that η 1 < 1 and η 2 < 1. Since the right side of (49) is a quadratic polynomial in p, it is easy to establish the following results.4) The minimum of R 0 (p) occurs atThis point is greater than one, since η 1 < 1 and η 2 < 1.

which is negative. From the results above and the fact that the right side of (39) is positive, it follows that (39) must have two real roots. It also follows from the results above that R 0 (p) is monotonically decreasing and (39) can not have two roots in [0, 1]. Since R 0 (p) is decreasing for p ∈ [0, 1], it follows that if (39) has exactly one root in [0, 1], the root is p . It is also quite clear that if (39) has a unique root in [0, 1], the root is the optimal solution of (37). Now we analyse cases, in which (39) does not have roots in [0, 1]. Since R 0 (p) is decreasing for p ∈ [0, 1], it follows that either both roots of (39) are greater than 1, or p < 0 and the other root is greater than 1. We now analyze these two cases separately. We first note that we have treated R 0 as an independent variable in (37). Thus, the derivative f (R 0 ) is with respect to R 0 . In the two cases in which Eq. (39) has no roots in [0, 1], we have to treat the objective function in (37) as a function of p. Specifically, the derivative of the objective function with respect to p is d dpwhere f (R 0 ) is given in (38). We claim that in the first case the objective function of (37 is increasing for all p ∈ [0, 1]. Thus, the optimal solution of (37) isp(t) = 1. We also claim that in the second case, the objective function of (37 is decreasing for all p ∈ [0, 1]. Thus, the optimal solution of (37) isp(t) = 0.We now prove the two claims. From (38) and (50),whereNote that h is a quadratic and concave function of R 0 with h(0) = y(t + 1) 2 /y(t) > 0. Function h has a zero at R 0 = c, where c is the right side of (39), i.e. c = −n + n 2 + 4y(t + 1) 2 2y(t) .If c < R 0 (1), both roots of (39) are greater than 1. This corresponds to the first case. If c > R 0 (0), then p < 0 and the other root is greater than 1. This corresponds to the second case. In the first case,for all p ∈ [0, 1]. From (51), it follows that df (R 0 (p))/dp > 0, since R 0 (p) < 0 for all p ∈ [0, 1]. This proves the first claim.In the second case, R 0 (p) < c for all p ∈ [0, 1]. Thus, h(R 0 (p)) > g(c) = 0 for all p ∈ [0, 1]. From (51), it follows that df (R 0 (p))/dp < 0. This proves the second claim.