Basic reproduction number of SEIRS model on regular lattice

: In this paper we give a basic reproduction number R 0 of an SEIRS model on regular lattice using next-generation matrix approach. Our result is straightforward but di ﬀ ers from the basic reproduction numbers for various fundamental epidemic models on the regular lattice which have been shown so far. Sometimes it is caused by the di ﬀ erence of derivation methods for R 0 although their threshold of infectious rates for epidemic outbreak remain the same. Then we compare our R 0 to the ones by these previous studies from several epidemic points of view.


Introduction
The basic reproduction number, usually expressed by R 0 , is one of the most important indices of epidemic models. This quantity gives the criterion for determining whether the infectious disease can spread into the whole population or not: when R 0 > 1 an epidemic will result from the introduction of infectious agent, whereas when R 0 < 1 the number of infected is expected to decline (e.g., [1]). Therefore it is known from the threshold principle, that R 0 = 1 gives its threshold.
This equation for theshold also coincides with the condition for the threshold between stability and instability of disease-free equilibrium, but we sometimes obtain different basic reproduction numbers even if the stability conditions for disease-free equilibrium are the same. When this happens, we can determine whether or not the infectious disease can expand, but it produces a possibility for us to make a wrong evaluation of its propagation speed.
In simplest cases, the basic reproduction number can be calculated as the ratio of newly reproduction rate of infectious individuals during the infected period of the focal infected individual. Unfortunately, however, for more complicated models, it requires some other contrivances to obtain the basic reproduction number in general.
In such situations, the most outstanding invention for obtaining R 0 is to use next-generation matirix approach [2], in which R 0 is defined as the spectral radius or dominant eigenvalue of the next-generation matrix. The computation of R 0 by the next-generation matrix is formulated by [3] and [4] (these methods by next-generation matrix are reviewed by [5]). Since then, a tremendous number of studies have adopted this method to obtain R 0 , and now it has become a standard approach to finding R 0 . On the other hand, we should notice that we can only find a few examples of obtaining R 0 in spatial epidemic models by using the next-generation model. In [6], R 0 was first obtained for their infectious model on the lattice by using next-generation matrix.
The rest of this paper is organized as follows. Section 2 reviews the basic reproduction number of SIR model, which gives us the most understandable derivation and result. In Section 3, we verify the procedure for calculating the basic reproduction number by using the next-generation matrix shown in the seminal work [3], and we apply it to SEIRS model on a regular lattice. In Section 4, we compare our result with the previous studies so far.

SIR model and its basic reproduction number R 0
First, we concentrate on the SIR model proposed by [7] as one of the most fundamental epidemic models. This SIR model can be written as follows: Here, ρ S , ρ I , and ρ R indicate the fractions of susceptible, infected and recovered individual, respectively. The infection occurs by the right-hand side (RHS) in Eq (1) or the first term of the RHS in Eq (2), which are assumed by the product of both fractions ρ S and ρ I , i.e., the law of mass action. The recovery from infection randomly occurs by the second term of the RHS in Eq (2) or the RHS in Eq (3). The rates β > 0 and σ > 0 indicate the unit time of speed of infection and recovery, respectively.
To check the behavior after introduction of a small amount of infectious individuals into susceptible population, or a small disturbance around the disease-free equilibrium (1, 0, 0), which is an equilibrium of the system (1)-(3), one has to determine whether the disease can invade the population or not. Thus, the linearization of Eq (2) around (1, 0, 0) becomes dρ I dt = βρ I − σρ I = (β − σ)ρ I and defining the ratio of infectious rate during the period of infected: then R 0 > 1 if and only if β > σ.
Compared with the above result, limiting transmission of the disease only within the nearest neighboring individuals reduces the value of R 0 . Therefore, we review the derivation of R 0 for SIR model by [8] (the notation in [8] slightly differs from the one presented below, but we can have an exact correspondence). First of all, Eq (2) corresponds to the following equation: where ρ S I is the fraction of the nearest-neighboring pair S and I, and q S /I is the conditional probability that the focal site has state S under the condition of state I in the nearest-neighbor site: ρ S I = ρ I q S /I . At the initial phase of an infection invading, we can write the rate of I in Eq (5) as Therefore, the basic reproduction number is defined as and then if R 0 > 1, then the disease can spread into population on a lattice. So we need to give the initial value of q S /I (0) at the neighborhood around the steady state without infectious individuals. After the invasion of an infectious individual, the spatial structure with the dynamical steady state would be soon reached, and so it can be used as an initial condition for the variable q S /I to determine whether or not an infection will succeed to invade. Then we consider the dynamics of q S /I : and so we obtain the steady state as q S /I = 1 − 2/z. Substituting this as an initial state into Eq (7) results as which indicates that the basic reproduction number for SIR model on lattice or network by Eq (10) is less than the one for the corresponding non-spatial model by Eq (4). Based on the above idea the basic reproduction numbers for basic epidemic models on regular lattice have been obtained by several studies [8,10,11]. In the next section we derive the basic reproduction number for SEIRS model, which is considered as more generalized model than SIR model, by using next-generation matrix introduced by [2] and studied extensively by [3].
3. Definition of R 0 for an epidemic model on the lattice by using next-generation matrix 3.1. Next-generation matrix and R 0 We can calculate R 0 through the next-generation matrix by [3] as follows. Consider the vector . Here x 1 , · · · , x m corresponds to m population densities of infected classes. Then disease free state is defined as The epidemic model is rewritten as i are the terms corresponding to the appearance of new infections in compartment i, transfer of individuals into compartment i by all other means, and transfer of individuals out of compartment i, respectively. Also, these functions satisfy the following conditions (A1)-(A5) with a disease free equilibrium x 0 : ..,n j=1,...,n have negative real parts.
Then R 0 is defined as the spectrum radius or dominant eigenvalue of next-generation matrix FV −1 : In the next subsections we use the above procedure to get R 0 's for SEIRS epidemic model in both cases of completely mixing and only the nearest neighboring infection on regular lattice.

R 0 for SEIRS epidemic model for completely mixing population
In addition to the system (1)-(3) we use the notation ρ E , ω, ν for the fraction of exposed individuals (who are in incubation period during which the individual has been infected but is not yet infectious), the rate of natural loss of immunity, transition rate to infectious, respectively. Then, we consider the following SEIRS model: Observe that dρ S /dt + dρ E /dt + dρ I /dt + dρ R /dt = 0 and each variable indicates the fraction of each state, we can assume the total sum of ρ S + ρ E + ρ I + ρ R = 1. Then we can choose only three variables, ρ S , ρ E , ρ I as independent, and so ρ R = 1 − (ρ S + ρ E + ρ I ). Therefore, we get the following system from Eqs (11)- (14): To express this system in non-dimensional terms, we choose new variables as follows: Then the above system (15)-(17) is transformed as We consider {ρ E , ρ I } as population densities of infected class. Hence, when we check conditions (A1)-(A5) for this system (18)-(20), we find three cases of vectors F and V as follows: Also, observe that the disease-free equilibrium is (ρ S , ρ E , ρ I ) = (1, 0, 0); the matrices F and V are given as,

respectively. Then
respectively. Therefore, basic reproduction numbers for these matrices become respectively. As pointed out by [3], we should choose F and V by appropriate epidemiological interpretation rather than mathematically. The first and the second choices give the same basic reproduction numbers, but both the second and the third include the terms of progression of infectiousness in the vector F . Therefore, we should adopt the first choice which considers only newly infected individuals.

R 0 for SEIRS epidemic model for the nearest-neighboring infection
Using the same notatons in the previous sections and others, such as ρ S S I , which indicates the fraction of the triplet S S I, the dynamics of SEIRS model on regular lattice is written as follows: where ρ σσ = ρ σ σ holds due to symmetry. Observe that by using (22)-(30), indicates that the system (31)-(33) including spatial structure corresponds to the non-spatial model (18)-(20) in the previous subsection. And so, the system of nine equations (22)-(30) with the following relation gives the dynamics of the nearest-neighboring pairs for spatial model: Unfortunately, however, the system (22)-(30) is not closed since triplet densities such as ρ S S I may change as time passes, whose dynamics also includes higher-order correlation functions (quartet densities) and these procedures continue forever. Therefore, we adopt an approximation called pair approximation [9] to this system. By pair approximation, conditional probabilities are defined: and the effect of the site with state σ from the site with state σ may be considered smaller compared with the site with state σ because the former is the next-nearest neighbor, but the latter is the nearest neighbor. Then we can approximate the conditional probability (34) as follows: Using Eqs (34) and (35), we obtain ρ σσ σ ≈ ρ σσ ρ σ σ ρ σ .
Pair approximation (36) with ρ S = ρ S S + ρ S E + ρ S I + ρ S R transforms Eqs (22)-(30) into the closed system as follows: We consider {ρ S E , ρ S I , ρ EE , ρ EI , ρ ER , ρ II , ρ IR } as population densities of infected class because these variables include E or I. When we check conditions (A1)-(A5) for the system (37)-(45) with the choice of F from the viewpoint of the reproduction of newly infected individuals, we can determine only one possible set of vectors F and V : Observe that the first and the fourth elements correspond to the changes of variables ρ S S and ρ S R , respectively, and that the disease-free equilibrium is Then, the next-generation matrix is The characteristic equation of this next-generation matrix (46) becomes Therefore the spectral radius of the next-generation matrix, or the basic reproduction number, for SEIRS model is When we set z → ∞ then we obtain which corresponds to R 0 for the standard SEIRS, SEIR, SIRS, SIR, and SIS models as Eq (4) and (21). Compared with Eq (48) and (47), epidemic breakout for the former only depends on the infectious ratẽ β, on the other hand, the latter depends not only on the infectious rate but also on other parameters, progression to infectiousnessν and rate of the natural loss of immunityω.
In Figure 1, we show the dependence of these three parameters on the threshold value of R 0 = 1. We can see that critical infectious rateβ increases as transition rate to infectiousν but decreases as rate of the natural loss of immunityω. This is reasonable because the period during infected 1/ν is shorter asν increases, then the strength of infection should be larger to maintain the value of R 0 . On the other hand, ifω increases, then the number of susceptibles, to some of which the focal infected individual has not yet transmitted disease, increases, and soβ can be reduced for the critical R 0 = 1.

R 0 for SEIR, SIRS, SIR, SIS epidemic models: simpler cases
We consider simpler epidemic models with reducing variables and parameters. When R never returns to S andω becomes zero, then we get SEIR model with pair approximation:   By similar calculation as SEIRS model for disease-free equilibrium, we obtain the next-generation matrix and its characteristic equation as Therefore the basic reproduction number for SEIR model is When the exposed period can be neglected and we don't need consider the state E, then we get SIRS model with pair approximation: We can consider {ρ S I , ρ II , ρ IR } as population densities of infected class, and noticing that the first and the third elements correspond to the changes of variables ρ S S and ρ S R , respectively, then the nextgeneration matrix and its characteristic equation become Therefore, the basic reproduction number for SIRS model is given as Next, we move to SIR model without the natural loss of immunity from R to S. SIR model with pair approximation is given as follows: By similar calculation as SIRS model for disease-free equilibrium, we obtain the next-generation matrix and its characteristic equation as Therefore, the basic reproduction number for SIR model is given by The last model SIS with pair approximation is given as follows: We can consider all the variables {ρ S I , ρ II } as population densities of infected class, then the nextgeneration matrix and its characteristic equation become Therefore, the basic reproduction number for SIS model is
SIR model can be obtained by taking the limit ofω → 0 for SIRS model (51). Therefore, R 0 for SIR model isβ which coincides with Eq (53). SIS model can be obtained by taking the limit ofω → ∞ for SIRS model (51) by regarding the transition from R to S as being instantaneous. So R 0 for SIS model is obtained as which coincides with Eq (54).

Comparison with Monte Carlo simulation
In this paper, we adopted an approximation to evaluate R 0 , so in order to check its accuracy, we would compare the result with the Monte Carlo simulation. To do this, we need an equation for some variable corresponding to R 0 , which increases for R 0 > 1 but decreases for R 0 < 1. Fortunately, we can find the variable for the SIR model, and so we explain it below.
SIR model gives the linearized equation of Eq (52) for ρ S I around disease free equilibrium gives The ratio of positive term to negative one (1 − 1/z)β/(1 +β/z) = (z − 1)β/(z +β) coincides with R 0 by Eq (53). We can expect the newly produced number of S-I pairs as this quantity in Monte Carlo simulation. Let us consider one infectious individual to invade into a susceptible population on two-dimensional square lattice space z = 4 ( Figure 3). There are four S-I pairs, then we trace how many S-I pairs are newly produced by the direct disease transmission of this infectious to the neighboring susceptibles before changing I to R.
(i) If the transition from I to R occurs, its probability is 1/(1+β) and no newly S-I pairs are produced.
(ii) If the transition from I to R does not occur and one of four Ss is infected, its probability is 1−1/(1+ β). Then if the transition from I to R occurs, the probability becomes 1 − 1/(1 +β) · 1/(1 + 3 4β ) and three S-I pairs are produced.
(iii) If the transition from I to R does not occur and one of the four Ss is infected. And then if the transition from I to R does not occur and one of the three Ss is infected. Then if the transition from I to R occurs, the probability becomes 1 − 1/(1 +β) · 1 − 1/(1 + 3 4β ) · 1/(1 + 2 4β ). 2 × 3 S-I pairs are produced.
(iv) If the transition from I to R does not occur and one of the four Ss is infected. And if the transition from I to R does not occur and one of the three Ss is infected. And if the transition from I to R does not occur and one of two Ss is infected. Then if the transition from I to R occurs, the ). 3 × 3 S-I pairs are produced.
(v) At the remaining probability 1 The total sum of (i)-(v) over 4 (averaged for one S-I pair) gives 3β/(4 +β), which is the case of z = 4 in Eq (53).
In the real simulation space there are more sites around these five sites and tertiary infections that can occur in these four Ss by the transmission from going around of secondary infections along the loop. So, the number of newly produced S-I pairs are expected to be somehow small. The results are shown in Figure 4.
R 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " e + 2 R p Y u s p i s Z k / 4 m A f 1 2 0 k k 8 V O s = " > A A A C t X i c j V E 9 S 8 N Q F D 2 N 3 9 9 V F 9 G l K I J T e X G p o + j i W D 9 q C 7 X U J L 5 q M E 1 C k g p S 6 i r 0 R z j p I K K 7 s 4 t / w M H F 3 d n B R d C T p A W 1 + H H D e 7 n v 3 n v O P Z e r u 5 b p B 0 I 8 J p S u 7 p 7 e v v 6 B w a H h k d G x 5 P j E t u / U P E P m D M d y v I K u + d I y b Z k L z M C S B d e T W l W 3 Z F 4 / X A 3 z + S P p + a Z j b w X H r i x V t X 3 b r J i G F j B U x Q W q 0 B D g A D 4 M e K h j A w 2 U I c r J O Z E W k a U 6 H b X l z C 1 P r z / t A s g 6 y T v s Y M Y J Y q V G S w j D V k k W O n E 5 z j G j d K R i k p e 0 o l L l U S L c w k v p j i f A B 7 G J m Z < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " e + 2 R p Y u s p i s Z k / 4 m A f 1 2 0 k k 8 V O s = " > A A A C t X i c j V E 9 S 8 N Q F D 2 N 3 9 9 V F 9 G l K I J T e X G p o + j i W D 9 q C 7 X U J L 5 q M E 1 C k g p S 6 i r 0 R z j p I K K 7 s 4 t / w M H F 3 d n B R d C T p A W 1 + H H D e 7 n v 3 n v O P Z e r u 5 b p B 0 I 8 J p S u 7 p 7 e v v 6 B w a H h k d G x 5 P j E t u / U P E P m D M d y v I K u + d I y b Z k L z M C S B d e T W l W 3 Z F 4 / X A 3 z + S P p + a Z j b w X H r i x V t X 3 b r J i G F j B U x Q W q 0 B D g A D 4 M e K h j A w 2 U I c r J O Z E W k a U 6 H b X l z C 1 P r z / t A s g 6 y T v s Y  Unfortunately we cannot succeed in finding a way of evaluating R 0 by Monte Carlo simulation for other models.

Discussion
We obtained the basic reproduction number for SEIRS model on the regular lattice as Eq (47) by using the next-generation matrix. Inaba [12] stated R 0 is determined uniquely, but from the viewpoint of [13], basic reproduction number by one method may not always agree with the one by another. Here we compare it to the results of other previous studies.
Ringa and Bauch [11] gives a basic reproduction number for SEIRS model on the lattice as find an appropriate quantity or variable to compare these results with Monte Carlo simulation except SIR model and then, we leave it as a future problem to estimate a basic reproduction number by Monte Carlo simulation for other basic epidemic models.
Here, we consider the basic reproduction number R 0 for fundamental epidemic models by using next-generation matrix only on a regular lattice. However, realistic epidemics may usually occur on heterogeneous networks, then we will study the effects of these spatial structures on R 0 by nextgeneration matrix approach in the future.