Computing R0 of dynamic models by a definition-based method

Objectives Computing the basic reproduction number (R0) in deterministic dynamical models is a hot topic and is frequently demanded by researchers in public health. The next-generation methods (NGM) are widely used for such computation, however, the results of NGM are usually not to be the true R0 but only a threshold quantity with little interpretation. In this paper, a definition-based method (DBM) is proposed to solve such a problem. Methods Start with the definition of R0, consider different states that one infected individual may develop into, and take expectations. A comparison with NGM has proceeded. Numerical verification is performed using parameters fitted by data of COVID-19 in Hunan Province. Results DBM and NGM give identical expressions for single-host models with single-group and interactive Rij of single-host models with multi-groups, while difference arises for models partitioned into subgroups. Numerical verification showed the consistencies and differences between DBM and NGM, which supports the conclusion that R0 derived by DBM with true epidemiological interpretations are better. Conclusions DBM is more suitable for single-host models, especially for models partitioned into subgroups. However, for multi-host dynamic models where the true R0 is failed to define, we may turn to the NGM for the threshold R0.


Introduction
The basic reproduction number (R 0 ) is 'one of the foremost and most valuable ideas that mathematical thinking has brought to epidemic theory' (Heesterbeek & Dietz, 1996). It is an epidemiology metric used to quantify the transmissibility of infectious diseases. R 0 is defined as the expectation of secondary cases produced by one infected individual in the entirely susceptible population, during its lifespan as infectious. The definition assumes that in the absence of any intervention (Becker et al., 2006), like quarantine, vaccination, etc., that is, all individuals are not immune. In its definition, R 0 is related to epidemiology characteristics, including the infectious period of a person who becomes infected, the frequency of contact during the infectious period, and the probability of transmission per contact between a susceptible individual and an infected person. Accordingly, the effectiveness of contact tracing and case isolation in controlling coronavirus disease 2019  depends partly on the precise estimation of R 0 , with R 0 of 1.5, outbreaks were controllable if 50% of contacts were traced. And with R 0 of 2.5 and 3.5, more than 70% and 90% of contacts had to be traced, respectively (Hellewell et al., 2020).
R 0 has the threshold property that if R 0 < 1, each infected individual produces, on average, less than one secondary case during its lifespan as infectious, and therefore the epidemic will gradually disappear, whereas if R 0 > 1, the disease may lead to a large number of new infections. Estimation of the correct R 0 is crucial in an epidemic emergency response because it informs later direction for control measures and can be compared to the effective reproduction number after interventions (Ma et al., 2022). Besides, it is also used to estimate the vaccine coverage needed for herd immunity (Garnett, 2005), the coverage should be greater than 1-1/R 0 (Desenclos, 2011;Fine et al., 2011), considering the effectiveness of vaccines is not 100%. And as transmissibility evolves, herd immunity thresholds can be estimated by modified R 0 (Gomes et al., 2022;Montalb an et al., 2020). Furthermore, following the results of the basic reproduction number, optimal control analyses were developed to reduce disease transmission, including co-infection (Tchoumi et al., 2021), as well as the most cost-efficient measures (Asamoah et al., 2022).
Since the concept of basic reproduction number was proposed by B€ ockh in 1886 till today (which George Macdonald applied it in modern epidemiology for the first time in 1952) (Heesterbeek et al., 2015;Heesterbeek, 2002;Heffernan et al., 2005), there are mainly two kinds of methods to calculate R 0 : driven by data and through dynamic models. Data-driven approaches, including intrinsic growth rate (Chowell et al., 2004), the final size equation (Arino et al., 2007), the average age of infection (Ajelli et al., 2008), etc. Through dynamic models, the methods separate into three general categories: Jacobian approach, the next-generation method (Martcheva, 2015), and estimate directly by R 0 's definition (Batabyal, 2021;Dearlove & Wilson, 2013;Getz et al., 2018;Hethcote, 2000). The next-generation method is likely the most frequently used method to calculate R 0 due to its generality, as R 0 is always computable if the compartment model is given. However, the definition, calculation, and interpretation of R 0 derive from NGM are anything but simple. Thoroughly grasping the mathematical method of R 0 calculation is still difficult for most practitioners and researchers in public health. Estimate R 0 directly through definition has been described in basic dynamic models without an explicit step-by-step derivation process (Batabyal, 2021;Dearlove & Wilson, 2013;Getz et al., 2018;Hethcote, 2000). Moreover, few of these methods agree with each other, even when based on the same dataset (Li & Blakeley, 2011). Therefore, this study formularizes and generalizes an intelligible and straightforward definition-based method (DBM) for computing R 0 of dynamic models of single host species for the first time, from which one can work out the expression of R 0 with just pen and paper. The result that this method produces is mutually coherent with the next-generation method (NGM) and is typically equivalent to the Van den Driessche and Watmough approach. Further, the result differs in computing R 0 for an entire population with multi-group models, while DBM gives the true epidemiological interpretation.

Modeling
We divide the compartmental dynamic models into four classes and associate each class with some applied models as examples. The assumptions, compartment definition, and parameter estimation to those models can be found in corresponding references. The flowchart, ordinary differential equations, and interpretations to variables and parameters are included in the supplementary materials. All parameters in compartmental models, except the transmission rate coefficient, which is usually obtained by data fitting, are obtained either from references or estimated using statistical inferences on original data.

Models of single host population divided by subgroups with single transmission route (SDS)
The multi-group SEIAR model of the previous studies ( Fig. 1)   is considered. For the universality and generality, we consider a modified model with a nonzero birth rate and death rate and assume the transition rate of A and I are nonequal.

Models of single host population with multi transmission routes (SM)
For single host models, there may be two or more routes of transmission. As the widely used Susceptible-Exposed-Infectious/Asymptomatic-Recovered-Water (SEIARW) model ( Fig. 2) (Chen, Leung, Zhou, et al., 2014), there are two routes of transmission: person to person and environment (water and food) to person.

Models of multi-host (MH)
The multi-host models involve two or more host populations of different species, like models for dengue fever (Yi et al., 2019), brucellosis , and severe fever with thrombocytopenia syndrome (Zhang et al., 2021), in which the transmission happened inside and between those different species. The vector-borne disease is a special case of multi-host models, that infections will never occur inside each one of host species, and can be formularized by setting the transmission coefficients inside each group be zero in the multi-host models.

Study design
We proposed generalized procedures of a definition-based method (DBM) for computing R 0 and compared it with the frequently used next-generation method (NGM) on models of SS, SDS, SM, MH. Numerical validation for SDS is proceeded using data of COVID-19 in Hunan Province. (Fig. 3). Flowchart of Multi Groups SEIAR model. The subscript i denotes the variable or parameter is specified to the i-th group. n denotes the total number of groups. Variables S i ; E i ; I i ; A i ; R i ; i ¼ 1; 2; /; n represent the susceptible, exposed, symptomatic infectious, asymptomatic infectious and recovered population in Nn is the total population size. Parameters L i are constant birth coefficients (for age-grouped model, only one L i of the youngest age group is nonzero); dr i denotes the mortality rate of group i; b ij is the coefficients describing the daily transmission rate from group i to group j; k is the relatively of transmission ability of asymptomatic cases compared with the symptomatic cases, it is assumed to be group irrelevant; f i is the case fatality rate of group i (if group irrelevant, then f ); p i is the probability that one infected individual in group i will developed into a asymptomatic case (if group irrelevant, then pÞ; u i is inverse of the average latent period of the symptomatic population in the i-th group, it is used to quantify the remove rate of compartment E i (the inverse of incubation period is usually adopted for practice); u 0 i is the inverse of the average latent period of the asymptomatic population in the i-th group; g i is the inverse of average infectious period for symptomatic cases; g 0 i is the inverse of average infectious period for asymptomatic. The solid arrows represent the transitions between variables. For an easier understanding of the proposed method, the procedure of a specific example of Susceptible-Acute infectious-Recovered/Chronic infectious (SIRC) model (Cui et al., 2020) (Fig. 4) will be proposed after each generalized step is stated.
Step 1. For any infected individual x, we start from the initial compartment of x [usually exposed (E) or infectious (I)] and compute the probability PðxÞ that x attains each infectious compartment. For the SIRC model, we have: Fig. 2. Flowchart of SEIARW model. Variables S; E; I; A; R, W represent population of susceptible, exposed, symptomatic infectious, asymptomatic infectious, recovered and pathogen in reservoir. Parameter L is a constant birth coefficient; b is the human-to-human transmission rate coefficient; b W is the reservoir-tohuman transmission rate coefficient; k is the relatively of transmission ability of asymptomatic cases compared with the symptomatic cases; dr is the natural death rate; f is the fatality; p represent the proportion of progressing to asymptomatic infectious stage; u is the inverse of the average latent period of the symptomatic cases; u 0 is the inverse of the average latent period of the asymptomatic cases; g is the inverse of the infectious period of symptomatic cases; g 0 is the inverse of average infectious period of asymptomatic cases. The solid arrows represent the transitions between variables; the dashed arrows represent pathogen shedding into environment. Fig. 3. Flowchart of study design. SS, SDS, SM, MH are models of four categories. Green checkmark arrows denote the method that works well for those models; yellow checkmark arrows denote the method is computable while it gives little epidemiological interpretation; the red cross mark arrow denotes the method is inadequate for those models.
Step 2. For all infectious compartments ½U 1 ;U 2 ;…;U r ] that x might attain, compute the newly infection rate for the case that x in U i ; i ¼ 1; 2; …; r; by letting U i equals 1, and other U j ; jsi equals zero in the term of newly infection (e:g: where f is a multi-variable function). For the SIRC model, by letting (I, C) ¼ (1, 0) and (0, 1) respectively, one obtains the secondary infection produced by x per unit time, for two cases that x in acute infectious (I) and x in chronic infectious (C): Step 3. Compute the length of time interval TðxÞ that x will stay in each infectious compartment using compartment removing rate (fatality of the disease, recovery rate, and natural death rate, usually set to be constants). For the SIRC model: Step 4. Assuming that the disease is currently at an initial stage of development, the number of infected takes only a small portion of the total population. This assumption implies Q ðxÞ can be treated as constant during the time interval TðxÞ.
Step 5. Taking expectations to eliminate the stochasticity of different infectious compartments x might develop into. The effective reproduction number is obtained: Step 6. The R 0 is obtained by assuming the entire population is susceptible and letting S ¼ N in the expression of R eff . For the SIRC model: Step 1, if the flowchart is bidirectional, then the computation of reachable probability may involve the assumption of Markov property and the use of the transition matrix.

Next-generation method (NGM)
The most frequently used Van den Driessche and Watmough approach in NGM is chosen as candidate method for comparison. The steps of NGM are listed as follows: Step 1. Divide all compartments into two categories: non-infected compartments (including susceptible, recovered, and fully immunized) and infected compartments Step 2. Write down the corresponding differential equations of the transition graph and divide the derivative vector of infected compartments into two parts: the first F is newly infected rate, and the second V is the rate of transition between different compartments.
Step 3. Take derivatives with respect to infected compartments for vectors F and V , and obtain the Jacobi matrices F and V: Step 4. Compute the next generation matrix FV À1 and its leading eigenvalue l max ðFV À1 Þ.
Step 5. Substitute the disease-free equilibrium into l max ðFV À1 Þ, and the R 0 defined by NGM is obtained.

Analysis and validation
The symbolic toolbox of MATLAB (R2021a) was used to perform the inverse and eigenvalue operations for symbolic matrices. Numerical experiments of DBM and NGM were performed using fitted transmission parameters and other parameters in a published age-specific multi-group SEIAR model . All codes of computation are available for download on GitHub and provided in supplementary materials (S1 Text: Supplementary Methods and Results).

R 0 for SS
For SS model, the results of DBM and NGM are mutually equivalent. Results of some popular models are listed as follows (derivations and detailed model frameworks are provided in S1 Text: Supplementary Methods and Results): SIR: SEIR: SEIAR: SIRC: 3.2. Interactive R ij and R 0 for SDS The interactive R ij between groups in a multi-group model, which is defined as the expectation of secondary infections that one infected individual in group i will produce in an entirely susceptible population of group j during its lifespan as infectious, is a quantity of great interest. All R ij forms a matrix that depicts transmission between groups.
We use a multi-group SEIAR model  as an example to illustrate the behaviors of DBM and NGM. The flowchart of this model is shown in (Fig. 1). The following derivation of DBM and NGM both leads to an identical result of R ij . And in constructing R 0 for the entire population via R ij matrix, we conclude that DBM are more suitable than NGM.

DBM for interactive R ij (SDS)
The derivation steps using DBM are as follows: Step 1. For any infected individual x i in group i, there are three possible states: E i , I i , and A i . The first step is to consider the initial state E i and compute the probability of x i passing other infectious compartments. It is shown in the transition graph ( Fig. 1), that for individual x i 2E i , it might develop into d r E i (natural death during exposure), I i and A i : Since Pðx i 2E i Þ ¼ 1; the probabilities are computed by ratio of transition rate given in the transition graph: where Pðx i 2E i Þ denotes the probability that individual x i will develop into the compartment E i ; and Pðx i 2I i jx i 2E i Þ denotes the conditional probability of x i will develop into the compartment I i with x i 2E i known.
Step 2. Letting one of three components of ðd r E i ; I i ; A i Þ equals to 1, and others equal 0 respectively and substitute into the term of newly infection rate b ij S j ðI i þ kA i Þ, one obtains the number of secondary infections in group j, that one infected individual x i will produce per unit time for the different compartments that x i might attain: Step 3. Take compartment I i as an example, if the death rates are omitted (d r ¼ f ¼ 0), then the infectious period of individual x i is 1=g (usually, the average course of infection is adopted). For the case that death rates are considered, we combined the terms ðd r þ f ÞI i and gI i , and hence the infectious period is 1=ðd r þ f þ gÞ. Similar procedures obtain the infectious period for other compartments of x i might attain: Step 4. At the beginning of disease transmission, SzN, and S=N is nearly a constant. As a good 1-order approximation, it is natural to assume that Q ðx i Þ is constant during the time interval Tðx i Þ.
Step 5. By taking expectations, the expectation of secondary infections in group j that one infected individual in group i will produce during its lifespan as infectious is given by: Step 6. Assuming group j is entirely susceptible, that by letting S j ¼ N j , one obtains the R ij defined in the beginning: 3.2.2. NGM for interactive R ij (SDS) The derivation steps using NGM are as follows: Step 1. Assume that the host population is divided into n groups. Firstly we divide all 5n compartments ðS i ; E i ; A i ; I i ; R i Þ k; i ¼ 1; 2; /; n into two categories: the first ðE i ; A i ; I i Þ k; i ¼ 1; 2; /; n are infected compartments, and the second ðS i ; R i Þ k; i ¼ 1; 2; /; n are non-infected compartments.
Step 2. Divide the derivation of ðE i ; A i ; I i Þ k; i ¼ 1; 2; /; n into two parts: the first part F denotes the rate of newly infection, and the second part V denotes the transition inside the infected compartments.
Step 3. Take derivatives of 3n-dimension vector-valued function F and V , respect variables ðE i ; A i ; I i Þ; i ¼ 1; 2; /; n. where d ij is the Kronecker Delta, i.e., d ij ¼ 1 for i ¼ j, and d ij ¼ 0 for isj.
We call F ij the newly infection matrix of group j to group i, call V ii be the transition matrix of group i. The inverse of V ii is further computed: V À1 ii ¼ 2 6 6 6 6 6 6 6 6 4 1 d r þ u À u p þ u 0 p 0 0 u 0 p ðd r þ g 0 Þ ðd r þ u À u p þ u 0 pÞ 1 d r þ g 0 0 u ð1 À pÞ ðd r þ f i þ gÞ ðd r þ u À u p þ u 0 pÞ 0 1 d r þ f i þ g 3 7 7 7 7 7 7 7 7 5 Step 4.  Step 5. We use the leading eigenvalue of M ij at the disease-free equilibrium to measure the infectivity of group j to group i: We can notice that equations (5) and (6) are mutually equivalent, that is, NGM is equivalent to DBM for the interactive R ij of this model.

DBM for R 0 (SDS)
Let Rðx i ; jÞ denote the random variable describe the secondary infections that individual x i in group i will produce in group j during its lifespan as infectious. We define group 0 as the entire population. For example, Rðx 0 ; 0Þ denotes the number of secondary infections that one infected individual x 0 (arbitrary in group 0, i.e., the entire population) will produce in the entire population during its lifespan as infectious. We use the tower property to compute the expectation of Rðx 0 ; 0Þ: E Rðx i ;jÞ; i¼1;2;/;n fEfRðx 0 ; jÞjRðx i ; jÞ; i ¼ 1; 2; /; ngg X. Guo, Y. Guo, Z. Zhao et al. Infectious Disease Modelling 7 (2022) 196e210 Here Pðx 0 2N i Þ denote the probability of an arbitrarily selected individual x 0 belongs to group i, and it is determined by the proportion of the size of group i. Random variable Rðx 0 ; 0Þ contains the stochasticity of x 0 may come from different groups and the stochasticity of x 0 may develop into different compartments. However, the random variable Rðx i ; jÞ contains only the stochasticity that x 0 may develop into different compartments.
We can see that the expectation results in a weighted averaging of interactive R ij .

NGM for R 0 (SDS)
In the NGM, R 0 is defined as the leading eigenvalue of the next-generation matrix M: In most cases, the leading eigenvalue cannot be formulated as eigenvalues of submatrices M ij , however, there are circumstances that l max ðMÞ still have a simple and comprehensive formulation.
identical to V 11 . Therefore: Using properties of Kronecker product (Petersen & Pedersen, 2008), one obtains: We can see that, with the simplification that the natural death rates of each group are equal, the R 0 derived by NGM can be formulated as the eigenvalue of matrix of the transmission coefficient b ij multiplies a group-irrelevant coefficient, that is, the eigenvalue of the R ij matrix. One may notice that NGM and DBM are both perform somehow 'averaging' to the interactive R ij . The weighted averaging is adopted for DBM, while NGM averages R ij by taking leading eigenvalue.

Computing by DBM (SM)
The secondary infections produced by one infected individual during its infectious period can be divided into two parts: the secondary infection generated by environmental media that this infected individual produced, and the secondary infections generated by direct contact (no environmental media involved). We compute those two parts respectively and then take summation. The steps are as follows: Step 1. Since Pðx 2EÞ ¼ 1, then for all infected compartments, the probability that one infected individual x will attain is given by: Step 2. Letting ðI; AÞ ¼ ð1; 0Þ and ð0; 1Þ in the transmission rate bSðI þkAÞ respectively, one obtains the number of secondary infections that individual x will produce per unit time by direct contact: Step 3. The infectious period of x is given by: Step 4. The secondary infections that individual x produced directly (without environmental media) during its lifespan as infectious are: Step 5. The producing rate of infectious environmental media that individual x produces is: And the amount of infectious environmental media that an individual x produced during its lifespan as infectious is: Step 6. Assuming that the current infectious environmental media is much greater than the amount that one individual can produce. This implies the infectivity of environmental media can be viewed as a constant for time interval TðxÞ. Hence the secondary infections that LðxÞ produces can be computed by letting W ¼ LðxÞ in b W SW, and multiplies the inverse of removing rate 1=ε: Step 7. Take the summation, the effective reproduction number is given by: Step 8. R 0 is obtained by assuming that all individuals are susceptible, and substitutes S ¼ N. Such a result is mutually equivalent to the result of NGM. (Derivation of NGM can be found in codes of supplementary materials) 3.5. Data validation for R ij and R 0 According to the interactive b ij of 4 time-segments fitted on the basis of real-world COVID-19 data, the matrices of R ij of different time segments are evaluated by DBM in Fig. 5 (mutually equivalent to NGM). The R 0 for the entire population is also evaluated by two methods. DBM gives R 0 ¼ 1:7476, while NGM gives R 0 ¼ 1:5032. Such a difference is explained by the property of NGM: the further consideration of the transmissibility of secondary cases for continuous transmission makes result of NGM serves only as a threshold for the stability of disease-free equilibrium, and it is hard to explain in cases involving subgroups.

Discussion
DBM is formatted and generalized in the study, and further compared with NGM. It is found that: (i) For single-host dynamic models with a single group, results of DBM are identical to NGM. (ii) For single-host dynamic models with multi groups, the results of two methods in computing interactive R ij are identical as well. However, in computing R 0 for an entire population, the results differ, while DBM gives the true epidemiological interpretation.
Some studies adopted calculated R 0 to estimate the transmissibility of infectious diseases, such as shigellosis (Z. Zhao, Chen, et al., 2021), HFMD , and COVID-19 (Zhao et al., 2021a). However, the detailed steps of deriving the R 0 expression are not illustrated. Therefore, we described the derivations of DBM in detail, and extended the procedure to more complicated models, making computation of R 0 accessible to the practitioners and researchers of public health. Also in these studies, only the parameters of the natural history of disease were included in dynamic models, whereas parameters of demographics were not considered. In this case, the calculated R 0 is more suitable for an outbreak of infectious disease. When it comes to estimating the transmissibility of infectious disease through surveillance data of continuous years, however, the lack of demographic parameters could cause overestimation of true R 0 . This overestimation can act as a barrier against our understanding towards the transmission features of disease and cause excessive investments on the interventions. Hence, we developed DBM to compute R 0 of single-host dynamic models, which parameters of demographic are sufficiently considered. The true R 0 is therefore achievable both in outbreaks and regular preventions.

Comparison between DBM and NGM
The difference occurs while calculating R 0 for the entire population with single-host dynamic models containing multi groups, by adopting DBM and NGM, respectively. Since the definition-based method uses the number of individuals in each subgroup as a proportion of the number of individuals in the entire host population to weight average R ij , we can arbitrarily Fig. 5. The R ij matrices for four segments. The real-world COVID-19 data we used for validation was collected by Hunan Provincial CDC from January 5 to February 19, 2020. The data included patient gender, age, inter-provincial travel history, case type (symptomatic/asymptomatic), exposure date, date of onset, and date of diagnosis. Berkeley Madonna was employed to do the least root-mean-square deviation curve fitting. And for the sensitivity analysis, six parameters were employed to examine the model's sensitivity in this study, each of which was divided into 1000 values based on its range, mean and standard deviation data were calculated. The R ij Matrices for 4-time segments (defined as stages in the disease transmission period). Parameters for calculating are obtained from published research. Subfigure A represents time segment 1, from Jan 5, 2021 to Jan 25, 2021; subfigure B represents time segment 2, from Jan 25, 2021 to Jan 31, 2021; subfigure C represents time segment 3, from Jan 31, 2021 to Feb 5, 2021; subfigure D represents time segment 4, from Feb 5, 2021 to Feb 19, 2021. select group i to obtain the R i (expectation of secondary infectious that one infectious individual in subgroup i produces in the whole host population). In contrast, in the NGM, the R 0 for the entire population is defined as the largest eigenvalue of the next-generation matrix and using the largest eigenvalue for 'average' has little epidemiological interpretation. Therefore, DBM is more suitable than NGM in the single-host dynamic models with multi groups.

Application in epidemics
In an epidemic outbreak, associated implications are urgently needed. If the data-driven approach of calculating R 0 is adopted, it can take a long time to collect enough data, and sometimes R 0 cannot be obtained until the end of the outbreak, for example, calculating R 0 by the final size equation requires the proportion of susceptible individuals in the whole population at the end of the outbreak (Brauer, 2008). Whilst calculating R 0 based on the dynamic compartment model, where some of the parameters can be obtained from a review of literature and others can be fitted to data from the early stages of an outbreak. This method allows for a rapid assessment of outbreak trends and informs the interventions required for control. Thus, the DBM of deriving R 0 applies to many infectious disease outbreaks in which single-host dynamic models can be established correspondingly, such as influenza (Chen et al., 2017), COVID-19, Ebola (Chen, Leung, Liu, et al., 2014), shigellosis, HFMD, etc. Besides, the derivation steps by DBM let practitioners gain a better understanding of the infectious disease dynamics. As we take apart the definition of R 0 , parameters that are essential in dynamic models may be plausible in interpreting public health interests. For example, the transmission rate and recovery rate determine overall speed and attack rate of epidemic, latent period determines duration of an epidemic. Interventions targeted at contact rate, infectious period and quarantine affect the transmission rate, recovery rate, and latent period, respectively (Ridenhour et al., 2018). The epidemiological application of R 0 will only become valuable with accurate and cautious interpretations.

Other methods for reproduction number
There are three kinds of reproduction numbers: basic reproduction number R 0 , effective reproduction number R eff, and time-varying or real-time reproduction number R(t). Both R 0 and R eff are functions of status (population status of entirely susceptible, non-intervention and social status of the certain economy and contact pattern for R 0 ; specific status for R eff ), while R(t) is a function of time. Once the status is expressed as functions of time (e.g., numerical solution of an ordinary differential equation), then we can immediately get R(t) from R eff by substitution. If the status at time 0 is simplified as entirely susceptible, non-intervention, then R 0 can be obtained by letting t ¼ 0 in R(t), i.e., R(0). Popular methods for time-varying R(t) use the distribution of generation time as a simplification of the natural history of disease and get free from the ordinary differential equations (Batista, 2021; Cori et al., 2013;Obadia et al., 2012;Thompson et al., 2019;Wallinga & Lipsitch, 2007). These methods are direct and fast implemented, while the resulting R(0) is only a reproduction number of initial time instance of an outbreak, and is not suitable to be interpreted as 'basic' reproduction number, especially with significant vaccination and NPIs (non-pharmacological interventions) at the initial time instance (the situation we faced in local outbreaks of . A more reasonable way for computing R 0 is to develop models with vaccination and intervention, and derive R 0 from the calibrated and degenerated model (degenerate by letting the intensity of vaccination and intervention be zero).

Limitations
There are still some limitations of calculating R 0 by DBM, which can only be applied to dynamic models with a single host population. For multi-host compartment models or co-dynamics models (Tchoumi et al., 2021), it is not able to establish a reasonable definition of R 0 , so the universal NGM is frequently used to get the threshold quantity of stability of disease-free equilibrium. In addition, in SIRC (C for chronic infection) models (Cui et al., 2020;Wang et al., 2020) for hepatitis C established by previous studies, the researchers used The Castillo-Chavez, Feng, and Huang approach of the next-generation matrix method to consider the changing of natural birth and death rates in susceptible populations over time, the result corresponds to making the susceptible population of R eff obtained by the definition-based method equal to the extreme of the total population, which is more applicable for chronic infectious diseases. Yet if we use Van den Driessche and Watmough approach of the next-generation method (NGM) to calculate R 0 in this case, the result is still equivalent to the definition-based method, if it is assumed that the natural birth and death rate are balanced, the results of the three methods are consistent. At the same time, because R 0 itself has its limitations, model builders who want to compare the transmissibility between diseases based on R 0 should also consider: Even for the same disease, there can be differences in population structure and social behavior in different periods and regions, adjustments should be made to the corresponding parameter assumptions (Delamater et al., 2019). Variables that affect contact patterns, such as birth rate, death rate, population density, cultural practices, etc., are usually similar within a region but vary in different areas (Guerra et al., 2017). Therefore, further studies might adopt more reasonable assumptions about the parameters in different regions at each stage if contact rate surveys are conducted in different regions at regular intervals.